Laboratorio di Problemi Inversi Esercitazione 2...
-
Upload
truongphuc -
Category
Documents
-
view
220 -
download
0
Transcript of Laboratorio di Problemi Inversi Esercitazione 2...
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
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 2 / 17
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).
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 3 / 17
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).
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 3 / 17
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).
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 3 / 17
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”!
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 4 / 17
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”!
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 4 / 17
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”!
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 4 / 17
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.
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 5 / 17
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.
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 6 / 17
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.
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 6 / 17
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!
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 7 / 17
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!
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 7 / 17
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!
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 7 / 17
Outline
1 Exploiting PSF structure for BCCB matrix computation
2 TSVD and Tikhonov implementation
3 Tikhonov parameter choice via GCV
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 8 / 17
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.
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 9 / 17
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.
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 9 / 17
The problem
x = B−1y = F∗ΦΣ−1Fy
Informally, when inverting with MATLAB:
Phi=speye(size(im));
phi=diag(Phi);
deb=ifft2((phi./s).* "fft2(B)")
Note the (possible) division by 0!
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 10 / 17
Truncated Singular Value Decomposition filtering
Idea: removing small singular values of σ by filtering them by multiplication by 0.
φi =
{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.
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 11 / 17
Truncated Singular Value Decomposition filtering
Idea: removing small singular values of σ by filtering them by multiplication by 0.
φi =
{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.
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 11 / 17
Truncated Singular Value Decomposition filtering
Idea: removing small singular values of σ by filtering them by multiplication by 0.
φi =
{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.
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 11 / 17
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?
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 12 / 17
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?
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 12 / 17
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 λ?
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 13 / 17
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 λ?
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 13 / 17
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 λ?
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 13 / 17
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 λ?
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 13 / 17
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).
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 14 / 17
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).
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 14 / 17
Outline
1 Exploiting PSF structure for BCCB matrix computation
2 TSVD and Tikhonov implementation
3 Tikhonov parameter choice via GCV
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 15 / 17
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).
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 16 / 17
GCV for Tikhonov filtering
GCV functional for Tikhonov filter matrix Φ:
G (λ) =
∑Ni=1
(y
σ2i +λ
)2
(∑Ni=1
1σ2i +λ
)2
It can be shown that the minimum of the functional G is unique and itcorresponds to the optimal choice of λ (leave-one-out principe).
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 16 / 17
GCV for Tikhonov filtering
GCV functional for Tikhonov filter matrix Φ:
G (λ) =
∑Ni=1
(y
σ2i +λ
)2
(∑Ni=1
1σ2i +λ
)2
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 λ.
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 16 / 17
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
Luca Calatroni (DIMA, Unige) Esercitazione 2, Lab. Prob. Inv. Aprile 2016 17 / 17