From a0af4c2e10bce17484005c2600113ff8cb81dc8d Mon Sep 17 00:00:00 2001
From: Nicola Vigano <nicola.vigano@esrf.fr>
Date: Thu, 17 Mar 2016 14:55:28 +0100
Subject: [PATCH] Grain-ODF: added another solver that also can use the uvw
shape functions
Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>
---
zUtil_Deformation/GtGrainODFuvwPatchSolver.m | 140 ++++++---------
zUtil_Deformation/GtGrainODFuvwSfSolver.m | 174 +++++++++++++++++++
zUtil_Deformation/GtGrainODFwSolver.m | 2 +-
3 files changed, 232 insertions(+), 84 deletions(-)
create mode 100644 zUtil_Deformation/GtGrainODFuvwSfSolver.m
diff --git a/zUtil_Deformation/GtGrainODFuvwPatchSolver.m b/zUtil_Deformation/GtGrainODFuvwPatchSolver.m
index 9c09bedd..2e762749 100644
--- a/zUtil_Deformation/GtGrainODFuvwPatchSolver.m
+++ b/zUtil_Deformation/GtGrainODFuvwPatchSolver.m
@@ -1,10 +1,11 @@
-classdef GtGrainODFuvwPatchSolver < GtGrainODFuvwSolver
+classdef GtGrainODFuvwPatchSolver < GtGrainODFuvwSfSolver
properties
+ shape_functions_extracted = {};
end
methods (Access = public)
function self = GtGrainODFuvwPatchSolver(parameters, varargin)
- self = self@GtGrainODFuvwSolver(parameters, varargin{:});
+ self = self@GtGrainODFuvwSfSolver(parameters, varargin{:});
end
end
@@ -16,35 +17,49 @@ classdef GtGrainODFuvwPatchSolver < GtGrainODFuvwSolver
end
function build_orientation_sampling(self, ref_gr, det_index, oversize, oversampling, shape_functions)
- self.build_orientation_sampling@GtGrainODFuvwSolver(ref_gr, det_index, oversize, oversampling, shape_functions);
-
- shape_functions_name = sprintf('tmp_shape_functions_uvw_%d_%d', ...
- size(self.voxel_centers, 1), self.sampler.get_total_num_orientations());
-
- if (evalin('base', ['exist(''' shape_functions_name ''', ''var'')']))
- fprintf('Loading shape functions..')
- c = tic();
- self.shape_functions = evalin('base', shape_functions_name);
- fprintf('\b\b: Done in %g seconds\n', toc(c));
- elseif (~exist([shape_functions_name '.mat'], 'file'))
- switch (self.shape_functions_type)
- case 'uvw'
- for ii_v = 1:size(self.voxel_centers, 1);
- disp_sampler = self.sampler.regenerate_displaced(self.voxel_centers(ii_v, :));
-
- self.shape_functions{ii_v} = gtDefCreateShapeFuncs(disp_sampler, 15, [], false);
- end
+ self.build_orientation_sampling@GtGrainODFuvwSfSolver(ref_gr, det_index, oversize, oversampling, shape_functions);
+
+ for ii_o = 1:self.sampler.get_total_num_orientations()
+ for ii_v = 1:size(self.voxel_centers, 1);
+ self.shape_functions_extracted{ii_o, ii_v} = {self.shape_functions{ii_v}{ii_o}(:).intm};
end
- sfs = self.shape_functions;
- save([shape_functions_name '.mat'], 'sfs', '-v7.3')
- assignin('base', shape_functions_name, sfs);
- else
- fprintf('Loading shape functions..')
- c = tic();
- sfs_file = load([shape_functions_name '.mat']);
- self.shape_functions = sfs_file.sfs;
- assignin('base', shape_functions_name, sfs_file.sfs);
- fprintf('\b\b: Done in %g seconds\n', toc(c));
+ end
+ end
+
+ function build_projection_matrices_sf_uvw(self)
+ ref_sel = self.sampler.selected;
+
+ grid_gr = self.sampler.get_orientations();
+ num_orients = numel(grid_gr);
+
+ num_voxels = size(self.voxel_centers, 1);
+ self.S = cell(num_voxels, 1);
+ self.St = cell(num_voxels, 1);
+
+ bls = self.sampler.bl(ref_sel);
+ bls_uvw_orig = cat(1, bls(:).bbuim, bls(:).bbvim, bls(:).bbwim);
+ bls_uvw_orig = reshape(bls_uvw_orig, [], 3, 2);
+ bls_uvw_orig = bls_uvw_orig(:, :, 1);
+ bls_uvw_orig(:, 3) = bls_uvw_orig(:, 3) - self.pre_paddings + 1;
+
+ for ii_v = 1:num_voxels
+ num_chars = fprintf('%03d/%03d', ii_v, num_voxels);
+
+ sf_uvw_shift = cell(num_orients, 1);
+ for ii_o = 1:num_orients
+ sf_uvw_shift{ii_o} = cat(1, ...
+ self.shape_functions{ii_v}{ii_o}(:).bbuim, ...
+ self.shape_functions{ii_v}{ii_o}(:).bbvim, ...
+ self.shape_functions{ii_v}{ii_o}(:).bbwim ...
+ );
+ sf_uvw_shift{ii_o} = reshape(sf_uvw_shift{ii_o}, [], 3, 2);
+ sf_uvw_shift{ii_o} = sf_uvw_shift{ii_o}(:, :, 1) - bls_uvw_orig;
+ end
+ sf_uvw_shift = cat(3, sf_uvw_shift{:});
+
+ self.S{ii_v} = sf_uvw_shift;
+
+ fprintf(repmat('\b', [1 num_chars]));
end
end
@@ -121,59 +136,22 @@ classdef GtGrainODFuvwPatchSolver < GtGrainODFuvwSolver
% fprintf(repmat('\b', [1 num_chars]));
% end
% end
-
- function build_projection_matrices_sf_uvw(self)
- ref_sel = self.sampler.selected;
-
- grid_gr = self.sampler.get_orientations();
- num_orients = numel(grid_gr);
-
- num_voxels = size(self.voxel_centers, 1);
- self.S = cell(num_voxels, 1);
- self.St = cell(num_voxels, 1);
-
- bls = self.sampler.bl(ref_sel);
- bls_uvw_orig = cat(1, bls(:).bbuim, bls(:).bbvim, bls(:).bbwim);
- bls_uvw_orig = reshape(bls_uvw_orig, [], 3, 2);
- bls_uvw_orig = bls_uvw_orig(:, :, 1);
- bls_uvw_orig(:, 3) = bls_uvw_orig(:, 3) - self.pre_paddings + 1;
-
- for ii_v = 1:num_voxels
- num_chars = fprintf('%03d/%03d', ii_v, num_voxels);
-
- sf_uvw_shift = cell(num_orients, 1);
- for ii_o = 1:num_orients
- sf_uvw_shift{ii_o} = cat(1, ...
- self.shape_functions{ii_v}{ii_o}(:).bbuim, ...
- self.shape_functions{ii_v}{ii_o}(:).bbvim, ...
- self.shape_functions{ii_v}{ii_o}(:).bbwim ...
- );
- sf_uvw_shift{ii_o} = reshape(sf_uvw_shift{ii_o}, [], 3, 2);
- sf_uvw_shift{ii_o} = sf_uvw_shift{ii_o}(:, :, 1) - bls_uvw_orig;
- end
- sf_uvw_shift = cat(3, sf_uvw_shift{:});
-
- self.S{ii_v} = sf_uvw_shift;
-
- fprintf(repmat('\b', [1 num_chars]));
- end
- end
end
methods (Access = protected)
function y = fp(self, x)
num_blobs = self.get_num_blobs();
+ num_orients = size(self.shape_functions_extracted, 1);
+ num_voxels = size(self.shape_functions_extracted, 2);
y = gtMathsGetSameSizeZeros(self.sino, self.data_type);
- for ii_v = 1:numel(self.S)
- vsfs = self.shape_functions{ii_v};
-
+ for ii_v = 1:num_voxels
uvw_shifts = self.S{ii_v};
- for ii_o = 1:numel(vsfs)
- sfs = vsfs{ii_o};
- bls = {sfs(:).intm};
+ for ii_o = 1:num_orients
+ bls = self.shape_functions_extracted{ii_o, ii_v};
+
for ii_b = 1:num_blobs
bls{ii_b} = bls{ii_b} .* x(ii_o, ii_v);
end
@@ -184,27 +162,23 @@ classdef GtGrainODFuvwPatchSolver < GtGrainODFuvwSolver
function y = bp(self, x)
num_blobs = self.get_num_blobs();
+ num_orients = size(self.shape_functions_extracted, 1);
+ num_voxels = size(self.shape_functions_extracted, 2);
- num_orients = numel(self.shape_functions{1});
- num_voxels = numel(self.S);
-
- y = zeros([num_orients, numel(self.S)], self.data_type);
+ y = zeros([num_orients, num_voxels], self.data_type);
for ii_v = 1:num_voxels
- vsfs = self.shape_functions{ii_v};
-
uvw_shifts = self.S{ii_v};
for ii_o = 1:num_orients
- sfs = vsfs{ii_o};
- bls = {sfs(:).intm};
+ bls = self.shape_functions_extracted{ii_o, ii_v};
- mask_sino = gtMathsGetSameSizeZeros(self.sino);
+ mask_sino = zeros(self.size_sino, self.data_type);
mask_sino = gtPlaceSubVolumes(mask_sino, bls, [uvw_shifts(:, :, ii_o), (0:num_blobs-1)']);
mask_sino = mask_sino .* x;
- y(ii_o, ii_v) = sum(sum(sum(sum(mask_sino, 1), 2), 3), 4);
+ y(ii_o, ii_v) = gtMathsSumNDVol(mask_sino);
end
end
end
diff --git a/zUtil_Deformation/GtGrainODFuvwSfSolver.m b/zUtil_Deformation/GtGrainODFuvwSfSolver.m
new file mode 100644
index 00000000..43f1a6c4
--- /dev/null
+++ b/zUtil_Deformation/GtGrainODFuvwSfSolver.m
@@ -0,0 +1,174 @@
+classdef GtGrainODFuvwSfSolver < GtGrainODFuvwSolver
+ properties
+ end
+
+ methods (Access = public)
+ function self = GtGrainODFuvwSfSolver(parameters, varargin)
+ self = self@GtGrainODFuvwSolver(parameters, varargin{:});
+ end
+ end
+
+ methods (Access = public)
+ function build_orientation_sampling(self, ref_gr, det_index, oversize, oversampling, shape_functions)
+ self.build_orientation_sampling@GtGrainODFuvwSolver(ref_gr, det_index, oversize, oversampling, shape_functions);
+
+ shape_functions_name = sprintf('tmp_shape_functions_uvw_%d_%d', ...
+ size(self.voxel_centers, 1), self.sampler.get_total_num_orientations());
+
+ if (evalin('base', ['exist(''' shape_functions_name ''', ''var'')']))
+ fprintf('Loading shape functions..')
+ c = tic();
+ self.shape_functions = evalin('base', shape_functions_name);
+ fprintf('\b\b: Done in %g seconds\n', toc(c));
+ elseif (~exist([shape_functions_name '.mat'], 'file'))
+ switch (self.shape_functions_type)
+ case 'uvw'
+ for ii_v = 1:size(self.voxel_centers, 1);
+ disp_sampler = self.sampler.regenerate_displaced(self.voxel_centers(ii_v, :));
+
+ self.shape_functions{ii_v} = gtDefCreateShapeFuncs(disp_sampler, 15, [], false);
+ end
+ end
+ sfs = self.shape_functions;
+ save([shape_functions_name '.mat'], 'sfs', '-v7.3')
+ assignin('base', shape_functions_name, sfs);
+ else
+ fprintf('Loading shape functions..')
+ c = tic();
+ sfs_file = load([shape_functions_name '.mat']);
+ self.shape_functions = sfs_file.sfs;
+ assignin('base', shape_functions_name, sfs_file.sfs);
+ fprintf('\b\b: Done in %g seconds\n', toc(c));
+ end
+ end
+
+ function build_projection_matrices_sf_uvw_nooo(self)
+ ref_sel = self.sampler.selected;
+
+ grid_gr = self.sampler.get_orientations();
+ num_orients = numel(grid_gr);
+
+ num_voxels = size(self.voxel_centers, 1);
+ self.S = cell(num_voxels, 1);
+ self.St = cell(num_voxels, 1);
+
+ bls = self.sampler.bl(ref_sel);
+ bls_uvw_orig = cat(1, bls(:).bbuim, bls(:).bbvim, bls(:).bbwim);
+ bls_uvw_orig = reshape(bls_uvw_orig, [], 3, 2);
+ bls_uvw_orig = bls_uvw_orig(:, :, 1);
+ bls_uvw_orig(:, 3) = bls_uvw_orig(:, 3) - self.pre_paddings + 1;
+
+ for ii_v = 1:num_voxels
+ num_chars = fprintf('%03d/%03d', ii_v, num_voxels);
+
+ sf_uvw_shift = cell(num_orients, 1);
+ for ii_o = 1:num_orients
+ sf_uvw_shift{ii_o} = cat(1, ...
+ self.shape_functions{ii_v}{ii_o}(:).bbuim, ...
+ self.shape_functions{ii_v}{ii_o}(:).bbvim, ...
+ self.shape_functions{ii_v}{ii_o}(:).bbwim ...
+ );
+ sf_uvw_shift{ii_o} = reshape(sf_uvw_shift{ii_o}, [], 3, 2);
+ sf_uvw_shift{ii_o} = sf_uvw_shift{ii_o}(:, :, 1) - bls_uvw_orig;
+ end
+ sf_uvw_shift = cat(3, sf_uvw_shift{:});
+
+ self.S{ii_v} = sf_uvw_shift;
+
+ fprintf(repmat('\b', [1 num_chars]));
+ end
+ end
+
+ function build_projection_matrices_sf_uvw(self)
+ ref_sel = self.sampler.selected;
+
+ bls = self.sampler.bl(ref_sel);
+ num_blobs = self.get_num_blobs();
+
+ num_ws = self.get_num_ws();
+ bls_bbws = cat(1, bls(:).bbwim);
+ min_w_conds = bls_bbws(:, 1) - self.pre_paddings + 1;
+ max_w_conds = min_w_conds + num_ws - 1;
+
+ bls_bbus = cat(1, bls(:).bbuim);
+ bls_bbvs = cat(1, bls(:).bbvim);
+
+ grid_gr = self.sampler.get_orientations();
+ num_orients = numel(grid_gr);
+
+ num_voxels = size(self.voxel_centers, 1);
+ self.S = cell(num_voxels, 1);
+ self.St = cell(num_voxels, 1);
+
+ bls_uvw_orig = cat(1, bls(:).bbuim, bls(:).bbvim, bls(:).bbwim);
+ bls_uvw_orig = reshape(bls_uvw_orig, [], 3, 2);
+ bls_uvw_orig = bls_uvw_orig(:, :, 1);
+ bls_uvw_orig(:, 3) = bls_uvw_orig(:, 3) - self.pre_paddings + 1;
+
+ for ii_v = 1:num_voxels
+ b_us = [];
+ b_vs = [];
+ b_ws = [];
+ b_cs = [];
+ b_is = [];
+ b_os = [];
+
+ num_chars = fprintf('%03d/%03d', ii_v, num_voxels);
+
+ for ii_o = 1:num_orients
+ sf = self.shape_functions{ii_v}{ii_o};
+ sf_uvw_shift = cat(1, sf(:).bbuim, sf(:).bbvim, sf(:).bbwim);
+ sf_uvw_shift = reshape(sf_uvw_shift, [], 3, 2);
+ sf_uvw_shift = sf_uvw_shift(:, :, 1) - bls_uvw_orig;
+
+ uvw_cs = cell(num_blobs, 1);
+ bl_is = cell(num_blobs, 1);
+ bl_shifts = cell(num_blobs, 1);
+
+ for ii_b = 1:num_blobs
+ inds = find(sf(ii_b).intm);
+ uvw_cs{ii_b} = sf(ii_b).intm(inds);
+
+% [inds_u, inds_v, inds_w] = ind2sub(sf(ii_b).bbsize, inds);
+ inds_u = mod(inds, sf(ii_b).bbsize(1));
+ inds_v = mod(floor(inds / sf(ii_b).bbsize(1)) + 1, sf(ii_b).bbsize(2));
+ inds_w = floor(inds / prod(sf(ii_b).bbsize(1:2))) + 1;
+
+ bl_is{ii_b} = ii_b(ones(numel(inds), 1));
+ bl_shifts{ii_b} = [inds_u, inds_v, inds_w] + sf_uvw_shift(bl_is{ii_b}, :);
+ end
+
+ uvw_cs = cat(1, uvw_cs{:});
+ bl_is = cat(1, bl_is{:});
+ bl_shifts = cat(1, bl_shifts{:});
+
+ b_us = [b_us; bl_shifts(:, 1)];
+ b_vs = [b_vs; bl_shifts(:, 2)];
+ b_ws = [b_ws; bl_shifts(:, 3)];
+
+ b_cs = [b_cs; uvw_cs];
+
+ b_is = [b_is; bl_is];
+ b_os = [b_os; ii_o(ones(numel(bl_is), 1))];
+ end
+
+% sino_indx = sub2ind(self.size_sino, b_us, b_vs, b_ws, b_is);
+ sino_indx = b_us ...
+ + self.size_sino(1) * (b_vs - 1 ...
+ + self.size_sino(2) * (b_ws - 1 ...
+ + self.size_sino(3) * (b_is - 1)));
+
+ self.S{ii_v} = sparse( ...
+ sino_indx, b_os, double(b_cs), ...
+ numel(self.sino), num_orients);
+
+ self.St{ii_v} = self.S{ii_v}';
+
+ fprintf(repmat('\b', [1 num_chars]));
+ end
+ end
+ end
+
+ methods (Access = protected)
+ end
+end
\ No newline at end of file
diff --git a/zUtil_Deformation/GtGrainODFwSolver.m b/zUtil_Deformation/GtGrainODFwSolver.m
index 71138aa0..42379bdd 100644
--- a/zUtil_Deformation/GtGrainODFwSolver.m
+++ b/zUtil_Deformation/GtGrainODFwSolver.m
@@ -24,7 +24,7 @@ classdef GtGrainODFwSolver < handle
verbose = false;
- data_type = 'single';
+ data_type = 'double';
end
methods (Access = public)
--
GitLab