function [int_profile, refs] = gtGrainComputeIncomingBeamIntensity(gr, parameters, det_index, refs)
% FUNCTION [gr, refs] = gtGrainComputeIncomingBeamIntensity(gr, parameters, det_index, refs)

    if (parameters.acq.interlaced_turns > 0)
        error('gtComputeIncomingBeamIntensity:wrong_conf', ...
            'Interlaced turns not supported, yet!')
    end

    if (isfield(gr.proj, 'ondet'))
        gr_included = gr.proj.ondet(gr.proj.included);
    else
        gr_included = gr.ondet(gr.included);
    end

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

    omegas = gr.allblobs.omega(gr_included);
    gc_det_pos = gtGeoGrainCenterSam2Det(gr.center, omegas, parameters);

    img_nums = omegas ./ gtGetOmegaStepDeg(parameters, det_index);
    [inds, ints] = gtMathsGetInterpolationIndices([gc_det_pos(:, [2 1]), img_nums / parameters.acq.refon]);

    valid = ints > 1e-4;
    inds = inds(valid, :);
    ints = ints(valid);

    % Because image positions in array start from 1?
    inds(:, 3) = inds(:, 3) + 1;

    if (~exist('refs', 'var'))
        refs = load_reference_imgs(parameters);
    end

    indexes = sub2ind(size(refs), inds(:, 1), inds(:, 2), inds(:, 3));
    vals_valid = refs(indexes) .* ints;

    int_profile = zeros(size(valid), 'like', vals_valid);
    int_profile(valid) = vals_valid;
    int_profile = reshape(int_profile, [], 8);
    int_profile = sum(int_profile, 2);
end

function refs = load_reference_imgs(parameters)
    img_nums = 0:parameters.acq.refon:(parameters.acq.nproj*2);
    base_dir = fullfile('0_rawdata', parameters.acq.name);
    img_0 = fullfile(base_dir, sprintf('refHST%04d.edf', 0));
    info = edf_info(img_0);
    fprintf(' - loading ref images: ')
    c = tic();
    tot_images = numel(img_nums);
    for ii = tot_images:-1:1
        num_chars = fprintf('%02d/%02d', tot_images - ii + 1, tot_images);
        img_name = fullfile(base_dir, sprintf('refHST%04d.edf', img_nums(ii)));
        refs(:, :, ii) = edf_read(img_name, [], false, info);
        fprintf(repmat('\b', [1, num_chars]))
    end
    fprintf('(%02d) in %f seconds.\n', tot_images, toc(c));
end