ETHZ: 227-0966-00L¶

Quantitative Big Imaging¶

Image Enhancement¶

February 27th, 2020¶

Anders Kaestner¶

Measurements are rarely perfect¶

Factors affecting the image quality¶

• Resolution (Imaging system transfer functions)
• Noise
• Contrast
• Inhomogeneous contrast
• Artifacts

No measurement system is perfect and therefore, you will also not get perfect images of the sample you measured. In the figure you can see the slice as it would be in a perfect world to the left and what you actually would measure with an imaging system. Some are of higher quality than others but there are still a collection of factors that have an impact on the image quality.

• Resolution Thanks to the optical system there will be some unsharpness caused by scintillator, lenses, and camera. this appears as blurring or the acquired images.
• Noise There are different sources of noise when an image is acquired. Most of them are Poisson/Binomial distributed point processes, but thanks to the light spread in the scintillator there can also be some spatially correlated noise.
• Contrast Sometimes there is a very low intensity difference between different regions of the sample. This is in particularly problematic when the signal to noise ratio is low.
• Inhomogeneous contrast The interaction between radiation and matter may introduce effects that appear as an inhomogeneous contrast. Typical such effects are beam hardening and scattering.
• Artifacts The artifacts are features in the image that are not present in the sample. Often this is caused by detector errors or interference of radiation background or seconday radiation events.

A typical processing chain¶

Today's lecture will focus on the enhancement

Traditionally, the processing of image data can be divided into a series of sub tasks that provide the final resul.

• Acquisition The data must obviously be acquired and stored. There are cases when simulated data is used. Then, the acquisition is replaced by the process to simulate the data.
• Enhancement The raw data is usually not ready to be processed in the form is comes from the acquisition. It usually has noise and artifacts as we saw on the previous slide. The enhancement step suppresses unwanted information in the data.
• Segmentation The segmenation identifies different regions based on different features such as intensity distribution and shape.
• Post processing After segmentation, there may be falsely identified regions. These are removed in a post processing step.
• Evaluation The last step of the process is to make conclusions based on the image data. It could be modelling material distrbutions, measuring shapes etc.

Noise and artifacts¶

The unwanted information¶

Noise types¶

• Spatially uncorrelated noise
• Event noise
• Stuctured noise

Noise examples¶

Gaussian noise¶

• Additive
• Easy to model
• Law of large numbers

Distribution function¶

$n(x)=\frac{1}{\sqrt{2\pi\sigma}}\exp{-\left(\frac{x-\mu}{2\sigma}\right)^2}$

Poisson noise¶

• Multiplicative
• Physically correct for event counting

Distribition function¶

$p(x)=\frac{\lambda^{k}}{k!} e^{-\lambda\,x}$

Compare Gaussian and Possion noise¶

Samples of two distributions¶

Noise models - Salt'n'pepper noise¶

• A type of outlier noise
• Noise frequency described as probability of outlier
• Can be additive, multiplicative, and independent replacement

Example model¶

