1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
|
% 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(:)));
%
|
