📄 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 Compressive Fourier 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 Compressive Fourier Measurements</h1> <introduction> <p>This numerical tour explores the reconstruction from randomized Fourier measurement with wavelet 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">Partial Randomized Fourier Measurements</a></li> <li><a href="#14">Wavelet 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>Partial Randomized Fourier Measurements<a name="8"></a></h2> <p>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. </p> <p>We load an image.</p><pre class="codeinput">n = 256;M = load_image(<span class="string">'cameraman'</span>, n);M = rescale(M);</pre><p>We compute a random mask that select a subset of frequencies. We enforce a few low frequency measurements that are not random.</p><pre class="codeinput"><span class="comment">% ratio of measurements</span>rho = .2;<span class="comment">% total number of measures</span>Q = round(rho*n^2);<span class="comment">% randomized values</span>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(:));<span class="keyword">if</span> v(1)<v(n*n) I = reverse(I);<span class="keyword">end</span><span class="comment">% mask</span>mask = zeros(n,n);mask(I(1:Q)) = 1;</pre><p>Partial measurements with noise added.</p><pre class="codeinput">F = fftshift(fft2(M)/n);<span class="comment">% noise level</span>sigma = .01;<span class="comment">% selection operator M->y</span>y = F(mask==1) + sigma*randn(Q,1);</pre><p>Pseudo inverse reconstruction.</p><pre class="codeinput">Mpseudo = zeros(n,n);Mpseudo(mask==1) = y;Mpseudo = real( ifft2( fftshift(Mpseudo)*n ) );</pre><p>We display the mask and the pseudo inverse reconstruction.</p><pre class="codeinput">clf;imageplot(mask, <span class="string">'Fourier mask'</span>, 1,2,1);imageplot(Mpseudo, <span class="string">'Pseudo inverse reconstruction'</span>, 1,2,2);</pre><img vspace="5" hspace="5" src="index_01.png"> <h2>Wavelet Reconstruction from Partial Tomoraphic Measurements<a name="14"></a></h2> <p>Total variation is well suited for cartoon image, but can be improved on complicated natural images by using a translation invariant wavelet transfrm. One can reconstruct the coefficients <tt>MW</tt> in this wavelet frame using L1 minimization. This is performed with iterative thresholding, that solves the following minimization. </p> <p><tt>min_{MW} norm(y - Tomography(Wavelet(MW)))^2 + lambda*norm(MW, 1)</tt> (<b>*</b>) </p> <p>where <tt>Wavelet(.)</tt> is the wavelet reconstruction operator. </p> <p>Initialization.</p><pre class="codeinput">lambda = .1;options.ti = 1;Jmin = 4;M1 = zeros(n,n);MW = perform_wavelet_transf(M1, Jmin, +1, options);</pre><p>Step 1: Gradient descent step. MW <- MW + Phi*( y - Phi MW);</p><pre class="codeinput">M1 = perform_wavelet_transf(MW, Jmin, -1, options);<span class="comment">% measures</span>F1 = fftshift(fft2(M1)/n); y1 = F1(mask==1);<span class="comment">% residual</span>RF = zeros(n,n);RF(mask==1) = y-y1;R = real( ifft2(fftshift(RF)*n) );<span class="comment">% descent</span>MW = MW + perform_wavelet_transf(R, Jmin, +1, options);</pre><p>Step 2: Soft thresholding step. Here we do not threshold the low frequency residual.</p><pre class="codeinput">MWl = MW(:,:,1);MW = perform_thresholding(MW, lambda, <span class="string">'soft'</span>);MW(:,:,1) = MWl;</pre><p>Step 3: Update the <tt>lambda</tt> parameter to fit the noise level. </p><pre class="codeinput">lambda = lambda * sqrt(Q)*sigma / norm( y-y1, <span class="string">'fro'</span> );</pre><p><i>Exercice 1:</i> (the solution is <a href="../private/inverse_cs_fourier/exo1.m">exo1.m</a>) Implement the iterative thresholding algorithm to minimize the wavelet sparsity. At each step, modify the <tt>lambda</tt> parameter to fit the noise level. </p><pre class="codeinput">exo1;</pre><img vspace="5" hspace="5" src="index_02.png"> <p class="footer"><br> Copyright ® 2008 Gabriel Peyre<br></p> </div> <!--##### SOURCE BEGIN #####%% Reconstruction from Compressive Fourier Measurements% This numerical tour explores the reconstruction from randomized % Fourier measurement with wavelet 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/');%% 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 measurementsrho = .2;% total number of measuresQ = round(rho*n^2);% randomized valuesA = 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% maskmask = zeros(n,n);mask(I(1:Q)) = 1;%%% Partial measurements with noise added.F = fftshift(fft2(M)/n);% noise levelsigma = .01;% selection operator M->yy = F(mask==1) + sigma*randn(Q,1);%%% Pseudo inverse reconstruction.Mpseudo = zeros(n,n);Mpseudo(mask==1) = y;Mpseudo = real( ifft2( fftshift(Mpseudo)*n ) );%%% 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% Total variation is well suited for cartoon image, but can be improved on% complicated natural images by using a translation invariant wavelet% transfrm. One can reconstruct the coefficients |MW| in this wavelet% frame using L1 minimization. This is performed with iterative% thresholding, that solves the following minimization.%%% |min_{MW} norm(y - Tomography(Wavelet(MW)))^2 + lambda*norm(MW, 1)| (***)%%% where |Wavelet(.)| is the wavelet reconstruction operator.%%% Initialization.lambda = .1;options.ti = 1;Jmin = 4;M1 = zeros(n,n);MW = perform_wavelet_transf(M1, Jmin, +1, options);%%% Step 1: Gradient descent step.% MW <- MW + Phi*( y - Phi MW);M1 = perform_wavelet_transf(MW, Jmin, -1, options);% measuresF1 = fftshift(fft2(M1)/n); y1 = F1(mask==1);% residualRF = zeros(n,n);RF(mask==1) = y-y1;R = real( ifft2(fftshift(RF)*n) );% descentMW = MW + perform_wavelet_transf(R, Jmin, +1, options);%%% Step 2: Soft thresholding step.% Here we do not threshold the low frequency residual.MWl = MW(:,:,1);MW = perform_thresholding(MW, lambda, 'soft');MW(:,:,1) = MWl;%%% Step 3: Update the |lambda| parameter to fit the noise level.lambda = lambda * sqrt(Q)*sigma / norm( y-y1, 'fro' ); %%% _Exercice 1:_ (the solution is <../private/inverse_cs_fourier/exo1.m exo1.m>)% Implement the iterative thresholding algorithm to minimize the wavelet sparsity. At% each step, modify the |lambda| parameter to fit the noise level.exo1;##### SOURCE END #####--> </body></html>⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -