Absorption volume reconstruction: added more freedom in the rotation of the geometry

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

5_reconstruction/gtAstraAbsorptionVolume.m

-function [abs, proj] = gtAstraAbsorptionVolume(parameters, proj)
+function [abs_struct, proj] = gtAstraAbsorptionVolume(parameters, proj, varargin)
 % GTASTRAABSORPTIONVOLUME reconstruct absorption volume using 3D SIRT or 2D FBP
 %     abs = gtAstraAbsorptionVolume(parameters)
 %     -----------------------------------------
@@ -11,6 +11,8 @@ function [abs, proj] = gtAstraAbsorptionVolume(parameters, proj)
 %     Version 002 26-09-2012 by LNervo
 %       Formatting, cleaning, commenting
 
+    conf = struct('save', false, 'filename', 'volume_absorption');
+    conf = parse_pv_pairs(conf, varargin);
 
     if ~exist('parameters','var') || isempty(parameters)
         parameters = gtLoadParameters();
@@ -48,26 +50,41 @@ function [abs, proj] = gtAstraAbsorptionVolume(parameters, proj)
     end
     proj.num_iter = rec_abs.num_iter;
 
+    if (~isfield(rec_abs, 'options') || isempty(rec_abs.options))
+        rec_abs_options = {};
+    else
+        rec_abs_options = [fieldnames(rec_abs.options), struct2cell(rec_abs.options)]';
+    end
+
     fprintf('Reconstructing..')
     c = tic();
     is_cone = isfield(parameters.labgeo, 'sourcepoint') && (~isempty(parameters.labgeo.sourcepoint));
     switch (rec_abs.algorithm)
         case {'3DART', 'SIRT'}
-            abs = gtAstra3D(proj, 'is_cone', is_cone);
+            vol = gtAstra3D(proj, 'is_cone', is_cone);
         case {'3DTV'}
-            if (~isfield(rec_abs, 'options') || isempty(rec_abs.options))
-                rec_abs_options = {};
-            else
-                rec_abs_options = [fieldnames(rec_abs.options), struct2cell(rec_abs.options)]';
-            end
-            abs = gtAstra3DTV(proj, 'is_cone', is_cone, rec_abs_options{:});
+            vol = gtAstra3DTV(proj, 'is_cone', is_cone, rec_abs_options{:});
         case '2DFBP'
-            abs = gtAstra_FBP2D(parameters, proj);
+            vol = gtAstra_FBP2D(parameters, proj);
         otherwise
             error('gtAstraAbsorptionVolume:wrong_algorithm', ...
                 'Algorithm for reconstruction not known. Please check parameters file');
     end
     fprintf('\b\b: Done. (%f seconds)\n', toc(c))
+
+    if (conf.save)
+        absorption_volume = fullfile(parameters.acq.dir, '5_reconstruction', [conf.filename '.mat']);
+        abs_struct = struct('abs', vol, 'rot_angle', [], 'rot_axis', []);
+        if (isfield(proj, 'rot_angle'))
+            abs_struct.rot_angle = proj.rot_angle;
+            if (isfield(proj, 'rot_axis') && ~isempty(proj.rot_axis))
+                abs_struct.rot_axis = proj.rot_axis;
+            else
+                abs_struct.rot_axis = [0 0 1];
+            end
+        end
+        save(absorption_volume, '-struct', 'abs_struct', '-v7.3');
+    end
 end
 
 

5_reconstruction/gtAstraPrepareAbsorptionStack.m

-function proj = gtAstraPrepareAbsorptionStack(parameters, rot_angle)
+function proj = gtAstraPrepareAbsorptionStack(parameters, varargin)
 % GTASTRAPREPAREABSORPTIONSTACK  Assembles the 3D projection stack of absorption
 %                                images and saves them to file absorption_stack.mat
 %
@@ -20,9 +20,12 @@ function proj = gtAstraPrepareAbsorptionStack(parameters, rot_angle)
 %     Version 002 26-09-2012 by LNervo
 %       Formatting, cleaning, commenting
 
-if (~exist('rot_angle', 'var') || isempty(rot_angle))
-    rot_angle = 0;
-end
+conf = struct( ...
+    'rot_angle', [], ...
+    'rot_axis', [0 0 1], ...
+    'save', true, ...
+    'filename', 'absorption_stack' );
+conf = parse_pv_pairs(conf, varargin);
 
 acq = parameters.acq;
 
@@ -64,7 +67,7 @@ end
 fprintf('Building geometry..')
 imgs_nums = 0 : rec_abs.interval : totproj-1;
 Omega_rad = angle_rad(imgs_nums + 1);
-Omega_deg = Omega_rad * 180 / pi + rot_angle;
+Omega_deg = Omega_rad * 180 / pi;
 
 bbpos = repmat(acq.bb, [numel(Omega_deg), 1]);
 % bbpos = [gtGeoGrainCenterSam2Det(samgeo.orig, 0, parameters) - (acq.bb(3:4) - 1) / 2, acq.bb(3:4)];
@@ -73,7 +76,16 @@ vol_center = repmat(samgeo.orig, [numel(Omega_deg), 1]);
 
 geom = gtGeoProjForReconstruction([], Omega_deg', vol_center, bbpos, [], ...
     detgeo, labgeo, samgeo, recgeo, 'ASTRA_absorption');
