diff --git a/5_reconstruction/gtReconstructGrainTwinClusters.m b/5_reconstruction/gtReconstructGrainTwinClusters.m
new file mode 100644
index 0000000000000000000000000000000000000000..b200c67edb73aaf3ab197411c352b69d1eea1bd5
--- /dev/null
+++ b/5_reconstruction/gtReconstructGrainTwinClusters.m
@@ -0,0 +1,185 @@
+function gtReconstructGrainTwinClusters(first, last, workingdirectory, phase_id, parameters, varargin)
+% gtReconstructGrainOrientation 6D reconstructions on a GPU machine
+% gtAstraReconstructGrainTwinClusters(first, last, workingdirectory, phaseID, [parameters], varargin)
+% -------------------------------------------------------
+ cd(workingdirectory);
+
+ if (~exist('parameters', 'var') || isempty(parameters))
+ parameters = gtLoadParameters();
+ end
+
+ parameters.fsim.check_spot = false;
+
+ conf = struct( ...
+ 'recompute_osbb', true, ...
+ 'ospace_resolution', [], ...
+ 'ospace_lims', [], ...
+ 'list', []);
+
+ [conf, ~] = parse_pv_pairs(conf, varargin);
+
+ if (isdeployed)
+ global GT_DB %#ok<NUSED,TLEV>
+ global GT_MATLAB_HOME %#ok<NUSED,TLEV>
+ load('workspaceGlobal.mat');
+ phase_id = str2double(phase_id);
+ first = str2double(first);
+ last = str2double(last);
+ end
+
+ sample = GtSample.loadFromFile;
+
+ if isempty(conf.list)
+ cluster_list = first : last;
+ else
+ cluster_list = conf.list;
+ % in case the code is deployed, we have to convert to strings
+ if ischar(cluster_list)
+ cluster_list = sscanf(cluster_list, '%d');
+ end
+ end
+
+ num_clusters = numel(cluster_list);
+
+ for ii_cl = 1 : num_clusters
+ cluster_id = cluster_list(ii_cl);
+ grain_ids = sample.phases{phase_id}.clusters(cluster_id).included_ids;
+
+ rec_opts = gtReconstruct6DGetParamenters(parameters);
+ if (isfield(parameters.rec, 'grains') ...
+ && isfield(parameters.rec.grains, 'options') ...
+ && ~isempty(parameters.rec.grains.options) ...
+ && (~isfield(parameters.rec.grains.options, 'use_predicted_scatter_ints') ...
+ || isempty(parameters.rec.grains.options.use_predicted_scatter_ints)) )
+ rec_opts.use_predicted_scatter_ints = true;
+ end
+
+ phase_dir = fullfile(parameters.acq.dir, '4_grains', ...
+ sprintf('phase_%02d', phase_id));
+ if (isstruct(grain_ids))
+ grs = grain_ids.samp_ors;
+ grains_str_ids = sprintf('_%04d', cat(2, grs.id));
+ else
+ grains_str_ids = sprintf('_%04d', grain_ids);
+ grain_cluster_file = fullfile(phase_dir, ...
+ sprintf('grain_cluster%s.mat', grains_str_ids));
+
+ fprintf('Loading the cluster file..')
+ grs = load(grain_cluster_file);
+ grs = grs.samp_ors;
+ fprintf('\b\b: Done.\n')
+ end
+ num_grains = numel(grs);
+
+ for ii_g = 1:num_grains
+ sampler(ii_g) = GtOrientationSampling(parameters, grs(ii_g)); %#ok<AGROW>
+
+ if (~isempty(conf.ospace_resolution))
+ if (numel(conf.ospace_resolution) == 1)
+ conf.ospace_resolution = conf.ospace_resolution(ones(num_grains, 1));
+ end
+ if (isempty(conf.ospace_lims) || (size(conf.ospace_lims, 1) < ii_g))
+ bb_gvdm = cat(1, grs(ii_g).bb_ors(:).R_vector);
+ else
+ diff_r_vecs = tand(conf.ospace_lims(ii_g, :) / 2);
+ bb_gvdm = [grs(ii_g).R_vector + diff_r_vecs(1:3); grs(ii_g).R_vector + diff_r_vecs(4:6)];
+ end
+ sampler(ii_g).make_res_simple_grid('cubic', conf.ospace_resolution(ii_g), bb_gvdm', rec_opts.ospace_oversize);
+ elseif (conf.recompute_osbb)
+ sampler(ii_g).make_simple_grid_estim_ODF('cubic', rec_opts.grid_edge, false, rec_opts.ospace_oversize);
+ else
+ bb_gvdm = cat(1, grs(ii_g).bb_ors(:).R_vector);
+ sampler(ii_g).make_even_simple_grid('cubic', rec_opts.grid_edge, bb_gvdm', rec_opts.ospace_oversize)
+ end
+ 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);
+
+ vols = algo.getCurrentSolution();
+
+ ranges = zeros(num_grains, 2);
+ ranges(1, 1) = 1;
+ for ii_g = 2:num_grains
+ ranges(ii_g, 1) = ranges(ii_g-1) + sampler(ii_g-1).get_total_num_orientations();
+ end
+ ranges(1:end-1, 2) = ranges(2:end, 1) - 1;
+ ranges(end, 2) = numel(vols);
+
+ for ii_g = num_grains:-1:1
+ or_sizes = sampler(ii_g).get_orientation_sampling_size();
+ gr_vols = vols(ranges(ii_g, 1):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{1});
+
+ [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', grs(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(grs(ii_g).proj, parameters, vol_size);
+
+ ODF6D(ii_g) = struct( ...
+ 'voxels_avg_R_vectors', {avg_R_vecs}, ...
+ 'intensity', {avg_R_vecs_int}, ...
+ 'shift', {shift}, ...
+ 'R_vectors', {sampler(ii_g).get_R_vectors()}, ...
+ 'voxels_stddev_R_vectors', {stddev_R_vecs}, ...
+ 'single_grain_ODF', {s_g_odf}, ...
+ 'single_grain_avg_R_vector', {avg_R_vec}, ...
+ 'kernel_average_misorientation', {kam}, ...
+ 'intra_granular_misorientation', {igm} );
+ end
+
+ cluster_rec_file = fullfile(phase_dir, ...
+ sprintf('grain_cluster_details%s.mat', grains_str_ids));
+ % Saving and cleaning at the same time
+ if (exist(cluster_rec_file, 'file'))
+ cl_rec = load(cluster_rec_file);
+ cl_rec.ODF6D = ODF6D;
+ else
+ cl_rec = struct('ODF6D', ODF6D); %#ok<NASGU>
+ end
+ save(cluster_rec_file, '-struct', 'cl_rec', '-v7.3');
+
+ if (algo.verbose)
+ clear('ODF6D')
+ [proj_blobs, proj_spots] = algo.getProjectionOfCurrentSolution();
+
+ % 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
+
+ for ii_g = num_grains:-1:1
+ gr_vols = vols(ranges(ii_g, 1):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)}, ...
+ 'compute_statistics', {algo.get_statistics()} );
+ end
+ 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
+ end
+end