$sp(x)=\left\{\begin{array}{ll} -1 & x\leq\lambda_1\\ 0 & \lambda_1< x \leq \lambda_2\\ 1 & \lambda_2<x \end{array}\right.\qquad \begin{array}{l}x\in\mathcal{U}(0,1)\\\lambda_1<\lambda_2\\ \lambda_1+\lambda_2 = \mbox{noise fraction} \end{array}$

Salt'n'pepper examples¶

Signal to noise ratio¶

The Signal to noise ratio measures the noise strengh in a signal

Definition¶

$SNR=\frac{mean(f)}{stddev(f)}$

Signal to noise ratio for Poisson noise¶

• For a Poisson distribution the SNR is :

$SNR \sim \sqrt{N}$

• $N$ is the number of particles $\sim$ exposure time

Useful python functions¶

Random number generators [numpy.random]¶

Generate an $m \times n$ random fields with different distributions:

• Gauss np.random.normal(mu,sigma, size=[rows,cols])
• Uniform np.random.uniform(low,high,size=[rows,cols])
• Poisson ```np.random.poisson(lambda, size=[rows,cols])

Statistics¶

• np.mean(f), np.var(f), np.std(f) Computes the mean, variance, and standard deviation of an image $f$.
• np.min(f),np.max(f) Finds minimum and maximum values in $f$.
• np.median(f), np.rank() Selects different values from the sorted data.

Basic filtering¶

What is a filter?¶

General definition¶

A filter is a processing unit that

• Enhances the wanted information
• Suppresses the unwanted information

Filter characteristics¶

Filters are characterized by the type of information they suppress

Low-pass filters¶

• Slow changes are enhanced
• Rapid changes are suppressed

High-pass filters¶

• Rapid changes are enhanced
• Slow changes are suppressed

Basic filters¶

Linear filters¶

Computed using the convolution operation

$$g(x)=h*f(x)=\int_{\Omega}f(x-\tau) h(\tau) d\tau$$

where

• $f$ is the image
• $h$ is the convolution kernel of the filter

Low-pass filter kernels¶

Mean or Box filter	Gauss filter
All weights have the same value.	$$G=\exp{-\frac{x^2+y^2}{2\,\sigma^2}}$$
Example: $$B=\frac{1}{25}\cdot\begin{array}{|c|c|c|c|c|} \hline 1 & 1 & 1 & 1& 1\\ \hline 1 & 1 & 1 & 1& 1\\ \hline 1 & 1 & 1 & 1& 1\\ \hline 1 & 1 & 1 & 1& 1\\ \hline 1 & 1 & 1 & 1& 1\\ \hline \end{array} $$	Example:

Different SNR using a Gauss filter¶

How is the convolution computed¶

Euclidean separability¶

The asociative and commutative laws apply to convoution

$$(a * b)*c=a*(b*c) \quad \mbox{ and } \quad a * b = b * a $$

A convolution kernel is called separable if it can be split in two or more parts:

$$\begin{array}{|c|c|c|} \hline \cdot & \cdot & \cdot\\ \hline \cdot & \cdot & \cdot\\ \hline \cdot & \cdot& \cdot\\ \hline \end{array}= \begin{array}{|c|} \hline \cdot \\ \hline \cdot \\ \hline \cdot \\ \hline \end{array} * \begin{array}{|c|c|c|} \hline \cdot & \cdot & \cdot\\ \hline \end{array}$$ $$\exp{-\frac{x^2+y^2}{2\,\sigma^2}}=\exp{-\frac{x^2}{2\,\sigma^2}}*\exp{-\frac{y^2}{2\,\sigma^2}}$$

Gain¶

Separability reduces the number of computations $\rightarrow$ faster processing

• 3$\times$3 $\rightarrow$ 9 mult and 8 add $\Leftrightarrow$ 6 mult and 4 add
• 3$\times$3$\times$3 $\rightarrow$ 27 mult and 26 add $\Leftrightarrow$ 9 mult and 6 add

The median filter¶

Comparing filters for different noise types¶

Filter example: Spot cleaning¶

Problem	Example	Possible solutions
Many neutron images are corrupted by spots that confuse following processing steps. The amount, size, and intensity varies with many factors.		Low pass filter Median filter Detect spots and replace by estimate

Spot cleaning algorithm¶

Parameters¶

• $N$ Width of median filter.
• $k$ Threshold level for outlier detection.

Spot cleaning - Compare performance¶

The ImageJ ways¶

• Despeckle Median ... please avoid this one!!!
• Remove outliers Similar to cleaning described algorithm

High-pass filters¶

High-pass filters enhance rapid changes $\rightarrow$ ideal for edge detection

Typical high-pass filters:¶

Gradients¶

$$\frac{\partial}{\partial\,x}=\frac{1}{2}\cdot\begin{array}{|c|c|} \hline -1 & 1\\ \hline \end{array}\qquad \frac{\partial}{\partial\,x}=\frac{1}{32}\cdot\begin{array}{|c|c|c|} \hline -3 & 0 & 3\\ \hline -10 & 0 & 10\\ \hline -3 & 0 & 3\\ \hline \end{array} $$

Laplacian¶

$$ \bigtriangleup=\frac{1}{2}\cdot\begin{array}{|c|c|c|} \hline 1 & 2 & 1\\ \hline 2 & -12 & 2\\ \hline 1 & 2 & 1\\ \hline \end{array} $$

Sobel¶

$$ G=|\nabla f|=\sqrt{\left(\frac{\partial}{\partial\,x}f\right)^2 + \left(\frac{\partial}{\partial\,y}f\right)^2} $$

Gradient example¶

Vertical edges¶

$\frac{\partial}{\partial x}=\frac{1}{32}\cdot\begin{array}{|c|c|c|} \hline -3 & 0 & 3\\ \hline -10 & 0 & 10\\ \hline -3 & 0 & 3\\ \hline \end{array}$

Horizontal egdges¶

$\frac{\partial}{\partial\,y}=\frac{1}{32}\cdot\begin{array}{|c|c|c|} \hline -3 & -10 & -3\\ \hline 0 & 0 & 0\\ \hline 3 & 10 & 3\\ \hline \end{array}$

Jaehne, 2005

Edge detection examples¶

Relevance of filters to machine learning¶

[NVIDIA Developer zone](https://developer.nvidia.com/discover/convolutional-neural-network)

Frequency space filters¶

Applications of the Fourier transform¶

The Fourier transform¶

Transform¶

$$G(\xi_1,\xi_2)=\mathcal{F}\{g\}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} g(x,y) \exp{-i(\xi_1\,x+\xi_2\,y)}\,dx\,dy$$

Inverse¶

$$ g(x,y)=\mathcal{F}^{-1}\{G\}=\frac{1}{(2\,\pi)^2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} G(\omega) \exp{i(\xi_1\,x+\xi_2\,y)}\,d\xi_1 \,d\xi_2$$

FFT (Fast Fourier Transform)¶

In practice - you never see the transform equations. The Fast Fourier Transform algorithm is available in numerical libraries and tools.

Jaehne, 2005

Some mathematical features of the FT¶

Additition¶

$\mathcal{F}\{a+b\} = \mathcal{F}\{a\}+\mathcal{F}\{b\}$

Convolution¶

$\mathcal{F}\{a * b\} = \mathcal{F}\{a\} \cdot \mathcal{F}\{b\} $ $\mathcal{F}\{a \cdot b\} = \mathcal{F}\{a\} * \mathcal{F}\{b\} $

Additive noise in Fourier space¶

Real space¶

Fourier space¶

Problem¶

How can we suppress noise without destroying relevant image features?

Spatial frequencies and orientation¶

Example - Stripe removal in Fourier space¶

• Transform the image to Fourier space

$$\mathcal{F_{2D}}\left\{a\right\} \Rightarrow A$$

• Multiply spectrum image by band pass filter

$$A_{filtered}=A \cdot H$$

• Compute the inverse transform to obtain the filtered image in real space

$$\mathcal{{F_{2D}}^{-1}}\left\{A_{filtered}\right\} \Rightarrow a_{filtered}$$

The effect of the stripe filter¶

Reconstructed CT slice before filter	Reconstructed CT slice after stripe filter

Intensity variations are suppressed using the stripe filter on all projections.

Technical details on Fourier space filters¶

When should you use convolution in Fourier space?¶

• Simplicity
• Kernel size
• Speed at repeated convolutions

Zero padding¶

The FFT is only working with data of size in $2^N$. If your data has a different length, you have to pad (fill with constant value) up the next $2^N$.

Python functions¶

Filters in the spatial domain¶

e.g. from scipy import ndimage

• ndimage.filters.convolve(f,h) Linear filter using kernel $h$ on image $f$.
• ndimage.filters.median_filter(f,\[n,m\]) Median filter using an $n \times m$ filter neighborhood

Fourier transform¶

• np.fft.fft2(f) Computes the 2D Fast Fourier Transform of image $f$
• np.fft.ifft2(F) Computes the inverse Fast Fourier Transform $F$.
• np.fft.fftshift() Rearranges the data to center the $\omega$=0. Works for 1D and 2D.

Complex numbers¶

• np.abs(f), np.angle(f) Computes amplitude and argument of a complex number.
• np.real(f), np.imag(f) Gives the real and imaginary parts of a complex number.

Scale spaces¶

Why scale spaces?¶

Motivation¶

Basic filters have problems to handle low SNR and textured noise.

Something new is required...

The solution¶

Filtering on different scales can take noise suppression one step further.

By Cmglee - Own work, CC BY-SA 3.0, Link

Wavelets - the basic idea¶

• The wavelet transform produces scales by decomposing a signal into two signals at a coarser scale containing \emph{trend} and \emph{details}.\
• The next scale is computed using the trend of the previous transform

$$WT\{s\}\rightarrow\{a_1,d_1\}, WT\{a_1\}\rightarrow\{a_2,d_2\}, \ldots, WT\{a_{N-1}\}\rightarrow\{a_N,d_N\} $$

• The inverse transform brings $s$ back using $\{a_N,d_1, \ldots,d_N\}$.
• Many wavelet bases exists, the choice depends on the application.

Applications of wavelets¶

• Noise reduction
• Analysis
• Segmentation
• Compression

Walker 2008 Mallat 2009

Wavelet transform of a 1D signal¶

Using symlet-4

Wavelet transform of an image¶

Wavelet transform of an image - example¶

Original	Wavelet transformed

Wavelet noise reduction¶

The noise is found in the detail part of the WT

• Make a WT of the signal to a level that corresponds to the scale of the unwanted information.
• Threshold the detail part $d_{\gamma}=|d|<\gamma \,?\, 0 : d$.
• Inverse WT back to normal scale $\rightarrow$ image is filtered.

Wavelet noise reduction - Image example¶

Example Filtered using two levels of Symlet-2 wavelet¶

Data: Neutron CT of a lead scroll

Python functions¶

dwt2/idtw2¶

Makes one level of the wavelet transform or its inverse using wavelet base specified by 'wn'.

wavedec2¶

Performs N levels of wavelet decomposition using a specified wavelet base.

wbmpen¶

Estimating threshold parameters for wavelet denoising.

wdencmp¶

Wavelet denoising and compression using information from wavedec2 and *wbmpen.

Parameterized scale spaces¶

PDE based scale space filters¶

These filters may work for applications where Linear and Rank filters fail.

Aubert 2002.

The starting point¶

The heat transport equation $$\frac{\partial T}{\partial t}=\kappa\,\nabla^2 T$$

• $T$ Image to filter (intensity $\equiv$ temperature)
• $\kappa$ Thermal conduction capacity

Original	Iterations

Controlling the diffusivity¶

We want to control the diffusion process...

• Near edges The Diffusivity $\rightarrow$ 0
• Flat regions The Diffusivity $\rightarrow$ 1

The contrast function $G$ is our control function $$G(x)=\frac{1}{1+\left(\frac{x}{\lambda}\right)^n}$$

• $\lambda$ Threshold level
• $n$ Steepness of the threshold function

Gradient controlled diffusivity¶

$$\frac{\partial u}{\partial t}=G(|\nabla u|)\,\nabla^2 u$$

Image	Diffusivity map

• $u$ Image to be filtered
• $G(\cdot)$ Non-linear function to control the diffusivity
• $\tau$ Time increment
• $N$ Number of iterations

The non-linear diffusion filter¶

A more robust filter is obtained with

$$\frac{\partial u}{\partial t}=G(|\nabla_{\sigma} u|)\,\nabla^2 u$$ - _$u$_ Image to be filtered - _$G(\cdot)$_ Non-linear function to control the contrast - _$\tau$_ Time increment per numerical iteration - _$N$_ Number of iterations - _$\nabla_{\sigma}$_ Gradient smoothed by a Gaussian filter, width $\sigma$

Diffusion filter example¶

Neutron CT slice from a real-time experiment observing the coalescence of cold mixed bitumen.

Original	Iterations of non-linear diffusion

Filtering as a regularization problem¶

The continued development¶

• 90's During the late 90's the diffusion filter was described in terms of a regularization problem.
• 00's Work toward regularization of total variation minimization.

TV-L1¶

$$u=\underset{u\in BV(\Omega)}{\operatorname{argmin}}\left\{\underbrace{|u|_{BV}}_{noise}+ \underbrace{\mbox{$\frac{\lambda}{2}$}\|f-u\|_{1}}_{fidelity}\right\}$$

Rudin-Osher-Fatemi model (ROF)¶

$$u=\underset{u\in BV(\Omega)}{\operatorname{argmin}}\left\{\underbrace{|u|_{BV}}_{noise}+ \underbrace{\mbox{$\frac{\lambda}{2}$}\|f-u\|^2_{2}}_{fidelity}\right\}$$

with $|u|_{BV}=\int_{\Omega}|\nabla u|^2$

The inverse scale space filter¶

The idea¶

We want smooth regions with sharp edges\ldots

• Turn the processing order of scale space filter upside down
• Start with an empty image
• Add large structures successively until an image with relevant features appears

The ISS filter - Some properties¶

• is an edge preserving filter for noise reduction.
• is defined by a partial differential equation.
• has a well defined termination point.

Burger et al. 2006

The ROF filter equation¶

The image $f$ is filtered by solving

$$\begin{eqnarray} \frac{\partial u}{\partial t}&=& \mathrm{div}\left(\frac{\nabla u}{|\nabla u|}\right)+\lambda\,(f-u+v)\nonumber\\ \frac{\partial v}{\partial t}&=& \alpha\, (f-u) \end{eqnarray}$$

Variables:¶

• $f$ Input image
• $u$ Filtered image
• $v$ Regularization term (feedback of previous iteration)

Filter parameters¶

• $\lambda$ Related to the scale of the features to suppress.
• $\alpha$ Quality refinement
• $N$ Number of iterations
• $\tau$ Time increment

Filter iterations¶

Neutron CT of dried lung filtered using 3D ISS filter

Original	Filter iterations

How to choose $\lambda$ and $\alpha$¶

The requirements varies between different data sets.

Initial conditions:¶

• Signal to noise ratio
• Image features (fine grained or wide spread)

Experiment:¶

• Scan $\lambda$ and $\alpha$
• Stop at $T=n\, \tau= \sigma$ use different $\tau$
• When does different effects occur, related to $\sigma$?

Solutions at different times¶

1	30	60	100	200	500	999

Solution time¶

The solution time ($T=N\cdot\tau$) is essential to the result

• $\tau$ large The solution is reached fast
• $\tau$ small The numerical accuracy is better

The choice of initial image¶

At some $T$ the solution with $u_0=f$ and $u_0=0$ converge.

Non-local means¶

Non-local smoothing¶

The idea¶

Smoothing normally consider information from the neighborhood like

• Local averages (convolution)
• Gradients and Curvatures (PDE filters)

Non-local smoothing average similiar intensities in a global sense.

• Every filtered pixel is a weighted average of all pixels.
• Weights computed using difference between pixel intensities.

Buades et al. 2005

Filter definition¶

The non-local means filter is defined as $$u(p)=\frac{1}{C(p)}\sum_{q\in\Omega}v(q)\,f(p,q)$$ where

• $v$ and $u$ input and result images.

• $C(p)$ is the sum of all pixel weights as

  $C(p)=\sum_{q\in\Omega}f(p,q)$

• $f(p,q)$ is the weighting function

  $f(p,q)=e^{-\frac{|B(q)-B(p)|^2}{h^2}}$

• B(x) is a neighborhood operator e.g. local average around $x$

Non-local means 2D - Example¶

Observations¶

• Good smoothing effect.
• Strong thin lines are preserved.
• Some patchiness related to filter parameter $t$, i.e. the size of $\Omega_i$.

Performance complications¶

Problem¶

The orignal filter compares all pixels with all pixels\ldots\

• Complexity $\mathcal{O}(N^2)$
• Not feasible for large images, and particular 3D images!

Solution¶

It has been shown that not all pixels have to be compared to achieve a good filter effect. i.e. $\Omega$ in the filter equations can be replaced by $\Omega_i<<\Omega$

Verification¶

How good is my filter?¶

Verify the correctness of the method¶

"Data massage"¶

Filtering manipulates the data, avoid too strong modifications otherwise you may invent new image features!!!

Verify the validity your method¶

• Visual inspection
• Difference images
• Use degraded phantom images in a "smoke test"

Verification using difference images¶

Compute pixel-wise difference between image $f$ and $g$

Difference images provide first diagnosis about processing performance

Performance testing - The smoke test¶

• Testing term from electronic hardware testing - drive the system until something fails due to overheating...
• In general: scan the parameter space for different SNR until the method fails to identify strength and weakness of the system.

Test strategy¶

1. Create a phantom image with relevant features.
2. Add noise for different SNR to the phantom.
3. Apply the processing method with different parameters.
4. Measure the difference between processed and phantom.
5. Repeat steps 2-4 $N$ times for better test statistics.
6. Plot the results and identify the range of SNR and parameters that produce acceptable results.

Data for evaluation - Phantom data¶

General purpose can be controlled

• Data with known features.
• Parameters can be changed.
  • Shape
  • Sharpness
  • Contrast
  • Noise (distribution and strength)

Data for evaluation - Labelled data¶

Often 'real' data

• Labeled by experts
• Used for training and validation
  • Training of model
  • Validation
  • Test

Evaluation metrics for images¶

An evaluation procedure need a metric to compare the performance

Mean squared error¶

$$MSE(f,g)=\sum_{p\in \Omega}(f(p)-g(p))^2$$

Structural similarity index¶

$$SSIM(f,g)=\frac{(2\mu_f\,\mu_g+C_1)(2\sigma_{fg}+C_2)}{(\mu_f^2+\mu_g^2+C_1)(\sigma_f^2+\sigma_g^2+C_2)}$$

• $\mu_f$, $\mu_g$ Local mean of $f$ and $g$.
• $\sigma{fg}$_ Local correlation between $f$ and $g$.
• $\sigma_f$, $\sigma_g$ Local standard deviation of $f$ and $g$.
• $C_1$, $C_2$ Constants based on the image dynamics (small numbers).

$$MSSIM(f,g)=E[SSIM(f,g)]$$

Wang 2009

Test run example¶

Phantom structures¶

Change SNR and contrast¶

Process¶

Plot results¶

Kaestner et al. 2006

Overview¶

Many filters¶

Details of filter performance¶

Take-home message¶

We have looked at different ways to suppress noise and artifacts:

• Convolution
• Median filters
• Wavelet denoising
• PDE filters

Which one you select depends on

• Purpose of the data
• Quality requirements
• Available time