gtDefSyntheticGrainCreate.m

function [grain, gv] = gtDefSyntheticGrainCreate(fedpars, grain, parameters, dmvol, intvol)
% [grain, gv] = gtDefSyntheticGrainCreate(fedpars, grain, parameters, dmvol, intvol)
%
% Sets up the framework, variables, random strain distribution for simulated data.
%  In the grain volume everything is defined in voxels, no dimensioning with pixel size.
%  In the projections and diff. blobs, everything is defined in detector pixels
%  (here comes the pixel size and setup geometry into play).
%    Thus the dimensioning/scaling happens in the forward projection
%    functions and their derivative functions, since they relate the
%    grain voxels to the detector pixels.
%
% INPUT
%   fedpars    - FED parameters
%                Some important ones:
%                  grainsize   - no. of voxels in grain
%                  grainenv    = grainsize
%                  ngv         = grainsize
%                  gvsize      = [1 1 1]
%                  loadbl      - blobs ID-s in grain.allblobs to be generated
%                  dcomps      - for rotations only use [1 1 1 0 0 0 0 0 0]
%                  blobsizeadd - blob volume padding
%                  vox000sam   - (X, Y, Z) coordinates of voxel (0, 0, 0) in the Sample reference frame
%   grain      - grain data as output from Indexter
%   parameters - parameters as in DCT paramaters file
%   dmvol      - all deformation components per voxel for the entire volume
%                (nx, ny, nz, ndef); if empty, random distribution is applied
%                with max limits as in fedpars.dapplim and max gradients
%                as in fedpars.dappgradlim
%   intvol     - Volume instensities

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Load and update geometry parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('Setting voxels up..')
c = tic();

labgeo = parameters.labgeo;
labgeo.omstep = gtGetOmegaStepDeg(parameters);
% Rotation matrix components:
labgeo.rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');

samgeo = parameters.samgeo;
recgeo = parameters.recgeo;
if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
    [parameters.detgeo, labgeo] = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
end
detgeo = parameters.detgeo;
% Update labgeo
parameters.labgeo = labgeo;

if (~isfield(parameters.fsim, 'oversizeVol') ...
        || isempty(parameters.fsim.oversizeVol))
    parameters.fsim.oversizeVol = 1.2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Load grain parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (~exist('dmvol', 'var') || isempty(dmvol))
    gr_size = fedpars.grainsize; % x, y, z number of voxels in grain
else
    gr_size = [size(dmvol, 1), size(dmvol, 2), size(dmvol, 3)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Volume sampling and volume elements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('\b\b: Done in %g seconds.\nAllocating memory for elements..', toc(c))
c = tic();

% total number of vol elements
nv = gr_size(1) * gr_size(2) * gr_size(3);
nv_ones = ones(nv, 1);

nd = sum(fedpars.dcomps);

% Element XYZ index
[gv.ind(:, 1), gv.ind(:, 2), gv.ind(:, 3)] = ind2sub(gr_size, 1:nv);
gv.ind = gv.ind';

% Element diffracting power
if (~exist('intvol', 'var') || isempty(intvol))
    intvol = ones(gr_size, fedpars.gvtype);
end

if (numel(intvol) ~= nv)
    error('gtFedTestLoadData_v2:wrong_argument', ...
        'The volume of voxel intensities should have the same size as dmvol''s first 3 coordinates')
end
gv.pow = reshape(intvol, 1, []);
if (strcmp(fedpars.gvtype, 'single'))
    gv.pow = single(gv.pow);
end

gv.used_ind = find(gv.pow > 0);

% Element actual strain state
gv.d = zeros(nd, nv, fedpars.gvtype);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Apply strain distribution to grain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('\b\b: Done in %g seconds.\n', toc(c))

disp('Active deformation components:')
disp(fedpars.dcomps)

disp('Applying random deformation..')
c = tic();

% Measured (true, real) strain distribution:
if (~exist('dmvol', 'var') || isempty(dmvol))
    [gvdm, dmvol, gradok] = gtFedTestGenerateRandomDef(fedpars, nd);
else
    [gvdm, dmvol_size] = gtDefDmvol2Gvdm(dmvol);
    if (strcmpi(fedpars.defmethod, 'small'))
        if (dmvol_size(4) == 3)
            mean_gvdm = mean(gvdm, 2);
            fprintf('Recentering GVDM and DMVOL to the mean orientation: (%f, %f, %f), from: (%f, %f, %f).\n', ...
                grain.R_vector, mean_gvdm)
%             gvdm_shift = grain.R_vector' - mean_gvdm;
            gvdm_shift = 0 - mean_gvdm;
            gvdm = gvdm + gvdm_shift(:, ones(size(gvdm, 2), 1));
            dmvol = gtDefGvdm2Dmvol(gvdm, dmvol_size);
        end
    end
    gradok = [];
end
gv.dm = gvdm;

fprintf('\b\b: Done in %g seconds.\nComputing voxel centers..', toc(c))
c = tic();

phantom_voxsize = recgeo(1).voxsize';
vox000sam = grain.center' - (gr_size' / 2) .* phantom_voxsize;

