diff --git a/zUtil_Deformation/GtGrainODFuvwSolver.m b/zUtil_Deformation/GtGrainODFuvwSolver.m
index f668980971291c02af7b646a3f19a8dffb13cff6..6b19033cf7a5c794d2eb1c1633f1375e3603f216 100644
--- a/zUtil_Deformation/GtGrainODFuvwSolver.m
+++ b/zUtil_Deformation/GtGrainODFuvwSolver.m
@@ -1,8 +1,11 @@
classdef GtGrainODFuvwSolver < GtGrainODFwSolver
properties
- voxel_centers = [];
- rspace_over_sampling = 1;
size_rspace_volume = [];
+
+ voxel_centers = [];
+
+ rspace_oversampling = 1;
+ rspace_downscaling = 1;
end
methods (Access = public)
@@ -12,22 +15,27 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
end
methods (Access = public) % Low Level API
- function build_orientation_sampling(self, ref_gr, det_index, oversize, oversampling, shape_functions)
- self.build_orientation_sampling@GtGrainODFwSolver(ref_gr, det_index, oversize, oversampling, shape_functions);
+ function build_orientation_sampling(self, ref_gr, det_index)
+ self.build_orientation_sampling@GtGrainODFwSolver(ref_gr, det_index);
self.chose_voxel_centers();
num_vox_centers = size(self.voxel_centers, 1);
switch (self.shape_functions_type)
- case 'nw'
+ case 'nw2uvw'
+ recgeo = self.parameters.recgeo(self.sampler.detector_index);
+
sf_nw = gtDefShapeFunctionsCreateNW(self.sampler, 2);
for ii_v = 1:num_vox_centers
fprintf('(%03d/%03d) ', ii_v, num_vox_centers);
disp_sampler = self.sampler.regenerate_displaced(self.voxel_centers(ii_v, :));
- self.shape_functions{ii_v} = gtDefShapeFunctionsNW2UVW(disp_sampler, sf_nw, false);
+ self.shape_functions{ii_v} = gtDefShapeFunctionsNW2UVW( ...
+ disp_sampler, sf_nw, 'recenter_sf', false, ...
+ 'rspace_oversampling', self.rspace_downscaling, ...
+ 'rspace_voxel_size', recgeo.voxsize * self.rspace_downscaling);
end
case 'uvw'
shape_functions_name = sprintf('tmp_shape_functions_uvw_%d_%d', ...
@@ -43,7 +51,8 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
fprintf('(%3d/%03d) ', ii_v, num_vox_centers);
disp_sampler = self.sampler.regenerate_displaced(self.voxel_centers(ii_v, :));
- self.shape_functions{ii_v} = gtDefShapeFunctionsFwdProjUVW(disp_sampler, 15, [], false);
+ self.shape_functions{ii_v} = gtDefShapeFunctionsFwdProjUVW( ...
+ disp_sampler, 'factor', 15, 'recenter_sf', false);
end
sfs = self.shape_functions;
save([shape_functions_name '.mat'], 'sfs', '-v7.3')
@@ -73,7 +82,7 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
with_shape_functions = ~isempty(self.shape_functions);
if (with_shape_functions)
- if (any(strcmpi(self.shape_functions_type, {'uvw', 'nw'})))
+ if (any(strcmpi(self.shape_functions_type, {'uvw', 'nw2uvw'})))
for ii_o = numel(self.shape_functions{1}):-1:1
proj_w_lims(:, :, ii_o) = cat(1, self.shape_functions{1}{ii_o}(:).bbwim);
end
@@ -147,10 +156,10 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
recgeo = self.parameters.recgeo(self.sampler.detector_index);
ref_gr = self.sampler.get_reference_grain();
proj = ref_gr.proj(self.sampler.detector_index);
- r_space_vol_size = [proj.vol_size_y, proj.vol_size_x, proj.vol_size_z];
+ rspace_vol_dims = ceil([proj.vol_size_y, proj.vol_size_x, proj.vol_size_z] / self.rspace_downscaling);
- max_dists_from_center = (r_space_vol_size - 1 ./ self.rspace_over_sampling) / 2;
- self.size_rspace_volume = r_space_vol_size .* self.rspace_over_sampling;
+ max_dists_from_center = (rspace_vol_dims - 1 ./ self.rspace_oversampling) / 2;
+ self.size_rspace_volume = rspace_vol_dims .* self.rspace_oversampling;
tot_voxels = prod(self.size_rspace_volume);
self.voxel_centers = zeros(tot_voxels, 3);
@@ -161,7 +170,7 @@ classdef GtGrainODFuvwSolver < GtGrainODFwSolver
for ss_y = linspace(-max_dists_from_center(2), max_dists_from_center(2), self.size_rspace_volume(2))
for ss_x = linspace(-max_dists_from_center(1), max_dists_from_center(1), self.size_rspace_volume(1))
- offset_voxel = [ss_x, ss_y, ss_z] .* recgeo.voxsize;
+ offset_voxel = [ss_x, ss_y, ss_z] .* recgeo.voxsize * self.rspace_downscaling;
self.voxel_centers(counter, :) = ref_gr.center + offset_voxel;
diff --git a/zUtil_Deformation/GtGrainODFwSolver.m b/zUtil_Deformation/GtGrainODFwSolver.m
index 0bc928b5874577e1797133e460f061e2f1ee7be7..fe9ae474f5e863b6fca7dab7eb1bdb72095ddf46 100644
--- a/zUtil_Deformation/GtGrainODFwSolver.m
+++ b/zUtil_Deformation/GtGrainODFwSolver.m
@@ -25,6 +25,9 @@ classdef GtGrainODFwSolver < handle
verbose = false;
data_type = 'double';
+
+ ospace_oversize = 1;
+ ospace_oversampling = 1;
end
methods (Access = public)
@@ -40,25 +43,16 @@ classdef GtGrainODFwSolver < handle
% INPUT (varargin):
% 'algorithm': {'cplsnn'} | 'sirt' | 'cplsl1nn' | 'cpl1nn'
% 'lambda': 1e-2
- % 'ospace_oversize': 1
- % 'ospace_oversampling': 1
- % 'shape_functions': {'none'} | 'w' | 'nw' | 'uvw'
% 'det_index': 1
%
conf = struct( ...
'algorithm', 'cplsnn', ...
'lambda', 1e-2, ...
- 'ospace_oversize', 1, ...
- 'ospace_oversampling', 1, ...
- 'rspace_oversampling', 1, ...
- 'shape_functions', 'none', ...
'det_index', 1 );
[conf, ~] = parse_pv_pairs(conf, varargin);
- self.build_orientation_sampling(ref_gr, conf.det_index, ...
- conf.ospace_oversize, conf.ospace_oversampling, ...
- conf.shape_functions );
+ self.build_orientation_sampling(ref_gr, conf.det_index);
self.build_sinogram();
self.build_projection_matrices();
@@ -88,29 +82,24 @@ classdef GtGrainODFwSolver < handle
end
methods (Access = public) % Low Level API
- function build_orientation_sampling(self, ref_gr, det_index, oversize, oversampling, shape_functions)
+ function build_orientation_sampling(self, ref_gr, det_index)
if (~exist('det_index', 'var') || isempty(det_index))
det_index = 1;
end
- if (~exist('oversampling', 'var') || isempty(oversampling))
- oversampling = 0;
- end
self.sampler = GtOrientationSampling(self.parameters, ref_gr, ...
'detector_index', det_index);
% 'verbose', self.verbose,
- self.sampler.make_grid_estim_ODF_resoluion('cubic', -oversampling, oversize);
+ self.sampler.make_grid_estim_ODF_resoluion('cubic', ...
+ -self.ospace_oversampling, self.ospace_oversize);
self.size_volume = size(self.sampler.lattice.gr);
- if (exist('shape_functions', 'var') && ~isempty(shape_functions))
- switch (shape_functions)
- case 'none'
- self.shape_functions = [];
- case 'w'
- self.shape_functions = gtDefShapeFunctionsCreateW(self.sampler);
- end
- self.shape_functions_type = shape_functions;
+ switch (self.shape_functions_type)
+ case 'none'
+ self.shape_functions = [];
+ case 'w'
+ self.shape_functions = gtDefShapeFunctionsCreateW(self.sampler);
end
end
@@ -197,7 +186,7 @@ classdef GtGrainODFwSolver < handle
self.build_projection_matrices_sf_none();
case 'w'
self.build_projection_matrices_sf_w();
- case {'nw', 'uvw'}
+ case {'nw2uvw', 'uvw'}
self.build_projection_matrices_sf_uvw();
end
fprintf('\b\b: Done in %f seconds.\n', toc(c));
diff --git a/zUtil_Deformation/gtDefShapeFunctionsFwdProjUVW.m b/zUtil_Deformation/gtDefShapeFunctionsFwdProjUVW.m
index 5f559e5cfde32a372e6d4fb8d6d61ec77aedd665..7d60a3f1ef5eb216deaa6be2069d90a4e30547fc 100644
--- a/zUtil_Deformation/gtDefShapeFunctionsFwdProjUVW.m
+++ b/zUtil_Deformation/gtDefShapeFunctionsFwdProjUVW.m
@@ -1,14 +1,13 @@
-function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, factor, dtype, recenter_sf)
+function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, varargin)
- if (~exist('dtype', 'var') || isempty(dtype))
- dtype = 'single';
- end
-
- if (~exist('recenter_sf', 'var') || isempty(recenter_sf))
- recenter_sf = true;
- end
-
- interpolated_voxel = false;
+ conf = struct( ...
+ 'recenter_sf', true, ...
+ 'dtype', [], ...
+ 'factor', 1, ...
+ 'rspace_oversampling', 1, ...
+ 'rspace_voxel_size', [1 1 1], ...
+ 'interpolated_voxel', false);
+ conf = parse_pv_pairs(conf, varargin);
num_detectors = numel(sampler.parameters.detgeo);
det_ind = sampler.detector_index;
@@ -18,14 +17,14 @@ function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, factor, dtype, rec
sampler.parameters.labgeo = labgeo;
fed_pars_detector = struct(...
- 'blobsizeadd', [0 0 0], 'psf', [], 'apply_uv_shift', recenter_sf);
+ 'blobsizeadd', [0 0 0], 'psf', [], 'apply_uv_shift', conf.recenter_sf);
fedpars = struct( ...
- 'bltype', dtype, ...
+ 'bltype', conf.dtype, ...
'dcomps', [1 1 1 0 0 0 0 0 0] == 1, ...
'defmethod', 'rod_rightstretch', ...
'detector', repmat( fed_pars_detector, [num_detectors, 1]));
- voxel_factor = factor([1 1 1]);
+ voxel_factor = conf.factor([1 1 1]);
num_sub_voxels = prod(voxel_factor);
ones_nsv = ones(num_sub_voxels, 1);
@@ -35,11 +34,11 @@ function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, factor, dtype, rec
gv.ind = gv.ind';
ref_gr = sampler.get_reference_grain();
- gv.cs = ref_gr.center(ones_nsv, :)';
+ gv.cs = compute_centers(ref_gr.center, conf.rspace_voxel_size, conf.rspace_oversampling)';
% computing the subsampling, used to deterine the shape functions
voxel_size = sampler.stats.sampling.gaps;
- if (interpolated_voxel)
+ if (conf.interpolated_voxel)
half_voxel_size = voxel_size; % Enlarging to neighbours
steps = voxel_size ./ voxel_factor * 2;
else
@@ -69,10 +68,10 @@ function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, factor, dtype, rec
end
end
- if (interpolated_voxel)
+ if (conf.interpolated_voxel)
gv.pow = prod(1 - abs(gv.d ./ half_voxel_size(ones_nsv, :)'), 1);
else
- gv.pow = ones(1, num_sub_voxels, dtype);
+ gv.pow = ones(1, num_sub_voxels, conf.dtype);
end
gv.pow = gv.pow / sum(gv.pow);
@@ -99,3 +98,26 @@ function shape_funcs = gtDefShapeFunctionsFwdProjUVW(sampler, factor, dtype, rec
fprintf('Done in %g seconds.\n', toc(c))
end
+
+function centers = compute_centers(ref_center, rspace_voxel_size, rspace_oversampling)
+ if (numel(rspace_oversampling) == 1)
+ rspace_oversampling = rspace_oversampling([1 1 1]);
+ end
+ max_dists_from_center = rspace_voxel_size .* (1 - 1 ./ rspace_oversampling) / 2;
+ tot_voxels = prod(rspace_oversampling);
+
+ centers = zeros(tot_voxels, 3);
+
+ counter = 1;
+
+ for ss_z = linspace(-max_dists_from_center(3), max_dists_from_center(3), rspace_oversampling(3))
+ for ss_y = linspace(-max_dists_from_center(2), max_dists_from_center(2), rspace_oversampling(2))
+ for ss_x = linspace(-max_dists_from_center(1), max_dists_from_center(1), rspace_oversampling(1))
+
+ centers(counter, :) = ref_center + [ss_x, ss_y, ss_z];
+
+ counter = counter + 1;
+ end
+ end
+ end
+end
diff --git a/zUtil_Deformation/gtDefShapeFunctionsNW2UVW.m b/zUtil_Deformation/gtDefShapeFunctionsNW2UVW.m
index 868b86107bcb861616d9ce4c58dac98d6ac3e0dc..03cfe4beb35bb2aaece08172fc70b1c9858113ce 100644
--- a/zUtil_Deformation/gtDefShapeFunctionsNW2UVW.m
+++ b/zUtil_Deformation/gtDefShapeFunctionsNW2UVW.m
@@ -1,10 +1,13 @@
-function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
-% FUNCTION sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw)
+function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, varargin)
+% FUNCTION sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, varargin)
%
- if (~exist('recenter_sf', 'var') || isempty(recenter_sf))
- recenter_sf = true;
- end
+ conf = struct( ...
+ 'recenter_sf', true, ...
+ 'dtype', [], ...
+ 'rspace_oversampling', 1, ...
+ 'rspace_voxel_size', [1 1 1] );
+ conf = parse_pv_pairs(conf, varargin);
det_ind = sampler.detector_index;
@@ -13,6 +16,7 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
samgeo = p.samgeo;
detgeo = p.detgeo(det_ind);
+ ref_gr = sampler.get_reference_grain();
ors = sampler.get_orientations();
ref_sel = sampler.selected;
@@ -25,6 +29,9 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
'Mismatching number of orientations and shape functions passed')
end
+ centers = compute_centers(ref_gr.center, conf.rspace_voxel_size, conf.rspace_oversampling);
+ tot_centers = size(centers, 1);
+
sfs_uvw = cell(tot_orientations, 1);
rotcomp_w = gtMathsRotationMatrixComp(labgeo.rotdir, 'row');
@@ -48,7 +55,7 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
blobs_nw_sizes = cat(1, sfs_in(:).bbsize);
or_uv = ors{ii_o}.allblobs.detector(det_ind).uvw(ref_sel, 1:2);
- if (recenter_sf)
+ if (conf.recenter_sf)
or_uv_shift = or_uv;
else
or_uv_shift = round(or_uv);
@@ -89,12 +96,9 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
% Calculating all the sample center positions at each omega inside
% the blobs
- ones_all_w = ones(numel(all_ws), 1);
- gcsam_v_all = ors{ii_o}.center(ones_all_w, :);
-
rot_s2l_all = gtMathsRotationTensor(all_ws, rotcomp_w);
- gclab_v_all = gtGeoSam2Lab(gcsam_v_all, rot_s2l_all, labgeo, samgeo, false);
- gclab_v_all = permute(gclab_v_all, [2 3 1]);
+ gclab_v_all = gtGeoSam2Lab(centers, rot_s2l_all, labgeo, samgeo, false, [], false);
+ gclab_v_all = reshape(gclab_v_all', 3, 1, tot_ws, []);
%%% Here we put the scattering vectors and center positions
%%% together
@@ -102,9 +106,9 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
uv_min = zeros(tot_blobs, 2);
uv_max = zeros(tot_blobs, 2);
- blobs_nw_area_sizes = prod(blobs_nw_sizes, 2);
+ blobs_nw_area_sizes = prod(blobs_nw_sizes, 2) * tot_centers;
- chunk = 32;
+ chunk = 16;
for ii_base_b = 0:chunk:tot_blobs-1
chunk_size = min(chunk, tot_blobs-ii_base_b);
@@ -120,7 +124,7 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
ii_b = ii_base_b + ii_var_b;
ones_n = ones(blobs_nw_sizes(ii_b, 1), 1);
- ones_w = ones(blobs_nw_sizes(ii_b, 2), 1);
+ ones_w = ones(blobs_nw_sizes(ii_b, 2) * tot_centers, 1);
inds_bl_chunk = inds_chunk(ii_var_b, 1):inds_chunk(ii_var_b, 2);
@@ -128,7 +132,8 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
dveclab = dveclab(:, :, ones_w);
dveclab_chunk(:, inds_bl_chunk) = reshape(dveclab, 3, []);
- gclab_v = gclab_v_all(:, :, inds_w(ii_b, 1):inds_w(ii_b, 2));
+ gclab_v = gclab_v_all(:, :, inds_w(ii_b, 1):inds_w(ii_b, 2), :);
+ gclab_v = reshape(gclab_v, size(gclab_v, 1), size(gclab_v, 2), []);
gclab_v = gclab_v(:, ones_n, :);
gclab_v_chunk(:, inds_bl_chunk) = reshape(gclab_v, 3, []);
end
@@ -161,10 +166,22 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
ws = 1:blobs_nw_sizes(ii_b, 2);
ws = ws(ones(blobs_nw_sizes(ii_b, 1), 1), :);
ws = reshape(ws, [], 1);
+ ws = reshape(ws(:, ones(tot_centers, 1)), [], 1);
- ones_nw = ones(prod(blobs_nw_sizes(ii_b, :)), 1);
+ ones_nw = ones(blobs_nw_area_sizes(ii_b), 1);
uv_bl = uv{ii_b} + uv_shift(ones_nw, :);
+ ints = reshape(sfs_in(ii_b).intm, [], 1);
+
+ if (isempty(conf.dtype))
+ uv_bl = cast(uv_bl, 'like', sfs_in(ii_b).intm);
+ else
+ uv_bl = cast(uv_bl, conf.dtype);
+ ints = cast(ints, conf.dtype);
+ end
+
+ ints = ints / tot_centers;
+ ints = reshape(ints(:, ones(tot_centers, 1)), [], 1);
uv_l = floor(uv_bl); % lower pixel indices
uv_h = uv_l + 1; % higher pixel indices
@@ -176,10 +193,9 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
[uv_h(:, 1), uv_l(:, 2), ws]; ...
[uv_h, ws] ];
- uv_l_c = cast(uv_h - uv_bl, 'like', sfs_in(ii_b).intm);
+ uv_l_c = uv_h - uv_bl;
uv_h_c = 1 - uv_l_c;
- ints = reshape(sfs_in(ii_b).intm, [], 1);
ints4 = [ ...
uv_l_c(:, 1) .* uv_l_c(:, 1) .* ints; ...
uv_l_c(:, 1) .* uv_h_c(:, 2) .* ints; ...
@@ -206,3 +222,27 @@ function sfs_uvw = gtDefShapeFunctionsNW2UVW(sampler, sfs_nw, recenter_sf)
fprintf(repmat('\b', [1, num_chars]));
fprintf('Done in %g seconds.\n', toc(c))
end
+
+function centers = compute_centers(ref_center, rspace_voxel_size, rspace_oversampling)
+ if (numel(rspace_oversampling) == 1)
+ rspace_oversampling = rspace_oversampling([1 1 1]);
+ end
+ max_dists_from_center = rspace_voxel_size .* (1 - 1 ./ rspace_oversampling) / 2;
+ tot_voxels = prod(rspace_oversampling);
+
+ centers = zeros(tot_voxels, 3);
+
+ counter = 1;
+
+ for ss_z = linspace(-max_dists_from_center(3), max_dists_from_center(3), rspace_oversampling(3))
+ for ss_y = linspace(-max_dists_from_center(2), max_dists_from_center(2), rspace_oversampling(2))
+ for ss_x = linspace(-max_dists_from_center(1), max_dists_from_center(1), rspace_oversampling(1))
+
+ centers(counter, :) = ref_center + [ss_x, ss_y, ss_z];
+
+ counter = counter + 1;
+ end
+ end
+ end
+end
+