Remove transparency constraint, allow saving absorption recon

RecGTVIC.m

@@ -1,7 +1,6 @@
 %% Tomographic Reconstruction with Gerchberg-Papoulis algorithm + Finite Object Support
 % Authors: Piotr Makowski, Wojciech Krauze, Paweł Ossowski
-% Warsaw University of Technology, 2023
-% Version: 0.7
+% Warsaw University of Technology, 2024
 
 clear; clc;
 close all
@@ -25,7 +24,7 @@
 resample_projections   = false; % default = false % resample projections to minimal safe resolution
 
 %% Save reconstruction?
-save_reconstruction    = 1;     % default = 0     % save the reconstruction REC(y,x,z)
+save_reconstruction    = 1;     % 0 - don't save; 1 - save RI; 2 - save RI and absorption (complex RI)
 
 %% Plots
 % show selected plots:
@@ -76,5 +75,5 @@
     save_recon;
 end
 
-if contains(plots,'1'); vis(RECON, [min(RECON(:)), max(RECON(:))]); colormap jet; end
+if contains(plots,'1'); vis(real(RECON), [min(real(RECON(:))), max(real(RECON(:)))]); colormap jet; end
 if contains(plots,'5'); Plots; end

lib/FDT.m

@@ -102,7 +102,7 @@
 % progress params
 Nc = 120;% cropped cube size for incremental difference evaluation
 % apply transparency constraint: imag(n_rec)=0
-do_TC = true;%true
+do_TC = false;% transparency constraint
 tukr = 0.25;%0.25  % Tukey window radius for projections (Approx=Born, interpFp=linc)
 % relaxed suppression of reconstruction padding  from Fourier oversampling (<=1, 0:disable)
 relaxLaR = 0;%1.00, 0.99(best with N_Kspace_z_padded_upsampled > N_Kspace_xy_padded), 0:disable !!!!
@@ -813,6 +813,6 @@
 KO(EWi) = EKO;
 
 % RECONSTRUCTION (y,x,z)
-RECON = permute(real(n_rec), [2 1 3]); % (y,x,z)
+RECON = permute((n_rec), [2 1 3]); % (y,x,z)
 % N_Kspace_z_padded_upsampled cropped to N_Kspace_xy_padded / alfa (cropping the Z axis to match sizes in X and Y directions):
 RECON = ndcrop(RECON, [N_Kspace_xy/projection_padding_xy N_Kspace_xy/projection_padding_xy N_Kspace_z_upsampled_cropped/Kspace_oversampling_z]);

lib/Plots.m

@@ -64,7 +64,7 @@
 % cross-section coords in [N_Kspace_xy_padded N_Kspace_xy_padded N_Kspace_xy_padded] cube
 cx0 = round(N_Kspace_xy_padded/projection_padding_xy/2)+1;
 cy0 = round(N_Kspace_xy_padded/projection_padding_xy/2)+1;
-cz0 = round(N_Kspace_xy_padded/projection_padding_xy/2)+0;
+cz0 = round(N_Kspace_z_padded_upsampled/2)+0;
 % % manual
 % cx0 = 202;%202(Meth)
 % cy0 = 194;%194(Meth)
@@ -75,12 +75,7 @@
 % Nc = 200;%100
 %%%%%%%%%%%%%%%%%%%%%%%
 
-
-if all(sino_params(5,:)==2) 
-    Kz_slice = 50;
-else; Kz_slice = size(KO,3)/2+1; 
-end
-
+Kz_slice = size(KO,3)/2+1; 
 
 % Collect reconstruction info (x,y,z)
 N_projection_padded = round(N_projection_padded/2)*2;
@@ -139,13 +134,13 @@
 if ~isempty(n_obj)		
     tmp_rec = OBJ(...
                 round(0.3*size(RECON,1)) : round(0.7*size(RECON,1)),...
-                round(0.3*size(RECON,1)) : round(0.7*size(RECON,1)),...
-                round(0.3*size(RECON,1)) : round(0.7*size(RECON,1)))-is_n0*n0;
+                round(0.3*size(RECON,2)) : round(0.7*size(RECON,2)),...
+                round(0.3*size(RECON,3)) : round(0.7*size(RECON,3)))-is_n0*n0;
 else
     tmp_rec = RECON(...
                 round(0.3*size(RECON,1)) : round(0.7*size(RECON,1)),...
-                round(0.3*size(RECON,1)) : round(0.7*size(RECON,1)),...
-                round(0.3*size(RECON,1)) : round(0.7*size(RECON,1)))-is_n0*n0;
+                round(0.3*size(RECON,2)) : round(0.7*size(RECON,2)),...
+                round(0.3*size(RECON,3)) : round(0.7*size(RECON,3)))-is_n0*n0;
 end
 nm = max(abs(tmp_rec(:)));
 clear tmp_rec
@@ -160,8 +155,8 @@
 
 %  slices
 corr_n0 = +(~is_n0&show_n0)*n0 -(~show_n0&is_n0)*n0;% correction for n0
-rec_zx = squeeze(RECON(cy,:,:) + corr_n0).';%(prim: make X horizontal)
-rec_yx = squeeze(RECON(:,:,cz) + corr_n0);
+rec_zx = real(squeeze(RECON(cy,:,:) + corr_n0).');%(prim: make X horizontal)
+rec_yx = real(squeeze(RECON(:,:,cz) + corr_n0));
 if ~isempty(OBJ) % calculate Quality Index for one slice
     obj_zx = squeeze(OBJ(cy,:,:) + corr_n0).'; %  one slice (prim: make X horizontal)
     diff_zx = rec_zx - obj_zx;
@@ -244,12 +239,12 @@
                   title 'RECON x-y'
     % 				  title(sprintf('cz=%d',cz)) 
                   caxis(clims)
-    if ~isempty(OBJ)
-                  title(sprintf('REC x-y; MAE=%.3fe-3, MSE=%.3fe-3',RMAEtab(nGPi+1)*1e3,RRMSEtab(nGPi+1)*1e3))
-    				  RMAE_REC = norm(RECON(:)-OBJ(:),1)/norm(OBJ(:),1);
-    				  RRMSE_REC = norm(RECON(:)-OBJ(:),2)/norm(OBJ(:),2);
-    				  title(sprintf('MAE=%.3fe-3, MSE=%.3fe-3',RMAE_REC*1e3,RRMSE_REC*1e3))
-    end
+    % if ~isempty(OBJ)
+    %               title(sprintf('REC x-y; MAE=%.3fe-3, MSE=%.3fe-3',RMAEtab(nGPi+1)*1e3,RRMSEtab(nGPi+1)*1e3))
+    % 				  RMAE_REC = norm(RECON(:)-OBJ(:),1)/norm(OBJ(:),1);
+    % 				  RRMSE_REC = norm(RECON(:)-OBJ(:),2)/norm(OBJ(:),2);
+    % 				  title(sprintf('MAE=%.3fe-3, MSE=%.3fe-3',RMAE_REC*1e3,RRMSE_REC*1e3))
+    % end
     subplot(248); hold on
                   if ~isempty(OBJ)
                   plot(real(obj_yx(cy,:)), 'b', 'LineWidth',1)

lib/save_recon.m

@@ -1,5 +1,9 @@
 disp('Saving reconstruction...')
 
+    if save_reconstruction == 1
+        RECON = real(RECON);
+    end
+
     FileName = 'REC[';
     if Masking
         FileName = strcat(FileName,'GPSC]');