Improved synthetic data volume super-sampling

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

zUtil_Deformation/gtDefFwdProjGvdm2UVP.m

@@ -88,7 +88,7 @@ function bl = gtDefFwdProjGvdm2UVP(grain, selectedph, gv, fedpars, parameters, d
         && strcmpi(fedpars.projector, 'nearest'));
 
     num_oversampling = size(gv.d, 3);
-    gvpow = gv.pow(1, uinds);
+    gvpow_uinds = gv.pow(1, uinds, :);
 
     % This usually happens for shape functions, where we only have one
     % center
@@ -102,6 +102,11 @@ function bl = gtDefFwdProjGvdm2UVP(grain, selectedph, gv, fedpars, parameters, d
         else
             gvpcs = gv.pcs(:, 1, ii_ss)';
         end
+        if (size(gv.pow, 3) > 1)
+            gvpow = gvpow_uinds(1, :, ii_ss);
+        else
+            gvpow = gvpow_uinds;
+        end
 
         gvd = gv.d(:, uinds, ii_ss);
 

zUtil_Deformation/gtDefFwdProjGvdm2UVW.m

@@ -72,7 +72,7 @@ function bl = gtDefFwdProjGvdm2UVW(grain, ref_sel, gv, fedpars, parameters, det_
         && strcmpi(fedpars.projector, 'nearest'));
 
     num_oversampling = size(gv.d, 3);
-    gvpow = gv.pow(1, uinds);
+    gvpow_uinds = gv.pow(1, uinds, :);
 
     % This usually happens for shape functions, where we only have one
     % center
@@ -86,6 +86,11 @@ function bl = gtDefFwdProjGvdm2UVW(grain, ref_sel, gv, fedpars, parameters, det_
         else
             gvpcs = gv.pcs(:, 1, ii_ss);
         end
+        if (size(gv.pow, 3) > 1)
+            gvpow = gvpow_uinds(1, :, ii_ss);
+        else
+            gvpow = gvpow_uinds;
+        end
 
         gvd = gv.d(:, uinds, ii_ss);
 

zUtil_Deformation/gtDefSyntheticGrainCreate.m

@@ -154,8 +154,14 @@ gv.cs = vox000sam(:, nv_ones) + (gv.ind - 0.5) .* phantom_voxsize(:, nv_ones);
 % Producing centers and deformation components for projection
 if (isfield(fedpars, 'volume_super_sampling') ...
         && (fedpars.volume_super_sampling > 1))
-    [gv.pcs, gv.d] = compute_oversampling_voxels(gv, dmvol, ...
-        phantom_voxsize, fedpars.volume_super_sampling, fedpars.gvtype);
+    if (~isfield(fedpars, 'volume_sampling_method') ...
+            || isempty(fedpars.volume_sampling_method))
+        [gv.pcs, gv.d] = compute_oversampling_voxels(gv, dmvol, ...
+            phantom_voxsize, fedpars.volume_super_sampling, fedpars.gvtype);
+    else
+        [gv.pcs, gv.d, gv.pow] = compute_oversampling_voxels_interpn(gv, dmvol, ...
+            fedpars.volume_super_sampling, fedpars.gvtype, fedpars.volume_sampling_method);
+    end
 else
     gv.pcs = gv.cs;
     gv.d = gv.dm;
@@ -392,4 +398,62 @@ function [gvpcs, gvd] = compute_oversampling_voxels(gv, dmvol, voxsize, vssampli
     end
 end
 
+function [gvpcs, gvd, gvpow] = compute_oversampling_voxels_interpn(gv, dmvol, vssampling, gvtype, method)
+    nv = size(gv.ind, 2);
+
+    num_subvoxels = vssampling ^ 3;
+    gvpcs = zeros(3, nv, num_subvoxels, gvtype);
+    gvd = zeros(3, nv, num_subvoxels, gvtype);
+
+    dmvol_size = [size(dmvol, 1), size(dmvol, 2), size(dmvol, 3)];
+    gvcs = reshape(gv.cs', [dmvol_size, 3]);
+
+    fprintf('\b\b: ')
+
+    for ii_c = 1:3
+        num_chars = fprintf('%03d/%03d', ii_c, 3);
+
+        gvpcs(ii_c, :, :) = expand_array(gvcs(:, :, :, ii_c), vssampling, true, 'linear');
+        gvd(ii_c, :, :) = expand_array(dmvol(:, :, :, ii_c), vssampling, false, method);
+
+        fprintf(repmat('\b', [1 num_chars]));
+    end
+
+    r_vec_magn = sqrt(sum(dmvol .^ 2, 4));
+    r_vec_magn = expand_array(r_vec_magn, vssampling, false, method);
+
+    r_vec_expand_magn = sqrt(sum(gvd .^ 2, 1));
+    r_vec_expand_magn = r_vec_expand_magn + (r_vec_expand_magn < eps('single'));
+
+    renorm_r_vec_magn = r_vec_magn ./ r_vec_expand_magn;
+    gvd = bsxfun(@times, gvd, renorm_r_vec_magn);
+
+    gvpow = reshape(gv.pow, dmvol_size);
+    gvpow = expand_array(gvpow, vssampling, false, method);
+end
+
+function arr = expand_array(arr, ssampling, is_spatial, method)
+    arr_size = [size(arr, 1), size(arr, 2), size(arr, 3)];
+
+    if (is_spatial)
+        arr = cat(1, 2 * arr(1, :, :) - arr(2, :, :), arr, 2 * arr(end, :, :) - arr(end-1, :, :));
+        arr = cat(2, 2 * arr(:, 1, :) - arr(:, 2, :), arr, 2 * arr(:, end, :) - arr(:, end-1, :));
+        arr = cat(3, 2 * arr(:, :, 1) - arr(:, :, 2), arr, 2 * arr(:, :, end) - arr(:, :, end-1));
+    else
+        arr = padarray(arr, [1 1 1], 'both');
+    end
+
+    corner = (1 - (1 / ssampling)) / 2;
+    xx = linspace(2 - corner, arr_size(1) + 1 + corner, ssampling * arr_size(1));
+    yy = linspace(2 - corner, arr_size(2) + 1 + corner, ssampling * arr_size(2));
+    zz = linspace(2 - corner, arr_size(3) + 1 + corner, ssampling * arr_size(3));
+    [xx, yy, zz] = meshgrid(xx, yy, zz);
+
+    arr = interp3(arr, xx, yy, zz, method);
+
+    arr = reshape(arr, [ssampling, arr_size(1), ssampling, arr_size(2), ssampling, arr_size(3)]);
+    arr = permute(arr, [2 4 6 1 3 5]);
+    arr = reshape(arr, 1, arr_size(1) * arr_size(2) * arr_size(3), ssampling ^ 3);
+end
+