# Optimizing distance calculation between vectors and pixels

2 views (last 30 days)

Show older comments

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 Comments

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.

### Accepted Answer

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

### More Answers (1)

Joss Knight
on 25 Sep 2023

`pdist2` is the function you're supposed to use to compute distances between points.

##### 1 Comment

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)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!