function varargout = gtReconstructGrainCluster(grain_ids, phase_id, parameters, varargin)
% FUNCTION gtReconstructGrainCluster cluster 6D reconstructions on a GPU machine
% gtReconstructGrainCluster(grainIDs, phaseID, [parameters])
% -------------------------------------------------------
if (~exist('parameters', 'var') || isempty(parameters))
parameters = gtLoadParameters();
end
conf = struct( ...
'ospace_resolution', [], ...
'ospace_lims', [], ...
'extra_output', false, ...
'save', true );
[conf, ~] = parse_pv_pairs(conf, varargin);
rec_opts = gtReconstruct6DGetParamenters(parameters);
if (~isempty(conf.ospace_resolution))
rec_opts.ospace_resolution = conf.ospace_resolution;
end
% handling input, and making sure that things are called by their name
if (isstruct(grain_ids))
cl = grain_ids;
else
fprintf('Loading the cluster file..')
c = tic();
cl = gtLoadCluster(phase_id, grain_ids);
fprintf('\b\b: Done in %g seconds.\n', toc(c))
end
if (isfield(cl, 'samp_ors'))
% this means it is a twin cluster, but this is the only place where
% it matters
if (isfield(cl, 'blobs'))
blobs = cl.blobs;
else
error([mfilename 'wrong_argument'], ...
'Old twin cluster reconstructions not supported any more! Please, re-produce the twin cluster file')
end
cl = cl.samp_ors;
grain_ids = [cl(:).id];
else
grain_ids = cl.gids;
blobs = [];
end
num_sub_ors = numel(cl);
if ((num_sub_ors > 1) && (numel(rec_opts.ospace_resolution) == 1))
rec_opts.ospace_resolution(2:num_sub_ors) = rec_opts.ospace_resolution;
end
for ii_g = 1:num_sub_ors
sampler(ii_g) = GtOrientationSampling(parameters, cl(ii_g)); %#ok<AGROW>
if (isempty(conf.ospace_lims) || (size(conf.ospace_lims, 1) < ii_g))
bb_gvdm = cat(1, cl(ii_g).bb_ors(:).R_vector);
ospace_oversize = rec_opts.ospace_oversize;
else
diff_r_vecs = tand(conf.ospace_lims(ii_g, :) / 2);
bb_gvdm = [cl(ii_g).R_vector + diff_r_vecs(1:3); cl(ii_g).R_vector + diff_r_vecs(4:6)];
ospace_oversize = 1;
end
sampler(ii_g).make_grid_resolution(rec_opts.ospace_resolution(ii_g), ...
rec_opts.max_grid_edge_points, bb_gvdm', 'oversize', ospace_oversize);
if (rec_opts.ospace_super_sampling > 1)
sampler(ii_g).make_supersampling_simple_grid([1 2 3], ...
rec_opts.ospace_super_sampling);
end
end
algo = gtReconstruct6DLaunchAlgorithm(sampler, rec_opts, parameters, 'blobs', blobs);
vols = algo.getCurrentSolution();
ors_ranges = algo.orientation_groups;
gam = zeros(num_sub_ors, 1);
gos = zeros(num_sub_ors, 1);
ODF6D = gtReconstructDataStructureDefinition('VOL6D', num_sub_ors);
fprintf('Producing output data-structures: ')
c = tic();
for ii_g = 1:num_sub_ors
num_chars = fprintf('%03d/%03d', ii_g, num_sub_ors);
ODF6D(ii_g).options = rec_opts;
ODF6D(ii_g).compute_statistics = algo.get_statistics();
or_sizes = sampler(ii_g).get_orientation_sampling_size();
gr_vols = vols(ors_ranges(ii_g, 1):ors_ranges(ii_g, 2));
[avg_R_vecs, avg_R_vecs_int, stddev_R_vecs] = sampler(ii_g).getAverageOrientations(gr_vols);
avg_R_vec = sampler(ii_g).getAverageOrientation(gr_vols);
s_g_odf = reshape(sampler(ii_g).getODF(gr_vols), or_sizes);
[kam, gam(ii_g)] = gtDefComputeKernelAverageMisorientation(avg_R_vecs, avg_R_vecs_int);
[igm, gos(ii_g)] = gtDefComputeIntraGranularMisorientation(avg_R_vecs, avg_R_vecs_int, 'R_vector', cl(ii_g).R_vector);
% Restoring initial volume size (depending on the rounding)
if (rec_opts.volume_downscaling > 1)
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);
end
vol_size = size(avg_R_vecs_int);
shift = gtFwdSimComputeVolumeShifts(cl(ii_g).proj, parameters, vol_size);
ODF6D(ii_g).voxels_avg_R_vectors = avg_R_vecs;
ODF6D(ii_g).intensity = avg_R_vecs_int;
ODF6D(ii_g).shift = shift;
ODF6D(ii_g).R_vectors = sampler(ii_g).get_R_vectors();
ODF6D(ii_g).voxels_stddev_R_vectors = stddev_R_vecs;
ODF6D(ii_g).single_grain_ODF = s_g_odf;
ODF6D(ii_g).single_grain_avg_R_vector = avg_R_vec;
ODF6D(ii_g).kernel_average_misorientation = kam;
ODF6D(ii_g).intra_granular_misorientation = igm;
fprintf(repmat('\b', [1, num_chars]));
end
fprintf('Done (%d) in %g seconds.\n', num_sub_ors, toc(c))
if (conf.save)
fprintf('Saving the cluster reconstruction file..')
c = tic();
% Saving and cleaning at the same time
cl_rec = struct('ODF6D', ODF6D);
gtSaveClusterRec(phase_id, grain_ids, cl_rec, 'clean', true);
fprintf('\b\b: Done in %g seconds.\n', toc(c))
end
if (nargout > 0)
varargout{1} = ODF6D;
end
if (conf.extra_output)
clear('ODF6D')
% Restoring initial volume size (depending on the rounding)
if (rec_opts.volume_downscaling > 1)
fprintf('Expanding volumes..')
c = tic();
vols = gtMathsUpsampleVolume(vols, rec_opts.volume_downscaling);
fprintf('\b\b: Done (%f seconds).\n', toc(c))
end
current_solution = algo.getCurrentSolution();
for ii_g = num_sub_ors:-1:1
algo.currentSolution = gtMathsGetSameSizeZeros(current_solution);
algo.currentSolution(ors_ranges(ii_g, 1):ors_ranges(ii_g, 2)) = current_solution(ors_ranges(ii_g, 1):ors_ranges(ii_g, 2));
[proj_blobs, proj_spots] = algo.getProjectionOfCurrentSolution();
gr_vols = vols(ors_ranges(ii_g, 1):ors_ranges(ii_g, 2));
ODF6D(ii_g) = struct( ...
'orientation_volumes', {gr_vols}, ...
'fwd_projected_blobs', {proj_blobs}, ...
'fwd_projected_spots', {proj_spots}, ...
'grain_average_misorientation', {gam(ii_g)}, ...
'grain_orientation_spread', {gos(ii_g)} );
end
if (conf.save)
phase_dir = fullfile(parameters.acq.dir, '4_grains', ...
sprintf('phase_%02d', phase_id));
grains_str_ids = sprintf('_%04d', sort(grain_ids));
grain_full_details_file = fullfile(phase_dir, ...
sprintf('grain_cluster_full_details%s.mat', grains_str_ids));
save(grain_full_details_file, 'ODF6D', '-v7.3');
end
if (nargout > 1)
varargout{2} = ODF6D;
end
end
end