Fixed and improved projection correction computation

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

zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m

@@ -499,7 +499,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                 fprintf('   * Using blobs'' intensities\n')
                 scatter_ints_avg = [bl(:).intensity]';
                 scatter_ints_avg = scatter_ints_avg ./ mean([bl(:).intensity]);
-                scatter_ints_ors = [scatter_ints_avg, scatter_ints_avg];
+                scatter_ints_ors = scatter_ints_avg(:, ones(1, numel(ors_allblobs)));
             end
 
             fprintf(' - Dealing with paddings:\n')
@@ -1311,25 +1311,33 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                 try
                     fprintf(' - Applying incoming beam intensity corrections\n')
                     [beam_intensities, refs] = gtGrainComputeIncomingBeamIntensity(ref_or, self.parameters, det_ind);
-                    beam_ints_avg = beam_intensities(sel);
-                    beam_ints_avg = beam_ints_avg / mean(beam_ints_avg);
+                    beam_in_ints_avg = beam_intensities(sel);
+                    beam_in_ints_avg = beam_in_ints_avg / mean(beam_in_ints_avg);
+                catch mexc
+                    gtPrintException(mexc, 'Reverting to non-corrected beam intensities')
+                    beam_in_ints_avg = 1;
+                end
 
+                try
                     fprintf(' - Applying beam attenuation corrections\n')
                     [atts_tot, ~, abs_vol] = gtGrainComputeBeamAttenuation(ref_or, self.parameters, det_ind);
 
-                    beam_ints_avg = beam_ints_avg .* atts_tot(sel);
+                    beam_out_ints_avg = atts_tot(sel);
                 catch mexc
-                    gtPrintException(mexc, 'Reverting to non-corrected beam intensities')
-                    beam_ints_avg = 1;
+                    gtPrintException(mexc, 'Reverting to non-corrected beam attenuations')
+                    beam_out_ints_avg = 1;
                 end
             else
-                beam_ints_avg = 1;
+                beam_in_ints_avg = 1;
+                beam_out_ints_avg = 1;
             end
-            scatter_ints_avg = scatter_ints .* ref_or.allblobs.lorentzfac(ab_sel) .* beam_ints_avg;
+
+            scatter_ints_avg = scatter_ints .* ref_or.allblobs.lorentzfac(ab_sel) .* beam_in_ints_avg .* beam_out_ints_avg;
 
             fprintf(' - Computing each orientation''s corrections..')
             c = tic();
-            if (self.use_predicted_scatter_ints > 2 && numel(beam_ints_avg) > 1)
+            if (self.use_predicted_scatter_ints > 2 ...
+                    && (numel(beam_in_ints_avg) > 1 || numel(beam_out_ints_avg) > 1))
                 num_ors = numel(ors_allbls);
                 scatter_ints_ors = scatter_ints(:, ones(num_ors, 1));
                 fprintf('\b\b: ')
@@ -1337,20 +1345,35 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                     if (mod(ii_o, 10) == 1)
                         num_chars = fprintf('%05d/%05d', ii_o, num_ors);
                     end
-                    try
-                        [beam_int_or, refs] = gtGrainComputeIncomingBeamIntensity(ref_or, self.parameters, det_ind, refs);
-                        beam_int_or = beam_int_or(sel);
-                        beam_int_or = beam_int_or / mean(beam_int_or);
+                    or = ref_or;
+                    or.allblobs = ors_allbls(ii_o);
+                    if (numel(beam_in_ints_avg) > 1)
+                        try
+                            [beam_in_int_or, refs] = gtGrainComputeIncomingBeamIntensity(or, self.parameters, det_ind, refs);
+                            beam_in_int_or = beam_in_int_or(sel);
+                            beam_in_int_or = beam_in_int_or / mean(beam_in_int_or);
+                        catch mexc
+                            gtPrintException(mexc, 'Reverting to non-corrected beam intensities')
+                            beam_in_int_or = 1;
+                        end
+                    else
+                        beam_in_int_or = 1;
+                    end
 
-                        atts_tot = gtGrainComputeBeamAttenuation(ref_or, self.parameters, det_ind, abs_vol);
+                    if (numel(beam_out_ints_avg) > 1)
+                        try
+                            atts_tot = gtGrainComputeBeamAttenuation(or, self.parameters, det_ind, abs_vol);
 
-                        beam_int_or = beam_int_or .* atts_tot(sel);
-                    catch mexc
-                        gtPrintException(mexc, 'Reverting to non-corrected beam intensities')
-                        beam_int_or = 1;
+                            beam_out_int_or = atts_tot(sel);
+                        catch mexc
+                            gtPrintException(mexc, 'Reverting to non-corrected beam intensities')
+                            beam_out_int_or = 1;
+                        end
+                    else
+                        beam_out_int_or = 1;
                     end
 
-                    scatter_ints_ors(:, ii_o) = scatter_ints_ors(:, ii_o) .* ors_allbls(ii_o).lorentzfac(ab_sel, :) .* beam_int_or;
+                    scatter_ints_ors(:, ii_o) = scatter_ints_ors(:, ii_o) .* ors_allbls(ii_o).lorentzfac(ab_sel, :) .* beam_in_int_or .* beam_out_int_or;
 
                     if (mod(ii_o, 10) == 1)
                         fprintf(repmat('\b', [1 num_chars]));
@@ -1358,7 +1381,8 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                 end
             else
                 lfacs = cat(2, ors_allbls(:).lorentzfac);
-                scatter_ints_ors = bsxfun(@times, scatter_ints .* beam_ints_avg, lfacs(ab_sel, :));
+                scatter_ints_ors = scatter_ints .* beam_in_ints_avg .* beam_out_ints_avg;
+                scatter_ints_ors = bsxfun(@times, scatter_ints_ors, lfacs(ab_sel, :));
             end
             fprintf('\b\b: Done in %g seconds.\n', toc(c))
         end