Optimizing nested loops to measure min distances in 3D

Thread view:

Enable EMail Alerts

Start New Thread

Thread rating:

Graham - 24 Jul 2008 18:44 GMT

I'm working on a block of code that measures the distance
between every point within a volume to a set of predefined
existant in the same volume.

Basically I'm interested in how far any given volume element
(i,j,k) is from a set of blood vessels (represented by sets
of Cartesian coordinates and a diameter). I've tried to do
this with nested loops but the run time is exceedingly long
as each run requires about 13E (300x300x50x3000) iterations.

I had read another thread that used a hankel matrix to
optimize a linear operation in a nested loop but I don't
understand how to apply it to my problem.
(http://www.mathworks.com/matlabcentral/newsreader/view_thread/165240)

Below is my current approach with a single dummy coordinate
set with only 3 points within the set.

xmax = 10;
ymax = 10;
zmax = 10;
vAVGF = [1 2 3 6;2 3 4 6;3 4 5 6];
mintissuedistance = [];
maxtissuedistance = [];
tissuedistance = [];
ijkdistance = [];
nvessels = 1;
for i = 1:xmax
for j = 1:ymax
for k = 1:zmax
for m = 1:nvessels
nvsegs = size(vAVGF);
nvsegs = nvsegs(1,1);
for n = 1:nvsegs
for n = 1:nvsegs
ijkd = sqrt((i-vAVGF(n,1))^2 +...
(j-vAVGF(n,1))^2 + ...
(k-vAVGF(n,3))^2);
if ijkd < round(vAVGF(n,4)/2)
else
ijkdistance = vertcat(ijkdistance,ijkd);
end
end
end
mintissuedistance =
vertcat(mintissuedistance,min(ijkdistance));
maxtissuedistance =
vertcat(maxtissuedistance,max(ijkdistance));
tissuedistance =
vertcat(tissuedistance,ijkdistance);

ijkdistance = [];
end
end
end

averagemintissuedistance = mean(mintissuedistance)
averagemaxtissuedistance = mean(maxtissuedistance)
averagetissuedistance = mean(tissuedistance)

Variables:
xmax, ymax, zmax - limits of the bounding volume
vAVGF - dummy vessel coordinates and diameter (x, y, z, d)
mintissuedistance, maxtissuedistance, tissuedistance - are
the preallocated arrays I'm interested in finding
ijkdistance - distance from volume element to vessel segment
nvessels - total number of vessels to compare against (1)
nvsegs - number of vessel segments (coordinate pairs)

Thank you for any help!

Reply to this Message

Graham - 24 Jul 2008 18:54 GMT

Sorry! I screwed up in the posted code, please use the
following instead:

xmax = 10;
ymax = 10;
zmax = 10;
vAVGF = [1 2 3 6;2 3 4 6;3 4 5 6];
mintissuedistance = [];
maxtissuedistance = [];
tissuedistance = [];
ijkdistance = [];
nvessels = 1;
for i = 1:xmax
for j = 1:ymax
for k = 1:zmax
for m = 1:nvessels
nvsegs = size(vAVGF);
nvsegs = nvsegs(1,1);
for n = 1:nvsegs
ijkd = sqrt((i-vAVGF(n,1))^2 +...
(j-vAVGF(n,1))^2 + ...
(k-vAVGF(n,3))^2);
if ijkd < round(vAVGF(n,4)/2)
else
ijkdistance = vertcat(ijkdistance,ijkd);
end
end
end
mintissuedistance =
vertcat(mintissuedistance,min(ijkdistance));
maxtissuedistance =
vertcat(maxtissuedistance,max(ijkdistance));
tissuedistance =
vertcat(tissuedistance,ijkdistance);

ijkdistance = [];
end
end
end

averagemintissuedistance = mean(mintissuedistance)
averagemaxtissuedistance = mean(maxtissuedistance)
averagetissuedistance = mean(tissuedistance)

Reply to this Message

Peter Boettcher - 24 Jul 2008 19:34 GMT

> Sorry! I screwed up in the posted code, please use the
> following instead:
[quoted text clipped - 39 lines]
> averagemaxtissuedistance = mean(maxtissuedistance)
> averagetissuedistance = mean(tissuedistance)

You could improve the code right away by preallocating your matrices and
using a counter to add to them.

But I'm not quite sure what kind of data structure you're trying to
create here. It looks like for each i/j/k index you want the distance
to the closest vAVGF entry, but only if it is within a certain
threshold. But then you want the average min/max/total distances for
the entire dataset...

Given that the i/j/k index set is evenly spaced and even integer, you
should really reverse the order of your search. Say I'm looking for the
closest integers to (3.2, 4.2, 1.7). Do I really need to loop from
(0,0,0) to (10,10,10) across all 3 dimensions, comparing each, to
realize that my closest is (3,4,2)? Take advantage of that even
spacing.

I realize actually want ALL the i/j/k indices that fall within the
threshold. Fine. You can at least limit your i/j/k search to a 3D
bounding box. In the above example, if the radius is 3, then you might
actually have to calculate the integer indices between:

0.2,1.2,0 to 6.2,7.2,4.7

Don't forget to round inward: (0,2,0) to (6,7,4) That's a MUCH smaller
set of indices to compute actual euclidean distance from.

The downside is you have to store a running minimum and a running
maximum for each i/j/k value. But that's easy. After you compute the
distances (over your smalling bounding box) into a 3D matrix, just take
the min of that whole thing with your running min matrix, and same with
the max.

Does that make sense? Shrinking the number of points you need to
compute is the most powerful way to reduce the runtime. If you can
express more specifically what you are trying to compute off of this,
you might be able to do better.

-Peter

Reply to this Message

Math Forum / Math Software / MATLAB / July 2008

Optimizing nested loops to measure min distances in 3D