index.html

信号处理系列导航

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 <../private/inverse_tomography/exo1.m 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) '.']));


%%
% 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(sum(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 <../private/inverse_tomography/exo2.m 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);

##### SOURCE END #####
-->
   </body>
</html>