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

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

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

    if (isfield(gr.proj, 'ondet'))
        gr_included = gr.proj(det_index).ondet(gr.proj(det_index).included);
    else
        gr_included = gr.ondet(gr.included);
    end
    
    if (~exist('radius', 'var') || isempty(radius))
        radius = 3;
    end

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

    img_nums = omegas ./ gtAcqGetOmegaStep(parameters, det_index);
    det_positions = [gc_det_pos, img_nums / parameters.acq(det_index).refon + 1];

    [disps_x, disps_y, disps_z] = ndgrid(-radius:radius, -radius:radius, 0);
    disps = [reshape(disps_x, 1, 1, []), reshape(disps_y, 1, 1, []), reshape(disps_z, 1, 1, [])];

    weights = exp( - (disps_x .^ 2 + disps_y .^ 2) ./ radius .^ 2 );
    weights = reshape(weights / (sum(weights(:))), 1, 1, []);

    det_positions_tot = bsxfun(@plus, det_positions, disps);

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

    int_profile = interp3(refs, det_positions_tot(:, 1, :), det_positions_tot(:, 2, :), det_positions_tot(:, 3, :), 'linear', 0);
    int_profile = sum(bsxfun(@times, int_profile, weights), 3);
end

function refs = load_reference_imgs(parameters, det_index)
    img_nums = 0:parameters.acq(det_index).refon:(parameters.acq(det_index).nproj*2);
    base_dir = fullfile('0_rawdata', parameters.acq(det_index).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));

    % At the moment, it seems to produce better results :|
    local_ints = sum(sum(refs, 1), 2) ./ prod([size(refs, 1), size(refs, 2)]);
    summed_refs = sum(refs, 3) ./ tot_images;
    refs = bsxfun(@times, summed_refs, local_ints);
end