Efficiently calculating the pixels crossed by a ray bundle

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

[Bruno Luong]

Your code won't run:

Unrecognized field name "z".

Error in xxx (line 22)

[array.X, array.Y, array.Z] = ndgrid(coordaxis.x,coordaxis.y,coordaxis.z);

[Matt J]

This is clearly an implementation of Siddon's algorithm. I hope this is not being done for the purposes of tomographic forward projection. There are already freely available GPU-accelerated libraries for that, e.g.,

https://github.com/CERN/TIGRE/blob/master/Frontispiece/MATLAB_installation.md

although not for the case of circular pixels.

[Bruno Luong]

I did not read your new code to see what changes from the old thread, but there is no assumption in the ols thread that the pixels ("voxel"?) are square. They can have random positions and the distance to rays are computed using euclidian distance. So when you threshold it it detect all the rays that cross the circle disk that is centered about the position.

Not sure what changes here.

[Matt J]

Since tx and ty are already sorted, using sort() in this line,

t = sort( [tx; ty], 1);

will be inefficient. You might use @Bruno Luong's FEX contribution Merge Sorted Arrays.

[Dominik Rhiem]

@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]

@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.

Answer 1

[Bruno Luong]

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

[Dominik Rhiem]

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]

@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.

Answer 2

[Matt J]

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.

[Dominik Rhiem]

Yes indeed, sorry for missing that piece of information. The radius of the circles is half of the edge length of a pixel. (And I just noticed that I neglected a factor of 1/2 in my code...)

This sounds promising, I will look into it.

[Dominik Rhiem]

I've tested it:

A = -v(2,:);
BB = v(1,:);
C = BB .* y_r + A .* x_r;
B2 = abs(A.' * vox_pos(1,:) + BB.' * vox_pos(2,:) - C.') <= vox.dx/2;

This yields the same results as the original code which calculated the minimum distance between the pixels and the rays, while being far quicker, so this is already an improvement, thank you. However, it is similarly brute force-ish in that it looks at every pixel-ray-combination, so I would still like to "narrow the list" with the code presented earlier in the hopes of speeding up calculations even more. This causes the same problem again, though; it is easy to solve with a loop:

for i = 1:N_rays
	idx_dummy = (B(:,i) == 1);
	B(idx_dummy,i) = abs(A(i) * vox_pos(1,idx_dummy) + BB(i) * vox_pos(2,idx_dummy) - C(i)) <= vox.dx/2;
end

But I don't see how to vectorise this.

[Matt J]

I don't understand everything in your code, but you should be able to find, using t, the actual (x,y) locations where the lines intersect the squares.

x = x_r + t(:).*vectors(1,:);
y = y_r + t(:).*vectors(2,:);

Assuming you have normalized your coordinates so that the grid lines area at integer coordinates,

i=floor(x);
j=floor(y);

should give the discrete (i,j) voxel coordinates of the line-square intersections.

[Dominik Rhiem]

@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]

@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.