Shape-functions: Added function to create UVW shape functions, by means of...

Shape-functions: Added function to create UVW shape functions, by means of fwd-project a selection of points in orientation-space

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

zUtil_Deformation/gtDefCreateShapeFuncs.m

+function shape_funcs = gtDefCreateShapeFuncs(sampler, factor, dtype, recenter_sf)
+
+    if (~exist('dtype', 'var') || isempty(dtype))
+        dtype = 'single';
+    end
+
+    if (~exist('recenter_sf', 'var') || isempty(recenter_sf))
+        recenter_sf = true;
+    end
+
+    interpolated_voxel = false;
+
+    det_ind = sampler.detector_index;
+
+    labgeo = sampler.parameters.labgeo;
+%     samgeo = sampler.parameters.samgeo;
+%     recgeo = sampler.parameters.recgeo(det_ind);
+%     detgeo = sampler.parameters.detgeo(det_ind);
+
+    labgeo.rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
+    sampler.parameters.labgeo = labgeo;
+
+    fed_pars_detector = struct(...
+        'blobsizeadd', [0 0 0], 'psf', [], 'apply_uv_shift', recenter_sf);
+    fedpars = struct( ...
+        'dcomps', [1 1 1 0 0 0 0 0 0] == 1, ...
+        'defmethod', 'rod_rightstretch', ...
+        'detector', repmat( fed_pars_detector, [numel(sampler.parameters.detgeo), 1]));
+
+    voxel_factor = factor([1 1 1]);
+    num_sub_voxels = prod(voxel_factor);
+    ones_nsv = ones(num_sub_voxels, 1);
+
+    gv = struct('cs', [], 'ind', [], 'dm', [], 'pow', []);
+
+    [gv.ind(:, 1), gv.ind(:, 2), gv.ind(:, 3)] = ind2sub(voxel_factor, 1:num_sub_voxels);
+    gv.ind = gv.ind';
+
+    ref_gr = sampler.get_reference_grain();
+    gv.cs = ref_gr.center(ones_nsv, :)';
+
+    % computing the subsampling, used to deterine the shape functions
+    voxel_size = sampler.stats.sampling.gaps;
+    if (interpolated_voxel)
+        half_voxel_size = voxel_size; % Enlarging to neighbours
+        steps = voxel_size ./ voxel_factor * 2;
+    else
+        half_voxel_size = voxel_size ./ 2;
+        steps = voxel_size ./ voxel_factor;
+    end
+    half_steps = steps ./ 2;
+    beg_pos = - half_voxel_size + half_steps;
+    end_pos = + half_voxel_size - half_steps;
+
+    x_steps = beg_pos(1):steps(1):end_pos(1);
+    y_steps = beg_pos(2):steps(2):end_pos(2);
+    z_steps = beg_pos(3):steps(3):end_pos(3);
+
+    num_steps = [numel(x_steps), numel(y_steps), numel(z_steps)];
+    ones_x_steps = ones(1, num_steps(1));
+
+    gv.d = zeros(3, prod(num_steps));
+
+    pos = 0;
+    for ii_z = z_steps
+        for ii_y = y_steps
+            gv.d(:, (pos+1):(pos+num_steps(1))) ...
+                = [ x_steps; ones_x_steps * ii_y; ones_x_steps * ii_z];
+
+            pos = pos + num_steps(1);
+        end
+    end
+
+    if (interpolated_voxel)
+        gv.pow = prod(1 - abs(gv.d ./ half_voxel_size(ones_nsv, :)'), 1);
+    else
+        gv.pow = ones(1, num_sub_voxels, dtype);
+    end
+    gv.pow = gv.pow / sum(gv.pow);
+
+    gv.used_ind = find(gv.pow > 0);
+    nuv = numel(gv.used_ind);
+    ones_nuv = ones(1, nuv);
+
+    sel_bls = find(sampler.selected);
+    inc_bls = sampler.ondet(sampler.included(sampler.selected));
+
+    % computing the size of the blobs
+    shape_funcs_size(:, 3) = ceil(2 * atand(sampler.stats.sampling.gaps(3)) ...
+        / labgeo.omstep .* ref_gr.allblobs.lorentzfac(inc_bls)) + 1;
+    shape_funcs_size(:, 1) = sampler.stats.proj.max_pixel_dist_u;
+    shape_funcs_size(:, 2) = sampler.stats.proj.max_pixel_dist_v;
+    if (interpolated_voxel)
+        shape_funcs_size = shape_funcs_size * 2;
+    end
+    shape_funcs_size = 2 * ceil(shape_funcs_size) + 1; % Imposing odd shape! -> for centering
+    half_sf_size = floor(shape_funcs_size / 2);
+
+    % Only one lattice allowed for the moment!
+    num_ors = numel(sampler.lattice(1).gr);
+    shape_funcs = cell(num_ors, 1);
+
+    nbl = numel(sel_bls);
+
+    fprintf('Computing UVW shape functions: ')
+    c = tic();
+    for ii_g = 1:num_ors
+        num_chars = fprintf('%03d/%03d', ii_g, num_ors);
+
+        or = sampler.lattice(1).gr{ii_g};
+
+        bl = repmat(gtFedBlobDefinition(), nbl, 1);
+
+        or_uvw = or.allblobs.detector(det_ind).uvw(sel_bls, :);
+
+        for ii_b = 1:nbl
+            % Load diffraction paramaters of blobs
+
+            % Initial strained non-signed plane normal in SAMPLE reference
+            bl(ii_b).plorig = or.allblobs.plorig(sel_bls(ii_b), :);
+            bl(ii_b).sinth0 = or.allblobs.sintheta(sel_bls(ii_b));  % initial sin(theta)
+            bl(ii_b).omega0 = or.allblobs.omega(sel_bls(ii_b));     % initial omega
+            bl(ii_b).omind  = or.allblobs.omind(sel_bls(ii_b));     % initial omega index (1..4)
+
+            % Initial centroid in image
+            bl(ii_b).u0im = or_uvw(ii_b, :);
+
+            % Projection information!!
+            bl(ii_b).bbuim = round(or_uvw(ii_b, 1) + [-half_sf_size(ii_b, 1), half_sf_size(ii_b, 1)]);
+            bl(ii_b).bbvim = round(or_uvw(ii_b, 2) + [-half_sf_size(ii_b, 2), half_sf_size(ii_b, 2)]);
+            bl(ii_b).bbwim = round(or_uvw(ii_b, 3) + [-half_sf_size(ii_b, 3), half_sf_size(ii_b, 3)]);
+            bl(ii_b).bbsize = [ ...
+                (bl(ii_b).bbuim(2) - bl(ii_b).bbuim(1) + 1), ...
+                (bl(ii_b).bbvim(2) - bl(ii_b).bbvim(1) + 1), ...
+                (bl(ii_b).bbwim(2) - bl(ii_b).bbwim(1) + 1) ];
+
+            % Blob origin and center [u, v, w]:
+            bl(ii_b).bbor   = [bl(ii_b).bbuim(1)-1, bl(ii_b).bbvim(1)-1, bl(ii_b).bbwim(1)-1];
+
+            bl(ii_b).bbcent = [ ...
+                round( (bl(ii_b).bbuim(1) + bl(ii_b).bbuim(2) ) / 2 ), ...
+                round( (bl(ii_b).bbvim(1) + bl(ii_b).bbvim(2) ) / 2 ), ...
+                round( (bl(ii_b).bbwim(1) + bl(ii_b).bbwim(2) ) / 2 ) ];
+
+            % u in blob = u in image minus u at blob origin
+            bl(ii_b).u0bl = bl(ii_b).u0im - bl(ii_b).bbor;
+            bl(ii_b).uv_shift = bl(ii_b).u0im(1:2) - bl(ii_b).bbcent(1:2);
+
+            bl(ii_b).int = zeros(bl(ii_b).bbsize, dtype);
+            bl(ii_b).intm = [];
+        end
+
+        for ii_b = nbl:-1:1
+            gvb(ii_b).uc      = zeros(3, nuv, dtype);
+            gvb(ii_b).ucd     = zeros(3, nuv, dtype);
+            gvb(ii_b).ucblmod = zeros(3, nuv, dtype);
+            gvb(ii_b).ucbl8   = zeros(1, nuv*8, dtype);
+            gvb(ii_b).ucbl8ini= zeros(1, nuv*8, dtype);
+            gvb(ii_b).ucbl8ok = 0;
+
+            % Initial u of gv centre relative to blob origin
+            gvb(ii_b).uc0bl = cast(bl(ii_b).u0im - bl(ii_b).bbor, dtype)';
+            gvb(ii_b).uc0bl = gvb(ii_b).uc0bl(:, ones_nuv);
+        end
+
+        [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, sampler.parameters, det_ind, false);
+
+        % Cutting borders
+        for ii_b = 1:nbl
+            diff_inds_w = [...
+                find(sum(sum(bl(ii_b).int, 1), 2), 1, 'first') - 1, ...
+                find(sum(sum(bl(ii_b).int, 1), 2), 1, 'last') - bl(ii_b).bbsize(3); ...
+                ];
+            bl(ii_b).bbwim = bl(ii_b).bbwim + diff_inds_w;
+
+            bl(ii_b).mbbu = bl(ii_b).bbuim;
+            bl(ii_b).mbbv = bl(ii_b).bbvim;
+            bl(ii_b).mbbw = bl(ii_b).bbwim;
+
+            inds_w = [1, bl(ii_b).bbsize(3)] + diff_inds_w;
+
+            bl(ii_b).intm = bl(ii_b).int(:, :, inds_w(1):inds_w(2));
+            bl(ii_b).int = [];
+
+            bl(ii_b).bbsize(3) = size(bl(ii_b).intm, 3);
+        end
+
+        shape_funcs{ii_g} = bl;
+
+        fprintf(repmat('\b', [1, num_chars]));
+    end
+
+    fprintf('Done in %g seconds.\n', toc(c))
+end