diff options
Diffstat (limited to 'demo/Demo2.m')
| -rw-r--r-- | demo/Demo2.m | 156 |
1 files changed, 0 insertions, 156 deletions
diff --git a/demo/Demo2.m b/demo/Demo2.m deleted file mode 100644 index 3c1592c..0000000 --- a/demo/Demo2.m +++ /dev/null @@ -1,156 +0,0 @@ -% Demonstration of tomographic reconstruction from noisy and corrupted by -% artifacts undersampled projection data using Students t penalty -% This is the missing wedge demo, run it after DemoFISTA_StudT - -% see ReadMe file for instructions -% clear all -% close all - -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.0); - -add_wedges % apply the missing wedge mask - -%% -fprintf('%s\n', 'Direct reconstruction using FBP...'); -FBP_1 = iradon(MW_sino_artifacts', theta, N); - -fprintf('%s %.4f\n', 'RMSE for FBP reconstruction:', RMSE(FBP_1(:), phantom(:))); - -figure(1); -% set(gcf, 'Position', get(0,'Screensize')); -subplot_tight(1,2,1, [0.05 0.05]); imshow(FBP_1,[-2 0.8]); 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-LS without regularization...'); -clear params -% define parameters -params.sino = MW_sino_artifacts; -params.N = N; % image size -params.angles = theta_rad; % angles in radians -params.iterFISTA = 132; %max number of outer iterations -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.show = 0; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -params.weights = Dweights; % statistical weighting -tic; [X_FISTA, error_FISTA, obj_FISTA, sinoFISTA] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS reconstruction:', min(error_FISTA(:))); - -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-LS 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-LS-TV...'); -clear params -% define parameters -params.sino = MW_sino_artifacts; -params.N = N; % image size -params.angles = theta_rad; % angles in radians -params.iterFISTA = 200; % max number of outer iterations -params.lambdaTV = 5.39e-05; % regularization parameter for TV problem -params.tol = 1.0e-04; % tolerance to terminate TV iterations -params.iterTV = 20; % the max number of TV iterations -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.weights = Dweights; % statistical weighting -params.show = 0; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -tic; [X_FISTA_TV, error_FISTA_TV, obj_FISTA_TV, sinoFISTA_TV] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS-TV reconstruction:', min(error_FISTA_TV(:))); - -figure(4); clf -subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-LS-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.sino = MW_sino_artifacts; -params.N = N; % image size -params.angles = theta_rad; % angles in radians -params.iterFISTA = 250; % max number of outer iterations -params.lambdaTV = 0.0019; % regularization parameter for TV problem -params.tol = 1.0e-04; % tolerance to terminate TV iterations -params.iterTV = 20; % the max number of TV iterations -params.X_ideal = phantom; % ideal phantom -params.ROI = ROI; % phantom region-of-interest -params.weights = Dweights; % statistical weighting -params.lambdaR_L1 = 0.002; % parameter to sparsify the "rings vector" -params.show = 0; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -tic; [X_FISTA_GH_TV, error_FISTA_GH_TV, obj_FISTA_GH_TV, sinoFISTA_GH_TV] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction:', min(error_FISTA_GH_TV(:))); - -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.sino = MW_sino_artifacts; -params.N = N; % image size -params.angles = theta_rad; % angles in radians -params.iterFISTA = 80; % max number of outer iterations -% params.L_const = 80000; % Lipshitz constant (can be chosen manually to accelerate convergence) -params.lambdaTV = 0.0016; % regularization parameter for TV problem -params.tol = 1.0e-04; % tolerance to terminate TV iterations -params.iterTV = 20; % the max number of TV iterations -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 = 0; % visualize reconstruction on each iteration -params.slice = 1; params.maxvalplot = 0.6; -tic; [X_FISTA_student_TV, error_FISTA_student_TV, obj_FISTA_student_TV, sinoFISTA_student_TV] = FISTA_REC(params); toc; - -fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction:', min(error_FISTA_student_TV(:))); - -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_2(:), phantom(:))); -fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS reconstruction:', min(error_FISTA(:))); -fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS-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 |