?? index.html
字號:
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"><html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. --> <title>Reconstruction from Partial Tomography Measurements</title> <meta name="generator" content="MATLAB 7.4"> <meta name="date" content="2008-10-17"> <meta name="m-file" content="index"> <LINK REL="stylesheet" HREF="style.css" TYPE="text/css"> </head> <body> <div class="content"> <h1>Reconstruction from Partial Tomography Measurements</h1> <introduction> <p>This numerical tour explores the reconstruction from tomographic measurement with TV regularization.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#1">Installing toolboxes and setting up the path.</a></li> <li><a href="#8">Tomography and Radon Transform</a></li> <li><a href="#15">Tomography Acquisition over the Fourier Domain</a></li> <li><a href="#21">TV Reconstruction from Partial Tomoraphic Measurements</a></li> </ul> </div> <h2>Installing toolboxes and setting up the path.<a name="1"></a></h2> <p>You need to download the <a href="../toolbox_general.zip">general purpose toolbox</a> and the <a href="../toolbox_signal.zip">signal toolbox</a>. </p> <p>You need to unzip these toolboxes in your working directory, so that you have <tt>toolbox_general/</tt> and <tt>toolbox_signal/</tt> in your directory. </p> <p><b>For Scilab user:</b> you must replace the Matlab comment '%' by its Scilab counterpart '//'. </p> <p><b>Recommandation:</b> You should create a text file named for instance <tt>numericaltour.sce</tt> (in Scilabe) or <tt>numericaltour.m</tt> to write all the Scilab/Matlab command you want to execute. Then, simply run <tt>exec('numericaltour.sce');</tt> (in Scilab) or <tt>numericaltour;</tt> (in Matlab) to run the commands. </p> <p>Execute this line only if you are using Matlab.</p><pre class="codeinput">getd = @(p)path(path,p); <span class="comment">% scilab users must *not* execute this</span></pre><p>Then you can add these toolboxes to the path.</p><pre class="codeinput"><span class="comment">% Add some directories to the path</span>getd(<span class="string">'toolbox_signal/'</span>);getd(<span class="string">'toolbox_general/'</span>);</pre><h2>Tomography and Radon Transform<a name="8"></a></h2> <p>The Radon transform compute the integral of the image along all the lines in the plane.</p> <p>Load the phantom image</p><pre class="codeinput">n = 256;M = load_image(<span class="string">'phantom'</span>,n);M = rescale(M,.08,.92);</pre><p>In practice, only a limited number of orientation and intecept are retained.</p><pre class="codeinput"><span class="comment">% number of orientation</span>ntheta = round(sqrt(2)*n);<span class="comment">% orientations in [0,pi]</span>theta = linspace(0,180,ntheta+1); theta(ntheta+1) = [];</pre><p>Now we compute the radon transform. Each <tt>R(i,j)</tt> measure the integral along a ray of angle <tt>theta(j)</tt>, with an intercept parameterized by <tt>i</tt>. </p><pre class="codeinput">R = radon(M,theta);<span class="comment">% display</span>clf;imageplot(M, <span class="string">'Image'</span>, 1,2,1);imageplot(R, <span class="string">'Radon transform'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_01.png"> <p>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. </p><pre class="codeinput"><span class="comment">% number of projections</span>nproj = 32;<span class="comment">% selected angles for inversion</span>sel = round(linspace(1,ntheta+1,nproj+1)); sel(nproj+1) = [];<span class="comment">% inversion with back projection</span>M1 = iradon(R(:,sel), theta(sel));<span class="comment">% display</span>clf;imageplot(M, <span class="string">'Image'</span>, 1,2,1);imageplot(clamp(M1), <span class="string">'Reconstruction'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_02.png"> <p>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. </p><pre class="codeinput"><span class="comment">% Fourier transform, orthogonal</span>F = fftshift(fft2(M))/n;F1 = fftshift(fft2(M1))/n;<span class="comment">% display</span>eta = 1e-1;clf;imageplot(log(eta+abs(F)), <span class="string">'Fourier transform of the image'</span>, 1,2,1);imageplot(log(eta+abs(F1)), <span class="string">'Fourier transform of the reconstruction'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_03.png"> <p><i>Exercice 1:</i> (the solution is <a href="../private/inverse_tomography/exo1.m">exo1.m</a>) Perform a reconstruction with an increasing number of angles, and display the reconstruction error. </p><pre class="codeinput">exo1;</pre><img vspace="5" hspace="5" src="index_04.png"> <h2>Tomography Acquisition over the Fourier Domain<a name="15"></a></h2> <p>The tomography measurement with a Randon transform is equivalent, over the Fourier domain, to a sub-sampling along rays.</p> <p>Set the number of rays used for the experiments.</p><pre class="codeinput">nrays = 18;</pre><p>For a given number of rays, we build the mask that perform the sub-sampling.</p><pre class="codeinput">Theta = linspace(0,pi,nrays+1); Theta(end) = [];mask = zeros(n);<span class="keyword">for</span> 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;<span class="keyword">end</span></pre><p>We display the mask and the masked Fourier transform</p><pre class="codeinput"><span class="comment">% Fourier transform, orthogonal</span>F = fftshift(fft2(M)/n);<span class="comment">% display</span>clf;imageplot(mask, <span class="string">'Fourier mask'</span>, 1,2,1);imageplot( log( eta + abs(mask.*F) ), <span class="string">'Masked frequencies'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_05.png"> <p>Tomographic measurement can thus be intepreted as a selection of a few Fourier frequencies.</p><pre class="codeinput"><span class="comment">% selection operator M->y</span>F = fftshift(fft2(M)/n);y = F(mask==1);<span class="comment">% number of measures</span>Q = length(y);disp(strcat([<span class="string">'Number of measurements Q='</span> num2str(Q) <span class="string">'.'</span>]));disp(strcat([<span class="string">'Sub-sampling Q/N='</span> num2str(length(y)/n^2,2) <span class="string">'.'</span>]));</pre><pre class="codeoutput">Number of measurements Q=5359.Sub-sampling Q/N=0.082.</pre><p>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). </p><pre class="codeinput">F1 = zeros(n);F1(mask==1) = y;M1 = real( ifft2( fftshift(F1)*n ) );<span class="comment">% display</span>clf;imageplot(M, <span class="string">'Original image'</span>, 1,2,1);imageplot(clamp(M1), <span class="string">'Pseudo inverse reconstruction'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_06.png"> <h2>TV Reconstruction from Partial Tomoraphic Measurements<a name="21"></a></h2> <p>Since the Phantom image is a cartoon image, it makes sense to perform the reconstruction while minimizing the TV norm of the image. </p> <p>Here we deal with noisy measurements.</p><pre class="codeinput">F = fftshift(fft2(M)/n);<span class="comment">% noise level</span>sigma = .03;<span class="comment">% selection operator M->y</span>y = F(mask==1) + sigma*randn(Q,1);</pre><p>The reconstruction is performed using the following TV minimization</p> <p><tt>min_{Mtv} norm(Tomography(Mtv)-y)^2 + lambda*normTV(Mtv)</tt> (*) </p><pre> Where |Tomography(.)| is the Fourier tomographic sub-sampling operator. The parameter |lambda| is optimized so that the solution satisfies</pre><p><tt>norm(Tomography(Mtv)-y)<sqrt(Q)*sigma</tt> (**) </p><pre> where |Q| is the number of measurements, and |sigma| is the noise level per measure.</pre><p>The minimization of (*) is performed by iterating between a gradient descent step of the energy <tt>norm(Tomography(Mtv)-y)^2</tt>, and Chambolle's algorithm to minimize the TV norm. </p> <p>For more details about Chambolle's algorithm, see the numerical tour on TV minimization.</p> <p>Initialization</p><pre class="codeinput">Mtv = zeros(n);<span class="comment">% regularization parameter</span>lambda = .1;<span class="comment">% vector field for Chambolle's algorithm</span>G = zeros(n,n,2);</pre><p>The first step corresponds to a projection on the measurement constraints. This corresponds in fact to the pseudo inverse solution. </p><pre class="codeinput">Ftv = fftshift(fft2(Mtv)/n);Ftv(mask==1) = y;M1 = real(ifft2( fftshift(Ftv*n) ));Mpseudo = M1;</pre><p>The second step corresponds to several steps of Chambolle's algorithm to minimize the TV norm of <tt>M1</tt>, using a fixed <tt>lambda</tt>. </p><pre class="codeinput"><span class="comment">% descent step in Chambolle's method</span>tau = 1/4;<span class="comment">% number of sub-iterations for Chambolle's method</span>niter_chambolle = 20;<span class="comment">% tolerance for early stop</span>tol = 1e-5;<span class="keyword">for</span> i=1:niter_chambolle <span class="comment">% gradient of the energy</span> dG = grad( div(G) - M1/lambda ); <span class="comment">% gradient descent</span> G = G + tau*dG; <span class="comment">% projection on Linfty constraints</span> d = repmat( sqrt(sum(G.^2,3)), [1 1 2] ); <span class="comment">% norm of the vectors</span> G = G ./ max(d,1); <span class="comment">% reconstruct</span> MtvOLD = Mtv; Mtv = M1 - lambda*div(G); <span class="comment">% check for early stop if no improvement is being made</span> <span class="keyword">if</span> norm(Mtv-MtvOLD, <span class="string">'fro'</span>)/norm(M1,<span class="string">'fro'</span>)<tol <span class="keyword">break</span>; <span class="keyword">end</span><span class="keyword">end</span></pre><p>The last step is an update of the <tt>lambda</tt> parameter to match the error constraints. </p><pre class="codeinput"><span class="comment">% measures</span>Ftv = fftshift(fft2(Mtv)/n);y1 = Ftv(mask==1);<span class="comment">% update</span>lambda = lambda * sqrt(Q)*sigma / norm( y-y1, <span class="string">'fro'</span> );</pre><p><i>Exercice 2:</i> (the solution is <a href="../private/inverse_tomography/exo2.m">exo2.m</a>) Put these three steps together in an iterative algorithm that reconstruction from the tomographic measurements by solving the TV minimization. </p><pre class="codeinput">exo2;</pre><img vspace="5" hspace="5" src="index_07.png"> <p>Display</p><pre class="codeinput">clf;imageplot(clamp(Mpseudo), strcat([<span class="string">'Pseudo inverse, SNR='</span> num2str(snr(M,Mpseudo),3) <span class="string">'dB'</span>]), 1,2,1);imageplot(clamp(Mtv), strcat([<span class="string">'TV, SNR='</span> num2str(snr(M,Mtv),3) <span class="string">'dB'</span>]), 1,2,2);</pre><img vspace="5" hspace="5" src="index_08.png"> <p class="footer"><br> Copyright ® 2008 Gabriel Peyre<br></p> </div> <!--##### SOURCE BEGIN #####%% Reconstruction from Partial Tomography Measurements% This numerical tour explores the reconstruction from tomographic% measurement with TV regularization.%% Installing toolboxes and setting up the path.%%% You need to download the % <../toolbox_general.zip general purpose toolbox>% and the <../toolbox_signal.zip 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 pathgetd('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 imagen = 256;M = load_image('phantom',n);?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -