# Optimizing distance calculation between vectors and pixels

4 vues (au cours des 30 derniers jours)
Dominik Rhiem le 22 Sep 2023
Modifié(e) : Bruno Luong le 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.
##### 7 commentairesAfficher 5 commentaires plus anciensMasquer 5 commentaires plus anciens
Bruno Luong le 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 le 25 Sep 2023
I've thought of something similar. I'll look into it. Thanks for the help!

Connectez-vous pour commenter.

### Réponse acceptée

Bruno Luong le 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
##### 1 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Dominik Rhiem le 25 Sep 2023
I like this a lot, code now takes only ~2/3 of runtime. Thanks!

Connectez-vous pour commenter.

### Plus de réponses (1)

Joss Knight le 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.

##### 1 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Bruno Luong le 25 Sep 2023
Modifié(e) : Bruno Luong le 25 Sep 2023
Dominik computes the distances between a set of points and a set of (half) lines (in 2D)

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Linear Algebra dans Help Center et File Exchange

### Community Treasure Hunt

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

Start Hunting!