-geom(:, 10:12) = round(geom(:, 10:12));
+if (~isempty(conf.rot_angle))
+    rot_comp = gtMathsRotationMatrixComp(conf.rot_axis', 'col');
+    rot_tensor = gtMathsRotationTensor(conf.rot_angle, rot_comp);
+
+    geom = reshape(geom, [], 3, 4);
+    geom = gtMathsMatrixProduct(geom, rot_tensor);
+    geom = reshape(geom, [], 12);
+end
+% This was ust a trick to correct for tilts, but it was probably a bad idea
+% geom(:, 10:12) = round(geom(:, 10:12));
 
 fprintf('\b\b: Done.\nBuilding sinogram: ');
 c = tic();
@@ -106,6 +118,7 @@ end
 
 fprintf('...done!\n')
 
+% Dealing with the padding
 new_stack_size = stack_size + [2, 0, 2] * rec_abs.padding;
 new_stack_shifts = [1, 0, 1] * rec_abs.padding;
 stack(stack < 0) = 0;
@@ -113,19 +126,24 @@ new_stack = gtPlaceSubVolume(zeros(new_stack_size, 'single'), stack, new_stack_s
 new_stack(1:rec_abs.padding, :) = new_stack((rec_abs.padding + 1) * ones(rec_abs.padding, 1), :);
 new_stack((end - rec_abs.padding + 1):end, :) = new_stack((end - rec_abs.padding) * ones(rec_abs.padding, 1), :);
 
-proj.padding    = rec_abs.padding;
-proj.stack      = new_stack;
-proj.Omega      = Omega_rad;
-proj.geom       = geom;
-proj.num_rows   = new_stack_size(3);
-proj.num_cols   = new_stack_size(1);
-proj.num_iter   = rec_abs.num_iter;
-proj.vol_size_x = new_stack_size(1);
-proj.vol_size_y = new_stack_size(1);
-proj.vol_size_z = new_stack_size(3);
-
-name = fullfile(acq.dir, 'absorption_stack.mat');
-disp(['Writing sinogram in ' name]);
-save(name, '-struct', 'proj');
-
-end % end of function
+proj = struct( ...
+    'padding', rec_abs.padding, ...
+    'rot_angle', conf.rot_angle, ... % Additional rotation
+    'rot_axis', conf.rot_axis, ...   % Additional rotation
+    'stack', new_stack, ...
+    'Omega', Omega_rad, ...
+    'geom', geom, ...
+    'num_rows', new_stack_size(3), ...
+    'num_cols', new_stack_size(1), ...
+    'num_iter', rec_abs.num_iter, ...
+    'vol_size_x', new_stack_size(1), ...
+    'vol_size_y', new_stack_size(1), ...
+    'vol_size_z', new_stack_size(3) );
+
+if (conf.save)
+    name = fullfile(acq.dir, [conf.filename '.mat']);
+    disp(['Writing sinogram in ' name]);
+    save(name, '-struct', 'proj');
+end
+
+end

5_reconstruction/gtSetupReconstruction.m

@@ -96,43 +96,39 @@ if (~exist(absorption_volume,'file'))
         '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', '-v7.3');
+        abs_struct = gtAstraAbsorptionVolume(parameters, [], 'save', true);
         sample.absVolFile = fullfile('5_reconstruction', 'volume_absorption.mat');
-        save(samplefile, 'sample');
+        sample.saveToFile();
     end
 else
     check = inputwdefault('Absorption volume found - would you like to overwrite it? [y/n]', 'n');
     if (strcmpi(check, 'y'))
-        abs = gtAstraAbsorptionVolume(parameters);
-        save(absorption_volume, 'abs', '-v7.3');
+        abs_struct = gtAstraAbsorptionVolume(parameters, [], 'save', true);
         sample.absVolFile = fullfile('5_reconstruction', 'volume_absorption.mat');
-        save(samplefile, 'sample');
+        sample.saveToFile();
     end
 end
 
 % add the possibility to redo volume threshold, even if existing
 absorption_mask = fullfile(parameters.acq.dir, '5_reconstruction', 'volume_mask.mat');
-if (~exist(absorption_mask,'file'))
+if (~exist(absorption_mask, 'file'))
     check = inputwdefault('Can not find absorption mask. Would you like to segment absorption volume? [y/n]', 'y');
     if (strcmpi(check, 'y'))
-        if (~exist('abs', 'var'))
-            abs = load(absorption_volume, 'abs');
-            abs = abs.abs;
+        if (~exist('abs_struct', 'var'))
+            abs_struct = load(absorption_volume);
         end
         disp([ 'Please create file volume_mask.mat from absorption volume ' ...
             '(right click in the volume to select a seed point and then press ''Segment'')' ]);
-        GtGuiThresholdVolume(abs);
+        GtGuiThresholdVolume(abs_struct.abs);
     end
 else
     check = inputwdefault('Absorption mask found - would you like to overwrite it? [y/n]', 'n');
     if (strcmpi(check, 'y'))
-        if (~exist('abs', 'var'))
-            abs = load(absorption_volume, 'abs');
-            abs = abs.abs;
+        if (~exist('abs_struct', 'var'))
+            abs_struct = load(absorption_volume);
         end
         disp('right click in the volume to select a seed point and then press ''Segment''');
-        GtGuiThresholdVolume(abs);
+        GtGuiThresholdVolume(abs_struct.abs);
     end
 end