change to FBP as default Reconstruction method for absorption volumes

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@430 4c865b51-4357-4376-afb4-474e03ccb993

5_reconstruction/gtAstraPrepareAbsorptionStack.m

@@ -38,7 +38,7 @@ rotcomp = gtMathsRotationMatrixComp(geo.rotdir);
 
 
 % stack = zeros(acq.bb(4),acq.bb(3),totproj/rec.absinterval);
-stack = zeros(acq.bb(3), totproj/rec.absinterval, acq.bb(4));
+stack = zeros(acq.bb(4), totproj/rec.absinterval, acq.bb(3));
 geom  = zeros(totproj/rec.absinterval, 12);
 
 jj = 1;
@@ -55,7 +55,7 @@ for ii = 0 : rec.absinterval : totproj-1
 
   Omega(jj) = angle_rad(ii + 1);   
   %%
-  stack(:, jj, :)  = -log(edf_read(absname))';
+  stack(:, jj, :)  = -log(edf_read(absname));
   Rottens = gtMathsRotationTensor(-Omega(jj)*180/pi, rotcomp); 
   detdiru = geo.detdiru * Rottens;
   detdirv = geo.detdirv * Rottens;

5_reconstruction/gtAstra_FBP2D.m

-function abs=gtAstra_FBP2D(parameters,proj)
+function abs=gtAstra_FBP2D(parameters, proj)
 
 acq=parameters.acq;
 rec=parameters.rec;
-geo=parameters.geo;
-
-
-% name=sprintf('%s/absorption_stack.mat',acq.dir);
-% if exist(name,'file')
-%     check=inputwdefault('Absorption stack already exists - do you want to (re)-create it ? y/n', 'n');
-%     if strcmpi(check,'n')
-%         disp('reading projections...');
-%         load(name);
-%     else
-%         proj=gtAstra_prepareAbsorptionstack(parameters);
-%     end
-%     
-% else
-%     
-%     proj=gtAstra_prepareAbsorptionstack(parameters);
-% end
-
+geo=parameters.labgeo;
+
+if ~exist('proj', 'var')
+    name=sprintf('%s/absorption_stack.mat',acq.dir);
+    if exist(name,'file')
+        check=inputwdefault('Absorption stack already exists - do you want to (re)-create it ? y/n', 'n');
+        if strcmpi(check,'n')
+            disp('reading projections...');
+            proj=load(name);
+        else
+            proj=gtAstraPrepareAbsorptionStack(parameters);
+        end
+
+    else
+
+        proj=gtAstraPrepareAbsorptionStack(parameters);
+    end
+end
 
     
 %% now start FBP reconstruction of absorption volume
 
-proj_data=permute(proj.stack,[3 2 1]);   % transforms stack of projection images to stack of sinograms
-
+%proj_data=permute(proj.stack,[3 2 1]);   % transforms stack of projection images to stack of sinograms
+proj_data = proj.stack;
 
 vsize=size(proj.stack,1);  % numer of slices to be reconstructed
-hsize=size(proj.stack,2);  % number of pixels in projection image
-nproj=size(proj.stack,3);  % number of projections
+hsize=size(proj.stack,3);  % number of pixels in projection image
+nproj=size(proj.stack,2);  % number of projections
 
 hsize_zp=ceil(hsize*sqrt(2));   % zeropad projection images
 offset=ceil((hsize_zp-hsize)/2);
@@ -68,7 +68,7 @@ for i=1:vsize
     alg_id = astra_mex_algorithm('create', cfg);
     
     % 'zeropadding' projections with the value
-    sino=gtPlaceSubImage(proj_data(:,:,i),sino,offset,1);
+    sino=gtPlaceSubImage(squeeze(proj_data(i,:,:)),sino,offset,1);
     
     sino(:,1:offset-1)=repmat(sino(:,offset),1,offset-1);
     sino(:,offset+hsize:end)=repmat(sino(:,offset+hsize-1),1,hsize_zp-offset-hsize+1);

5_reconstruction/gtSetupReconstruction.m

@@ -32,17 +32,22 @@ end
 
 absorption_volume = fullfile(acq.dir, '5_reconstruction', 'volume_absorption.mat');
 if ~exist(absorption_volume,'file')
-    disp('Can not find absorption volume - launching reconstruction of absorption volume');
-    abs = gtAstraAbsorptionVolume(parameters);
-    save(absorption_volume, 'abs');
-    sample.absVolFile = fullfile('5_reconstruction', 'volume_absorption.mat');
-    save(samplefile, 'sample'); 
+    check = inputwdefault('Can not find absorption volume - would you like to launch reconstruction of absorption volume (y/n)', 'y');
+    if strcmpi(check, 'y')
+        abs = gtAstraAbsorptionVolume(parameters);
+        save(absorption_volume, 'abs');
+        sample.absVolFile = fullfile('5_reconstruction', 'volume_absorption.mat');
+        save(samplefile, 'sample'); 
+    end
 end
 
 absorption_mask = fullfile(acq.dir, '5_reconstruction', 'volume_mask.mat');
 if ~exist(absorption_mask,'file')
-    disp('Please create file volume_mask.mat from absorption volume (right click in the volume to select a seed point)'); 
-    GtGuiThresholdVolume(abs);
+    check = inputwdefault('Can not find absorption mask - would you like to segment absorption volume (y/n)', 'y');
+    if strcmpi(check, 'y')
+        disp('Please create file volume_mask.mat from absorption volume (right click in the volume to select a seed point)'); 
+        GtGuiThresholdVolume(abs);
+    end
 end
 
 

zUtil_DataStructures/gtCreateGtSample.m

@@ -53,8 +53,10 @@ function sample = gtCreateGtSample()
             phase.setCompleteness(ii, grainInfo.completeness);
             difspotID = grainInfo.difspotID;
             for kk = 1 : length(grainInfo.selected)
-                spot2grain{difspotID(grainInfo.selected(kk))} = [spot2grain{difspotID(grainInfo.selected(kk))}, ii]; 
-                spot2phase{difspotID(grainInfo.selected(kk))} = [spot2phase{difspotID(grainInfo.selected(kk))}, jj];    
+                if grainInfo.selected(kk)
+                    spot2grain{difspotID(kk)} = [spot2grain{difspotID(kk)}, ii]; 
+                    spot2phase{difspotID(kk)} = [spot2phase{difspotID(kk)}, jj];    
+                end
             end
             
             fprintf('\b\b\b\b\b\b\b\b\b');