Optimizing distance calculation between vectors and pixels

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]

What are typical n_vex and n_pix ?

The better way is perhaps not compute all pair distances, is there any structure of voxels and ray that we use to eliminate pairs that does not need to compute (distance larger than pixel_size) ?

[Bruno Luong]

instead of working with distance may be you can work with distance squared, which save sqrt call. Also here are some faster ways to compute (squared) distance (on CPU)

dummy = rand(2,10000,10000);

tic; d=sqrt(dummy(1,:,:).^2+dummy(2,:,:).^2); toc

Elapsed time is 1.028496 seconds.

tic; d2=(dummy(1,:,:).^2+dummy(2,:,:).^2); toc

Elapsed time is 0.828537 seconds.

tic; d2=sum(dummy.^2,1); toc

Elapsed time is 0.635611 seconds.

tic; d2=sum(dummy.*dummy,1); toc

Elapsed time is 0.638204 seconds.

[Bruno Luong]

Also I wouldn' use vecnorm either

vectors = rand(2,10000,10000);;
tic; the_norm = vecnorm(vectors); the_norm_2 = the_norm.^2; toc

Elapsed time is 1.582805 seconds.

tic; the_norm_2 = sum(vectors.*vectors,1); toc

Elapsed time is 0.640674 seconds.

[Dominik Rhiem]

Thanks for your many replies, I will investigate them today!

n_vec starts at around 2000, and n_pix starts at around 20.000. However, both of these values fluctuate. The latter will decrease with later iterations; as soon as I have found a ray bundle that crosses the pixel, I throw out the pixel. The former is a bit more erratic; I throw out rays that do not hit any object, but rays that hit a surface will be split into a reflected and a refracted ray, i.e. they will double.

[Dominik Rhiem]

@Bruno Luong

Also, indeed we throw out those rays that do not approach a pixel close enough. I don't see an easy way to ignore them for distance calculation, though. The one rough idea I have in mind is that the distance to a pixel is a function that is piecewise smooth. Imagine one ray bundle: the outmost rays will not hit the pixel, but the innermost will. There will be two "border rays" that barely hit the pixel; every ray between those will also hit the pixel, while all rays outside of them will not hit the pixel. We could probably make use of this criterion somehow. However, in later iterations, we will also have several ray bundles (recall my other thread about "grouping" the rays).

[Bruno Luong]

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]

@Bruno Luong

I've thought of something similar. I'll look into it. Thanks for the help!

Answer 1

[Bruno Luong]

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]

I like this a lot, code now takes only ~2/3 of runtime. Thanks!

Answer 2

[Joss Knight]

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]

Dominik computes the distances between a set of points and a set of (half) lines (in 2D)