Outline Filters in Image Processing Image Restoration David Svoboda, Marek Kasfk email: svoboda@fi.muni.cz Centre for Biomedical Image Analysis Faculty of Informatics, Masaryk University, Brno, CZ CBIA December 6, 2019 David Svoboda, Marek Kasfk (CBIASFI) Filters in Image Processing autumn 2019 1/42 Motivation Images in optical microscopy are affected by blur and by noise. This blur is almost visible in z-axis. Images also suffer from noise due to low light intensities in confocal imaging. Cell at the start of optical path Cell after acquisition David Svoboda, Marek Kasfk (CBIASFI) Filters in Image Processing autumn 2019 4 / 42 Motivation » Blur • Noise ^ Restoration » Fundamentals a Unconstrained o Constrained • Iterative David Svoboda, Marek Kasfk (CBIASFI) Filters in Image Processing autumn 2019 2 / 42 Motivation The task of image restoration = transforming the acquired image to its original form. Acquired cell Restored cell image David Svoboda, Marek Kasfk (CBIASFI) Filters in Image Processing autumn 2019 5 / 42 Motivation Image Blur The blur can be described by Point Spread Function (PSF) response of optical setup to an infinitely small point source to the input. All the points are influenced by this function. PSF is of light placed Blur can be caused by different sources: © Move of the camera during acquisition 0 Defoe us (D Physical limits David Svoboda, Marek Kasfk (CBIAOFI) Filters in linage Processing autumn 2019 6 / 42 Motivation Noise Noise is present in almost every real image. It can be caused by, for example: o Environmental conditions during image acquisition (temperature) • Quality of the sensing elements (hot pixels) • Interference during image/data transmission Noise is "a random change" of pixel values. Our interest is usually focused on the three basic types of noise: • Additive noise (amplifiers) o Impulse noise (hot/cold pixels in CCD) o Poisson noise (photon-shot noise) David Svoboda, Marek Kasfk (CBIASFI) Filters in Image Processing autumn 2019 8 / 42 Motivation Image Blur Examples line - move disc - defocused camera Airy disc - physical limit David Svoboda, Marek Kasfk (CBIASFI) Filters in Image Processing autumn 2019 7 / 42 Motivation Noise Examples Additive noise Impulse noise David Svoboda, Marek Kasfk (CBIASFI) Filters in Image Processing autumn 2019 9 / 42 Motivation Noise Additive noise The most common type of noise. Gray values and noise are independent: g = f + n where f is the original image, n is the noise, and g is the noisy image Noise n may have different distributions: • Gaussian distribution (amplifiers) • Rayleigh distribution (radar) • Exponential distribution (laser imaging) • Gamma distribution (laser imaging) David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 10 / 42 Restoration Fundamentals Digital image restoration tries to restore original image from acquired image WITH the knowledge of characteristics of degradations. Examples of restorations: » intensity correction • chromatic aberration correction o deconvolution David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 13 / 42 Motivation Noise Poisson Noise o Important type of noise in CCD imaging (photon-shot noise, thermal noise) • Poisson noise is not additive and depends on the signal. • The noisy image f of the original one g is given by random Poisson process, which describes photon collection for each pixel position o Let photons occur at CCD pixel (/,_/) with the average rate g,j (photons/s) and let the original image is observed for t seconds. Then the probability that exactly X photons have occurred at pixel is given by Axe"A P(X) X! where A = g,j ■ t and X = 0,1, 2, David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 11 / 42 Deconvolution Deconvolution is an inverse process to convolution. It tries to remove blur from image. This inverse process has to deal with noise as well. Classification of deconvolution methods: =>- According to linearity of image processing: • Linear - linear filtering is performed • Non-linear - non-linear filtering is performed =>- According to knowledge of PSF: • Blind Deconvolution - PSF is unknown o Non-blind Deconvolution - PSF is known David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 14 / 42 Deconvolution Image degradation: image image RESTORATION estimate Ulb 1 UK 1 IUIM f f g i noise n Deconvolution tries to invert degradation of an image. Such process is possible only in some cases. Sometimes, the image is irrecoverably damaged and we cannot restore most of details. Notice: Blur removes some frequencies from the image. In this case, we cannot restore these frequencies. David Svoboda, Marek Kašík (CBIAOFI) Filters in linage Processing autumn 2019 15 / 42 Deconvolution Example Point-spread function used in the previous slide: 1 • 1 M 1 1 ■ Point-spread function of confocal microscope David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 17 / 42 Deconvolution Example An example of blurred image may be an image acquired in confocal microscopy. Picture is blurred with the point-spread function (PSF) of microscope. Cell at the start of optical path Cell after acquisition (blur + noise) David Svoboda, Marek Kašík (CBIAOFI) Filters in Image Processing Acquired image after deconvolution autumn 2019 16 / 42 Convolution x Matrix multiplication Block circulant matrix 1 2 3 4 -1 1 -2 2 A* B -1 2 -3 8 -2 8 BfAp 1 2 -1 1 Aext = 3 4 Bext = -2 2 0 0 0 0 -i 0 1 0 0 0 -2 0 2 i -1 0 0 0 0 2 -2 0 0 1 -1 0 0 0 0 2 -2 -2 0 2 -1 0 1 0 0 0 2 -2 0 1 -1 0 0 0 0 0 2 -2 0 1 - 1 0 0 0 0 0 0 -2 0 2 -1 0 1 0 0 0 2 -2 0 1 -1 0 0 0 0 0 2 - 2 0 1 -1 - - - 1 ■ - -1 - 2 -1 0 2 3 -5 4 = -3 0 8 0 -6 0 -2 . 0 . . 8 . o Cp, Ap ... linearized versions of Cext and Aext 9 Matrix Bt, . . . block circulant version of Bext David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 18 / 42 Convolution x Matrix multiplication Block circulant matrix Any convolution g = h * f can be written as a matrix multiplication. The above equation can be written as g=Hf where g and f are understood as linearized matrices, i.e. vectors • H is block circulant version of matrix h Notice: The math background and complexity is the same as in the case of the convolution. It is only a notation. David Svoboda, Marek Kašík (CBIAOFI) Filters in linage Processing autumn 2019 19 / 42 Unconstrained Restoration Setting the derivative of 1/1/(0 to zero with respect to f produces: <9l/l/(0 df and solving for f yields: -2HT(g-Hf) = 0 f = H-xg which can be written (in the case of space-invariant PSF by using of convolution theorem): f = FT -i ( FT(g) ] I FT(h) j Disadvantages: 9 Problems with zero-amplitude frequencies. • The approach doesn't deal with noise, o Inverse matrix may not exists. David Svoboda, Marek Kasfk (CBIAOFI) Filters in Image Processing autumn 2019 21 / 42 Unconstrained Restoration Unconstrained restoration is a base method for image deconvolution. It is very fast and applicable to data without noise. This method assumes following image formation process: g = Hf g . . . acquired image (including degradation) H... point spread function f ... ideal image we are searching for f ... our estimate of f We seek to minimize the function: 1/1/(0 = ||e(0||2 = \\g - Hf\\2 = (g- Hf)T(g - Hf) where ||a|| = Va1a is the Euclidean norm of vector and e(0 = g — Hf is a vector of residuals. David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 20 / 42 Wiener Filtering Cross-correlation (revision) Cross-correlation Cross-correlation function Rfg indicates the relative degree to which two functions f and g agree for various amounts of misalignment. It is given by Rfg(r)= J f(t)g(t + r)dt Auto-correlation Auto-correlation is special case of cross-correlation: David Svoboda, Marek Kašík (CBIASFI) OO Rf(r)= J f(t)f(t + r)dt Filters in Image Processing autumn 2019 22 / 42 Wiener Filtering Fourier analysis (revision) Power spectrum Power spectrum V of signal f is the FT of autocorrelation of f. FT[Rf(r)} = FT[f(t) * f(-t)} = F{k)F{-k) = F(k)T*(k) = \F(k)\2 = Vf(k) Lena Power spectra of Lena David Svoboda, Marek Kašík (CBIAOFI) Filters in Image Processing autumn 2019 23 / 42 Wiener Filtering Derivation MSE(W) (f(t) - at) dt -OO oo f2(t)-2f(t)f(t) + f2(t) dt = Ti + T2 + T3 f(t)2dt = Rf(0) where Tl = T2 = 73 = David Svoboda, Marek Kašík (CBIASFI) oo oo oo 2 J f(t)f(t)dt = -2 J f(t)[(W * s)(t)]dt =-2 J W(t)Rfs(t)dt — oo —oo —oo 30 oo oo oc 2{t)dt= J (W * s){t){W * s)(t)dt = j J W(t)W{T)R5(T - t)drdt — oo —oo Filters in Image Processing autumn 2019 25 / 42 Wiener Filtering Derivation We are searching for filter W that suppresses noise and keeps the signal quality: f=W*(f + n)=W*s o f is noise-free input o 1/1/ is unknown filter • n is additive noise o f is signal s = f + n filtered with W We try to minimize: MSE(W) = j e2(t)dt= J (7(t) - f(t))2 dt David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 24 / 42 Wiener Filtering Derivation oo oo oc MSE{W) = /?f(0)-2 J W{t)Rf5{t)dt+ J J W{t)W{T)R5{r - t)drdt We try to minimize MSE(W), i.e. dMSEjW) dW Solution Vf Rf5 = W*R5 —>• FT —^ Vfs = FT(W)-Vs -> FT(W) = David Svoboda, Marek Kašík (CBIAOFI) Filters in Image Processing autumn 2019 26 / 42 Wiener Filtering Provided noise and signal are NOT correlated the transfer function of filter 'W minimizing MSE is defined: FT{W) where: • Pf is power spectrum of noise-free signal » Pn is power spectrum of noise Notice: 'W is called Wiener filter. Vfs Vf mm2 Vs Vf + Vn FUO\2 + FT(")|2 David Svoboda, Marek Kašík (CBIAOFI) Filters in linage Processing autumn 2019 27 / 42 Wiener Deconvolution An example Lena to Deconvolved Lena David Svoboda, Marek Kašík (CBIAOFI) Filters in Image Processing autumn 2019 29 / 42 Wiener Deconvolution Unconstrained restoration followed by Wiener filtering we get optimal deconvolution respecting the MSE condition: f = FT -i (FT(g) Vf \FT(h) Vf + Vk w here Vk FT{n) FT(h) is a cross power spectrum of noise and PSF. There exists more practical version of Wiener deconvolution: -i (FT{g)Vg-Vn f = FJ- \FT(h) Vg David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 28 / 42 Constrained Least Square Restoration Idea To suppress noise in the result of deconvolution we should impose some constraints to solution. Such constraint can be the convolution of result with some convolution kernel and the minimization of power of this image. The task is to minimize: W(f) = \\Qf\\2 + X(\\g - Hf\\2 - \\n\\2) g-Hf = n where Q .. . convolution kernel, which imposes some constraint on resulting image f • A ... Lagrange multiplier (magic factor) Example: Concerning Q, we can use Laplace filter which boosts high frequencies (noise). David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 30 / 42 Constrained Restoration Now we need to solve equation <9l/l/(0 df 2QT Qf — 2\HT{g — Hf) = 0 Solving for f then yields: f = (\HTH + QTQ)-1XHTg David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 31 / 42 Iterative Methods In the most cases it is not possible to find solution of problem of restoration by simple linear filtering. Hence there exist different methods for restoration, which are iterative. This methods iteratively improve initial estimate until it is acceptable. Most common iterative methods: • EM-MLE (Expectation Maximization - Maximum Likelihood Estimation) • ICTM (Iterative Constrained Tikhonov-Miller algorithm) Used criteria: • number of iterations 9 relative change between two iterations • difference between blurred estimate and input David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 33 / 42 Constrained Least Square Deconvolution An example Lena David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing Deconvolved Lena autumn 2019 32 / 42 EM-MLE EM-MLE is one of the best methods for restoration of images degraded by noise with Poisson statistics. Its principle consists in maximizing certain functional. Such maximization is performed by Expectation-Maximization method, which is iterative numerical method. The functional which is maximized is likelihood functional and is expressed by: l(0 = UpM = U J J [li(uj)]NJ e-v(uj) where o p(x) is a probability density function following Poisson distribution • Nj is j-th point value of acquired image » fi(uj) is blurred version of our estimate (Hf). Notice: EM-MLE method is also known as Richardson-Lucy method. David Svoboda, Marek Kašík (CBIAOFI) Filters in linage Processing autumn 2019 34 / 42 EM-MLE Illustration L(í)=P(x0)P(x1)P(x2)P(x3)P(x4) 1 g(x) Hř(x) CL 4 - < I 3 3 cl cl 7 i N Q- 4< Xo Xi X2 X3 X4 X Xrj Xi X; X3 X4 X On the left is acquired image with Poisson distributions on values of g(x). On the right, you can see the same distributions, but used probabilities values (green spots) are from positions of (H/")(x) (green lines). Product of this probabilities is likelihood value. David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 35 / 42 EM-MLE Example Acquired cell Restored cell image David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 37 / 42 EM-MLE After derivation of this method we get simple expression for one iteration: fk+l = fkHT g Hfk where b is the known background • Hfk and HT() are matrix-vector multiplications • all other operations are point-wise David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 36 / 42 ICTM ICTM method considers additive noise with Gaussian distribution. Due to iterative form, we can impose non-linear constraints to result of deconvolution. One of the base constraints is non-negativity of resulting intensities. This cannot be achieved by non-iterative method. In this method we seek to minimize functional of the form: mf) = ^\\Hf-g\\2+7\\Qf\\2)- Notice: The functional of this form is solvable by algorithm of conjugate gradients (see numerical methods - FI:PV027 Optimization). David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 38 / 42 ICTM Algorithm Computation consists in the following steps: = P(fk + akdk), where dk is direction of k — th step, ak is size of the k — th step, fk is result of the last iteration and P() is projection operator, which clips values to non-negative values. Direction dk+1 is computed in this way: jk+l |f 11 :dk-rk. rk-l\\2 where rk = Afk-b,A = HTH + jQTQ and b = HTg. Termination of algorithm: threshold > W i+i W W+1 where 1/1/' is value of l/l/() at iteration /. David Svoboda, Marek Kasfk (CBIAOFI) Filters in linage Processing autumn 2019 39 / 42 Bibliography o Gonzalez, R. C, Woods, R. E., Digital image processing / 2nd ed. Upper Saddle River: Prentice Hall, c2002, pages 793, ISBN 0201180758 • Castleman, R. C, Digital Image Processing, Prentice Hall, 1996 9 Verveer, P. J., Computational and Optical Methods for Improving Resolution and Signal Quality in Fluorescence Microscopy, PhD thesis, Technical University Delft, 1998 David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 41 / 42 ICTM An example Acquired cell Restored cell image David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 40 / 42 You should know the answers ... • Explain the difference between a constrained and unconstrained image restoration. • Explain the equation g = H * f + n. • Explain why we cannot invert the convolution theorem to eliminate the consequences of convolution with known PSF. « Explain the meaning of the multiplication (F|) in EM-MLE algorithm. 9 What is a difference between Wiener filtering and Wiener deconvolution? o How do we get Vn in order we could perform the Wiener deconvolution? 9 Are you able to implement Constrained least square restoration in your favourite programming language? How do we stop iterative restoration methods? David Svoboda, Marek Kašík (CBIASFI) Filters in Image Processing autumn 2019 42 / 42