gtTwinForwardSimulate.m

function grain = gtTwinForwardSimulate(phaseID, grain, twin_center, varargin)
% GTTWINFORWARDSIMULATE  Finds diffraction spots in the database relative to the twin variant
%     grain = gtTwinForwardSimulate(phaseID, grain, twin_center, varargin)
%     --------------------------------------------------------------------------------------
%
%     INPUT:
%       phaseID          = <int>    phase number {1}
%       grain            = <struct> grain structure of interest
%
%     OPTIONAL INPUT (varargin):
%       save_grain      = <logical> flag to save or not the twins_####.mat
%       variantIDs      = <double>  array of ID for variants {all}
%
%     OUTPUT:
%       grain  = cell structure with fields for each variant:
%                - difstack       <double XxYxN>  stack of normalized difspots
%                - proj_geom      <double Nx12>   input for grain reconstruction / forward projection using roi (hsize, vsize)
%                - proj_geom_full <double Nx12>   input for grain reconstruction / forward projection (full images)
%                - proj_geom_abs  <double Nx12>   input for grain reconstruction / forward projection of extinction spots
%                - ondet          <double N>      indices of those spots (from the list in allblobs) which are expected on the detector
%                - uvw            <double Nx3>    spot positions on the detector
%                - hklsp          <double Nx3(4)> signed hk(i)l of reflection
%                - flag           <int    Nx3>    status flag indicating result of forward simulation
%                - allblobs       <struct>        forward simulation output produced by gtFedGenerateRandomGrain
%                - status         <string 5>      explanation string corresponding to flag
%
%     Version 001 16-10-2012 by LNervo
%
%     SUB-FUNCTIONS:
%[sub]- sfGetRoi

if isdeployed
    global GT_DB
    global GT_MATLAB_HOME
    load('workspaceGlobal.mat');
end

gtDBConnect();

parameters = [];
load('parameters.mat');
acq    = parameters.acq;
samgeo = parameters.samgeo;
labgeo = parameters.labgeo;
fsim   = parameters.fsim;

app.save_grain = false;
app.variantIDs = [];

if ~exist('phaseID','var') || isempty(phaseID)
    phaseID = 1;
end
if ~exist('twin_center','var')
    twin_center = [];
end

if isdeployed
    phaseID = str2double(phaseID);
    fsim.display_figure = false;
else
    [app,rej_pars] = parse_pv_pairs(app, varargin);
end
fsim.check_spot = true;
fsim.display_figure = true;

difspotfile = fullfile(parameters.acq.dir, '2_difspot', '00000', 'difspot00001.edf');
if exist(difspotfile, 'file');
    info = edf_info(difspotfile);
else
    info = [];
end

%%% this section should be replaced with a function dealing with different
% versions of parameters names:  fed, v1, v2...
detrefu      = labgeo.detrefu;        % detector u (line direction) reference point (center) in lab
detrefv      = labgeo.detrefv;        % detector v (column direction )reference point (center) in lab
detdiru0     = labgeo.detdiru;        % will use detdiru0 to distinguish from rotated version of this vector
detdirv0     = labgeo.detdirv;
rotdir       = labgeo.rotdir;
beamdir      = labgeo.beamdir;

pixelsize    = mean([labgeo.pixelsizeu, labgeo.pixelsizev]);
detrefpos    = labgeo.detrefpos ./ pixelsize;   % detector reference position in pixels
detposcorner = gtGeoDetOrig(labgeo) ./ pixelsize;  % the (0,0) corner of the detector in Lab system (edge of detector pixel)
detrefpos_pix = parameters.labgeo.detrefpos ./ pixelsize;

rotcomp      = gtMathsRotationMatrixComp(rotdir,'row');
omstep       = 180 / acq.nproj;
imagerange   = round(fsim.omegarange/omstep);   % number of images summed up when searching raw data for a missing spot

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load info from database
fprintf('\nLoading reflections from database...\n');
segmentedSpots = gtDBLoadTable([parameters.acq.name 'difspot'], 'difspotID');
fprintf('Done!\n')

