function atts = gtGrainComputeBeamAttenuation(gr, p, det_index, abs_vol, sampling_factor)

    if (~exist('sampling_factor', 'var') || isempty(sampling_factor))
        sampling_factor = 1;
    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

    rot_l2s = gr.allblobs.srot(:, :, gr_included);
    % They are reversed because we also use them to find the intersaction
    % with the volume border
    beam_dirs_in = - p.labgeo.beamdir * reshape(rot_l2s, 3, []);
    beam_dirs_in = permute(reshape(beam_dirs_in, 1, 3, []), [3 2 1]);

    beam_dirs_out = gr.allblobs.dvecsam(gr_included, :);

    gr_center = gr.proj(det_index).centerpix + size(abs_vol) / 2 + 1;

    % ray tracing
    dists_borders = cat(3, [1 1 1], size(abs_vol)) - gr_center(1, :, [1 1]);

    dists_renorm_in = compute_renormalized_distances(dists_borders, beam_dirs_in);
    dists_renorm_out = compute_renormalized_distances(dists_borders, beam_dirs_out);

    num_refl = numel(gr_included);
    atts = zeros(num_refl, 2);

    for ii = 1:num_refl
%         dist_in = dists_borders(d_inds_in_dir(d_inds_in_pm(ii)));
%         samp_points_in = linspace(0, dists_renorm_in(ii), ceil(abs(dist_in) * sampling_factor) );

        samp_factor_in = ceil(dists_renorm_in(ii) * sampling_factor);

        samp_points_in = linspace(0, dists_renorm_in(ii), samp_factor_in);
        samp_points_in = bsxfun(@times, beam_dirs_in(ii, :), samp_points_in');
        samp_points_in = bsxfun(@plus, samp_points_in, gr_center);

%         dist_out = dists_borders(d_inds_out_dir(d_inds_out_pm(ii)));
%         samp_points_out = linspace(0, dists_renorm_out(ii), ceil(abs(dist_out) * sampling_factor) );

        samp_factor_out = ceil(dists_renorm_out(ii) * sampling_factor);

        samp_points_out = linspace(0, dists_renorm_out(ii), samp_factor_out);
        samp_points_out = bsxfun(@times, beam_dirs_out(ii, :), samp_points_out');
        samp_points_out = bsxfun(@plus, samp_points_out, gr_center);

        [inds_in, ints_in] = gtMathsGetInterpolationIndices(samp_points_in);
        valid_in = ints_in > eps('single');
        inds_in = inds_in(valid_in, :);
        ints_in = ints_in(valid_in);

        [inds_out, ints_out] = gtMathsGetInterpolationIndices(samp_points_out);
        valid_out = ints_out > eps('single');
        inds_out = inds_out(valid_out, :);
        ints_out = ints_out(valid_out);

        inds_in = sub2ind(size(abs_vol), inds_in(:, 1), inds_in(:, 2), inds_in(:, 3));
        vals_in = abs_vol(inds_in) .* ints_in;

        inds_out = sub2ind(size(abs_vol), inds_out(:, 1), inds_out(:, 2), inds_out(:, 3));
        vals_out = abs_vol(inds_out) .* ints_out;

        abs_in = sum(vals_in) / samp_factor_in * dists_renorm_in(ii);
        abs_out = sum(vals_out) / samp_factor_out * dists_renorm_out(ii);

        atts(ii, :) = [abs_in, abs_out];

%         f = figure();
%         ax = axes('parent', f);
%         hold(ax, 'on')
%         scatter3(ax, samp_points_in(:, 1), samp_points_in(:, 2), samp_points_in(:, 3));
%         scatter3(ax, samp_points_out(:, 1), samp_points_out(:, 2), samp_points_out(:, 3));
%         show_bbox(ax, cat(1, [1 1 1], size(abs_vol)))
%         hold(ax, 'off')
%         pause
    end

    atts = exp(-atts);

%     figure, semilogy(ints)
%     figure, semilogy(prod(ints, 2))
end

function [dists_renorm, d_inds_pm, d_inds_dir] = compute_renormalized_distances(dists_borders, beam_dirs)
    dists_renorm = bsxfun(@times, dists_borders, 1 ./ beam_dirs);
    [dists_renorm, d_inds_pm] = max(dists_renorm, [], 3); % Takes negative out
    [dists_renorm, d_inds_dir] = min(dists_renorm, [], 2); % Takes shortest, and removes last remaining nans
end

function show_bbox(ax, borders)
    min_sampled_R_vecs = min(borders, [], 1);
    max_sampled_R_vecs = max(borders, [], 1);
    bbox_R_vecs = [ ...
        min_sampled_R_vecs(1), min_sampled_R_vecs(2), min_sampled_R_vecs(3); ...
        min_sampled_R_vecs(1), min_sampled_R_vecs(2), max_sampled_R_vecs(3); ...
        min_sampled_R_vecs(1), max_sampled_R_vecs(2), min_sampled_R_vecs(3); ...
        min_sampled_R_vecs(1), max_sampled_R_vecs(2), max_sampled_R_vecs(3); ...

        max_sampled_R_vecs(1), min_sampled_R_vecs(2), min_sampled_R_vecs(3); ...
        max_sampled_R_vecs(1), min_sampled_R_vecs(2), max_sampled_R_vecs(3); ...
        max_sampled_R_vecs(1), max_sampled_R_vecs(2), min_sampled_R_vecs(3); ...
        max_sampled_R_vecs(1), max_sampled_R_vecs(2), max_sampled_R_vecs(3); ...
        ];
    faces = [ ...
        1 5; 2 6; 3 7; 4 8; ...
        1 3; 2 4; 5 7; 6 8; ...
        1 2; 3 4; 5 6; 7 8; ...
        ];

    scatter3(ax, bbox_R_vecs(:, 1), bbox_R_vecs(:, 2), bbox_R_vecs(:, 3), 30, 'y', 'filled');
    patch('parent', ax, 'Faces', faces, 'Vertices', bbox_R_vecs, 'FaceColor', 'w');
end