6D-Reconstruction: added Lorentz factor correction for slice interpolation

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

5_reconstruction/gtRecGrainsDefaultParameters.m

@@ -13,6 +13,6 @@ function par_rec = gtRecGrainsDefaultParameters(algo)
         case {'6DL1', '6DTV', '6DTVL1'}
             par_rec = struct(...
                 'algorithm', upper(algo), 'num_iter', 50, 'list', [], ...
-                'options', struct('grid_edge', 7, 'super_sampling', 1, 'num_interp', 4, 'lambda', 1e-3) );
+                'options', struct('grid_edge', 7, 'super_sampling', 1, 'num_interp', 1, 'lambda', 1e-3) );
     end
 end

5_reconstruction/gtReconstructGrainOrientation.m

@@ -71,11 +71,8 @@ function [vols, sampling, gr, algo, sampler] = gtReconstructGrainOrientation(gra
     end
 
     rec_factory = Gt6DReconstructionAlgorithmFactory();
-    if (num_interp > 1)
-        [algo, good_gr] = rec_factory.getSubBlobReconstructionAlgo(gr.bl, gr.selected, sampler, '3D', num_interp);
-    else
-        [algo, good_gr] = rec_factory.getBlobReconstructionAlgo(gr.bl, gr.selected, sampler, '3D');
-    end
+    [algo, good_gr] = rec_factory.getSubBlobReconstructionAlgo(gr.bl, gr.selected, sampler, '3D', num_interp);
+
     switch (upper(algorithm))
         case '6DL1'
             algo.cp(num_iter, lambda);

zUtil_Deformation/Gt6DReconstructionAlgorithmFactory.m

@@ -118,6 +118,9 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             fprintf('\b\b (%f s).\n', toc(c));
             tot_blobs = numel(bl);
 
+            ref_gr = sampler.get_reference_grain();
+            etas = ref_gr.allblobs.eta(ref_gr.ondet(ref_gr.included(ref_gr.selected)));
+
             [grains, grains_ss] = sampler.get_orientations();
             grains = [grains{:}];
             % Array of structures that contain all the info relative to
@@ -126,11 +129,11 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
 
             omega_int_tab = Gt6DReconstructionAlgorithmFactory.getOmegaIntegerTab(grains_props, sampler.parameters);
 
-            [sub_blob_slices, proj_coeffs] = self.getSubBlobSlices(bl, num_interp);
+            [sub_blob_slices, proj_coeffs, slices_interp] = self.getSubBlobSlices(bl, num_interp, etas);
             blobs = cell(size(sub_blob_slices));
 
             [slice_lims, lims, abs_lims, ~, extremesOfBlobs, extremesOmIntTab] = ...
-                Gt6DReconstructionAlgorithmFactory.getSubBlobLims(omega_int_tab, bl, num_interp);
+                Gt6DReconstructionAlgorithmFactory.getSubBlobLims(omega_int_tab, bl, slices_interp);
 
             [overflow_b, overflow_o] = Gt6DReconstructionAlgorithmFactory.compute_overflow(omega_int_tab, extremesOfBlobs);
 
@@ -146,7 +149,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                 grains_props_ss = Gt6DReconstructionAlgorithmFactory.get_selected_props(grains_ss, selected_bls);
                 omega_int_tab_ss = Gt6DReconstructionAlgorithmFactory.getOmegaIntegerTab(grains_props_ss, sampler.parameters);
                 [slice_lims_ss, lims_ss, abs_lims_ss, ~, extremesOfBlobs_ss, extremesOmIntTab_ss] ...
-                    = Gt6DReconstructionAlgorithmFactory.getSubBlobLims(omega_int_tab_ss, bl, num_interp);
+                    = Gt6DReconstructionAlgorithmFactory.getSubBlobLims(omega_int_tab_ss, bl, slices_interp);
 
                 [overflow_ss_b, overflow_ss_o] = Gt6DReconstructionAlgorithmFactory.compute_overflow(omega_int_tab_ss, extremesOfBlobs_ss);
 
@@ -183,7 +186,8 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                 if (slice_lims_tot(ii, 2) > size(sub_blob_slices{ii}, 3))
                     depth_blob = size(bl(ii).intm, 3);
                     fprintf('\n AHIA! numInterp: %d, lims: %d:%d, blob size: 1:%d, virtual slices 1:%d\n\n', ...
-                        num_interp, lims_tot(ii, :), depth_blob, depth_blob + mod(num_interp - mod(depth_blob, num_interp), num_interp) + 1)
+                        slices_interp(ii), lims_tot(ii, :), depth_blob, ...
+                        depth_blob + mod(slices_interp(ii) - mod(depth_blob, slices_interp(ii)), slices_interp(ii)) + 1)
                 end
 
                 fprintf(' - ExremesBlobs %4d:%4d, ExtremesProjs %4d:%4d, Lims %3d:%3d, SliceLims %d:%d, Size 1:%d, ', ...
@@ -226,67 +230,67 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
     end
 
     methods (Access = protected)
-        function [subBlobSlices, projCoeffs] = getSubBlobSlices(self, blobs, numInterp)
-            subBlobSlices = cell(size(blobs));
-            subBlobCoeffs = cell(size(blobs));
-            projCoeffs = cell(size(blobs));
-
-            blocky_groups = false;
+        function [sub_blob_slices, proj_coeffs, slices_interp] = getSubBlobSlices(self, blobs, num_interp, etas)
+            num_blobs = numel(blobs);
+            sub_blob_slices = cell(num_blobs, 1);
+            sub_blob_coeffs = cell(num_blobs, 1);
+            proj_coeffs = cell(num_blobs, 1);
+
+            slices_interp = num_interp(ones(num_blobs, 1));
+            if (exist('etas', 'var'))
+                lorentz_factor = 1 ./ abs(sind(etas));
+                % Let's try to not penalize rounding errors and being a bit
+                % tollerant with small deviations from integer numbers
+                % (rounding in [-0.75 0.25), instead of [-0.5 0.5) ), but
+                % still trying to avoid gaps in the projections
+                slices_interp = slices_interp .* ceil(lorentz_factor - 0.25);
+%                 slices_interp = slices_interp .* round(lorentz_factor);
+            end
 
             fprintf('Creating sub-blob volumes: ')
             c = tic();
-            for ii = 1:numel(subBlobSlices)
-                num_chars = fprintf('%03d/%03d', ii, numel(subBlobSlices));
+            for ii = 1:numel(sub_blob_slices)
+                num_chars = fprintf('%03d/%03d', ii, numel(sub_blob_slices));
+
+                slice_interp = slices_interp(ii);
 
-                [blobSize(1), blobSize(2), blobSize(3)] = size(blobs(ii).intm);
+                [blob_size(1), blob_size(2), blob_size(3)] = size(blobs(ii).intm);
                 % Number of interpolated slices should be at least enough
                 % to accomodate all the slices of the blob
-                numSlices = ceil((blobSize(3) -1) / numInterp) +1;
+                num_slices = ceil((blob_size(3) -1) / slice_interp) +1;
 
-                covered_slices_size = (numSlices-1) * numInterp + 1;
+                covered_slices_size = (num_slices-1) * slice_interp + 1;
 
 %                 num_chars = num_chars + fprintf(' Blob size %d, numSlices %d, convertedSliceSize %d', ...
 %                     blobSize(3), numSlices, covered_slices_size);
 
-                subBlobSlices{ii} = zeros(blobSize(1), blobSize(2), numSlices, self.data_type);
-                for n = numSlices:-1:1
-                    coeffs = abs( (numInterp * (n-1) +1) - (1:covered_slices_size));
-                    coeffs = 1 - coeffs / numInterp;
+                sub_blob_slices{ii} = zeros(blob_size(1), blob_size(2), num_slices, self.data_type);
+                for n = num_slices:-1:1
+                    coeffs = abs( (slice_interp * (n-1) +1) - (1:covered_slices_size));
+                    coeffs = 1 - coeffs / slice_interp;
 
                     % these will be the blob slices that relate to the
                     % given ii interpolated slice
-                    if (blocky_groups)
-                        goodCoeffs = find(coeffs >= 0.5);
-                    else
-                        goodCoeffs = find(coeffs > 0);
-                    end
+                    goodCoeffs = find(coeffs > 0);
                     % The corresponding coefficients
                     coeffs = coeffs(goodCoeffs);
 
-                    if (blocky_groups)
-                        coeffs(1:end) = 1;
-                        if (mod(numInterp, 2) == 0)
-                            coeffs(1) = 0.5;
-                            coeffs(end) = 0.5;
-                        end
-                    end
-
-                    safe_coeffs_pos = goodCoeffs <= blobSize(3);
+                    safe_coeffs_pos = goodCoeffs <= blob_size(3);
                     safe_good_coeffs = goodCoeffs(safe_coeffs_pos);
                     safe_coeffs = coeffs(safe_coeffs_pos);
                     % Let's build the interpolated slices
                     blob = blobs(ii).intm;
                     blob = blob ./ sum(blob(:));
                     for w = 1:numel(safe_coeffs)
-                        subBlobSlices{ii}(:, :, n) = ...
-                            subBlobSlices{ii}(:, :, n) + safe_coeffs(w) .* blob(:, :, safe_good_coeffs(w));
+                        sub_blob_slices{ii}(:, :, n) = ...
+                            sub_blob_slices{ii}(:, :, n) + safe_coeffs(w) .* blob(:, :, safe_good_coeffs(w));
                     end
                     % Save the metadata
-                    subBlobCoeffs{ii}(n) = struct('indexes', goodCoeffs, 'coeffs', coeffs);
+                    sub_blob_coeffs{ii}(n) = struct('indexes', goodCoeffs, 'coeffs', coeffs);
                 end
 
                 % Perpare next metadata structure
-                projCoeffs{ii} = struct( ...
+                proj_coeffs{ii} = struct( ...
                     'slices', arrayfun(@(x){[]}, 1:covered_slices_size), ...
                     'coeffs', arrayfun(@(x){[]}, 1:covered_slices_size) );
 
@@ -301,17 +305,17 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             % Per each blob we want to easily find out what the projections
             % project to. The metadata from subBlobCoeffs is the other way
             % around.
-            for ii = 1:numel(subBlobSlices)
-                numSlices = numel(subBlobCoeffs{ii});
+            for ii = 1:numel(sub_blob_slices)
+                num_slices = numel(sub_blob_coeffs{ii});
 
-                for n = numSlices:-1:1
-                    indexes = subBlobCoeffs{ii}(n).indexes;
-                    coeffs = subBlobCoeffs{ii}(n).coeffs;
+                for n = num_slices:-1:1
+                    indexes = sub_blob_coeffs{ii}(n).indexes;
+                    coeffs = sub_blob_coeffs{ii}(n).coeffs;
                     for m = numel(indexes):-1:1
-                        projCoeffs{ii}(indexes(m)).slices = ...
-                            [ projCoeffs{ii}(indexes(m)).slices, n ];
-                        projCoeffs{ii}(indexes(m)).coeffs = ...
-                            [ projCoeffs{ii}(indexes(m)).coeffs, coeffs(m) ];
+                        proj_coeffs{ii}(indexes(m)).slices = ...
+                            [ proj_coeffs{ii}(indexes(m)).slices, n ];
+                        proj_coeffs{ii}(indexes(m)).coeffs = ...
+                            [ proj_coeffs{ii}(indexes(m)).coeffs, coeffs(m) ];
                     end
                 end
             end
@@ -401,13 +405,13 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
                 | (extremesOmIntTab(:, 2) > extremesOfBlobs(:, 2));
         end
 
-        function [sliceLims, lims, absLims, overflow, extremesOfBlobs, extremesOmIntTab] = getSubBlobLims(omegaIntegerTab, bl, numInterp)
+        function [sliceLims, lims, absLims, overflow, extremesOfBlobs, extremesOmIntTab] = getSubBlobLims(omegaIntegerTab, bl, slices_interp)
             [~, ~, ~, extremesOfBlobs, extremesOmIntTab] = Gt6DReconstructionAlgorithmFactory.getBlobLims(omegaIntegerTab, bl);
 
             blobs_sizes = extremesOfBlobs(:, 2) - extremesOfBlobs(:, 1) +1;
 
             extremesOfBlobs(:, 2) = extremesOfBlobs(:, 1) + (blobs_sizes-1) ...
-                + mod(numInterp - mod((blobs_sizes-1), numInterp), numInterp);
+                + mod(slices_interp - mod((blobs_sizes-1), slices_interp), slices_interp);
 
             % Blob limits in the full omega stack
             absLims = [ ...
@@ -416,7 +420,7 @@ classdef Gt6DReconstructionAlgorithmFactory < handle
             % blob limits in its volume size
             lims = absLims - extremesOfBlobs(:, [1 1]) + 1;
 
-            sliceLims = (lims -1) / numInterp;
+            sliceLims = (lims -1) ./ slices_interp(:, [1 1]);
             sliceLims = [floor(sliceLims(:, 1)), ceil(sliceLims(:, 2))] + 1;
 
             overflow = ...

zUtil_Deformation/GtOrientationSampling.m

@@ -292,6 +292,10 @@ classdef GtOrientationSampling < handle
             or_ss = self.gr_ss(:);
         end
 
+        function ref_gr = get_reference_grain(self)
+            ref_gr = self.ref_gr;
+        end
+
         function coeffs_ss = get_ss_coeffs(self)
             coeffs_ss = self.coeffs_gr_ss;
         end