Newer
Older
Nicola Vigano
committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
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