diff options
Diffstat (limited to 'supp')
| -rw-r--r-- | supp/RMSE.m | 12 | ||||
| -rw-r--r-- | supp/add_wedges.m | 7 | ||||
| -rw-r--r-- | supp/zing_rings_add.m | 45 |
3 files changed, 38 insertions, 26 deletions
diff --git a/supp/RMSE.m b/supp/RMSE.m index 734e4b8..002f776 100644 --- a/supp/RMSE.m +++ b/supp/RMSE.m @@ -1,7 +1,7 @@ -function err = RMSE(signal1, signal2) -%RMSE Root Mean Squared Error - -err = sum((signal1 - signal2).^2)/length(signal1); % MSE -err = sqrt(err); % RMSE - +function err = RMSE(signal1, signal2)
+%RMSE Root Mean Squared Error
+
+err = sum((signal1 - signal2).^2)/length(signal1); % MSE
+err = sqrt(err); % RMSE
+
end
\ No newline at end of file diff --git a/supp/add_wedges.m b/supp/add_wedges.m index 8b8f2a7..5bc215c 100644 --- a/supp/add_wedges.m +++ b/supp/add_wedges.m @@ -27,4 +27,9 @@ grayImage(173:end,:) = 0; %figure; imshow(grayImage, [0 1]); MW_sino_artifacts = sino_zing_rings.*double(grayImage); -%Dweights = Dweights.*double(grayImage);
\ No newline at end of file +% !!! +% note that we do not re-calculate Dweights for MW_sino_artifacts +% if one does: Dweights = Dweights.*double(grayImage); +% then the MW artifacts will be reduced substantially, +% through weighting. However we would like to explore +% the effect of the penalty instead.
\ No newline at end of file diff --git a/supp/zing_rings_add.m b/supp/zing_rings_add.m index 023ac27..d197b1f 100644 --- a/supp/zing_rings_add.m +++ b/supp/zing_rings_add.m @@ -1,14 +1,23 @@ % uncomment this part of script to generate data with different noise characterisitcs fprintf('%s\n', 'Generating Projection Data...'); -multfactor = 1000; + % Creating RHS (b) - the sinogram (using a strip projection model) -vol_geom = astra_create_vol_geom(N, N); -proj_geom = astra_create_proj_geom('parallel', 1.0, P, theta_rad); -proj_id_temp = astra_create_projector('strip', proj_geom, vol_geom); -[sinogram_id, sinogramIdeal] = astra_create_sino(phantom./multfactor, proj_id_temp); -astra_mex_data2d('delete',sinogram_id); -astra_mex_algorithm('delete',proj_id_temp); +% vol_geom = astra_create_vol_geom(N, N); +% proj_geom = astra_create_proj_geom('parallel', 1.0, P, theta_rad); +% proj_id_temp = astra_create_projector('strip', proj_geom, vol_geom); +% [sinogram_id, sinogramIdeal] = astra_create_sino(phantom, proj_id_temp); +% astra_mex_data2d('delete',sinogram_id); +% astra_mex_algorithm('delete',proj_id_temp); + +%% +% inverse crime data generation +[sino_id, sinogramIdeal] = astra_create_sino3d_cuda(phantom, proj_geom, vol_geom); +astra_mex_data3d('delete', sino_id); + +% [id,x] = astra_create_backprojection3d_cuda(sinogramIdeal, proj_geom, vol_geom); +% astra_mex_data3d('delete', id); +%% % % % adding Gaussian noise % eta = 0.04; % Relative noise level @@ -20,7 +29,7 @@ astra_mex_algorithm('delete',proj_id_temp); %% % adding zingers val_offset = 0; -sino_zing = sinogramIdeal; +sino_zing = sinogramIdeal'; vec1 = [60, 80, 80, 70, 70, 90, 90, 40, 130, 145, 155, 125]; vec2 = [350, 450, 190, 500, 250, 530, 330, 230, 550, 250, 450, 195]; for jj = 1:length(vec1) @@ -58,19 +67,17 @@ sino_zing_rings(sino_zing_rings < 0) = 0; % adding Poisson noise dose = 50000; -dataExp = dose.*exp(-sino_zing_rings); % noiseless raw data -dataPnoise = astra_add_noise_to_sino(dataExp,2*dose); % pre-log noisy raw data (weights) -Dweights = dataPnoise; -sinogram = log(dose./dataPnoise); %log corrected data -> sinogram -sinogram = abs(sinogram); -clear dataPnoise dataExp +multifactor = 0.002; -% normalizing -sinogram = sinogram.*multfactor; -sino_zing_rings = sinogram; -Dweights = multfactor./Dweights; +dataExp = dose.*exp(-sino_zing_rings*multifactor); % noiseless raw data +dataPnoise = astra_add_noise_to_sino(dataExp, dose); % pre-log noisy raw data (weights) +sino_zing_rings = log(dose./max(dataPnoise,1))/multifactor; %log corrected data -> sinogram +Dweights = dataPnoise'; % statistical weights +sino_zing_rings = sino_zing_rings'; +clear dataPnoise dataExp + +% w = dose./exp(sinogram*multifactor); % getting back raw data from log-cor -% % figure(1); % set(gcf, 'Position', get(0,'Screensize')); % subplot(1,2,1); imshow(phantom,[0 0.6]); title('Ideal Phantom'); colorbar; |