Geo: reverted modularization of transform functions

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

7_fed2/gtFedFwdProjExact.m

@@ -20,8 +20,6 @@ function [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, parameters, det_num
     labgeo = parameters.labgeo;
     samgeo = parameters.samgeo;
 
-    sam2lab_rescaled = gtGeoSam2Lab_GetTransformRescaled(samgeo);
-
     dcomps = fedpars.dcomps;
     defmethod = fedpars.defmethod;
 
@@ -85,7 +83,7 @@ function [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, parameters, det_num
             % reference state)
             plsam = plorig(ii, :)';  % ! use plorig and not plsam to keep omega order below!
 
-            [ucbl, ~] = predict_uvw(plsam, gvcs, defT, bl(ii), labgeo, samgeo, detgeo, sam2lab_rescaled);
+            [ucbl, ~] = predict_uvw(plsam, gvcs, defT, bl(ii), labgeo, samgeo, detgeo);
 
             if (apply_uv_shifts)
                 ucbl(1:2, :) = ucbl(1:2, :) - bl(ii).uv_shift(ones(1, nv), :)';
@@ -244,7 +242,7 @@ function [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, parameters, det_num
     end
 end
 
-function [uvw_bl, rot_s2l] = predict_uvw(plsam, gvcs, defT, bl, labgeo, samgeo, detgeo, sam2lab_rescaled)
+function [uvw_bl, rot_s2l] = predict_uvw(plsam, gvcs, defT, bl, labgeo, samgeo, detgeo)
 %
     [pl_samd, drel] = gtStrainPlaneNormals(plsam, defT); % unit column vectors
 
@@ -268,7 +266,7 @@ function [uvw_bl, rot_s2l] = predict_uvw(plsam, gvcs, defT, bl, labgeo, samgeo,
     dvec_lab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');
 
     rot_s2l = permute(rot_l2s, [2 1 3]);
-    gvcs_lab = gtGeoSam2Lab_Position(gvcs', sam2lab_rescaled, rot_s2l, labgeo, samgeo);
+    gvcs_lab = gtGeoSam2Lab(gvcs', rot_s2l, labgeo, samgeo, false);
 
     uvw_bl = gtFedPredictUVWMultiple([], dvec_lab, gvcs_lab', ...
         detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...

7_fed2/gtFedFwdProjExactAlt.m

-function [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, parameters, det_num, verbose)
-% FUNCTION [gvb, bl] = gtFedFwdProjExact(gv, gvb, bl, fedpars, parameters, det_num, verbose)
-
-    if (~exist('det_num', 'var'))
-        det_num = 1;
-    end
-    if (~exist('verbose', 'var'))
-        verbose = true;
-    end
-
-    o = GtConditionalOutput(verbose);
-
-    o.fprintf('Forward projection (%s):\n', mfilename)
-
-    detgeo = parameters.detgeo(det_num);
-    labgeo = parameters.labgeo;
-    samgeo = parameters.samgeo;
-
-    nbl = length(bl);
-    uinds = gv.used_ind;
-    nv = numel(uinds);
-    gvpow = gv.pow(1, uinds);
-
-    % Rescaled sam -> lab coordinate trasform matrix
-    sam2lab = gtGeoSam2Lab_GetTransformRescaled(samgeo);
-
-    dcomps = fedpars.dcomps;
-    defmethod = fedpars.defmethod;
-
-    psf = get_conf_psf(fedpars, det_num, nbl);
-    apply_uv_shifts = get_conf_apply_uv_shifts(fedpars, det_num);
-
-    plorig = vertcat(bl(:).plorig);
-    bbor   = vertcat(bl(:).bbor);
-    bbsize = vertcat(bl(:).bbsize);
-
-    % The plane normals need to be brought in the Lab reference where the
-    % beam direction and rotation axis are defined.
-    % Use the Sample -> Lab orientation transformation assuming omega=0;
-    % (vector length preserved for free vectors)
-    plorig = gtGeoSam2Lab(plorig, eye(3), labgeo, samgeo, true, false);
-
-    num_super_sampling = size(gv.d, 3);
-
-    o.fprintf(' - Supersampling: ')
-
-    for ii_ss = 1:num_super_sampling
-        o.fprintf('%03d/%03d:\n', ii_ss, num_super_sampling);
-
-        gvd   = gv.d(:, uinds, ii_ss);
-
-        % Deformation tensors (relative to reference state) from their components
-        defT = gtFedDefTensorFromComps(gvd, dcomps, defmethod, 0);
-
-        gvcs = gv.cs(:, uinds, ii_ss);
-        uvw = cell(1, nbl);
-
-        o.fprintf('   * Computing indices and bbsizes: ')
-        t = tic();
-
-        for ii = 1:nbl
-            num_chars = o.fprintf('%03d/%03d', ii, nbl);
-            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-            % Calculate new detector coordinates
-            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-            % New deformed plane normals and relative elongations (relative to
-            % reference state)
-            plsam = plorig(ii, :)';  % ! use plorig and not plsam to keep omega order below!
-
-            uvw_bl = predict_uvw(plsam, gvcs, defT, bl(ii), labgeo, samgeo, detgeo, sam2lab);
-
-            if (apply_uv_shifts)
-                uvw_bl(1:2, :) = uvw_bl(1:2, :) - bl(ii).uv_shift(ones(1, nv), :)';
-            end
-
-            % Detector coordinates U,V in blob
-            bbort = bbor(ii, :)';
-            uvw_bl = uvw_bl - bbort(:, ones(1, nv));
-
-            % U,V shift relative to reference state
-            gvb(ii).ucd = uvw_bl - gvb(ii).uc0bl - gvb(ii).uc;
-            gvb(ii).uc  = uvw_bl - gvb(ii).uc0bl;
-
-            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-            % Sum intensities
-            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-            % Distribute value over 8 pixels
-
-            % Modulus of uvw position
-            gvb(ii).ucblmod = mod(uvw_bl, 1);
-
-            % uvw indices of 8 nearest neighbours (nx3)
-            [uvw{ii}, gvb(ii).ucbl8] = compute_projected_positions(uvw_bl, bbsize(ii, :));
-
-            o.fprintf(repmat('\b', [1, num_chars]));
-        end
-        o.fprintf('Done in %f s\n', toc(t));
-
-        verify_uvw_limits(o, uvw, bbsize, psf, fedpars, det_num);
-
-        o.fprintf('   * Projecting volumes: ')
-        t = tic();
-        for ii = 1:nbl
-            num_chars = o.fprintf('%03d/%03d', ii, nbl);
-
-            ints = compute_projected_intensities(gvb(ii), gvpow, bl(ii), nv);
-
-            % Accumulate all intensities in the blob voxel-wise
-            %   'uvw' needs to be nx3; 'ints' is now 1xnx8
-            try
-                bl(ii).int = bl(ii).int + accumarray(uvw{ii}, ints, bbsize(ii,:));
-            catch mexc
-                disp(' ')
-                disp('ERROR')
-                disp(['The error is probably caused by the projected intensities' ...
-                    ' falling outside the blob volume. Try increasing the blob' ...
-                    ' volume padding: fedpars.detector(det_ind).blobsizeadd'])
-                disp('Blob ID:')
-                disp(ii)
-                disp('Blob size:')
-                disp(bbsize(ii,:))
-                disp('Min projected U,V,W coordinates:')
-                disp(min(uvw{ii},[],1))
-                disp('Max projected U,V,W coordinates:')
-                disp(max(uvw{ii},[],1))
-                new_exc = GtFedExceptionFwdProj(mexc, det_num, minimum_blob_size_add);
-                throw(new_exc)
-            end
-
-            o.fprintf(repmat('\b', [1, num_chars]));
-        end
-    end
-
-    if (~isempty(psf))
-        o.fprintf('Done in %f s\n', toc(t));
-        o.fprintf(' - Applying PSF:')
-
-        for ii = 1:nbl
-            num_chars = o.fprintf('%03d/%03d', ii, nbl);
-
-            bl(ii).int = convn(bl(ii).int, psf{ii}, 'same');
-
-            o.fprintf(repmat('\b', [1, num_chars]));
-        end
-    end
-
-    o.fprintf('Done in %f s\n', toc(t));
-end
-
-function [uvw_bl, rot_s2l] = predict_uvw(plsam, gvcs, defT, bl, labgeo, samgeo, detgeo, sam2lab_rescaled)
-%
-    [pl_samd, drel] = gtStrainPlaneNormals(plsam, defT); % unit column vectors
-
-    % New sin(theta)
-    sinth = bl.sinth0 ./ drel;
-
-    % Predict omega angles: 4 for each plane normal
-    [om, pllab, ~, rot_l2s] = gtFedPredictOmegaMultiple(pl_samd, ...
-        sinth, labgeo.beamdir', labgeo.rotdir', labgeo.rotcomp, bl.omind);
-
-    % Delete those where no reflection occurs
-    if any(isnan(om));
-        inds_bad = find(isnan(om));
-        gvd(:, inds_bad(1:10))
-        error('gtFedFwdProjExact:bad_R_vectors', ...
-            'No diffraction from some elements after deformation (%d over %d).', ...
-            numel(inds_bad), numel(om))
-    end
-
-    % Diffraction vector
-    dvec_lab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');
-
-    rot_s2l = permute(rot_l2s, [2 1 3]);
-    gvcs_lab = gtGeoSam2Lab_Position(gvcs', sam2lab_rescaled, rot_s2l, labgeo, samgeo);
-
-    uvw_bl = gtFedPredictUVWMultiple([], dvec_lab, gvcs_lab', ...
-        detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...
-        [detgeo.detrefu, detgeo.detrefv]', om, labgeo.omstep);
-end
-
-function verify_uvw_limits(o, uvw, bbsize, psf, fedpars, det_num)
-%
-    nbl = numel(uvw);
-    o.fprintf('   * Computing max BBox size and feasibilities:\n')
-    max_bbox_proj_size = [inf, inf, inf; 0, 0, 0];
-    max_bbsize_blob = [1, 1, 1; 0, 0, 0];
-    psf_size = zeros(nbl, 2);
-    for ii = 1:nbl
-        bbsize_blob = bbsize(ii,:);
-        min_ind_blob = min(uvw{ii}, [], 1);
-        max_ind_blob = max(uvw{ii}, [], 1);
-
-        max_bbsize_blob(2, :) = max(max_bbsize_blob(2, :), bbsize_blob);
-
-        if (~isempty(psf))
-            psf_size(ii, :) = (size(psf{ii}) - 1) / 2;
-        end
-
-        max_bbox_proj_size(1, :) = min(max_bbox_proj_size(1, :), min_ind_blob);
-        max_bbox_proj_size(2, :) = max(max_bbox_proj_size(2, :), max_ind_blob);
-
-        if (any(min_ind_blob < 1) || any(max_ind_blob > bbsize_blob))
-            fprintf('     + Blob: %d, Blobsize: (%d, %d, %d), bbox projection [(%d, %d, %d), (%d, %d, %d)]\n', ...
-                ii, bbsize_blob, min_ind_blob, max_ind_blob);
-        elseif (any(min_ind_blob(1:2) < psf_size(ii, :)) || any(max_ind_blob(1:2) > (bbsize_blob(1:2) - psf_size(ii, :))))
-            fprintf('     + Blob: %d, Blobsize: (%d, %d), bbox projection [(%d, %d), (%d, %d)], psf (%d, %d)\n', ...
-                ii, bbsize_blob(1:2), min_ind_blob(1:2), max_ind_blob(1:2), psf_size(ii, :));
-        end
-    end
-
-    if (~isempty(psf))
-        max_bbox_proj_size(1, 1:2) = max_bbox_proj_size(1, 1:2) - psf_size;
-        max_bbox_proj_size(2, 1:2) = max_bbox_proj_size(2, 1:2) + psf_size;
-        o.fprintf('     Maximum BBox size: [(%d, %d, %d), (%d, %d, %d)]\n', max_bbox_proj_size(1, :), max_bbox_proj_size(2, :));
-    else
-        o.fprintf('     Maximum BBox size: [(%d, %d, %d), (%d, %d, %d)]\n', max_bbox_proj_size(1, :), max_bbox_proj_size(2, :));
-    end
-    o.fprintf('     Maximum Blob size: [(%d, %d, %d), (%d, %d, %d)]\n', max_bbsize_blob(1, :), max_bbsize_blob(2, :));
-
-    blob_size_add = get_conf_blob_size_add(fedpars, det_num);
-
-    minimum_blob_size_add = blob_size_add + max( ...
-        [max_bbsize_blob - max_bbox_proj_size; max_bbox_proj_size - max_bbsize_blob], [], 1);
-    o.fprintf('     Current Blob size add: (%d, %d, %d) -> minimum suggested: (%d, %d, %d)\n', ...
-        blob_size_add(det_num, :), minimum_blob_size_add);
-end
-
-function ints = compute_projected_intensities(gvb, gvpow, bl, nv)
-%
-    mu = gvb.ucblmod;
-    mu1 = mu(1, :);
-    mu2 = mu(2, :);
-    mu3 = mu(3, :);
-
-    mu12  = mu1 .* mu2;
-    mu13  = mu1 .* mu3;
-    mu23  = mu2 .* mu3;
-    mu123 = mu12 .* mu3;
-
-    ints = zeros(1, numel(gvpow) * 8, 'like', bl.int);
-
-    ints(     1 :   nv) = gvpow.*(-mu1 - mu2 - mu3 + mu12 + mu13 + mu23 - mu123 + 1);
-    ints(  nv+1 : 2*nv) = gvpow.*(mu3 - mu13 - mu23 + mu123);
-    ints(2*nv+1 : 3*nv) = gvpow.*(mu2 - mu12 - mu23 + mu123);
-    ints(3*nv+1 : 4*nv) = gvpow.*(mu23 - mu123);
-    ints(4*nv+1 : 5*nv) = gvpow.*(mu1 - mu12 - mu13 + mu123);
-    ints(5*nv+1 : 6*nv) = gvpow.*(mu13 - mu123);
-    ints(6*nv+1 : 7*nv) = gvpow.*(mu12 - mu123);
-    ints(7*nv+1 : 8*nv) = gvpow.*mu123;
-end
-
-function [uvw8, uvw8_ind] = compute_projected_positions(uvw, bbsize)
-%
-    ul = floor(uvw); % lower pixel indices
-    uh = ul + 1;     % higher pixel indices
-
-    % uvw indices of 8 nearest neighbours (3xnx8)
-    uvw8 = [ ul, ...
-        [ul([1,2],:); uh(3,:)], ...
-        [ul(1,:); uh(2,:); ul(3,:)], ...
-        [ul(1,:); uh([2,3],:)], ...
-        [uh(1,:); ul([2,3],:)], ...
-        [uh(1,:); ul(2,:); uh(3,:)], ...
-        [uh([1,2],:); ul(3,:)], ...
-        uh ];
-
-    if (nargout > 1)
-        % Don't use sub2ind, it's slow
-        uvw8_ind = uvw8(1, :) ...
-            + bbsize(1) .* (uvw8(2, :)-1) ...
-            + bbsize(1) .* bbsize(2) .* (uvw8(3, :)-1);
-    end
-
-    uvw8 = uvw8';
-end
-
-function blob_size_add = get_conf_blob_size_add(fedpars, det_num)
-    if (isfield(fedpars, 'detector'))
-        blob_size_add = fedpars.detector(det_num).blobsizeadd;
-    else
-        blob_size_add = fedpars.blobsizeadd(det_num, :);
-    end
-end
-
-function apply_uv_shifts = get_conf_apply_uv_shifts(fedpars, det_num)
-    if (isfield(fedpars, 'detector') ...
-            && isfield(fedpars.detector, 'apply_uv_shift') ...
-            && ~isempty(fedpars.detector(det_num).apply_uv_shift) ...
-            && fedpars.detector(det_num).apply_uv_shift)
-        apply_uv_shifts = true;
-    else
-        apply_uv_shifts = false;
-    end
-end
-
-function psf = get_conf_psf(fedpars, det_num, nbl)
-    if (isfield(fedpars, 'detector') ...
-            && isfield(fedpars.detector, 'psf') ...
-            && ~isempty(fedpars.detector(det_num).psf))
-        psf = fedpars.detector(det_num).psf;
-        if (~iscell(psf))
-            psf = { psf };
-        end
-        if (numel(psf) == 1)
-            ones_bl = ones(1, nbl);
-            psf = psf(ones_bl);
-        end
-    else
-        psf = {};
-    end
-end
-
-

zUtil_Deformation/gtDefFwdProjGvdm.m

@@ -20,8 +20,6 @@ function bl = gtDefFwdProjGvdm(grain, ref_sel, gv, fedpars, parameters, det_ind,
     labgeo = parameters.labgeo;
     samgeo = parameters.samgeo;
 
-    sam2lab_rescaled = gtGeoSam2Lab_GetTransformRescaled(samgeo);
-
     g = gtMathsRod2OriMat(grain.R_vector);
 
     pls_orig = grain.allblobs.plorig(ref_sel, :);
@@ -120,7 +118,7 @@ function bl = gtDefFwdProjGvdm(grain, ref_sel, gv, fedpars, parameters, det_ind,
             dvec_lab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');
 
             rot_s2l = permute(rot_l2s, [2 1 3]);
-            gvcs_lab = gtGeoSam2Lab_Position(gvcs', sam2lab_rescaled, rot_s2l, labgeo, samgeo);
+            gvcs_lab = gtGeoSam2Lab(gvcs', rot_s2l, labgeo, samgeo, false);
 
             uvw_bl = gtFedPredictUVWMultiple([], dvec_lab, gvcs_lab', ...
                 detgeo.detrefpos', detgeo.detnorm', detgeo.Qdet, ...

zUtil_Geo/gtGeoLab2Sam.m

@@ -54,7 +54,10 @@ if (~isfield(samgeo,'lab2sam') || (isfield(samgeo,'lab2sam') && isempty(samgeo.l
         % the transformation matrix 
         sam2lab = [samgeo.dirx; samgeo.diry; samgeo.dirz];
     else
-        sam2lab = gtGeoSam2Lab_GetTransformRescaled(samgeo);
+        sam2lab = [ ...
+            samgeo.dirx * samgeo.voxsize(1); ...
+            samgeo.diry * samgeo.voxsize(2); ...
+            samgeo.dirz * samgeo.voxsize(3) ];
     end
 
     lab2sam = inv(sam2lab);
@@ -70,13 +73,14 @@ end
 %%%%%%%%%%%%%%%%%%%%%
 
 n = size(labXYZ, 1);
+ones_n = ones(n, 1);
 
 [size_omega(1), size_omega(2), size_omega(3)] = size(omega);
 % Tensor which rotates with omega from Lab back to Sample:
 if (all(size_omega(1:2) == [3 3]))
     % This means that the rotation tensors were passed already!
     if (size_omega(3) == 1)
-        rottensors = omega(:, :, ones(1, n));
+        rottensors = omega(:, :, ones_n);
     else
         rottensors = omega;
     end
@@ -85,12 +89,25 @@ else
     rottensors = gtMathsRotationTensor(omega, rotcomp);
 end
 
-if (freevec)
-    % free vectors
-    samXYZ = gtGeoLab2Sam_Direction(labXYZ, lab2sam, rottensors);
-else
-    % offset vectors
-    samXYZ = gtGeoLab2Sam_Position(labXYZ, lab2sam, rottensors, labgeo, samgeo);
+if (~freevec)
+    % Vectors relative to rotation axis position
+    labXYZ = labXYZ - labgeo.rotpos(ones_n, :);
+end
+
+% Apply rottensor element-wise
+labXYZ_t = [labXYZ, labXYZ, labXYZ]';
+labXYZ_t = reshape(labXYZ_t, 3, []);
+rottensors_t = reshape(rottensors, 3, []);
+pro = labXYZ_t .* rottensors_t;
+
+% Sum, reshape, transpose to get rotated vector
+labXYZ_rot = reshape(sum(pro, 1), 3, [])';
+
+if (~freevec)
+    % Set offset back to get rotated vector
+    labXYZ_rot = labXYZ_rot + labgeo.rotpos(ones_n, :) - samgeo.orig(ones_n, :);
 end
 
+samXYZ = labXYZ_rot * lab2sam;
+
 end % of function

zUtil_Geo/gtGeoLab2Sam_Direction.m

-function samXYZ = gtGeoLab2Sam_Direction(labXYZ, lab2sam, rottensors)
-% FUNCTION samXYZ = gtGeoLab2Sam_Direction(labXYZ, lab2sam, rottensors)
-
-    n = size(labXYZ, 1);
-
-    % Multiply element-wise
-    pro = reshape([labXYZ, labXYZ, labXYZ]', 3, n * 3) .* reshape(rottensors, 3, n * 3);
-
-    % Sum, reshape, transpose to get rotated vector
-    labXYZ_rot = reshape(sum(pro, 1), 3, n)';
-
-    % Change ref. base
-    samXYZ = labXYZ_rot * lab2sam;
-end

zUtil_Geo/gtGeoLab2Sam_Position.m

-function samXYZ = gtGeoLab2Sam_Position(labXYZ, lab2sam, rottensors, labgeo, samgeo)
-% FUNCTION samXYZ = gtGeoLab2Sam_Position(labXYZ, lab2sam, rottensors, labgeo, samgeo)
-
-    n = size(labXYZ, 1);
-    ones_n = ones(n, 1);
-
-    % Vectors relative to rotation axis position
-    labXYZrel = labXYZ - labgeo.rotpos(ones_n, :);
-
-    % Apply rottensor element-wise
-    labXYZrel_t = [labXYZrel, labXYZrel, labXYZrel]';
-    labXYZrel_t = reshape(labXYZrel_t, 3, n * 3);
-    rottensors_t = reshape(rottensors, 3, n * 3);
-    pro = labXYZrel_t .* rottensors_t;
-
-    % Sum, reshape, transpose, set offset back to get rotated vector
-    labXYZ_rot = reshape(sum(pro, 1), 3, n)' + labgeo.rotpos(ones_n, :);
-
-    % Change ref. base
-    samXYZ = (labXYZ_rot - samgeo.orig(ones_n, :)) * lab2sam;
-end

zUtil_Geo/gtGeoSam2Lab.m

@@ -53,7 +53,10 @@ if (~isfield(samgeo,'sam2lab') || (isfield(samgeo,'sam2lab') && isempty(samgeo.s
         % the transformation matrix
         sam2lab = [samgeo.dirx; samgeo.diry; samgeo.dirz];
     else
-        sam2lab = gtGeoSam2Lab_GetTransformRescaled(samgeo);
+        sam2lab = [ ...
+            samgeo.dirx * samgeo.voxsize(1); ...
+            samgeo.diry * samgeo.voxsize(2); ...
+            samgeo.dirz * samgeo.voxsize(3) ];
     end
     
 elseif (rot_flag)
@@ -67,13 +70,14 @@ end
 %%%%%%%%%%%%%%%%%%%%%
 
 n = size(samXYZ,1);
+ones_n = ones(n, 1);
 
 [size_omega(1), size_omega(2), size_omega(3)] = size(omega);
 % Tensor which rotates with omega from Sample to Lab
 if (all(size_omega(1:2) == [3 3]))
     % This means that the rotation tensors were passed already!
     if (size_omega(3) == 1)
-        rottensors = omega(:, :, ones(1, n));
+        rottensors = omega(:, :, ones_n);
     else
         rottensors = omega;
     end
@@ -82,12 +86,27 @@ else
     rottensors = gtMathsRotationTensor(-omega, rotcomp);
 end
 
-if (freevec)
-    % free vectors
-    labXYZ = gtGeoSam2Lab_Direction(samXYZ, sam2lab, rottensors);
-else
-    % offset vectors
-    labXYZ = gtGeoSam2Lab_Position(samXYZ, sam2lab, rottensors, labgeo, samgeo);
+% Vectors in Lab reference at the origin before rotation
+labXYZprerot = samXYZ * sam2lab;
+
+if (~freevec)
+    % Vectors relative to rot. axis position in Lab reference
+    vv           = samgeo.orig - labgeo.rotpos;
+    labXYZprerot = labXYZprerot + vv(ones_n,:);
+end
+
+% Multiply element-wise to avoid loop
+labXYZprerot_t = [labXYZprerot, labXYZprerot, labXYZprerot]';
+labXYZprerot_t = reshape(labXYZprerot_t, 3, []);
+rottensors_t = reshape(rottensors, 3, []);
+pro = labXYZprerot_t .* rottensors_t;
+
+% Sum, reshape, transpose
+labXYZ = reshape(sum(pro, 1), 3, [])';
+
+if (~freevec)
+    % Set offset back
+    labXYZ = labXYZ + labgeo.rotpos(ones_n, :);
 end
 
 end % of function

zUtil_Geo/gtGeoSam2Lab_Direction.m

-function labXYZ = gtGeoSam2Lab_Direction(samXYZ, sam2lab, rottensors)
-% FUNCTION labXYZ = gtGeoSam2Lab_Direction(samXYZ, sam2lab, rottensors)
-
-    n = size(samXYZ, 1);
-
-    % Vectors in Lab reference at the origin before rotation
-    labXYZprerot = samXYZ*sam2lab;
-
-    % Multiply element-wise to avoid loop
-    pro = reshape([labXYZprerot, labXYZprerot, labXYZprerot]',3,n*3).*reshape(rottensors,3,n*3);
-
-    % Sum, reshape, transpose
-    labXYZ = reshape(sum(pro,1),3,n)';
-end

zUtil_Geo/gtGeoSam2Lab_GetTransformRescaled.m

-function sam2lab = gtGeoSam2Lab_GetTransformRescaled(samgeo)
-% FUNCTION sam2lab = gtGeoSam2Lab_GetTransformRescaled(samgeo)
-
-    sam2lab = [ ...
-        samgeo.dirx * samgeo.voxsize(1); ...
-        samgeo.diry * samgeo.voxsize(2); ...
-        samgeo.dirz * samgeo.voxsize(3) ];
-end

zUtil_Geo/gtGeoSam2Lab_Position.m

-function labXYZ = gtGeoSam2Lab_Position(samXYZ, sam2lab, rottensors, labgeo, samgeo)
-% FUNCTION labXYZ = gtGeoSam2Lab_Position(samXYZ, sam2lab, rottensors, labgeo, samgeo)
-
-    n = size(samXYZ, 1);
-
-    % Vectors relative to rot. axis position in Lab reference
-    vv           = samgeo.orig - labgeo.rotpos;
-    labXYZprerot = samXYZ*sam2lab + vv(ones(n,1),:);
-
-    % Apply rottensor element-wise to avoid loop
-    labXYZprerot_t = [labXYZprerot, labXYZprerot, labXYZprerot]';
-    labXYZprerot_t = reshape(labXYZprerot_t, 3, n*3);
-    rottensors_t = reshape(rottensors, 3, n*3);
-    pro = labXYZprerot_t .* rottensors_t;
-
-    % Sum, reshape, transpose, set offset back
-    labXYZ = reshape(sum(pro,1),3,n)' + labgeo.rotpos(ones(n,1),:);
-end

zUtil_Geo/gtGeoSam2Sam.m

@@ -52,7 +52,10 @@ function [xyzB, samA2samB] = gtGeoSam2Sam(xyzA, samgeoA, samgeoB, freevec, ...
             % the transformation matrix
             sam2labA = [samgeoA.dirx; samgeoA.diry; samgeoA.dirz];
         else
-            sam2labA = gtGeoSam2Lab_GetTransformRescaled(samgeoA);
+            sam2labA = [ ...
+                samgeoA.dirx * samgeoA.voxsize(1); ...
+                samgeoA.diry * samgeoA.voxsize(2); ...
+                samgeoA.dirz * samgeoA.voxsize(3) ];
         end
     elseif (rot_flag)
         error('Transformation matrix already defined. Set ''rot_flag'' to false.')
@@ -85,7 +88,10 @@ function [xyzB, samA2samB] = gtGeoSam2Sam(xyzA, samgeoA, samgeoB, freevec, ...
             % the transformation matrix
             sam2labB = [samgeoB.dirx; samgeoB.diry; samgeoB.dirz];
         else
-            sam2labB = gtGeoSam2Lab_GetTransformRescaled(samgeoB);
+            sam2labB = [ ...
+                samgeoB.dirx * samgeoB.voxsize(1); ...
+                samgeoB.diry * samgeoB.voxsize(2); ...
+                samgeoB.dirz * samgeoB.voxsize(3) ];
         end
         lab2samB = inv(sam2labB);