6D-Reconstruction: improved blobs/projections renormalization

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

5_reconstruction/gtReconstruct6DGetParamenters.m

@@ -45,10 +45,6 @@ function rec_opts = gtReconstruct6DGetParamenters(parameters)
             || isempty(rec_opts.rspace_oversize))
         rec_opts.rspace_oversize = 1.2;
     end
-    if (~isfield(rec_opts, 'use_predicted_scatter_ints') ...
-            || isempty(rec_opts.use_predicted_scatter_ints))
-        rec_opts.use_predicted_scatter_ints = false;
-    end
     if (~isfield(rec_opts, 'shape_functions_type') ...
             || isempty(rec_opts.shape_functions_type))
         rec_opts.shape_functions_type = 'none';
@@ -65,4 +61,12 @@ function rec_opts = gtReconstruct6DGetParamenters(parameters)
             || isempty(rec_opts.tv_strategy))
         rec_opts.tv_strategy = 'groups';
     end
+    if (~isfield(rec_opts, 'use_predicted_scatter_ints') ...
+            || isempty(rec_opts.use_predicted_scatter_ints))
+        rec_opts.use_predicted_scatter_ints = false;
+    end
+    if (~isfield(rec_opts, 'use_matrix_row_rescaling') ...
+            || isempty(rec_opts.use_matrix_row_rescaling))
+        rec_opts.use_matrix_row_rescaling = false;
+    end
 end

5_reconstruction/gtReconstruct6DLaunchAlgorithm.m

@@ -10,6 +10,7 @@ function algo = gtReconstruct6DLaunchAlgorithm(sampler, rec_opts, parameters, va
         'rspace_super_sampling', rec_opts.rspace_super_sampling, ...
         'rspace_oversize', rec_opts.rspace_oversize, ...
         'use_predicted_scatter_ints', rec_opts.use_predicted_scatter_ints, ...
+        'use_matrix_row_rescaling', rec_opts.use_matrix_row_rescaling, ...
         'shape_functions_type', rec_opts.shape_functions_type, ...
         'det_index', conf.det_index );
 

zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m

@@ -9,6 +9,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
         rspace_oversize = 1.2;
 
         use_predicted_scatter_ints = false;
+        use_matrix_row_rescaling = false;
 
         % Different types of shape functions:
         %   'none'   : no shape functions are used
@@ -427,20 +428,43 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             [blobs_w_lims_pad, projs_w_lims, paddings] = ...
                 self.get_blob_lims(w_tab_lims, bl, ref_ws, blobs_w_interp);
 
-            blobs = self.pad_blobs(bl, paddings, blobs_w_interp);
+            fprintf('Dealing with blobs/scattering intensities:\n')
             if (self.use_predicted_scatter_ints)
-                blobs = self.renormalize_blobs(sampler, blobs);
+                fprintf(' - Computing families'' scattering intensities\n')
+                scatter_ints = self.get_scattering_intensities(sampler);
+            else
+                fprintf(' - Using blobs'' intensities\n')
+                scatter_ints = [bl(:).intensity]';
+            end
+
+            blobs = self.pad_blobs(bl, paddings, blobs_w_interp);
+            if (self.use_matrix_row_rescaling)
+                renorm_int = [bl(:).intensity]';
+            else
+                % If both 'use_matrix_row_rescaling' and
+                % 'use_predicted_scatter_ints' are false, this will be a
+                % list of ones
+                renorm_int = [bl(:).intensity]' ./ scatter_ints;
+            end
+            % They were renormalized to 1 in pad_blobs
+            for ii_b = 1:numel(blobs)
+                blobs{ii_b} = blobs{ii_b} .* renorm_int(ii_b);
             end
 
             % Let's restrict blobs to the reached regions
+            fprintf('Blob (n.) | Blob w lims            | Projs w lims | Padding | n-interp | Size | Blob Int | Scatter Int\n')
+            fprintf('----------+------------------------+--------------+---------+----------+------+----------+-------------\n')
             for ii_b = 1:tot_blobs
