6D-reconstructor: introduced specific projection for each W in the W-shape-functions

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

zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m

@@ -347,7 +347,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             fprintf('Computing blobs <-> sinograms coefficients: ')
             c = tic();
             % ASTRA Geometry
-            geometries = self.compute_projection_geometry_numinterp(bl, ref_gr, ors_allblobs);
+            geometries = self.compute_proj_geom_numinterp(bl, ref_gr, ors_allblobs);
 
             [geometries, offsets] = self.makeSubBlobGeometries(geometries, ...
                 proj_coeffs, extreemes_blobs_w_tot, w_tab, abs_lims_tot);
@@ -369,7 +369,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                     w_tab_oss = self.get_w_tab(orient_props_ss);
                     w_tab_oss = gtGrainAnglesTabularFix360deg(w_tab_oss, ref_ws, self.parameters);
 
-                    geometries_ss{ii_o} = self.compute_projection_geometry_numinterp(bl, ref_gr, orient_props_ss);
+                    geometries_ss{ii_o} = self.compute_proj_geom_numinterp(bl, ref_gr, orient_props_ss);
 
                     [geometries_ss{ii_o}, offsets_ss{ii_o}] = self.makeSubBlobGeometries( ...
                         geometries_ss{ii_o}, proj_coeffs, ...
@@ -408,7 +408,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
         end
 
         function algo_params = get_algorithm_params_shape_funcs(self, sampler, shape_funcs, sf_type)
-            fprintf('Extracting blobs on detector: ')
+            fprintf('- Extracting blobs on detector: ')
             c = tic();
             sel_bls = sampler.selected;
             bl = sampler.bl(sel_bls);
@@ -422,19 +422,6 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             sel_allb_idxs = sampler.ondet(sampler.included(sel_bls));
             omegas = ref_gr.allblobs.omega(sel_allb_idxs);
 
-            [orients, orients_ss] = sampler.get_orientations();
-
-            num_ospace_oversampling = prod(sampler.ss_factor);
-            if (num_ospace_oversampling == 1)
-                ors_struct = [orients{:}];
-            else
-                ors_struct = cat(4, orients_ss{:});
-                ors_struct = [ors_struct{:}];
-            end
-            % Array of structures that contain all the info relative to
-            % each orientation, for each blob
-            ors_allblobs = self.get_selected_props(ors_struct, sel_bls);
-
             [w_tab, ref_ws] = self.get_shape_funcs_w_tab(shape_funcs, omegas);
             [lims_blobs_w, lims_projs_w, paddings, ext_blobs_w_orig] = self.get_blob_lims(w_tab, bl, ref_ws);
 
@@ -450,28 +437,32 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                 blobs{ii} = permute(blobs{ii}, [1 3 2]);
             end
 
-            fprintf('Computing blobs <-> sinograms coefficients: ')
+            fprintf('- Computing projection geometries: ')
             c = tic();
+
+            [orients, orients_ss] = sampler.get_orientations();
+
+            num_ospace_oversampling = prod(sampler.ss_factor);
+            if (num_ospace_oversampling == 1)
+                ors_struct = [orients{:}];
+            else
+                ors_struct = cat(4, orients_ss{:});
+                ors_struct = [ors_struct{:}];
+            end
+            % Array of structures that contain all the info relative to
+            % each orientation, for each blob
+            ors_allblobs = self.get_selected_props(ors_struct, sel_bls);
+
+            % ASTRA Geometry
             switch (lower(sf_type))
                 case 'w'
-                    offsets = self.make_proj_geom_shape_funcs_w(...
-                        shape_funcs, lims_blobs_w, w_tab);
+                    [geometries, offsets] = self.compute_proj_geom_shapefuncs_w( ...
+                        bl, ref_gr, ors_allblobs, shape_funcs, lims_blobs_w);
                 case {'nw2uvw', 'uvw'}
+                    geometries = self.compute_proj_geom_numinterp(bl, ref_gr, ors_allblobs);
                     offsets = self.make_proj_geom_shape_funcs_uvw(...
                         shape_funcs, lims_blobs_w, w_tab);
             end
-            fprintf('Done (%f s).\n', toc(c));
-
-            % ASTRA Geometry
-            geometries = self.compute_projection_geometry_numinterp(bl, ref_gr, ors_allblobs);
-
-            if (isfield(bl, 'psf'))
-                psf = arrayfun(@(x){ permute(x.psf, [1 3 2]) }, bl);
-            elseif (isfield(proj, 'psf'))
-                psf = { permute(proj.psf, [1 3 2]) };
-            else
-                psf = {};
-            end
 
             tot_orient = sampler.get_total_num_orientations();
             if (num_ospace_oversampling == 1)
@@ -496,6 +487,15 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                 geometries = {};
                 offsets = {};
             end
+            fprintf('Done (%f s).\n', toc(c));
+
+            if (isfield(bl, 'psf'))
+                psf = arrayfun(@(x){ permute(x.psf, [1 3 2]) }, bl);
+            elseif (isfield(proj, 'psf'))
+                psf = { permute(proj.psf, [1 3 2]) };
+            else
+                psf = {};
+            end
 
             algo_params = struct( ...
                 'blobs', {blobs}, ...
@@ -754,7 +754,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             sliceLims = [floor(sliceLims(:, 1)), ceil(sliceLims(:, 2))] + 1;
         end
 
-        function geoms = compute_projection_geometry_numinterp(self, bl, ref_gr, ors_allblobs)
+        function geoms = compute_proj_geom_numinterp(self, bl, ref_gr, ors_allblobs)
             bls_bbuim = cat(1, bl.bbuim);
             bls_bbvim = cat(1, bl.bbvim);
             bls_bbsize = cat(1, bl.bbsize);
@@ -788,7 +788,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
         end
 
         function [geometries, offsets] = makeSubBlobGeometries( self, ...
-                geometries, proj_coeffs, extreemes_blobs_w, w_tab, abs_lims)
+                geometries, proj_coeffs, lims_blobs_w, w_tab, abs_lims)
             num_geoms = numel(geometries);
 
             w_int_tab(:, :, 1) = floor(w_tab);
@@ -796,7 +796,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             w_int_coeffs(:, :, 1) = 1 - (w_int_tab(:, :, 2) - w_tab);
             w_int_coeffs(:, :, 2) = 1 - w_int_coeffs(:, :, 1);
 
-            slice_pos_in_blob = w_int_tab - repmat(extreemes_blobs_w(:, 1), [1 num_geoms 2]) + 1;
+            slice_pos_in_blob = w_int_tab - repmat(lims_blobs_w(:, 1), [1 num_geoms 2]) + 1;
 
             % Assign offsets
             offsets = cell(num_geoms, 1);
@@ -863,48 +863,75 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             end
         end
 
-        function offsets = make_proj_geom_shape_funcs_w(self, shape_funcs, extreemes_blobs_w, w_tab)
-            num_geoms = size(w_tab, 2);
-            num_blobs = size(w_tab, 1);
+        function [geometries, offsets] = compute_proj_geom_shapefuncs_w(self, bl, ref_gr, ors_allblobs, sf_w, lims_blobs_w)
+            bls_bbuim = cat(1, bl.bbuim);
+            bls_bbvim = cat(1, bl.bbvim);
+            bls_bbsize = cat(1, bl.bbsize);
+            bb_bls_uv = [bls_bbuim(:, 1), bls_bbvim(:, 1), bls_bbsize(:, [1 2])];
 
-            slice_pos_in_blob = w_tab - repmat(extreemes_blobs_w(:, 1), [1 num_geoms 2]) + 1;
+            num_blobs = numel(bl);
+            num_orients = numel(ors_allblobs);
 
-            % Assign offsets
-            offsets = cell(num_geoms, 1);
-            for ii = 1:num_geoms
-                sino_offsets = cell(1, num_blobs);
-                blob_offsets = cell(1, num_blobs);
-                proj_offsets = cell(1, num_blobs);
-                proj_coeffs = cell(1, num_blobs);
+            rot_comp_w = gtMathsRotationMatrixComp(self.parameters.labgeo.rotdir', 'col');
 
-                for n = 1:num_blobs
-                    shape_func = shape_funcs{ii}(n);
-                    ones_num_slices = ones(1, shape_func.bbsize);
+            om_step = gtGetOmegaStepDeg(self.parameters);
 
-                    slices = slice_pos_in_blob(n, ii, 1):slice_pos_in_blob(n, ii, 2);
+            geometries = cell(num_orients, 1);
+            offsets = cell(num_orients, 1);
+            for ii_o = 1:num_orients
+                sf = sf_w{ii_o};
 
-                    sino_offsets{n} = n * ones_num_slices;
-                    blob_offsets{n} = n * ones_num_slices;
-                    proj_offsets{n} = slices;
-                    proj_coeffs{n} = reshape(shape_func.intm, 1, []);
+                % Assign offsets
+                blob_offsets = cell(1, num_blobs);
+                ws = cell(1, num_blobs);
+                for ii_b = 1:num_blobs
+                    ws{ii_b} = sf(ii_b).bbwim(1):sf(ii_b).bbwim(2);
+                    blob_offsets{ii_b} = ii_b * ones(1, sf(ii_b).bbsize);
                 end
-
-                offsets{ii}.('sino_offsets') = [sino_offsets{:}];
-                offsets{ii}.('blob_offsets') = [blob_offsets{:}];
-                offsets{ii}.('proj_offsets') = [proj_offsets{:}];
-                offsets{ii}.('proj_coeffs')  = [proj_coeffs{:}];
-
-                offsets{ii}.proj = struct( ...
+                ws = [ws{:}];
+                sino_offsets = 1:numel(ws);
+                blob_offsets = [blob_offsets{:}];
+                proj_offsets = ws - lims_blobs_w(blob_offsets, 1)' + 1;
+                proj_coeffs = cat(1, sf.intm)';
+
+                offsets{ii_o}.('sino_offsets') = sino_offsets;
+                offsets{ii_o}.('blob_offsets') = blob_offsets;
+                offsets{ii_o}.('proj_offsets') = proj_offsets;
+                offsets{ii_o}.('proj_coeffs')  = proj_coeffs;
+
+                offsets{ii_o}.proj = struct( ...
                     'sino_offset', sino_offsets, 'blob_offset', blob_offsets, ...
                     'proj_offsets', proj_offsets, 'proj_coeffs', proj_coeffs );
+
+                % Computing geometries
+                rot_w_l2s = gtMathsRotationTensor(ws' .* om_step, rot_comp_w);
+
+                dveclab = ors_allblobs(ii_o).dvec(blob_offsets, :);
+                dvecsam = gtGeoLab2Sam(dveclab, rot_w_l2s, ...
+                    self.parameters.labgeo, self.parameters.samgeo, true);
+
+                proj_geom = gtGeoProjForReconstruction(...
+                    dvecsam, rot_w_l2s, ref_gr.center, ...
+                    bb_bls_uv(blob_offsets, :), [], ...
+                    self.parameters.detgeo(self.det_index), ...
+                    self.parameters.labgeo, ...
+                    self.parameters.samgeo, ...
+                    self.parameters.recgeo(self.det_index), ...
+                    'ASTRA_grain');
+
+                if (self.volume_downscaling > 1)
+                    proj_geom(:, 4:12) = proj_geom(:, 4:12) / self.volume_downscaling;
+                end
+
+                geometries{ii_o} = double(proj_geom);
             end
         end
 
-        function offsets = make_proj_geom_shape_funcs_uvw(self, shape_funcs, extreemes_blobs_w, w_tab)
+        function offsets = make_proj_geom_shape_funcs_uvw(self, shape_funcs, lims_blobs_w, w_tab)
             num_geoms = size(w_tab, 2);
             num_blobs = size(w_tab, 1);
 
-            slice_pos_in_blob = w_tab - repmat(extreemes_blobs_w(:, 1), [1 num_geoms 2]) + 1;
+            slice_pos_in_blob = w_tab - repmat(lims_blobs_w(:, 1), [1 num_geoms 2]) + 1;
 
             % Assign offsets
             offsets = cell(num_geoms, 1);