Laboratorio di Problemi Inversi Esercitazione 2...

transcript

Laboratorio di Problemi InversiEsercitazione 2: filtraggio spettrale

Luca Calatroni

Dipartimento di Matematica, Universita degli studi di Genova

Aprile 2016.

Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 1 / 17

Outline

1 Exploiting PSF structure for BCCB matrix computation

2 TSVD and Tikhonov implementation

3 Tikhonov parameter choice via GCV

A structural way to build blur operators

I Load the Shepp-Logan phantom

im1=phantom(’Modified Shepp-Logan’, 256);

I Generate Gaussian PSF of the same size, σ2 = 3 to blur the image.

We want to exploit the structure of the PSF and the type of boundary conditionsconsidered to build efficiently the blur matrix B (similar to how imfilter works).

Easy case: generic PSF, periodic boundary conditions → BCCB matrix: startingfrom PSF, elements are circularly shifted with respect to the centre of thePSF. . . see theory!

I Determine the centre of Gaussian PSF.

I Use circshift to build first column of the blur matrix B (still a matrix).

Why all this? Can’t we use MATLAB imfilter?

Need for inversion!No blur operator provided by MATLAB when using fspecial, only PSF andfiltered results! From an inverse problem perspective, the operator B (needed forinversion) is not provided by MATLAB.

The spectral properties of B can further be exploited for quick implementation!

. . . and actually the computation of B is never directly required!

but. . . “only the first column is needed”!

Notation

For what follows: denote by

I σi , i = 1, . . . ,N the (complex) eigenvalues of B → diagonal matrix Σ;

I F the 2-dimensional Fourier transform matrix;

I F∗ the 2-dimensional anti-Fourier transform matrix.

MATLAB: use fft2 and ifft2 as (rescaled) matrix multiplication by F and F∗,respectively.

. . . Vectorisation! y=Bx. . .

Warning: fft2 and ifft2 act on 2D arrays. So, from here onwards

1 First apply the functions to matrices;

2 Then, reshape.

Useful properties of BBCB matrices

1 The set of eigenvectors of every BCCB matrix can be expressed in terms ofthe 2D-Fourier transform:

B = F∗ΣF ,

which entails:FB = ΣF → Fb1 = Σf1

where b1 and f1 are the first columns of B and F matrix, respectively. Wehave:

f1 =1√N

11...1

2 ⇒ to compute σ, the eigenvalues of B, apply fft2 to first column of blurmatrix (which comes, unintuitively, as a matrix), then vectorise the result.

Useful properties of BBCB matrices

1 The set of eigenvectors of every BCCB matrix can be expressed in terms ofthe 2D-Fourier transform:

B = F∗ΣF ,

which entails:

FB = ΣF → Fb1 = Σf1=1√Nσ

where b1 and f1 are the first columns of B and F matrix, respectively. Wehave:

f1 =1√N

11...1

2 ⇒ to compute σ, the eigenvalues of B, apply fft2 to first column of blur

matrix (which comes, unintuitively, as a matrix), then vectorise the result.

Why all this?

We will. . .. . . create MATLAB functions that will serve as regularised solvers for thedeblurring problem. Input? We will focus on Tikhonov and TSVD, bothimplemented by means of the σi .

Do it wisely! Using this strategy based on the use of fft2 reducescomputational costs with respect to using the eig function.

I Step 0: use all this to create the function naive deblur to computedeblurred version of im1 using:

x = B−1y = F∗Σ−1Fy

- compute eigenvalues as above.- standard Fourier inversion (take the real part).- we know it won’t work!

Why all this?

x = B−1y = F∗Σ−1Fy

Why all this?

x = B−1y = F∗Σ−1Fy

Outline

Let’s take a look at the SVD. . .

Let B bw BCCB and let us introduce Φ, filters (like the ones we saw last time forthresholding/regularisation), we have:

x = B−1y = F∗ΦΣ−1Fy

Towards filtering. . .

I refine naive deb giving Φ as input. In the naive case, what is Φ?

I vectorise Σ and blurred image as σ, y , respectively.

I introduce auxiliary filtering vector s filt = . . . and “invert by multiplying”.

I remember always to reshape before applying fft2, ifft2 and to take thereal part.

I Output: deblurred image.

Let’s take a look at the SVD. . .

Let B bw BCCB and let us introduce Φ, filters (like the ones we saw last time forthresholding/regularisation), we have:

Towards filtering. . .

I refine naive deb giving Φ as input. In the naive case, what is Φ?

I vectorise Σ and blurred image as σ, y , respectively.

I introduce auxiliary filtering vector s filt = . . . and “invert by multiplying”.

I remember always to reshape before applying fft2, ifft2 and to take thereal part.

I Output: deblurred image.

The problem

Informally, when inverting with MATLAB:

Phi=speye(size(im));

phi=diag(Phi);

deb=ifft2((phi./s).* "fft2(B)")

Note the (possible) division by 0!