graincenter = grain.center ./ pixelsize;    % graincenter contains COM coordinates in pixels

spotsCommProps = getDiffractionSpotsCommonProperties(grain, parameters.fsim, detrefpos_pix, omstep);

% calc radius of u,v coordinate
getradius = @(u,v) sqrt((u - labgeo.detrefu).^2 + (v - labgeo.detrefv).^2);

if ~isempty(twin_center) && isRowVectorWithLength(twin_center,3)
    grain.center_parent = grain.center;
    grain.center = twin_center;
    fprintf('Using twin center (%f %f %f) for simulation of grain %d\n',grain.center,grain.id)
    fprintf('parent grain center (%f %f %f)\n',grain.center_parent)
end

if ~isfield(grain, 'g_twins') || ( isfield(grain, 'g_twins') )
    for kk = 1:length(grain.g_twins)
        grain = gtCalculateGrain(grain, parameters, kk, 'variants', true, varargin{:});
    end
end

if isempty(app.variantIDs)
    app.variantIDs = 1:length(grain.g_twins);
end

gr = grain;


% loop over the variants
for kk = app.variantIDs
    fprintf('\nRunning forward simulation for twin variant %d for grain %d ...\n',kk,grain.id);

    variants = grain.g_twins(kk);
    % assume same center for parent and twin
    graincenter = grain.center ./ pixelsize;    % graincenter contains COM coordinates in pixels
    out.g        = variants.g;
    out.R_vector = variants.R_vector;

    %%% Get the average difspot X and Y sizes for the current spots in the grain,
    %  define theshold values which will be used during the database search later

    dif_Xsize    = mean(segmentedSpots.BoundingBoxXsize);
    dif_Ysize    = mean(segmentedSpots.BoundingBoxYsize);

    % take the biggest spot and add a relative extra margin (multiplicative factor oversize)
    hsize = round(max(segmentedSpots.BoundingBoxXsize)*fsim.oversize);
    % ...same for vertical direction
    vsize = round(max(segmentedSpots.BoundingBoxYsize)*fsim.oversize);

    u = variants.allblobs.detector(1).uvw(:, 1);
    v = variants.allblobs.detector(1).uvw(:, 2);
    uvw = variants.allblobs.detector(1).uvw;
    radiuses = getradius(uvw(:,1),uvw(:,2));
    ondet = findSpotsOnDetector(spotsCommProps, uvw(:,1), uvw(:,2), radiuses, parameters.acq, parameters.seg.bbox);

    uvw = uvw(ondet, :);
    % number of exspected diffraction spots
    numspots = numel(ondet);
    grain.spotid = zeros(numspots,1);
    % We will classify spots with a flag  [1...5]:
    % 1 = 'No intensity at expected spot position'
    % 2 = 'Found intensity in FULL images at predicted position'
    % 3 = 'Found segmented difspot at expected position - but position / size do not match strict criteria'

    flag      = zeros(numspots, 1);                % see above
    difstack  = zeros(hsize, numspots, vsize);     % stack of zero padded diffraction spots (projections used for 3D reconstruction)
    bb        = zeros(numspots, 4);                % difspot BoundingBox from database
    status    = cell(numspots, 1);                 % string corresponding to flag (summarizing result of forward simulation)
    spotid    = zeros(numspots, 1);                % Mysql difspottable difspotID
    check     = zeros(numspots, 2);
    hklsp     = zeros(numspots,size(variants.allblobs.hklsp,2));    % signed hk(i)l indices of reflection
    omega     = zeros(numspots, 1);                % rotation angle in degrees
    proj_geom = zeros(numspots, 12);               % Vector used for ROI grain reconstruction in ASTRA (grain shifted to center of roi volume)
    proj_geom_full = zeros(numspots, 12);          % Vector describing full projection geometry (full images, grain at nominal position in sample volume)
    proj_geom_abs  = zeros(numspots, 12);          % Vector describung absorption image projection geometry

    r         = zeros(numspots, 3);                % direction of diffracted beam used in ASTRA
    r_th      = zeros(numspots, 3);                % theoretically predicted direction of diffracted beam (based on R_vector)
    r_ex      = zeros(numspots, 3);                % estimated direction of diffracted beam, based on grain & difspot COM positions
    U         = zeros(numspots, 1);                % COM u positions of difspots from difspottable
    V         = zeros(numspots, 1);                % COM v positions of difspots from difspottable
    full      = zeros(acq.xdet, acq.ydet);


    %%% Main loop:  examine all hklsp and search corresponding spots / intensity on the detector..
    for ii = 1 : numspots
        currentSpotID = ondet(ii);
        pl      = variants.allblobs.pl(currentSpotID, :);
        theta   = variants.allblobs.theta(currentSpotID);
        om_th   = variants.allblobs.omega(currentSpotID);
        image_n = round(om_th/omstep);
        eta     = variants.allblobs.eta(currentSpotID);
        hklsp(ii,:)  = variants.allblobs.hklsp(currentSpotID, :);

        % Maximum offset (in number of images) from the predicted omega position
        % allow for a 1/sin(eta) dependency of this value - limited to a
        % maximum angular offset of fsim.MaxOmegaOffset degrees: this could be
        % improved...
        thr_max_offset = round( fsim.omegarange/(abs(sind(eta))) /omstep );

        % Look for difspots appying loose search criteria - is there any intensity at the predicted spot
        % position ?  Choose the one which is closest to the expected size

        %if isempty(candidateIDs)
        % no candidates have been found: segmentation / overlap / intensity / ... problem ?

        flag(ii) = 1; % No intensity at expected spot position

        if fsim.check_spot
            % in case bb_size for search BoundingBox (used for check of
            % completeness) is empty, use the same size as the difspot stack
            if isempty(fsim.bb_size)
                bbhsize = hsize;
                bbvsize = vsize;
            else
                bbhsize = fsim.bb_size(1);
                bbvsize = fsim.bb_size(2);
            end

            % Adjust search bbox in order to avoid reading information
            % stemming from regions outside the image border
            bbox_ustart = round(u(currentSpotID) - bbhsize/2);
            bbox_vstart = round(v(currentSpotID) - bbvsize/2);
            bbox_uend   = min(acq.xdet, bbox_ustart + bbhsize);
            bbox_vend   = min(acq.ydet, bbox_vstart + bbvsize);
            bbox(1) = max(1, bbox_ustart);
            bbox(2) = max(1, bbox_vstart);
            bbox(3) = (bbox_uend - bbox_ustart + 1);
            bbox(4) = (bbox_vend - bbox_vstart + 1);
            start_image = max(0, image_n - imagerange);
            end_image   = min(image_n + imagerange, 2*acq.nproj - 1);

            % read in small image stack at predicted position read ROI from
            % raw images - m is the mean value for each image in the stack
            [roi, m] = sfGetRoi(start_image, end_image, acq, bbox);
            [intensity, ~] = max(m);
            %test = sum(roi,3)/sum(roi(:))*dif_Xsize*dif_Xsize*dif_Ysize;
            test = sum(roi,3);
            if (intensity > 0.2)
                flag(ii) = 2; % Found intensity in FULL images at predicted position
            end

            if (fsim.assemble_figure)
                full = gtPlaceSubImage2(test, full, bbox(1), bbox(2), 'summed');
            end % app.assemble_figure
        end % app.check_spot

        % use theoretically predited directions (default)
        omega(ii)  = om_th;
        Rottens    = gtMathsRotationTensor(-omega(ii), rotcomp);
        % theoretically predicted beam direction
        r(ii,:)    = beamdir*Rottens + 2*sind(theta)*pl;

    end % loop over predicted spot positions

    found_reflections  = length(find(flag > 1));
    completeness       = found_reflections / numspots;

    newspots = length(find(flag == 2));
    included = find(flag == 2);
    selected = find(flag(included) == 3);  % indices of the selected spots in list of included spots

    if (~fsim.check_spot)
        fprintf(['Completeness estimates may be too pessimistic - activate ' ...
            'option ''check_spot'' for more reliable results\n'])
    end

    fprintf('\n Grain number %d:\n', grain.id);

    fprintf(['Forward simulation found intensity at %d out of %d expected ' ...
        'spot positions \n'], found_reflections, numspots);
    fprintf('Completeness estimate:  %f\n', completeness);
    fprintf('%d spots will be used for reconstruction\n', length(selected));


    %%% Now fill in the information for the output
    %  remove reflections for which no segmented spot has been found

    out.included  = flag(included) == 3;
    out.selected  = flag(included) == 4;%>= fsim.selected;  % vector of logicals with pre-selected spots active
    out.id        = grain.id;
    out.parent_R_vector  = grain.R_vector;
    out.phaseid   = phaseID;
    out.center    = graincenter*pixelsize;
    out.ondet     = ondet;
    out.allblobs  = variants.allblobs;
    out.full      = full;
    out.flag      = flag;
    out.status    = status;
    out.check     = check;
    out.spotid = grain.spotid;

    out.difspotID = spotid(included);
    out.omega     = omega(included);
    out.hklsp     = hklsp(included,:);
    out.uvw       = variants.allblobs.detector(1).uvw(ondet(included),:);
    out.uv_exp    = [U(included), V(included)];
    out.bb        = bb(included,:);
    out.r_th      = r_th(included,:);

    out.completeness    = completeness;
    out.proj.stack      = difstack(:, included, :);
    out.proj.geom       = proj_geom(included,:);
    out.proj.geom_full  = proj_geom_full(included,:);
    out.proj.geom_abs   = proj_geom_abs(included,:);
    out.proj.num_rows   = vsize;
    out.proj.num_cols   = hsize;
    out.proj.vol_size_x = round(hsize / sqrt(2));
    out.proj.vol_size_y = round(hsize / sqrt(2));
    out.proj.vol_size_z = round(vsize / sqrt(2));
    out.proj.num_iter   = parameters.rec.num_iter;
    out.vol             = [];
    out.segbb           = zeros(1,6);
    out.threshold       = 0;
    out.shift(1) = round((acq.bb(3)-hsize)/2 + graincenter(1)) + 1;
    out.shift(2) = round((acq.bb(3)-hsize)/2 + graincenter(2)) + 1;
    out.shift(3) = round((acq.bb(4)-vsize)/2 + graincenter(3)) + 1;

    if fsim.display_figure
        out.full = full;
        gtShowFsim(out,out.phaseid,'clims',[-200 1000],'f_title',sprintf('twin variant %d for grain %d ...\n',kk,grain.id))
    end

    grain.g_twins(kk).out = out;

