diff options
author | Daniil Kazantsev <dkazanc@hotmail.com> | 2017-07-03 22:35:23 +0100 |
---|---|---|
committer | Daniil Kazantsev <dkazanc@hotmail.com> | 2017-07-03 22:35:23 +0100 |
commit | 329a104d4cb5ba50a59fb80e58de0453ba49f075 (patch) | |
tree | 278a879fb4000c488b3e07dbd6cac6bb9d9aeb7e /demos/Demo1.m | |
parent | e55c200119ebf9fd42755cb2fea7c3d286ffe96b (diff) | |
download | regularization-329a104d4cb5ba50a59fb80e58de0453ba49f075.tar.gz regularization-329a104d4cb5ba50a59fb80e58de0453ba49f075.tar.bz2 regularization-329a104d4cb5ba50a59fb80e58de0453ba49f075.tar.xz regularization-329a104d4cb5ba50a59fb80e58de0453ba49f075.zip |
Major reorganization, updated routines
Diffstat (limited to 'demos/Demo1.m')
-rw-r--r-- | demos/Demo1.m | 174 |
1 files changed, 174 insertions, 0 deletions
diff --git a/demos/Demo1.m b/demos/Demo1.m new file mode 100644 index 0000000..486b97c --- /dev/null +++ b/demos/Demo1.m @@ -0,0 +1,174 @@ +% Demonstration of tomographic reconstruction from noisy and corrupted by
+% artifacts undersampled projection data using Students't penalty
+% Optimisation problem is solved using FISTA algorithm (see Beck & Teboulle)
+
+% see Readme file for instructions
+%%
+% compile MEX-files ones
+% cd ..
+% cd main_func
+% compile_mex
+% cd ..
+% cd demos
+%%
+
+close all;clc;clear all;
+% adding paths
+addpath('../data/');
+addpath('../main_func/');
+addpath('../supp/');
+
+load phantom_bone512.mat % load the phantom
+load my_red_yellowMAP.mat % load the colormap
+% load sino1.mat; % load noisy sinogram
+
+N = 512; % the size of the tomographic image NxN
+theta = 1:1:180; % acquisition angles (in parallel beam from 0 to Pi)
+theta_rad = theta*(pi/180); % conversion to radians
+P = 2*ceil(N/sqrt(2))+1; % the size of the detector array
+ROI = find(phantom > 0);
+
+% using ASTRA to set the projection geometry
+% potentially parallel geometry can be replaced with a divergent one
+Z_slices = 1;
+det_row_count = Z_slices;
+proj_geom = astra_create_proj_geom('parallel3d', 1, 1, det_row_count, P, theta_rad);
+vol_geom = astra_create_vol_geom(N,N,Z_slices);
+
+zing_rings_add; % generating data, adding zingers and stripes
+%%
+fprintf('%s\n', 'Direct reconstruction using FBP...');
+FBP_1 = iradon(sino_zing_rings', theta, N);
+
+fprintf('%s %.4f\n', 'RMSE for FBP reconstruction:', RMSE(FBP_1(:), phantom(:)));
+
+figure(1);
+subplot_tight(1,2,1, [0.05 0.05]); imshow(FBP_1,[0 0.6]); title('FBP reconstruction of noisy and corrupted by artifacts sinogram'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - FBP_1).^2,[0 0.1]); title('residual: (ideal phantom - FBP)^2'); colorbar;
+colormap(cmapnew);
+
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-PWLS without regularization...');
+clear params
+% define parameters
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = sino_zing_rings; % sinogram
+params.iterFISTA = 45; %max number of outer iterations
+params.X_ideal = phantom; % ideal phantom
+params.ROI = ROI; % phantom region-of-interest
+params.show = 1; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6;
+params.weights = Dweights; % statistical weighting
+tic; [X_FISTA, output] = FISTA_REC(params); toc;
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction is:', min(error_FISTA(:)));
+error_FISTA = output.Resid_error; obj_FISTA = output.objective;
+
+figure(2); clf
+%set(gcf, 'Position', get(0,'Screensize'));
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA,[0 0.6]); title('FISTA-PWLS reconstruction'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA).^2,[0 0.1]); title('residual'); colorbar;
+colormap(cmapnew);
+figure(3); clf
+subplot_tight(1,2,1, [0.05 0.05]); plot(error_FISTA); title('RMSE plot'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); plot(obj_FISTA); title('Objective plot'); colorbar;
+colormap(cmapnew);
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-PWLS-TV...');
+clear params
+% define parameters
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = sino_zing_rings;
+params.iterFISTA = 45; % max number of outer iterations
+params.Regul_LambdaTV = 0.0015; % regularization parameter for TV problem
+params.X_ideal = phantom; % ideal phantom
+params.ROI = ROI; % phantom region-of-interest
+params.weights = Dweights; % statistical weighting
+params.show = 1; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6;
+tic; [X_FISTA_TV, output] = FISTA_REC(params); toc;
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS-TV reconstruction is:', min(error_FISTA_TV(:)));
+error_FISTA_TV = output.Resid_error; obj_FISTA_TV = output.objective;
+
+figure(4); clf
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-PWLS-TV reconstruction'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_TV).^2,[0 0.1]); title('residual'); colorbar;
+colormap(cmapnew);
+figure(5); clf
+subplot_tight(1,2,1, [0.05 0.05]); plot(error_FISTA_TV); title('RMSE plot'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); plot(obj_FISTA_TV); title('Objective plot'); colorbar;
+colormap(cmapnew);
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-GH-TV...');
+clear params
+% define parameters
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = sino_zing_rings;
+params.iterFISTA = 50; % max number of outer iterations
+params.Regul_LambdaTV = 0.0015; % regularization parameter for TV problem
+params.X_ideal = phantom; % ideal phantom
+params.ROI = ROI; % phantom region-of-interest
+params.weights = Dweights; % statistical weighting
+params.Ring_LambdaR_L1 = 0.002; % parameter to sparsify the "rings vector"
+params.Ring_Alpha = 20; % to accelerate ring-removal procedure
+params.show = 0; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6;
+tic; [X_FISTA_GH_TV, output] = FISTA_REC(params); toc;
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction is:', min(error_FISTA_GH_TV(:)));
+error_FISTA_GH_TV = output.Resid_error; obj_FISTA_GH_TV = output.objective;
+
+figure(6); clf
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_GH_TV,[0 0.6]); title('FISTA-GH-TV reconstruction'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]);imshow((phantom - X_FISTA_GH_TV).^2,[0 0.1]); title('residual'); colorbar;
+colormap(cmapnew);
+
+figure(7); clf
+subplot_tight(1,2,1, [0.05 0.05]); plot(error_FISTA_GH_TV); title('RMSE plot'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); plot(obj_FISTA_GH_TV); title('Objective plot'); colorbar;
+colormap(cmapnew);
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-Student-TV...');
+clear params
+% define parameters
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = sino_zing_rings;
+params.iterFISTA = 55; % max number of outer iterations
+params.L_const = 0.1; % Lipshitz constant (can be chosen manually to accelerate convergence)
+params.Regul_LambdaTV = 0.00152; % regularization parameter for TV problem
+params.X_ideal = phantom; % ideal phantom
+params.ROI = ROI; % phantom region-of-interest
+params.weights = Dweights; % statistical weighting
+params.fidelity = 'student'; % selecting students t fidelity
+params.show = 1; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6;
+params.initilize = 1; % warm start with SIRT
+tic; [X_FISTA_student_TV, output] = FISTA_REC(params); toc;
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction is:', min(error_FISTA_student_TV(:)));
+error_FISTA_student_TV = output.Resid_error; obj_FISTA_student_TV = output.objective;
+
+figure(8);
+set(gcf, 'Position', get(0,'Screensize'));
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_student_TV,[0 0.6]); title('FISTA-Student-TV reconstruction'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_student_TV).^2,[0 0.1]); title('residual'); colorbar;
+colormap(cmapnew);
+
+figure(9);
+subplot_tight(1,2,1, [0.05 0.05]); plot(error_FISTA_student_TV); title('RMSE plot'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); plot(obj_FISTA_student_TV); title('Objective plot'); colorbar;
+colormap(cmapnew);
+%%
+% print all RMSE's
+fprintf('%s\n', '--------------------------------------------');
+fprintf('%s %.4f\n', 'RMSE for FBP reconstruction:', RMSE(FBP_1(:), phantom(:)));
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction:', min(error_FISTA(:)));
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS-TV reconstruction:', min(error_FISTA_TV(:)));
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction:', min(error_FISTA_GH_TV(:)));
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction:', min(error_FISTA_student_TV(:)));
+%
\ No newline at end of file |