Optimizing distance calculation between vectors and pixels

2 次查看（过去 30 天）
Dominik Rhiem
Dominik Rhiem 2023-9-22
编辑： Bruno Luong 2023-9-25
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 个评论
显示 5更早的评论
Bruno Luong
Bruno Luong 2023-9-25
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bruno Luong
Bruno Luong 2023-9-22
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 个）

Joss Knight
Joss Knight 2023-9-25

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.

类别

在 Help Center 和 File Exchange 中查找有关 Image Segmentation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by