Skip to content
Snippets Groups Projects
Commit a0af4c2e authored by Nicola Vigano's avatar Nicola Vigano
Browse files

Grain-ODF: added another solver that also can use the uvw shape functions 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent e060a5b4
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
zUtil_Deformation/GtGrainODFuvwPatchSolver.m
+ 57
− 83
View file @ a0af4c2e
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
......
zUtil_Deformation/GtGrainODFuvwSfSolver.m 0 → 100644
+ 174
− 0
View file @ a0af4c2e
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
zUtil_Deformation/GtGrainODFwSolver.m
+ 1
− 1
View file @ a0af4c2e
......@@ -24,7 +24,7 @@ classdef GtGrainODFwSolver < handle
verbose = false;
data_type = 'single';
data_type = 'double';
end
methods (Access = public)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment