Reconstruction from Compressive Fourier Measurements
This numerical tour explores the reconstruction from randomized Fourier measurements with orthogonal and translation invariant sparse wavelet regularization.
Contents
Installing toolboxes and setting up the path.
You need to download the general purpose toolbox and the signal toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_general/ and toolbox_signal/ in your directory.
For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.
Recommandation: You should create a text file named for instance numericaltour.sce (in Scilabe) or numericaltour.m to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.
Execute this line only if you are using Matlab.
getd = @(p)path(path,p); % scilab users must *not* execute this
Then you can add these toolboxes to the path.
% Add some directories to the path getd('toolbox_signal/'); getd('toolbox_general/');
Partial Randomized Fourier Measurements
Compressive sampling / Compressed sensing corresponds to the random acquisition of linear measurement about a signal or an image. We explore a specific instantiation of this method, by taking randomized Fourier measurements. This provides a better incoherence between the measurements and the signal than tomography, which sub-samples the signal along lines in the Fourier domain.
We load an image.
n = 256;
M = load_image('cameraman', n);
M = rescale(M);
We compute a random mask that select a subset of frequencies. We enforce a few low frequency measurements that are not random.
% ratio of measurements rho = .2; % total number of measures Q = round(rho*n^2); % randomized values A = rand(n,n); [Y,X] = meshgrid(-n/2+1:n/2, -n/2+1:n/2); A = A + (sqrt(X.^2+Y.^2)<20); [v,I] = sort(A(:)); if v(1)<v(n*n) I = reverse(I); end % mask mask = zeros(n,n); mask(I(1:Q)) = 1;
We defines centered, normalized, Fourier transforms. With Matlab, we use inlined callback function, and with Scilab, the functions are implemented in a separated file.
if using_matlab() Fourier = @(M)fftshift(fft2(M)/n); iFourier = @(F)real( ifft2( fftshift(F)*n ) ); end
Build the compressed sensing acquisition function, CS(M) computes the CS measurements from M.
if using_matlab() S.type = '()'; S.subs = {find(mask==1)}; CS = @(M)subsref(Fourier(M),S); else global Ics; Ics = find(mask==1); end
The dual operator maps the measures to an image, by zero padding.
if using_matlab() iCS = @(y)iFourier(subsasgn(zeros(n,n),S,y)); end
Set the noise level that perturbates the measures.
sigma = .01;
Partial measurements with noise added.
y = CS(M) + sigma*randn(Q,1);
Pseudo inverse reconstruction. Since the measurements is an orthogonal projection, its pseudo inverse is simply the dual operator.
% Mpseudo = real( ifft2( fftshift(Mpseudo)*n ) );
Mpseudo = iCS(y);
We display the mask and the pseudo inverse reconstruction.
clf; imageplot(mask, 'Fourier mask', 1,2,1); imageplot(Mpseudo, 'Pseudo inverse reconstruction', 1,2,2);
Wavelet Reconstruction from Partial Tomoraphic Measurements
The compressed sensing recovery can be performed using L1 minimization in a wavelet basis. We start here with an orthogonal basis, which is fast, but gives poor denoising results.
This corresponds to the following optimization on the coefficients MW:
min_{MW} norm(y - CS o Wavelet(MW))^2 + lambda*norm(MW, 1) (*)
where Wavelet(.) is the wavelet reconstruction operator. Note that MW has a much larger number of entries than the original image M. The recovered image M1 is obtained as M1 = Wavelet(MW).
This L1 optimization can be solved with iterative thresholding, also known as Foward Backward splitting. This is a generalization of the projected gradient descent where the projection is replaced by a soft thresholding. Note that faster algorithms exists, but iterative thresholding are very simple to implement.
We set the regularization parameter lambda. The larger, the more agressive the regularization (denoising effect) is.
lambda = .05;
Parameters for the orthogonal wavelet transform.
options.ti = 0; Jmin = 4;
We extend the stack size (for Scilab users).
extend_stack_size(4);
Initialization of the coefficients.
MW = zeros(n,n);
Step 1: Gradient descent step:
MW <- MW + Phi'*( y - Phi*MW)
Where Phi is the operator CS o Wavelet, and thus Phi' = iWavelet o iCS, where iWavelet is the inverse wavelet transform.
Compute the gradient descent. Since Wavelet(MW) is goind to be usefull as it is the restored image, we store it appart in a variable M1.
M1 = perform_wavelet_transf(MW, Jmin, -1, options); MW = MW + perform_wavelet_transf( iCS(y-CS(M1)), Jmin, +1, options);
Step 2: Soft thresholding step:
Here we do not threshold the low frequency residual because it gives better denoising result.
a = MW(1:2:2^Jmin,1:2:2^Jmin);
MW = perform_thresholding(MW, lambda, 'soft');
MW(1:2:2^Jmin,1:2:2^Jmin) = a;
Exercice 1: (the solution is exo1.m) Implement the iterative thresholding algorithm to minimize the wavelet sparsity. At each step, modify the lambda parameter to fit the noise level. Display the evolution the energy that is minimized by the algorithm.
exo1;
Display the result.
clf; imageplot(clamp(Mpseudo), strcat(['Pseudo inverse, SNR=' num2str(snr(M,Mpseudo),3) 'dB']), 1,2,1); imageplot(clamp(M1), strcat(['Wavelets, SNR=' num2str(snr(M,M1),3) 'dB']), 1,2,2);
Update of the lambda Parameter
The value of lambda is difficult to choose. One can instead update this parameter so that the error norm(y-y1) fit the noise level sqrt(Q)*sigma.
After step 1 and 2, we can update the lambda parameter to fit the noise.
lambda = lambda * sqrt(Q)*sigma / norm( y-y1 );
Intuitively, if the error norm(y-y1) is too large, then lambda should be decayed.
Exercice 2: (the solution is exo2.m) Implement the iterative thresholding algorithm with the update of the lambda parameter. Here we update the value of lambda only each 3 iteration to stabilize the adaptation process. Display the evolution of lambda, and the evolution of the eror norm(y-y1)/(sqrt(Q)*sigma).
exo2;
Display the result.
clf; imageplot(clamp(Mpseudo), strcat(['Pseudo inverse, SNR=' num2str(snr(M,Mpseudo),3) 'dB']), 1,2,1); imageplot(clamp(M1), strcat(['Wavelets, SNR=' num2str(snr(M,M1),3) 'dB']), 1,2,2);
Translation Invariant Wavelet Regularization
To obtain better result, we optimize the wavelet coefficient of the image in a translation invariant (i.e. redundant) wavelet tight frame.
Set up the translation invariant transform.
options.ti = 1;
Exercice 3: (the solution is exo3.m) Implement the same iterative thresholding algorithm with the update of the lambda parameter. Store the result in an image M2
exo3;
Display the result.
clf; imageplot(clamp(M1), strcat(['Orthogonal, SNR=' num2str(snr(M,M1),3) 'dB']), 1,2,1); imageplot(clamp(M2), strcat(['Invariant, SNR=' num2str(snr(M,M2),3) 'dB']), 1,2,2);