-                fprintf('Blob: %03d (n. %03d) --> Blob w lims %4d:%4d -> %4d:%4d, Projs w lims %4d:%4d, Paddings %2d:%2d, Size 1:%2d, Numinterp %2d -> Intensity: %3.3f\n', ...
+                fprintf('%03d (%03d) | %4d:%4d -> %4d:%4d |    %4d:%4d |   %2d:%2d |       %2d | 1:%2d | %8s | %s\n', ...
                     ii_b, sel_incl_indx(ii_b), ...
                     bl(ii_b).bbwim, blobs_w_lims_pad(ii_b, :),...
                     round(projs_w_lims(ii_b, :)), ...
-                    paddings(ii_b, :), size(blobs{ii_b}, 2), ...
-                    blobs_w_interp(ii_b), gtMathsSumNDVol(blobs{ii_b}));
+                    paddings(ii_b, :), blobs_w_interp(ii_b), ...
+                    size(blobs{ii_b}, 2), ...
+                    sprintf('%3.3g', gtMathsSumNDVol(blobs{ii_b})), ...
+                    sprintf('%3.3g', scatter_ints(ii_b)) );
             end
+            fprintf('----------+------------------------+--------------+---------+----------+------+----------+-------------\n')
 
             fprintf('Computing ASTRA proj geometry and blobs <-> sinograms coefficients: ')
             c = tic();
@@ -464,6 +488,18 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             end
             fprintf('Done (%f s).\n', toc(c));
 
+            fprintf('Rescaling projection coefficients, based on scattering intensity: ')
+            c = tic();
+            if (self.use_matrix_row_rescaling)
+                for ii_o = 1:numel(offsets)
+                    for ii_b = 1:numel(scatter_ints)
+                        inds = offsets{ii_o}.blob_offsets == ii_b;
+                        offsets{ii_o}.proj_coeffs(inds) = offsets{ii_o}.proj_coeffs(inds) .* scatter_ints(ii_b);
+                    end
+                end
+            end
+            fprintf('Done (%f s).\n', toc(c));
+
             % Taking care of the orientation-space super-sampling
             tot_orient = sampler.get_total_num_orientations();
             if (num_ospace_oversampling == 1)
@@ -562,9 +598,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                 num_chars = fprintf('%03d/%03d', ii, num_blobs);
 
                 blob = blobs(ii).intm;
-                if (self.use_predicted_scatter_ints)
-                    blob = blob ./ sum(blob(:)) .* blobs(ii).intensity;
-                end
+                blob = blob ./ gtMathsSumNDVol(blob);
                 blob(blob < 0) = 0;
 
                 pad_interp_blobs{ii} = zeros(blobs_sizes(ii, :), self.data_type);
@@ -909,21 +943,47 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             end
         end
 
-        function blobs = renormalize_blobs(self, sampler, blobs)
+        function scatter_ints = get_scattering_intensities(self, sampler)
             ref_gr = sampler.get_reference_grain();
+
+            sel = sampler.ondet(sampler.included(sampler.selected));
+
             if (isfield(self.parameters.cryst(ref_gr.phaseid), 'int_exp'))
+                fprintf(' - Using experimental intensities for each family (from: parameters.cryst.int_exp)\n')
                 cryst_scatter_ints = self.parameters.cryst(ref_gr.phaseid).int_exp;
             else
+                fprintf(' - Using theoretically predicted intensities for each family\n')
                 cryst_scatter_ints = self.parameters.cryst(ref_gr.phaseid).int;
             end
+            cryst_scatter_ints = cryst_scatter_ints ./ mean(cryst_scatter_ints);
 
-            sel = sampler.ondet(sampler.included(sampler.selected));
             thetatypes = ref_gr.allblobs.thetatype(sel);
             scatter_ints = cryst_scatter_ints(thetatypes);
+            scatter_ints = reshape(scatter_ints, [], 1);
+
             lorentz_fact = ref_gr.allblobs.lorentzfac(sel);
-            for ii_b = 1:numel(blobs)
-                blobs{ii_b} = blobs{ii_b} ./ scatter_ints(ii_b) ./ lorentz_fact(ii_b);
+
+            if (self.use_predicted_scatter_ints > 1)
+                fprintf(' - Applying beam corrections\n')
+                try
+                    beam_intensities = gtGrainComputeIncomingBeamIntensity(ref_gr, self.parameters, self.det_index);
+                    beam_ints = beam_intensities(sampler.selected);
+                    beam_ints = beam_ints / mean(beam_ints);
+
+                    sample = GtSample.loadFromFile();
+                    vol_abs = load(sample.absVolFile);
+                    atts = gtGrainComputeBeamAttenuation(ref_gr, self.parameters, self.det_index, vol_abs.abs);
+
+                    beam_ints = beam_ints .* atts(sampler.selected);
+                catch mexc
+                    gtPrintException(mexc, 'Reverting to non-corrected beam intensities')
+                    beam_ints = 1;
+                end
+            else
+                beam_ints = 1;
             end
+
+            scatter_ints = scatter_ints .* lorentz_fact .* beam_ints;
         end
     end