# Optimizing distance calculation between vectors and pixels

2 views (last 30 days)
Dominik Rhiem on 22 Sep 2023
Edited: Bruno Luong on 25 Sep 2023
I have written a Ray Tracer which I am currently trying to optimise, as there are 1-2 functions that take up around 90% of runtime. Take note that I have already vectorised those functions, and I am computing them on a GPU. This is the slowest function:
function distances = dist_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.
vox_pos = gpuArray(vox_pos);
vectors = gpuArray(vectors);
ray_pos = gpuArray(ray_pos);
vox_pos = reshape(vox_pos, 2, 1, size(vox_pos,2));
the_norm = vecnorm(vectors);
vox_vec = vox_pos - ray_pos;
%project the vectors onto the vectors connecting the starting positions and
%the voxels
%these three lines are about 25% computation time (10, 10, 5)
projections = (sum(vectors .* vox_vec,1)) ./ the_norm.^2;
%include safeguard: we only want rays in the positive directions
projections(projections<0) = Inf;
closest_points = ray_pos + projections .* vectors;
%this is a 2 x n_vec x n_pix array
dummy = vox_pos - closest_points;
%this is around 50% of computation time
distances = squeeze(sqrt(dummy(1,:,:).^2 + dummy(2,:,:).^2));
%now we test whether a ray that hits an object hits the object or the pixel
%first
if exist('t_hit_min', 'var')
t_hit_min = gpuArray(t_hit_min);
projections = squeeze(projections);
%we need a failsafe here; if we only have 1 voxel, projections is 1*x,
%while the right side is x*1 in size
%this is a simple but effective method
dummy = repmat(t_hit_min, 1, size(vox_pos,3));
if size(vox_pos,3) == 1
dummy = dummy.';
end
idx_dummy = projections < dummy;
distances(~idx_dummy) = Inf;
end
Note that the line where I calculate the minimum distances between all pixels and all rays is taking the longest time (I used the profiler to check the code.)
I am "brute force calculating" the distances between all pixels and rays (at their closest points) because I can then simply check for valid rays (i.e. those that come close enough to the pixels) by
valid_rays_dir = dist_dir < pixel_size;
valid_rays_hit = dist_hit < pixel_size;
This is given to another function. Let me first explain a phenomenon why this is needed in the first place.
Imagine a pixel that is submerged in an object. Due to a combination of object geometry and refraction, that pixel may be hit by two or more separate rays (or rather sub-bundles of rays) that come from different sides of the object. My goal here is to select the "most important" of those sub-bundles (for an in-depth look at the current solution, see this thread: https://de.mathworks.com/matlabcentral/answers/2022537-identifying-regions-in-matrix-rows?s_tid=srchtitle ), which I currently do by selecting which bundle contains the highest amount of rays. In other words, the "valid_rays" form bundles which I need to identify.
Any suggestions on how to optimize the code/calculations shown here? For large calculations (~100.000 antennae) it still takes too long for my taste (~24h). If there is anything unclear here, please let me know and I will provide more information.
Bruno Luong on 25 Sep 2023
A simple idea of reduction pair is:
• compute the bounding box for each ray then
• compute only the distance of the pixels falling ising the bounding box +/- 1 pixel.
That will disrupt your 3D array calculation where all the pixels are in the third dimension.
Dominik Rhiem on 25 Sep 2023
I've thought of something similar. I'll look into it. Thanks for the help!

Bruno Luong on 22 Sep 2023
This is different function that return the squared of the distance.
I simplify the code and use instructions that I think it's faster.
On my PC:
• dist2_ray_pix_GPU: 0.0316 second
• dist_ray_pix_GPU: 0.0825 second
So it looks like I save about half of the time.
If you not need the sqrt on top it would save even more.
clear
vox_pos = randn(2,10000);
vectors = randn(2,1000);
ray_pos = randn(2,1);
timeit(@() gather(sqrt(dist2_ray_pix_GPU(vox_pos, vectors, ray_pos))), 1)
timeit(@() gather(dist_ray_pix_GPU(vox_pos, vectors, ray_pos)), 1)
function distances = dist_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.
vox_pos = gpuArray(vox_pos);
vectors = gpuArray(vectors);
ray_pos = gpuArray(ray_pos);
vox_pos = reshape(vox_pos, 2, 1, size(vox_pos,2));
the_norm = vecnorm(vectors);
vox_vec = vox_pos - ray_pos;
%project the vectors onto the vectors connecting the starting positions and
%the voxels
%these three lines are about 25% computation time (10, 10, 5)
projections = (sum(vectors .* vox_vec,1)) ./ the_norm.^2;
%include safeguard: we only want rays in the positive directions
projections(projections<0) = Inf;
closest_points = ray_pos + projections .* vectors;
%this is a 2 x n_vec x n_pix array
dummy = vox_pos - closest_points;
%this is around 50% of computation time
distances = squeeze(sqrt(dummy(1,:,:).^2 + dummy(2,:,:).^2));
%now we test whether a ray that hits an object hits the object or the pixel
%first
if exist('t_hit_min', 'var')
t_hit_min = gpuArray(t_hit_min);
projections = squeeze(projections);
%we need a failsafe here; if we only have 1 voxel, projections is 1*x,
%while the right side is x*1 in size
%this is a simple but effective method
dummy = repmat(t_hit_min, 1, size(vox_pos,3));
if size(vox_pos,3) == 1
dummy = dummy.';
end
idx_dummy = projections < dummy;
distances(~idx_dummy) = Inf;
end
end
function d2 = 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 on 25 Sep 2023
I like this a lot, code now takes only ~2/3 of runtime. Thanks!

### More Answers (1)

Joss Knight on 25 Sep 2023

I feel like I haven't fully understood what you're after here, but pdist2 is the function you're supposed to use to compute distances between points.

Bruno Luong on 25 Sep 2023
Edited: Bruno Luong on 25 Sep 2023
Dominik computes the distances between a set of points and a set of (half) lines (in 2D)