function [result_viewer, result, post_result, algo] = gtReconstructTestCase(test_data, det_index, varargin)
if (~exist('det_index', 'var') || isempty(det_index))
det_index = 1;
end
conf = struct( ...
'do_plot', false, ...
'single_orientation', false, ...
'rspace_super_resolution', 1);
conf = parse_pv_pairs(conf, varargin);
rec_opts = gtReconstruct6DGetParamenters(test_data.parameters);
if (~isfield(test_data, 'gv'))
test_data.gv.dm = gtDefDmvol2Gvdm(test_data.dmvol);
end
gvdm = test_data.gv.dm(:, test_data.gv.used_ind);
% gvcs = test_data.gv.cs(:, test_data.gv.used_ind);
% if (~isfield(test_data, 'intvol'))
% size_int_vol = [ ...
% size(test_data.dmvol, 1), ...
% size(test_data.dmvol, 2), ...
% size(test_data.dmvol, 3) ];
% test_data.intvol = ones(size_int_vol, 'like', test_data.dmvol);
% end
% gvpow = test_data.intvol(test_data.gv.used_ind);
cryst = test_data.parameters.cryst(test_data.grain.phaseid);
symm = gtCrystGetSymmetryOperators(cryst.crystal_system, cryst.spacegroup);
gtGetMaxDisorientation(gvdm, symm, 'diameter_average');
gtGetMaxDisorientation(gvdm, symm, 'diameter_zero');
if (conf.do_plot)
f = figure();
ax = axes('parent', f);
hold(ax, 'on')
scatter3(ax, gvdm(1, 1:10:end)', gvdm(2, 1:10:end)', gvdm(3, 1:10:end)', 20);
end
if (conf.rspace_super_resolution > 1)
test_data.parameters.recgeo(det_index(1)).voxsize = ...
test_data.parameters.recgeo(det_index(1)).voxsize / conf.rspace_super_resolution;
test_data.grain.proj(det_index(1)).vol_size_x = ...
test_data.grain.proj(det_index(1)).vol_size_x * conf.rspace_super_resolution;
test_data.grain.proj(det_index(1)).vol_size_y = ...
test_data.grain.proj(det_index(1)).vol_size_y * conf.rspace_super_resolution;
test_data.grain.proj(det_index(1)).vol_size_z = ...
test_data.grain.proj(det_index(1)).vol_size_z * conf.rspace_super_resolution;
end
sampler = GtOrientationSampling(test_data.parameters, test_data.grain, 'detector_index', det_index(1));
if (conf.do_plot)
% Testing ability of GtOrientationSampling to detect the correct
% orientation-space boundingbox
sampler.make_grid_resolution(rec_opts.ospace_resolution, rec_opts.max_grid_edge_points);
scatter3(ax, sampler.R_vectors(:, 1), sampler.R_vectors(:, 2), sampler.R_vectors(:, 3), 20, 'r');
drawnow();
sampler = GtOrientationSampling(test_data.parameters, test_data.grain, 'detector_index', det_index(1));
end
is_topotomo = (numel(det_index) == 1) && strcmpi(test_data.parameters.acq(det_index).type, 'topotomo');
gvdm_lims = [min(gvdm, [], 2), max(gvdm, [], 2)];
gvdm_lims(:, 1) = gvdm_lims(:, 1) - tand(rec_opts.ospace_resolution / 4);
gvdm_lims(:, 2) = gvdm_lims(:, 2) + tand(rec_opts.ospace_resolution / 4);
if (is_topotomo)
gvdm_lims(3, :) = [-1, 1] * tand(rec_opts.ospace_resolution / 4);
rec_opts.max_grid_edge_points = [rec_opts.max_grid_edge_points([1 1]), 1];
% error_axes = 1:2;
else
% error_axes = 1:3;
end
for ii_d = 1:numel(det_index)
% Using the optimal orientation-space boundingbox
if (conf.single_orientation)
sampler.make_undeformed_sampling('det_ind', det_index(ii_d));
else
sampler.make_grid_resolution(rec_opts.ospace_resolution, ...
rec_opts.max_grid_edge_points, gvdm_lims, ...
'oversize', 1, 'det_ind', det_index(ii_d));
end
if (rec_opts.ospace_super_sampling > 1)
sampler.make_supersampling_simple_grid([1 2 3], rec_opts.ospace_super_sampling, det_index(ii_d));
end
end
if (conf.do_plot)
% ospace_bb = [min(gvdm, [], 2) - 0.002, max(gvdm, [], 2) + 0.002];
% sampler.make_grid('cubic', 3, ospace_bb, rec_opts.ospace_oversize);
scatter3(ax, sampler.R_vectors(:, 1), sampler.R_vectors(:, 2), sampler.R_vectors(:, 3), 20, 'g');
hold(ax, 'off')
drawnow();
end
if (conf.single_orientation)
detgeo = test_data.parameters.detgeo(det_index(1));
recgeo = test_data.parameters.recgeo(det_index(1));
samgeo = test_data.parameters.samgeo;
labgeo = test_data.parameters.labgeo;
proj = test_data.grain.proj(det_index(1));
ab = test_data.grain.allblobs(det_index(1));
proj.num_iter = rec_opts.num_iter;
proj.stack = proj.stack(:, proj.selected, :);
proj.num_cols = size(proj.stack, 1);
proj.num_rows = size(proj.stack, 3);
ab_sel = proj.ondet(proj.included(proj.selected));
dvecsam = ab.dvecsam(ab_sel, :);
rot_w_l2s = ab.srot(:, :, ab_sel);
bbuim = cat(1, proj.bl(proj.selected).bbuim);
bbvim = cat(1, proj.bl(proj.selected).bbvim);
bb_pos = [bbuim(:, 1), bbvim(:, 1), ...
bbuim(:, 2) - bbuim(:, 1) + 1, ...
bbvim(:, 2) - bbvim(:, 1) + 1, ];
volume_size = [proj.vol_size_y, proj.vol_size_x, proj.vol_size_z];
spacing = mean(volume_size) * (rec_opts.rspace_oversize - 1);
volume_size = ceil(volume_size + spacing);
proj.vol_size_y = volume_size(1);
proj.vol_size_x = volume_size(2);
proj.vol_size_z = volume_size(3);
proj.geom = gtGeoProjForReconstruction(dvecsam, ...
rot_w_l2s, test_data.grain.center, bb_pos, [], ...
detgeo, labgeo, samgeo, recgeo, 'ASTRA_grain');
vols = { gtAstra3D(proj, 'is_cone', false) };
else
algo = gtReconstruct6DLaunchAlgorithm(sampler, rec_opts, ...
test_data.parameters, 'det_index', det_index);
vols = algo.getCurrentSolution();
end
or_sizes = sampler.get_orientation_sampling_size();
[avg_R_vecs, avg_R_vecs_int, stddev_R_vecs] = sampler.getAverageOrientations(vols);
avg_R_vec = sampler.getAverageOrientation(vols);
s_g_odf = reshape(sampler.getODF(vols), or_sizes);
vol_size = size(avg_R_vecs_int);
shift = gtFwdSimComputeVolumeShifts(test_data.grain.proj, test_data.parameters, vol_size, det_index(1));
[kam, gam] = gtDefComputeKernelAverageMisorientation(avg_R_vecs, avg_R_vecs_int);
[igm, gos] = gtDefComputeIntraGranularMisorientation(avg_R_vecs, avg_R_vecs_int, 'R_vector', test_data.grain.R_vector);
% Restoring initial volume size (depending on the rounding)
if (rec_opts.volume_downscaling > 1)
fprintf('Expanding volumes..')
c = tic();
avg_R_vecs_int = gtMathsUpsampleVolume(avg_R_vecs_int, rec_opts.volume_downscaling);
avg_R_vecs = gtMathsUpsampleVolume(avg_R_vecs, rec_opts.volume_downscaling);
stddev_R_vecs = gtMathsUpsampleVolume(stddev_R_vecs, rec_opts.volume_downscaling);
kam = gtMathsUpsampleVolume(kam, rec_opts.volume_downscaling);
igm = gtMathsUpsampleVolume(igm, rec_opts.volume_downscaling);
vols = gtMathsUpsampleVolume(vols, rec_opts.volume_downscaling);
fprintf('\b\b: Done (%f seconds).\n', toc(c))
end
ODF6D = gtReconstructDataStructureDefinition('VOL6D');
ODF6D.options = rec_opts;
if (~conf.single_orientation)
ODF6D.compute_statistics = algo.get_statistics();
end
ODF6D.voxels_avg_R_vectors = avg_R_vecs;
ODF6D.intensity = avg_R_vecs_int;
ODF6D.shift = shift;
ODF6D.R_vectors = sampler.get_R_vectors();
ODF6D.voxels_stddev_R_vectors = stddev_R_vecs;
ODF6D.single_grain_ODF = s_g_odf;
ODF6D.single_grain_avg_R_vector = avg_R_vec;
ODF6D.kernel_average_misorientation = kam;
ODF6D.intra_granular_misorientation = igm;
result = gt6DAssembleResult(test_data, sampler.get_orientations(), vols, det_index);
result.ODF6D = ODF6D;
t = GtThreshold(test_data.parameters);
t.param.rec.thresholding.do_region_prop = false;
t.param.rec.thresholding.use_levelsets = true;
t.param.rec.thresholding.mask_border_voxels = 1;
gr_rec = struct('vol', avg_R_vecs_int, 'shift', shift);
result.SEG = t.thresholdAutoSingleGrain(gr_rec);
result.voxsize = test_data.parameters.recgeo(det_index(1)).voxsize;
post_result = gt6DPostProcessOrientationSpread(test_data, result, true);
% post_result = gt6DPostProcessOrientationSpread(test_data, result, true, error_axes);
if (~conf.single_orientation)
[proj_blobs, proj_spots] = algo.getProjectionOfCurrentSolution();
post_result.ODF6D = struct( ...
'orientation_volumes', {vols}, ...
'fwd_projected_blobs', {proj_blobs}, ...
'fwd_projected_spots', {proj_spots}, ...
'grain_average_misorientation', {gam}, ...
'grain_orientation_spread', {gos}, ...
'compute_statistics', {algo.get_statistics()} );
end
result_viewer = GtOrientationResultView(test_data, result, post_result);
end