bpdq_demo_2d

 bpdq_demo_2d : Demo for BPDQ-TV reconstruction in 2-d

 This demo creates a synthetic angiogram image, takes randomly subsampled
 Fourier measurements, and quantizes the measurements.

 Then the BPDQ-TV algorithm is run to reconstruct the image from the 
 quantized measurements.
 
 This file is part of BPDQ Toolbox (Basis Pursuit DeQuantizer)
 Copyright (C) 2009, the BPDQ Team (see the file AUTHORS distributed with
 this library) (See the notice at the end of the file.)

0001 % bpdq_demo_2d : Demo for BPDQ-TV reconstruction in 2-d
0002 %
0003 % This demo creates a synthetic angiogram image, takes randomly subsampled
0004 % Fourier measurements, and quantizes the measurements.
0005 %
0006 % Then the BPDQ-TV algorithm is run to reconstruct the image from the
0007 % quantized measurements.
0008 %
0009 % This file is part of BPDQ Toolbox (Basis Pursuit DeQuantizer)
0010 % Copyright (C) 2009, the BPDQ Team (see the file AUTHORS distributed with
0011 % this library) (See the notice at the end of the file.)
0012 
0013 % With parameters as set in this demo, execution time was 120 seconds
0014 % on 2.4 GHz Intel Core 2 Duo processor.
0015 
0016 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0017 % Demo Parameters :
0018 % N - image dimension (image will be NxN)
0019 % nellipse - # of ellipses in syntheric angiogram
0020 % fsubfactor - Fourier subsampling factor
0021 % alpha - quantization bin width
0022 % p - BPDQ moment (p=2 corresponds to standard BPDN)
0023 
0024 % (Quantization threshold alpha is computed depending on nbins and on
0025 %  the signal instance)
0026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0027 N=256;
0028 fsubfactor=6;
0029 p=10;
0030 alpha=50;
0031 nellipse=10;
0032 aseed=6;
0033 seed=1;
0034 
0035 bpdq_opts={'dr_gamma',.1,'dr_lambda',1,'dr_verbosity',0,'dr_maxiter',50};
0036 
0037 dim=[N,N];
0038 fprintf('Generating %g x %g synthetic angiogram\n',N,N);
0039 im=bpdq_generate_angiogram(N,N,nellipse,aseed);
0040 
0041 K=.5*floor(N^2/fsubfactor);
0042 fprintf(['Selecting %g random fourier locations, i.e. 1/%i of the %g ' ...
0043          'available\n'],2*K, fsubfactor, N^2);
0044 rand('seed',seed);
0045 [useind,loc]=bpdq_random_fourier_locations(dim,K,'forcecenter',1);
0046 
0047 % fourier subsampling operators
0048 samp=@(x) bpdq_fft_subsample(x,useind,dim);
0049 samp_t=@(x) bpdq_fft_subsample_t(x,useind,dim); % transpose of samp
0050 
0051 %% utility functions
0052 % unpack complex vector into vector of reals
0053 c2r = @(x) [real(x(:));imag(x(:))];
0054 r2c = @(x) x(1:end/2)+sqrt(-1)*x(end/2+1:end);
0055 
0056 % "complete" measurement operators (from sparsity domain to measurements)
0057 A=@(x) samp((x(:)));
0058 At=@(x) vec(samp_t(x));
0059 
0060 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0061 % make measurements
0062 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0063 fprintf('Computing quantized measurements alpha=%g\n',alpha);
0064 
0065 y=samp(im);
0066 % quantize
0067 yq=bpdq_quantize(y,alpha);
0068 Nbins=2*max(abs(c2r(y)))/alpha; % equivalent # of bins
0069 % fprintf('%g equivalent number of quantization bins\n',Nbins);
0070 
0071 % compute epsilon for BPDQ program
0072 dM=numel(useind);
0073 M=2*dM;
0074 
0075 fprintf('Running BPDQ-TV reconstruction with p=%g\n',p);
0076 epsilon = bpdq_err_p(p,alpha,M);
0077 [xstar,D_bpdq]=bpdq_2d_tv(yq,A,At,dim,epsilon,p,bpdq_opts{:});
0078 im_bpdq=real(reshape(xstar,dim));
0079 SNR_bpdq= bpdq_compute_snr(im,im_bpdq);
0080 
0081 fprintf('Running BPDN-TV (p=2) reconstruction\n');
0082 epsilon2 = bpdq_err_p(2,alpha,M);
0083 [xstar2,D_bpdn]=bpdq_2d_tv(yq,A,At,dim,epsilon2,2,bpdq_opts{:});
0084 im_bpdn=real(reshape(xstar2,dim));
0085 SNR_bpdn=bpdq_compute_snr(im,im_bpdn);
0086 
0087 rmin=min([im(:);im_bpdn(:);im_bpdq(:)]);
0088 rmax=max([im(:);im_bpdn(:);im_bpdq(:)]);
0089 drange=[rmin,rmax];
0090 
0091 fprintf('Displaying results and convergence curves\n');
0092 figure(1)
0093 fprintf('Original Image in Figure 1\n');
0094 bpdq_show_im(im,drange,1);
0095 title('Original Image');
0096 
0097 figure(2)
0098 fprintf('BPDQ_%g Reconstruction in Figure 2\n', p);
0099 bpdq_show_im(im_bpdq,drange,1);
0100 title(sprintf('bpdq-tv (p=%g) reconstruction, SNR %g dB',p,SNR_bpdq));
0101 drawnow;
0102 
0103 figure(3)
0104 fprintf('BPDN Reconstruction in Figure 3\n', p);
0105 bpdq_show_im(im_bpdn,drange,1);
0106 title(sprintf('bpdn-tv reconstruction, SNR %g dB',SNR_bpdn));
0107 drawnow;
0108 
0109 figure(4)
0110 fprintf('Convergence Analysis in Figure 4\n');
0111 subplot(1,2,1)
0112 plot(log10(D_bpdq.dr_relerr_save));
0113 title('DR convergence (for bpdq)');
0114 xlabel('iteration');
0115 ylabel('log relative error');
0116 
0117 subplot(1,2,2)
0118 plot(log10(D_bpdn.dr_relerr_save));
0119 title('DR convergence (for bpdn)');
0120 xlabel('iteration')
0121 ylabel('log relative error');
0122 
0123 % The BPDQ Toolbox is free software: you can redistribute it and/or modify
0124 % it under the terms of the GNU General Public License as published by
0125 % the Free Software Foundation, either version 3 of the License, or
0126 % (at your option) any later version.
0127 %
0128 % The BPDQ Toolbox is distributed in the hope that it will be useful,
0129 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0130 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0131 % GNU General Public License for more details.
0132 %
0133 % You should have received a copy of the GNU General Public License
0134 % along with The BPDQ Toolbox.  If not, see <http://www.gnu.org/licenses/>.