gtForwardSimulate_v2.m

function out = gtForwardSimulate_v2(first, last, workingdirectory, phaseID, varargin)
% GTFORWARDSIMULATE_V2  Finds missing diffraction spots and updates grain structure
%     out = gtForwardSimulate_v2(first, last, workingdirectory, phaseID, varargin)
%     ----------------------------------------------------------------------------
%
%     INPUT:
%       first            = <int>    first grainid
%       last             = <int>    last grainid
%       workingdirectory = <string> working directory
%       phaseID          = <int>    phase number {1}
%
%     OPTIONAL INPUT (varargin):
%       bb_size         = <int>     size of search area  ([] will use default projection size)
%       assemble_figure = <logical> collects all difspots and assembles them into a fullimage
%       save_grain      = <logical> flag to save grain_###.mat {true}
%
%     OUTPUT:
%       out  = updated grain structure <struct>
%              fields:
%                - 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 003 14-06-2012 by LNervo
%       Use gtFsimDefaultParameters
%
%     Version 004 21/6/2012 by AKing
%       Trivial change to make segbb zeros(1,6)
%
%
%     SUB-FUNCTIONS:
%[sub]- sfGetRoi


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

cd(workingdirectory);
gtDBConnect();

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

if isfield(parameters,'fsim')
    fsim = parameters.fsim;
else
    fsim = gtFsimDefaultParameters();
end

if ~isfield(fsim, 'save_grain')
    fsim.save_grain = true;
end

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

if isdeployed
    first   = str2double(first);
    last    = str2double(last);
    phaseID = str2double(phaseID);
    fsim.display_figure = false;
else
    % in interactive mode we can modify parameters
    fsim = parse_pv_pairs(fsim, varargin);
end

phaseID_str = sprintf('phase_%02d',phaseID);
phases_num    = length(parameters.cryst);    % number of phases in the material

cryst         = parameters.cryst(phaseID);
grainfilename = fullfile(acq.dir, '4_grains', phaseID_str, 'index.mat');
fprintf('\nLoading indexing results...\n');
grains        = load(grainfilename);
grain_num     = length(grains.grain);

