Efficiently calculating the pixels crossed by a ray bundle

5 次查看（过去 30 天）
Dominik Rhiem
Dominik Rhiem 2023-10-24
评论： Dominik Rhiem 2023-10-25
I have created a script that calculates which pixels on a regular square grid are crossed by which ray(s) of a bundle. Up to the line "B(idx_lin) = 1;", the code is quite efficient, I would say. It's calculated on the GPU for the simple reason that there are easily thousands of rays, and these calculations are sped up quite a lot on the GPU.
Recently, however, I have tried to go from square pixels to circles (with the circle diameter being the pixel edge length and the circles lying in the pixels' centers), and I have unfortunately not found a solution for vectorisation yet, instead solving it with a loop, which is highly inefficient, especially on a GPU. I have tested it and the results are correct, i.e. the solution works, its performance is just bad. The basis for the function min_dist2_ray_pix_GPU can be found here, for the interested reader. The function simply calculates the minimum distance between all rays and all pixels supplied to it.
Without a preselection of pixels, this operation is on the order of n_rays * n_pix_x * n_pix_y. The code shown here performs an operation on the order n_rays * (n_pix_x + n_pix_y), which is far quicker, creating a preselection of pixels which is then supplied to the minimum distance calculation, making that calculation far quicker in theory, if it was vectorised.
In principle, the issue is that each ray now needs a different selection of pixels it passes through, in most cases not even the same number of pixels, which I don't know how to supply to the function all at once. This is quite easily solved in a loop, as shown, but it makes the code slow. So how can this be vectorised? Or is there maybe a wholly different approach yielding the same results which I overlooked?
%all of this is usually supplied to the function, but I explicitly write it
%here so that calculations are possible
x_r = 0;
y_r = -0.3;
x_r = repmat(x_r, 1, 1001);
y_r = repmat(y_r, 1, 1001);
vectors_dummy = [0;1];
vectors = zeros(2,1001);
angles = linspace(-30,30,1001) * pi/180;
vectors = rot_matrix_2d(vectors_dummy,angles);
vectors = squeeze(vectors).';
vox.dx = 0.002;
vox.dy = 0.002;
vox.xstart = -0.2;
vox.xend = 0.2;
vox.ystart = -0.2;
vox.yend = 0.2;
vox.Nx = 201;
vox.Ny = 201;
coordaxis.x = linspace(vox.xstart, vox.xend, vox.Nx);
coordaxis.y = linspace(vox.ystart, vox.yend, vox.Ny);
coordaxis.z = 0;
[array.X, array.Y, array.Z] = ndgrid(coordaxis.x,coordaxis.y,coordaxis.z);
vox_pos = zeros(3, length(array.X(:)));
vox_pos(1,:) = array.X(:);
vox_pos(2,:) = array.Y(:);
N_rays = size(vectors, 1);
v = gpuArray(vectors.');
x_r = gpuArray(x_r);
y_r = gpuArray(y_r);
%define outer grid around pixels
x2 = gpuArray(vox.xstart - vox.dx/2:vox.dx:vox.xend + vox.dx/2);
y2 = gpuArray(vox.ystart - vox.dy/2:vox.dy:vox.yend + vox.dy/2);
% Calculate t-values for the intersection points with all the grid lines
tx = ( x2(:) - x_r ) ./ v(1,:);
ty = ( y2(:) - y_r ) ./ v(2,:);
% Sort the t-values for each ray
t = sort( [tx; ty], 1);
M = size(t,1);
if exist('t_hit_min', 'var')
t(t>t_hit_min.') = NaN;
end
% Calculate the distance between subsequent intersection points
dist_t = diff(t);
% Define mid point between two intersection points
p = reshape([x_r; y_r], 2, 1, N_rays) + ...
reshape(t(1:end-1,:) + 1/2*dist_t, 1, M - 1, N_rays) .* reshape(v, 2, 1, N_rays);
% Define nearest grid points to mid points between two intersection points
idx_x = (round((p(1,:,:) - vox.xstart)/ vox.dx) + 1);
idx_y = (round((p(2,:,:) - vox.ystart)/ vox.dy) + 1);
idx_x(idx_x < 1 | idx_x > vox.Nx) = NaN;
idx_y(idx_y < 1 | idx_y > vox.Ny) = NaN;
%% Calculate matrix
mask = isfinite(idx_x) & isfinite(idx_y);
idx_ray = repmat(reshape(1:N_rays, 1, 1, []), 1, size(idx_x, 2), 1);
idx_lin = sub2ind([vox.Nx, vox.Ny, N_rays], idx_x(mask), idx_y(mask), idx_ray(mask));
%this operation is necessary, as in some cases, the same linear index can
%occur multiple times (which I don't fully understand)
idx_lin = unique(idx_lin);
%boolean is more efficient and also sufficient
B = false(vox.Nx*vox.Ny, N_rays, 'gpuArray');
B(idx_lin) = 1;
for i = 1:N_rays
idx_dummy = (B(:,i) == 1);
d2 = min_dist2_ray_pix_GPU(vox_pos(1:2,B(:,i)), vectors(i,:).', [x_r(i); y_r(i)]);
B(idx_dummy,i) = (d2 < (vox.dx/2)^2);
end
function rot_vec = rot_matrix_2d(vec, theta)
if length(theta)>1
theta = reshape(theta, 1, 1, length(theta));
end
R = [cos(theta) sin(theta) ; -sin(theta) cos(theta)];
rot_vec = pagemtimes(R, vec);
end
  6 个评论
显示 4更早的评论
Dominik Rhiem
Dominik Rhiem 2023-10-24
@Bruno Luong
Fixed the error there.
As you said, the old code used circles. However, the old code was also a "brute force" method, as it calculates the distance between every pixel and ray, i.e. this is still pretty inefficient overall. This is what I wanted to improve upon, hence I wrote the code you see in this thread. And the code is far more efficient in fact, up to B(idx_lin) = 1! The issue that arised was that the square pixels, which are considered in the first part of the new code, result in worse signal to noise ratio in the end result (as this code is embedded in a larger code) - this is before I introduced the loop at the end of the code here, which goes back to considering the circles. Basically, I want to keep the efficiency of this new code and at the same time use circles again. I hope this explains the process.
Dominik Rhiem
Dominik Rhiem 2023-10-24
@Matt J Tigre is specifically for tomographic reconstruction, not forward projection, as far as I see - the code you see here is not really either of those. And from a quick view, it seems their algorithms are inapplicable to our specific use case, as their use case is CT. While we borrow from there, our specifics are quite different, as we work with radar; we already have our own algorithm stack, which is probably at this point easier to work with and expand than modifying an unfamiliar toolbox for a case it was not designed for nor intended to be used with.
I'm already aware of the merge sorted arrays code, I have simply not yet gotten around to take care of it.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bruno Luong
Bruno Luong 2023-10-24
编辑：Bruno Luong 2023-10-25
I only add an altervative to the for-loop at the end.
On my PC the runtime if about 10 time faster.
%all of this is usually supplied to the function, but I explicitly write it
%here so that calculations are possible
x_r = 0;
y_r = -0.3;
x_r = repmat(x_r, 1, 1001);
y_r = repmat(y_r, 1, 1001);
vectors_dummy = [0;1];
% vectors = zeros(2,1001);
angles = linspace(-30,30,1001) * pi/180;
vectors = rot_matrix_2d(vectors_dummy,angles);
vectors = squeeze(vectors).'; % 1001 x 2
vox.dx = 0.002;
vox.dy = 0.002;
vox.xstart = -0.2;
vox.xend = 0.2;
vox.ystart = -0.2;
vox.yend = 0.2;
vox.Nx = 201;
vox.Ny = 201;
coordaxis.x = linspace(vox.xstart, vox.xend, vox.Nx);
coordaxis.y = linspace(vox.ystart, vox.yend, vox.Ny);
coordaxis.z = 0;
[array.X, array.Y, array.Z] = ndgrid(coordaxis.x,coordaxis.y,coordaxis.z);
vox_pos = zeros(3, length(array.X(:)));
vox_pos(1,:) = array.X(:);
vox_pos(2,:) = array.Y(:); % 3 x nvoxels
N_rays = size(vectors, 1);
v = gpuArray(vectors.'); % 2 x 1001
x_r = gpuArray(x_r); % 1 x 1001
y_r = gpuArray(y_r); % 1 x 1001
tic
%define outer grid around pixels
x2 = gpuArray(vox.xstart - vox.dx/2:vox.dx:vox.xend + vox.dx/2); % 1 x 202
y2 = gpuArray(vox.ystart - vox.dy/2:vox.dy:vox.yend + vox.dy/2); % 1 x 202
% Calculate t-values for the intersection points with all the grid lines
tx = ( x2(:) - x_r ) ./ v(1,:); % 202 x 1001
ty = ( y2(:) - y_r ) ./ v(2,:); % 202 x 1001
% Sort the t-values for each ray
t = sort( [tx; ty], 1); % 404 x 1001
M = size(t,1); % == 404
if exist('t_hit_min', 'var')
t(t>t_hit_min.') = NaN;
end
% Calculate the distance between subsequent intersection points
dist_t = diff(t); % 403 x 1001
% Define mid point between two intersection points
p = reshape([x_r; y_r], 2, 1, N_rays) + ...
reshape(t(1:end-1,:) + 1/2*dist_t, 1, M - 1, N_rays) .* reshape(v, 2, 1, N_rays); % 3 x 403 x 1001
% Define nearest grid points to mid points between two intersection points
idx_x = (round((p(1,:,:) - vox.xstart)/ vox.dx) + 1); % 1 x 403 x 1001
idx_y = (round((p(2,:,:) - vox.ystart)/ vox.dy) + 1); % 1 x 403 x 1001
idx_x(idx_x < 1 | idx_x > vox.Nx) = NaN;
idx_y(idx_y < 1 | idx_y > vox.Ny) = NaN;
%% Calculate matrix
mask = isfinite(idx_x) & isfinite(idx_y);
idx_ray = repmat(reshape(1:N_rays, 1, 1, []), 1, size(idx_x, 2), 1); % 1 x 403 x 1001
idx_lin = sub2ind([vox.Nx, vox.Ny, N_rays], idx_x(mask), idx_y(mask), idx_ray(mask)); %long column vector
%this operation is necessary, as in some cases, the same linear index can
%occur multiple times (which I don't fully understand)
idx_lin = unique(idx_lin);
%boolean is more efficient and also sufficient
B = false(vox.Nx*vox.Ny, N_rays, 'gpuArray');
B(idx_lin) = 1;
toc % Elapsed time is 0.086350 seconds.
BORG = B; % Just for debug
tic
for i = 1:N_rays
idx_dummy = (B(:,i) == 1);
d2 = min_dist2_ray_pix_GPU(vox_pos(1:2,B(:,i)), vectors(i,:).', [x_r(i); y_r(i)]);
B(idx_dummy,i) = (d2 < (vox.dx/2)^2);
end
toc % Elapsed time is 0.907424 seconds.
B1 = B; % Save result, Just for debug
B = BORG; % restore original array, Just for debug
tic
[idx_v,idx_r] = find(B);
VOX_POS = vox_pos(1:2, idx_v);
R_VEC = vectors(idx_r,:).';
R_POS = [x_r(idx_r); y_r(idx_r)];
d2 = nc_min_dist2_ray_pix_GPU(VOX_POS, R_VEC, R_POS);
B(sub2ind(size(B),idx_v,idx_r)) = (d2 < (vox.dx/2)^2);
toc % Elapsed time is 0.091884 seconds.
B2 = B; % Save result, Just for debug
% Check correctness
isequal(B1, B2)
%%
function rot_vec = rot_matrix_2d(vec, theta)
if length(theta)>1
theta = reshape(theta, 1, 1, length(theta));
end
R = [cos(theta) sin(theta) ; -sin(theta) cos(theta)];
rot_vec = pagemtimes(R, vec);
end
%%
% (Non cross) distance
function d2 = nc_min_dist2_ray_pix_GPU(VOX_POS, R_VEC, R_POS, t_hit_min)
% this function calculates the distance between a bundle of
% vectors/rays, which start at some point(s) ray_pos, and a bunch of
% pixels/voxels given by their positions vox_pos
% vectors: 2 x n
% vox_pos: 2 x n
% ray_pos: 2 x n
% if vectors which hit an object's edge are considered, the parameter t is
% given in t_hit_min, where edge = ray_pos + t * vector.
vox_pos = gpuArray(VOX_POS);
vectors = gpuArray(R_VEC);
ray_pos = gpuArray(R_POS);
vox_vec = vox_pos - ray_pos;
the_norm2 = sum(vectors.*vectors,1);
projections = sum(vectors./the_norm2 .* vox_vec,1); % divide on 2D array rather than3D array
%include safeguard: we only want rays in the positive directions
projections(projections<0) = Inf;
% this is a 2 x n_vec x n_pix array
dummy = vox_vec - projections .* vectors;
d2 = sum(dummy.*dummy,1);
%now we test whether a ray that hits an object hits the object or the pixel
%first
if nargin >= 4
t_hit_min = gpuArray(t_hit_min);
idx_dummy = projections < t_hit_min;
d2(~idx_dummy) = Inf;
end
end
%%
function d2 = min_dist2_ray_pix_GPU(vox_pos, vectors, ray_pos, t_hit_min)
% this function calculates the distance between a bundle of
% vectors/rays, which start at some point(s) ray_pos, and a bunch of
% pixels/voxels given by their positions vox_pos
% vectors: 2 x n_vec
% vox_pos: 2 x n_pix
% ray_pos: 2 x n_vec
% if vectors which hit an object's edge are considered, the parameter t is
% given in t_hit_min, where edge = ray_pos + t * vector.
n_vec = size(vectors,2);
n_pix = size(vox_pos,2);
vox_pos = gpuArray(vox_pos);
vectors = gpuArray(vectors);
ray_pos = gpuArray(ray_pos);
vox_vec = reshape(vox_pos, 2, 1, n_pix) - ray_pos;
the_norm2 = sum(vectors.*vectors,1);
projections = sum(vectors./the_norm2 .* vox_vec,1); % divide on 2D array rather than3D array
%include safeguard: we only want rays in the positive directions
projections(projections<0) = Inf;
% this is a 2 x n_vec x n_pix array
dummy = vox_vec - projections .* vectors;
d2 = sum(dummy.*dummy,1);
d2 = reshape(d2, n_vec, n_pix);
%now we test whether a ray that hits an object hits the object or the pixel
%first
if nargin >= 4
t_hit_min = gpuArray(t_hit_min);
projections = reshape(projections, n_vec, n_pix);
idx_dummy = projections < t_hit_min;
d2(~idx_dummy) = Inf;
end
end
  2 个评论
Dominik Rhiem
Dominik Rhiem 2023-10-25
Thank you, that vectorisation is pretty much what I have been looking for! I have combined it with the approach described by @Matt J and managed to improve my code with that.
Dominik Rhiem
Dominik Rhiem 2023-10-25
@Bruno Luong
This is the final version:
N_rays = size(vectors, 1);
v = gpuArray(vectors.');
x_r = gpuArray(x_r);
y_r = gpuArray(y_r);
x2 = gpuArray(vox.xstart - vox.dx/2:vox.dx:vox.xend + vox.dx/2);
y2 = gpuArray(vox.ystart - vox.dy/2:vox.dy:vox.yend + vox.dy/2);
tx = ( x2(:) - x_r ) ./ v(1,:);
ty = ( y2(:) - y_r ) ./ v(2,:);
if exist('t_hit_min', 'var')
tx = [tx; t_hit_min.'];
end
% Sort the t-values for each ray
t = sort( [tx; ty], 1);
M = size(t,1);
if exist('t_hit_min', 'var')
t(t>t_hit_min.') = NaN;
end
t(t<0) = 0;
% Calculate the distance between subsequent intersection points
dist_t = diff(t);
% Define mid point between two intersection points
p = reshape([x_r; y_r], 2, 1, N_rays) + ...
reshape(t(1:end-1,:) + 1/2*dist_t, 1, M - 1, N_rays) .* reshape(v, 2, 1, N_rays);
% Define nearest grid points to mid points between two intersection points
idx_x = (round((p(1,:,:) - vox.xstart)/ vox.dx) + 1);
idx_y = (round((p(2,:,:) - vox.ystart)/ vox.dy) + 1);
idx_x(idx_x < 1 | idx_x > vox.Nx) = NaN;
idx_y(idx_y < 1 | idx_y > vox.Ny) = NaN;
%% Calculate matrix
mask = isfinite(idx_x) & isfinite(idx_y);
idx_ray = repmat(reshape(1:N_rays, 1, 1, []), 1, size(idx_x, 2), 1);
idx_lin = sub2ind([vox.Nx, vox.Ny, N_rays], idx_x(mask), idx_y(mask), idx_ray(mask));
idx_lin = unique(idx_lin);
B = false(vox.Nx*vox.Ny, N_rays, 'gpuArray');
%this doesnt work for sparse gpu arrays
B(idx_lin) = 1;
%using the pre-selection:
[idx_v,idx_r] = find(B);
VOX_POS = vox_pos(1:2, idx_v);
R_VEC = v(:,idx_r);
R_POS = [x_r(idx_r); y_r(idx_r)];
A = -R_VEC(2,:);
BB = R_VEC(1,:);
C = BB .* R_POS(2,:) + A .* R_POS(1,:);
B(sub2ind(size(B),idx_v,idx_r)) = abs(A .* VOX_POS(1,:) + BB .* VOX_POS(2,:) - C) <= vox.dx/2;
HOWEVER... sadly I noticed that the issue I had was not exactly due to the square pixels...
The issue was the parameter t_hit_min, which indicates at which value for the parameter t a vector hits a surface, where the vector is parametrized by:
x = x0 + t * vx;
y = y0 + t * vy;
The issue is that some rays may start within the grid, and they are not supposed to go backwards. Therefore t(t<0) = 0. This was the issue that caused my images to look worse. Unfortunately (and weirdly tbh), this fix makes the code a lot slower. I will have to keep looking into this.

请先登录，再进行评论。

更多回答（1 个）

Matt J
Matt J 2023-10-24
编辑：Matt J 2023-10-24
I have tried to go from square pixels to circles,
I'm assuming the circles are inscribed inside the square pixels? If so, I would first use the code you already have to detect intersections with the square pixels. In order for a line to intersect a circle, it must first intersect the square that circumscribes it, which narrows down the list of candidate intersections.
Once you've narrowed the list (or even if you haven't), it is easy to test intersections of lines with a circles in a vectorized manner. Given a line whose equation is in the form A*x+B*y=C where norm([A,B])=1 and C>=0 and given a circle of radius R and centered at (P,Q), the line intersects the circle if,
abs(A*P+B*Q-C)<=R
which can easily be vectorized over multiple lines and circles.
  5 个评论
显示 3更早的评论
Dominik Rhiem
Dominik Rhiem 2023-10-25
@Matt J finding the intersection points is not the goal, but the idea you provided already helped plenty. I have improved my code by combining your approach with the vectorisation provided by @Bruno Luong in his answer. I would like to accept both answers, but sadly, I can't.
Bruno Luong
Bruno Luong 2023-10-25
编辑：Bruno Luong 2023-10-25
@Dominik Rhiem Can you please post the final version of all the modifs integrated? I lost track of what ideas you use, and curious to see them.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Matrices and Arrays 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by