Reconstruction from Partial Tomography Measurements
This numerical tour explores the reconstruction from tomographic measurement with TV 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 Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab commands 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/');
Tomography and Radon Transform
The Radon transform compute the integral of the image along all the lines in the plane.
Load the phantom image
n = 256;
M = load_image('phantom',n);
M = rescale(M,.08,.92);
In practice, only a limited number of orientation and intecept are retained.
% number of orientation ntheta = round(sqrt(2)*n); % orientations in [0,pi] theta = linspace(0,180,ntheta+1); theta(ntheta+1) = [];
Now we compute the radon transform. Each R(i,j) measure the integral along a ray of angle theta(j), with an intercept parameterized by i.
R = radon(M,theta); % display clf; imageplot(M, 'Image', 1,2,1); imageplot(R, 'Radon transform', 1,2,2);
In practice, the reconstruction is performed using only a limitted number of projections. Since there is many possible inverse for these partial measurements, the reconstruction is performed with a pseudo inverse. This pseudo-inverse is implemented in a stable and fast maner with filtered back-projection.
% number of projections nproj = 32; % selected angles for inversion sel = round(linspace(1,ntheta+1,nproj+1)); sel(nproj+1) = []; % inversion with back projection M1 = iradon(R(:,sel), theta(sel)); % display clf; imageplot(M, 'Image', 1,2,1); imageplot(clamp(M1), 'Reconstruction', 1,2,2);
The pseudo inverse corresponds to imposing the value of the Fourier transform of the reconstruction only about 1D lines that are orthogonal to the angle used for the Radon transform.
% Fourier transform, orthogonal F = fftshift(fft2(M))/n; F1 = fftshift(fft2(M1))/n; % display eta = 1e-1; clf; imageplot(log(eta+abs(F)), 'Fourier transform of the image', 1,2,1); imageplot(log(eta+abs(F1)), 'Fourier transform of the reconstruction', 1,2,2);
Exercice 1: (the solution is exo1.m) Perform a reconstruction with an increasing number of angles, and display the reconstruction error.
exo1;
Tomography Acquisition over the Fourier Domain
The tomography measurement with a Randon transform is equivalent, over the Fourier domain, to a sub-sampling along rays.
Set the number of rays used for the experiments.
nrays = 18;
For a given number of rays, we build the mask that perform the sub-sampling.
Theta = linspace(0,pi,nrays+1); Theta(end) = []; mask = zeros(n); for theta = Theta t = linspace(-1,1,3*n)*n; x = round(t.*cos(theta)) + n/2+1; y = round(t.*sin(theta)) + n/2+1; I = find(x>0 & x<=n & y>0 & y<=n); x = x(I); y = y(I); mask(x+(y-1)*n) = 1; end
We display the mask and the masked Fourier transform
% Fourier transform, orthogonal F = fftshift(fft2(M)/n); % display clf; imageplot(mask, 'Fourier mask', 1,2,1); imageplot( log( eta + abs(mask.*F) ), 'Masked frequencies', 1,2,2);
Tomographic measurement can thus be intepreted as a selection of a few Fourier frequencies.
% selection operator M->y F = fftshift(fft2(M)/n); y = F(mask==1); % number of measures Q = length(y); disp(strcat(['Number of measurements Q=' num2str(Q) '.'])); disp(strcat(['Sub-sampling Q/N=' num2str(length(y)/n^2,2) '.']));
Number of measurements Q=5359. Sub-sampling Q/N=0.082.
The transposed operator corresponds to the pseudo inverse reconstruction (because the measurement operator is in fact an orthogonal projection). It is similar to the filtered back-projection (excepted that the Fourier sub-sampling is now on a discrete grid, which is not really faithful to the geometry of tomographic acquisition).
F1 = zeros(n); F1(mask==1) = y; M1 = real( ifft2( fftshift(F1)*n ) ); % display clf; imageplot(M, 'Original image', 1,2,1); imageplot(clamp(M1), 'Pseudo inverse reconstruction', 1,2,2);
TV Reconstruction from Partial Tomoraphic Measurements
Since the Phantom image is a cartoon image, it makes sense to perform the reconstruction while minimizing the TV norm of the image.
Here we deal with noisy measurements.
F = fftshift(fft2(M)/n); % noise level sigma = .03; % selection operator M->y y = F(mask==1) + sigma*randn(Q,1);
The reconstruction is performed using the following TV minimization
min_{Mtv} norm(Tomography(Mtv)-y)^2 + lambda*normTV(Mtv) (*)
Where |Tomography(.)| is the Fourier tomographic sub-sampling operator. The parameter |lambda| is optimized so that the solution satisfies
norm(Tomography(Mtv)-y)<sqrt(Q)*sigma (**)
where |Q| is the number of measurements, and |sigma| is the noise level per measure.
The minimization of (*) is performed by iterating between a gradient descent step of the energy norm(Tomography(Mtv)-y)^2, and Chambolle's algorithm to minimize the TV norm.
For more details about Chambolle's algorithm, see the numerical tour on TV minimization.
Initialization
Mtv = zeros(n); % regularization parameter lambda = .1; % vector field for Chambolle's algorithm G = zeros(n,n,2);
The first step corresponds to a projection on the measurement constraints. This corresponds in fact to the pseudo inverse solution.
Ftv = fftshift(fft2(Mtv)/n); Ftv(mask==1) = y; M1 = real(ifft2( fftshift(Ftv*n) )); Mpseudo = M1;
The second step corresponds to several steps of Chambolle's algorithm to minimize the TV norm of M1, using a fixed lambda.
% descent step in Chambolle's method tau = 1/4; % number of sub-iterations for Chambolle's method niter_chambolle = 20; % tolerance for early stop tol = 1e-5; for i=1:niter_chambolle % gradient of the energy dG = grad( div(G) - M1/lambda ); % gradient descent G = G + tau*dG; % projection on Linfty constraints d = repmat( sqrt(sum3(G.^2,3)), [1 1 2] ); % norm of the vectors G = G ./ max(d,1); % reconstruct MtvOLD = Mtv; Mtv = M1 - lambda*div(G); % check for early stop if no improvement is being made if norm(Mtv-MtvOLD, 'fro')/norm(M1,'fro')<tol break; end end
The last step is an update of the lambda parameter to match the error constraints.
% measures Ftv = fftshift(fft2(Mtv)/n); y1 = Ftv(mask==1); % update lambda = lambda * sqrt(Q)*sigma / norm( y-y1, 'fro' );
Exercice 2: (the solution is exo2.m) Put these three steps together in an iterative algorithm that reconstruction from the tomographic measurements by solving the TV minimization.
exo2;
Display
clf; imageplot(clamp(Mpseudo), strcat(['Pseudo inverse, SNR=' num2str(snr(M,Mpseudo),3) 'dB']), 1,2,1); imageplot(clamp(Mtv), strcat(['TV, SNR=' num2str(snr(M,Mtv),3) 'dB']), 1,2,2);
