function [proj, verts] = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, omegas, bb, parameters, stackUSize, stackVSize, selected, det_ind)
% FUNCTION proj = gtFwdSimBuildProjGeometry(diff_beam_dirs, gr_center, omegas, bb, parameters, stackUSize, stackVSize, selected, det_ind)
%
    if (~exist('det_ind', 'var') || isempty(det_ind))
        det_ind = 1;
    end

    acq = parameters.acq;
    fsim = parameters.fsim;
    labgeo = parameters.labgeo;
    samgeo = parameters.samgeo;
    recgeo = parameters.recgeo(det_ind);
    detgeo = parameters.detgeo(det_ind);
    if (isfield(parameters, 'clone') & parameters.clone.active)
        cloning = true;
    else
        cloning = false;
    end
    num_projs = size(diff_beam_dirs, 1);

    spot_shifts = gtFwdSimGetStackShifts(stackUSize, stackVSize, bb, true);

    % Bouning boxes in detector coordinates!
    bbpos_det_grain = [ ...
        (bb(:, 1) - spot_shifts.u), (bb(:, 2) - spot_shifts.v), ...
        stackUSize(ones(num_projs, 1)), stackVSize(ones(num_projs, 1)) ...
        ];

    % Bouning boxes in detector coordinates!
    bbpos_det_full = repmat([1 1 acq.xdet acq.ydet], num_projs, 1);
    bbpos_det_abs = repmat(acq.bb, num_projs, 1);

    %%% Now we actually build the geometries

    % Calculate the projection geometry for the spots:
    %  - Vector used for ROI grain reconstruction in ASTRA (grain shifted
    %    to center of roi volume)
    [gr_shift_lab, gr_shift_sam] = gtMatchGetSampleShifts(parameters, omegas);
    proj_geom = gtGeoProjForReconstruction(diff_beam_dirs, ...
        omegas, gr_center, bbpos_det_grain, [], ...
        detgeo, labgeo, samgeo, recgeo, 'ASTRA_grain', gr_shift_sam);

    % Calculate the projection geometry for the full images:
    %  - Vector describing full projection geometry (full images, grain at
    %    nominal position in sample volume)
    proj_geom_full = gtGeoProjForReconstruction( ...
        diff_beam_dirs, omegas, [], bbpos_det_full, ...
        [], detgeo, labgeo, samgeo, recgeo, 'ASTRA_full', gr_shift_sam);

    % Geometry for extinction spots
    proj_geom_abs = gtGeoProjForReconstruction([], omegas, [], ...
        bbpos_det_abs, [], detgeo, labgeo, samgeo, recgeo,  ...
        'ASTRA_absorption', gr_shift_lab);

    proj = gtFwdSimProjDefinition('fwd_sim');

    % diffraction geometry
    proj.geom      = proj_geom;
    proj.geom_full = proj_geom_full;
    proj.geom_abs  = proj_geom_abs;

    proj.num_rows  = stackVSize;
    proj.num_cols  = stackUSize;

    proj.centerpix = gtGeoSam2Sam(gr_center, samgeo, recgeo, true);

    use_polyhedron = ~cloning & (numel(find(selected)) >= 3);
    if (use_polyhedron)
        % This should behave better with vertical detector
        % (if no bug is introduced with it :D)
        verts = gtFwdSimComputeCircumscribingPolyhedron(...
            gr_center, diff_beam_dirs(selected, :), omegas(selected), ...
            bb(selected, :), parameters, det_ind);

        vol_size = 2 * max(abs(max(verts, [], 1)), abs(min(verts, [], 1)));
    else
        verts = [];
        % We should in the future handle properly vertical detector
        % (general geometry) maybe determining a convex shape of the grain!
        vol_size = [stackUSize, stackUSize, stackVSize] / fsim.oversize;
    end

    vol_size = round(vol_size * fsim.oversizeVol);

    proj.vol_size_x = vol_size(2);
    proj.vol_size_y = vol_size(1);
    proj.vol_size_z = vol_size(3);
end