difspotfile = fullfile(acq.dir, '2_difspot', '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)

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


fprintf('\nLoading reflections from database...\n');
difspottable = [parameters.acq.name 'difspot'];

% Load difspot data in 'allspots'  /  spotpair data in 'allpairs'
mymcmd = ['SELECT difspotID, CentroidX, CentroidY, CentroidImage, '...
          'BoundingBoxXsize, BoundingBoxYsize, BoundingBoxXorigin, BoundingBoxYorigin '...
          'FROM ' difspottable];

[difspotIDs, CentroidXs, CentroidYs, CentroidImages, BoundingBoxXsizes, ...
    BoundingBoxYsizes, BoundingBoxXorigins, BoundingBoxYorigins] = mym(mymcmd);

numspots           = max(difspotIDs);
difspotID          = NaN(numspots, 1);
CentroidX          = NaN(numspots, 1);
CentroidY          = NaN(numspots, 1);
CentroidImage      = NaN(numspots, 1);
BoundingBoxXsize   = NaN(numspots, 1);
BoundingBoxYsize   = NaN(numspots, 1);
BoundingBoxXorigin = NaN(numspots, 1);
BoundingBoxYorigin = NaN(numspots, 1);

for ii = 1: length(difspotIDs)
    jj = difspotIDs(ii);
    difspotID(jj)          = jj; 
    CentroidX(jj)          = CentroidXs(ii);
    CentroidY(jj)          = CentroidYs(ii);
    CentroidImage(jj)      = CentroidImages(ii);
    BoundingBoxXsize(jj)   = BoundingBoxXsizes(ii);
    BoundingBoxYsize(jj)   = BoundingBoxYsizes(ii);
    BoundingBoxXorigin(jj) = BoundingBoxXorigins(ii);
    BoundingBoxYorigin(jj) = BoundingBoxYorigins(ii);
end

    
mymcmd    = ['SELECT ldirX, ldirY, ldirZ, omega, omegaB, difAID, difBID FROM ' acq.pair_tablename];
[allpairs.ldirX, allpairs.ldirY, allpairs.ldirZ, allpairs.omegaA, allpairs.omegaB, allpairs.difAID, allpairs.difBID] = mym(mymcmd);

difspot_num = length(difspotIDs);   % number of difspots in database
spot2grain  = cell(difspot_num, 1); % cell array linking difspots with grainids
spot2phase  = cell(difspot_num, 1);

filetable   = [acq.name '_filetable'];
filename    = fullfile(acq.dir, '4_grains', 'spot2grain.mat');
spotTable   = GtMergeTable(filetable, filename, {'spot2grain', 'spot2phase'});

samplefile  = fullfile(acq.dir, '4_grains', 'sample.mat');
sample      = GtSample.loadFromFile(samplefile);

% initialise sample structure
phase = GtPhase(parameters.cryst(phaseID).name, grain_num);
sample.phases{phaseID} = phase;


for n = first : last
    fprintf('\nRunning forward simulation for grain %d...\n',n);
    
    gr = grains.grain{n};     %  output from INDEXTER
    
    graincenter = gr.center ./ pixelsize;    % graincenter contains COM coordinates in pixels 

    %%% Predict spot positions
    gr     = gtCalculateGrain(gr, parameters);
    u      = gr.allblobs.uv(:, 1);
    v      = gr.allblobs.uv(:, 2);
    % radial distance of spot center from detector center
    radius = sqrt((u-detrefu).^2 + (v-detrefv).^2);

    %%% 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    = gr.stat.bbxsmean; % average X Size of difspots
    dif_Ysize    = gr.stat.bbysmean; % average Y Size of difspots

    dif_XsizeMin = max(5, dif_Xsize - fsim.bbsize_factor * gr.stat.bbxsstd);  % limit to 5 pixels minimum width
    dif_XsizeMax = dif_Xsize + fsim.bbsize_factor * gr.stat.bbxsstd;
    dif_YsizeMin = max(5, dif_Ysize - fsim.bbsize_factor * gr.stat.bbysstd);
    dif_YsizeMax = dif_Ysize + fsim.bbsize_factor * gr.stat.bbysstd;

    % allow for a lateral deviation of spot postions based on grain orientation statics
    %dif_Offset   = fsim.Rdist_factor * gr.stat.Rdiststd * norm(detrefpos)
    dif_Offset    = fsim.Rdist_factor * deg2rad(gr.stat.dangstd) * norm(detrefpos);

    grainBBox = [max(gr.bbxs), max(gr.bbys)];
    % take the biggest spot and add a relative extra margin (multiplicative factor oversize)
    hsize = round(grainBBox(1) * fsim.oversize);
    % ...same for vertical direction
    vsize = round(grainBBox(2) * fsim.oversize);

    %%% Get the indices of all spots which fully intersect the active region of the detector
    ondet = find(all([u > dif_Xsize/2, u < acq.xdet - dif_Xsize/2 , ...
                      v > dif_Ysize/2 , v < acq.ydet - dif_Ysize/2, ...
                      radius < parameters.acq.maxradius], 2));

    % Exclude spots which overlap with segmentation boundingbox
    seg = find(all([( (u + dif_Xsize/2 > parameters.seg.bbox(1)) ...
                    & (u - dif_Xsize/2 < parameters.seg.bbox(1) + parameters.seg.bbox(3))), ...
                    ( (v + dif_Ysize/2 > parameters.seg.bbox(2)) ...
                    & (v - dif_Ysize/2 < parameters.seg.bbox(2) + parameters.seg.bbox(4)))], 2));

    ondet    = setdiff(ondet, seg); % indices of spots intersection active area of detector
    numspots = length(ondet);       % number of exspected diffraction spots


    % 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'
    % 4 = 'Found segmented difspot in database matching the strict search criteria' - we will use theoretically predicted angles for reconstruction
    % 5 = 'Difspot is part of a pair identified by INDEXTER' - use Friedelpair information for reconstruction

    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(gr.hkl,1));    % 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      = gr.allblobs.pl(currentSpotID, :);
        theta   = gr.allblobs.theta(currentSpotID);
        om_th   = gr.allblobs.omega(currentSpotID);
        image_n = round(om_th/omstep);
        eta     = gr.allblobs.eta(currentSpotID);
        hklsp(ii,:)  = gr.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(min(fsim.MaxOmegaOffset, (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
        
        
        
        image_ok      = find((CentroidImage >= image_n - thr_max_offset) & (CentroidImage <= image_n + thr_max_offset));
        intensity_ok  = find(abs(CentroidX(image_ok) - u(currentSpotID)) < max(dif_Xsize/2, BoundingBoxXsize(image_ok)/2 - dif_Xsize/2) & ...
                        abs(CentroidY(image_ok) - v(currentSpotID)) < max(dif_Ysize, BoundingBoxYsize(image_ok)/2 - dif_Ysize/2));
            
        candidateIDs   = image_ok(intensity_ok);
       
        
        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;

                if (intensity > parameters.seg.thr_grow_low)
                    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

        else % isempty(candidateIDs)
            % in case there are several candidates, take the one which has
            % the Area close to the expected Area.  This selection "rule" could be
            % improved...

            [area, area_index]   = sort(abs(BoundingBoxXsize(candidateIDs) .* BoundingBoxYsize(candidateIDs) - dif_Xsize*dif_Ysize));
            [image, image_index] = sort(abs(CentroidImage(candidateIDs) - image_n));
            candidateID  = candidateIDs(image_index(1));
            
            U(ii)  = CentroidX(candidateID);
            V(ii)  = CentroidY(candidateID);
            Z      = CentroidImage(candidateID);
            bb(ii, :)  = [BoundingBoxXorigin(candidateID), BoundingBoxYorigin(candidateID), ...
                            BoundingBoxXsize(candidateID), BoundingBoxYsize(candidateID)];
            spotid(ii) = candidateID;
            
            flag(ii) = 3; % Found a segmented difspot at expected position - we will check if it does match strict criteria
            info.dim_1 = bb(ii, 3);
            info.dim_2 = bb(ii, 4);
            spot = gtGetSummedDifSpot(candidateID, parameters, 1, info);
            spot = spot / sum(spot(:)) * dif_Xsize * dif_Xsize * dif_Ysize;   % here we force the intensity to be the same in all spots -

            if (fsim.assemble_figure)
                full = gtPlaceSubImage2(spot, full, bb(ii, 1), bb(ii, 2), 'summed');
            end
             
            cent_ok    = sqrt((U(ii) - u(currentSpotID)).^2 + (V(ii) - v(currentSpotID)).^2) <= dif_Offset;
            size_ok    = all([between(bb(ii,3), dif_XsizeMin, dif_XsizeMax), between(bb(ii,4), dif_YsizeMin, dif_YsizeMax)]);
            
            
            check(ii,:)  = [cent_ok, size_ok];
            
       
            if all(check(ii,:)), flag(ii) = 4; end % otherwise we keep flag == 3

            if find(gr.difspotidA == candidateID)
                % If available, we want to use the angles as determined from the
                % Friedel pairs - interrogate spotpair table
                flag(ii) = 5; % This spot is part of a pair and has been used by INDEXTER              
                index = find(allpairs.difAID == candidateID);
                omega(ii) = allpairs.omegaA(index);
                r(ii,:) = [allpairs.ldirX(index), allpairs.ldirY(index), allpairs.ldirZ(index)];
                Rottens = gtMathsRotationTensor(-omega(ii), rotcomp);
            elseif find(gr.difspotidB == candidateID)
                flag(ii) = 5; % This spot is part of a pair and has been used by INDEXTER
                index = find(allpairs.difBID == candidateID);
                omega(ii) = allpairs.omegaB(index);  % in this case we have to reverse the line direction
                r(ii,:)   = -[allpairs.ldirX(index), allpairs.ldirY(index), allpairs.ldirZ(index)];
                Rottens = gtMathsRotationTensor(-omega(ii), rotcomp);
            else
                % we found a segmentedd spot -  can use the theoretically
                % predicted direction (r_th) or apparent direction (based on
                % COM of grain and spot) for backprojection

                % theoretically prediceted omega angle
                 
                
                if (fsim.use_th)
                    % 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;
                else
                    % use experimentally determined omega angles
                    omega(ii)  = Z*omstep;
                    Rottens    = gtMathsRotationTensor(-omega(ii), rotcomp);
                    % projection direction as inferred from grain / spot
                    % center positions
                    r(ii,:)    = gtGeoDiffVecInSample([U(ii), V(ii)], omega(ii), gr.center, labgeo, samgeo);
                end
                
           
            end

            
            hoffset = round((hsize - bb(ii,3)) / 2);
            voffset = round((vsize - bb(ii,4)) / 2);

            spot = gtPlaceSubImage2(spot, zeros(vsize, hsize), hoffset+1, ...
                                    voffset+1, 'summed');
            difstack(:, ii, :) = spot';

            % position of Detector
            detpos  = detposcorner ...
                      + (bb(ii,1) - hoffset - 1 + hsize/2) * detdiru0 ...
                      + (bb(ii,2) - voffset - 1 + vsize/2) * detdirv0;
            
            
            det_abs = detposcorner ...
                      + (acq.bb(1) - 1 + acq.bb(3)/2) * detdiru0 ...
                      + (acq.bb(2) - 1 + acq.bb(4)/2) * detdirv0;
                  
            detdiru = detdiru0 * Rottens;
            detdirv = detdirv0 * Rottens;
            % take into account the coordinate shift corresponding to the offset
            % of the reconstructed volume
            detpos  = detpos * Rottens - graincenter;
            det_abs = det_abs * Rottens;
            r_abs   = beamdir * Rottens;

            % Compact version, since they're all <1x3> vectors, so we can read
            % them easily
            proj_geom(ii, :)     = [r(ii,:), detpos, detdiru, detdirv];
            proj_geom_abs(ii, :) = [r_abs, det_abs, detdiru, detdirv];

            % calculate the projection geometry for the full images as well
            detposFull = detposcorner + detrefu*detdiru0 + detrefv*detdirv0;
            proj_geom_full(ii, :)   = proj_geom(ii,:);
            proj_geom_full(ii, 4:6) = (detposFull * Rottens);
        end % if isempty(candidateID)
    end % loop over predicted spot positions

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

    newspots = length(find(flag == 4));
    included = find(flag >= fsim.included);
    selected = find(flag(included) >= fsim.selected);  % 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', n);

    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));
    fprintf(['%d new spots have been detected in addition to the %d ' ...
             'spotpairs identified by INDEXTER\n'], newspots, gr.nof_pairs);

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

    
    out.selected  = flag(included) >= fsim.selected;  % vector of logicals with pre-selected spots active
    out.id        = gr.id;
    out.phaseid   = phaseID;
    out.center    = graincenter;
    out.ondet     = ondet;
    out.allblobs  = gr.allblobs;
    out.full      = full;
    out.flag      = flag;
    out.status    = status;
    out.check     = check;
    out.pairid    = gr.pairid;
    
    out.difspotID = spotid(included);
    out.omega     = omega(included);
    out.hklsp     = hklsp(included,:);
    out.uvw       = gr.allblobs.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;
    % We should in the future handle properly vertical detector (general geometry)
    out.proj.vol_size_x = round(grainBBox(1) * fsim.oversizeVol)
    out.proj.vol_size_y = round(grainBBox(1) * fsim.oversizeVol);
    out.proj.vol_size_z = round(grainBBox(2) * fsim.oversizeVol);
    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)-out.proj.vol_size_x)/2 + graincenter(1)) + 1;
    out.shift(2) = round((acq.bb(3)-out.proj.vol_size_y)/2 + graincenter(2)) + 1;
    out.shift(3) = round((acq.bb(4)-out.proj.vol_size_z)/2 + graincenter(3)) + 1;

    if fsim.display_figure
        out.full = full;
        gtShowFsim(out,fsim.Grayscale)
    end


    if fsim.save_grain
        filename = fullfile(acq.dir, '4_grains', phaseID_str, ...
                            sprintf('grain_%04d.mat', gr.id));
        save(filename, '-struct', 'out');
    end
    
    spot2grain(spotid(included)) = {n};
    spot2phase(spotid(included)) = {phaseID};

    sample.phases{phaseID}.setR_vector(n, gr.R_vector);
    sample.phases{phaseID}.setCompleteness(n, completeness);
    sample.phases{phaseID}.setCenter(n, graincenter);
    sample.phases{phaseID}.setSelectedDiffspots(n, out.selected);
    sample.phases{phaseID}.setDifspotID(n, spotid(included)); % fix bug 05/10/2012 LNervo
   
end

spotTable.fillTable(spot2grain, 'spot2grain');
spotTable.fillTable(spot2phase, 'spot2phase');

sample.mergeToFile(filetable, samplefile, 'sample');

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