% Elements' positions in the sample volume in mm!
gv.cs = vox000sam(:, nv_ones) + (gv.ind - 0.5) .* phantom_voxsize(:, nv_ones);

% Producing centers and deformation components for projection
if (isfield(fedpars, 'volume_super_sampling') ...
        && (fedpars.volume_super_sampling > 1))
    [gv.pcs, gv.d] = compute_oversampling_voxels(gv, dmvol, fedpars, phantom_voxsize);
else
    gv.pcs = gv.cs;
    gv.d = gv.dm;
end

if (strcmp(fedpars.gvtype, 'single'))
    fprintf('\b\b: Done in %g seconds.\nConverting type to single..', toc(c))
    c = tic();
    fn = fieldnames(gv);
    for ii = 1:length(fn)
        gv.(fn{ii}) = single(gv.(fn{ii}));
    end
end

fprintf('\b\b: Done in %g seconds.\n', toc(c))

if (~all(gradok))
    disp('Gradients are higher than the limits for components:')
    disp(find(~gradok)')
end

disp('Largest deformation component:')
disp(max(dmvol(:)))
disp('Smallest deformation component:')
disp(min(dmvol(:)))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Computing allblobs info
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Cleaning grain struct
if (isfield(grain, 'id'))
    grain.id = 1;
end
grain = struct('id', grain.id, 'phaseid', grain.phaseid, ...
    'center', grain.center, 'R_vector', grain.R_vector);

grain = gtCalculateGrain(grain, parameters);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Iinitial blob parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

num_dets = numel(detgeo);
grain.proj(1:num_dets) = struct( ...
    'stack', [], 'num_cols', [], 'num_rows', [], 'bl', [], ...
    'vol_size_x', [], 'vol_size_y', [], 'vol_size_z', [], 'shift', [], ...
    'ondet', [], 'included', [], 'selected', [], 'centerpix', []);

for n = 1:num_dets
    % Checking if we should just project all of them, or wether the user
    % has something to say abou it.
    if (isfield(fedpars, 'detector') && isfield(fedpars.detector(n), 'ondet'))
        ondet = fedpars.detector(n).ondet;
        included = fedpars.detector(n).included;
        selected = fedpars.detector(n).selected;
    else
        ondet = find(grain.allblobs.detector(n).ondet);
        included = (1:numel(ondet))';
        selected = true(size(included));
    end

    ref_inc = ondet(included);
    nbl = numel(ref_inc);

    is_nearfield = norm(detgeo(n).detrefpos) < 90;
    % If it is farfield, we should adjust the voxel size to be 1/4 of the
    % pixel size
%     if (~is_nearfield)
%         voxsize = mean([detgeo(n).pixelsizeu, detgeo(n).pixelsizev]) / 4;
%         recgeo(n).voxsize = voxsize([1 1 1]);
%         parameters.recgeo(n) = recgeo(n);
%     end
    fprintf( ['Generating blob data (for detector: %d, blobs on the' ...
        ' detector: %d):\n'], n, nbl);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Simulate "measured" diffraction patterns
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    fprintf(' - Simulating true ("measured") blob intensities..')
    c = tic();
    bl = gtDefFwdProjGvdm2UVW(grain, ref_inc, gv, fedpars, parameters, n);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% Assembling the final proj structure
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    fprintf('\b\b (%f seconds)\n - Producing data-structure..\n', toc(c))
    c = tic();

    stack = arrayfun(@(x){sum(x.intm, 3)}, bl);
    stack = cat(3, stack{:});

    vol_size = gtGeoSam2Sam(gr_size, recgeo(1), recgeo(n), 1);
    vol_size = ceil(vol_size * parameters.fsim.oversizeVol);
    if (is_nearfield && size(intvol, 3) == 2) % This handles the 2D case
        vol_size(3) = max(vol_size(3), 8);
    end

    proj = gtFwdSimProjDefinition();

    proj.stack = permute(stack, [1 3 2]);
    proj.num_cols = size(proj.stack, 3);
    proj.num_rows = size(proj.stack, 1);
    proj.bl = bl;
    proj.vol_size_x = vol_size(2);
    proj.vol_size_y = vol_size(1);
    proj.vol_size_z = vol_size(3);
    proj.centerpix = gtGeoSam2Sam(grain.center, samgeo, recgeo(n), 0);
    proj.shift = gtFwdSimComputeVolumeShifts(grain.proj(n), parameters, vol_size);
    proj.ondet = ondet;
    proj.included = included;
    proj.selected = selected;

    grain.proj(n) = proj;

    fprintf('\b\b(%f seconds)\n - True blob dimensions:\n', toc(c))
    fprintf(' %3d) %4d %4d %4d\n', [(1:nbl)', vertcat(bl(:).mbbsize)]')
end

end

function [gvpcs, gvd] = compute_oversampling_voxels(gv, dmvol, fedpars, voxsize)
    nv = size(gv.ind, 2);
    nv_ones = ones(nv, 1);

    vssampling = fedpars.volume_super_sampling;

    num_subvoxels = vssampling ^ 3;
    gvpcs = zeros(3, nv, num_subvoxels, fedpars.gvtype);
    gvd = zeros(3, nv, num_subvoxels, fedpars.gvtype);

    dmvol_size = [size(dmvol, 1), size(dmvol, 2), size(dmvol, 3)]';
    dmvol_size = dmvol_size(:, nv_ones);

    orig_gvdm = gv.dm;

    counter = 0;
    corner = (1 - (1 / vssampling)) / 2 .* voxsize;

    fprintf('\b\b: ')

    for ss_z = linspace(-corner(3), corner(3), vssampling)
        for ss_y = linspace(-corner(2), corner(2), vssampling)
            for ss_x = linspace(-corner(1), corner(1), vssampling)

                num_chars = fprintf('%03d/%03d', counter, num_subvoxels);

                offset_subvoxel = [ss_x; ss_y; ss_z];

                gvpcs(:, :, counter) = gv.cs + offset_subvoxel(:, nv_ones);

                offset_pos = offset_subvoxel ./ voxsize;
                offset_pos = gv.ind + offset_pos(:, nv_ones);

                res_r_vecs = zeros(3, nv);
                res_ints = zeros(3, nv);

                for z_r = 0:1
                    for y_r = 0:1
                        for x_r = 0:1
                            rounding = [x_r, y_r, z_r]';
                            pos = floor(offset_pos) + rounding(:, nv_ones);
                            int_coeff = 1 - abs(offset_pos - pos);

                            valid = all(int_coeff > 0, 1) ...
                                & all(pos > 0, 1) ...
                                & all(pos <= dmvol_size, 1);

                            valid_x = pos(1, valid);
                            valid_y = pos(2, valid);
                            valid_z = pos(3, valid);

                            ind_valid_pos = valid_x ...
                                + (valid_y - 1) * size(dmvol, 1) ...
                                + (valid_z - 1) * size(dmvol, 1) * size(dmvol, 2);

                            % Filtering the positions that have no
                            % intensity
                            valid(valid) = valid(valid) & (gv.pow(ind_valid_pos) > 0);

                            valid_x = pos(1, valid);
                            valid_y = pos(2, valid);
                            valid_z = pos(3, valid);

                            ind_valid_pos = valid_x ...
                                + (valid_y - 1) * size(dmvol, 1) ...
                                + (valid_z - 1) * size(dmvol, 1) * size(dmvol, 2);

                            res_ints(:, valid) = res_ints(:, valid) + int_coeff(:, valid);

                            res_r_vecs(:, valid) = res_r_vecs(:, valid) ...
                                + int_coeff(:, valid) .* orig_gvdm(:, ind_valid_pos);
                        end
                    end
                end

                renorm_res_ints = res_ints + (res_ints == 0);
                gvd(:, :, counter) = res_r_vecs ./ renorm_res_ints;

                counter = counter + 1;

                fprintf(repmat('\b', [1 num_chars]));
            end
        end
    end
end