Major reorganization, updated routines

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;