Truncated Singular Value Decomposition filtering

Idea: removing small singular values of σ by filtering them by multiplication by 0.

{1, i = 1, . . . , k ,

0, i = k + 1, . . . ,N.

Task: modify Φ to eliminate “small” singular values.

I Write function tsvd fft2 with the same input as before with additionaltol=0.01 value to threshold.

I Set small elements of σ to 0.

I Use find to compute the vectors of the non-null elements of Φ.

I For these elements, compute the filtering, reshape, use ifft2 to compute theoutput deblurred image.

{1, i = 1, . . . , k ,

0, i = k + 1, . . . ,N.

{1, i = 1, . . . , k ,

0, i = k + 1, . . . ,N.

TSVD: take-home message and variations on the theme

Problem solved!!Observe that the problems of the inversion are solved by the multiplication offilters which never blow up (after thresholding)!

1 1st variation: verify the robustness of the approach after adding noise?

2 2nd variation: how does the choice of tol affect the deblurring result?

TSVD: take-home message and variations on the theme

Problem solved!!Observe that the problems of the inversion are solved by the multiplication offilters which never blow up (after thresholding)!

1 1st variation: verify the robustness of the approach after adding noise?

2 2nd variation: how does the choice of tol affect the deblurring result?

Tikhonov filtering

Idea: “correct” Φ by imposing that the denominator never vanishes.

φi =σ2i

σ2i + λ

, λ > 0,

where λ is the Tikhonov regularisation parameter.

Task: write Tikhonov filtering function.

I Write function tik fft2 with the same input as before with additionalλ=0.01 regularisation parameter..

I Write the filters appropriately. . . careful: do not divide by zero!

I Compute the filtering, reshape, use ifft2 to compute the output deblurredimage.

+ Variations: noise? Choice of λ?

Tikhonov filtering

φi =σ2i

σ2i + λ

, λ > 0,

Tikhonov filtering

φi =σ2i

σ2i + λ

, λ > 0,

Tikhonov filtering

φi =σ2i

σ2i + λ

, λ > 0,

The choice of λ

Tikhonov functional:

min‖y − Bx‖22 + λ‖x‖2

Heuristically:

I small λ: low regularisation, high data-fitting. Reconstruction is still noisysince oscillations are preserved.

I large λ: high regularisation, low data-fitting. Reconstruction isover-smoothed.

How to choose the optimal λ?

The optimal λ is somewhere in between. Is there an automatic method to chooseit? Which prior information do we need?

I prior information on the noise level σ2: discrepancy principle.

I no prior knowledge on the noise level: Generalised Cross Validation(GCV).

The choice of λ

Tikhonov functional:

min‖y − Bx‖22 + λ‖x‖2

Heuristically:

I small λ: low regularisation, high data-fitting. Reconstruction is still noisysince oscillations are preserved.

I large λ: high regularisation, low data-fitting. Reconstruction isover-smoothed.

How to choose the optimal λ?

The optimal λ is somewhere in between. Is there an automatic method to chooseit? Which prior information do we need?

I prior information on the noise level σ2: discrepancy principle.

I no prior knowledge on the noise level: Generalised Cross Validation(GCV).

Outline

GCV for Tikhonov filtering

GCV functional for general filter matrix Φ:

G (λ) =‖(Id −Φ)y‖2

(tr(Id −Φ))2,

so its form depends on the filter used. Note that this does not make sense whenΦ = Id (non-filtered case).

GCV functional for Tikhonov filter matrix Φ:

G (λ) =

∑Ni=1

σ2i +λ

(∑Ni=1

1σ2i +λ

It can be shown that the minimum of the functional G is unique and itcorresponds to the optimal choice of λ (leave-one-out principe).

GCV functional for Tikhonov filter matrix Φ:

G (λ) =

∑Ni=1

σ2i +λ

(∑Ni=1

1σ2i +λ

It can be shown that the minimum of the functional G is unique and itcorresponds to the optimal choice of λ (leave-one-out principe).

Task: implement GCV function for Tikohonov filter. In particular:

I Input: PSF, blurred image. Output: optimal λ.

I Write the GCV functional as a @ GCV (function handle). Input?

I → use fminbound MATLAB inbuilt function to get the minimum of GCVfunctional. Use box constraints to bound the solution:

λ=fminbnd(@ GCV, min(abs(s)), max(abs(s)), [], s, bhat);

I Show result with optimal λ.

References

Inspiration and references are taken from:

- P. C. Hansen, J. G. Nagy, D. P. O’Leary, Deblurring Images, Matrices,Spectra and Filtering, SIAM, Philadelphia, 2006.

- Online material and inspiring MATLAB toolboxes:http://www.imm.dtu.dk/~pcha/HNO/

- Slides on blur matrices:http:

//www.mathcs.emory.edu/~nagy/courses/fall06/ID_lecture1.pdf

Laboratorio di Problemi Inversi Esercitazione 2...

Documents