end % loop over the variants
%
% if app.save_grain
%     filename = fullfile(acq.dir, '4_grains', sprintf('phase_%02d',phaseID), ...
%     sprintf('twins_%04d.mat', grain.id));
%     save(filename, 'grain');
% end

end % end of function

%%
% sub-functions

function [spot, m] = sfGetRoi(start_image, end_image, acq, bb)
% GTFORWARDSIMULATE/SFGETROI Sums a roi in raw images

nimages = end_image - start_image + 1;
spot    = zeros(bb(4), bb(3), nimages);
m       = zeros(nimages, 1);
centerx = round((bb(3)/2 -5) : (bb(3)/2 +5));
centery = round((bb(4)/2 -5) : (bb(4)/2 +5));

if all([bb(1) > 0, bb(2) > 0, bb(1)+bb(3) < acq.xdet, bb(2)+bb(4) < acq.ydet])
    fullImgsDir = fullfile(acq.dir, '1_preprocessing', 'full');
    indexes = start_image : end_image;
    filename = fullfile(fullImgsDir, sprintf('full%04d.edf', indexes(1)));
    info     = edf_info(filename);


    for img_n = 1:nimages
        fullImgName = fullfile(fullImgsDir, sprintf('full%04d.edf', indexes(img_n)));
        spot(:, :, img_n) = edf_read(fullImgName, bb, [], info);
        m(img_n) = gtImgMeanValue(spot(centery, centerx, img_n));
    end
end

end