Trying to avoid a simple triple for loop.

1 次查看(过去 30 天)
Dearest Matlab users,
I'm trying to simplify the following simple code, because as you may observe is running really slow when called multiple (~1M) times. What I would like to do is simply check the surroundings (cubes of increasing size 3x3x3,5x5x5,...) of a voxel on a 3D volume, counting the number of voxels in the surrounding that are 1,0 and -1.
Thanks so much for any advice you may suggest to speed up this code :)
cubedims = [3 5 9 15];
conts = zeros(3,4);
%Given a point [xi yi zi]
for ic = 1:ncubs
lad = cubedims(ic);
for xi = 1:lad
for yi=1:lad
for zi=1:lad
if sum([xi yi zi] > 0) == 3 && sum([xi yi zi] <= lasizev) == 3
switch testmask(xi,yi,zi) %this is a 3D volume
case 1
conts(1,ic) = conts(1,ic)+1;
case 0
conts(2,ic) = conts(2,ic)+1;
case -1
conts(3,ic) = conts(3,ic)+1;
end
end
end
end
end
conts(:,ic) = conts(:,ic) ./ (lad*lad*lad);
end

回答(1 个)

Daniel Shub
Daniel Shub 2013-4-3
编辑:Daniel Shub 2013-4-3
It looks like you can eliminate the first part of the if statement since xi, yi and zi are always greater than 0. Instead of looping all the way to lad, you can stop at lasizev. You can then eliminate the if statement. Finally, you can also replace the switch statement with conts(2-testmask(...), ic).
This then leads to an vectorised solution which may or may not be faster:
conts(1, ic) = sum(testmask(1:lasizev, 1:lasizev, 1:lasizev))==1;
  1 个评论
fred sampedro
fred sampedro 2013-4-3
Thanks so much!
Sorry I forgot to mention that lasizev is the size of the volume (size(testmask) = 144x144x213), so that conditional statement simply checks that the point xi,yi,zi is within the volume limits to avoid an exception. Maybe I could eliminate this if statement if before the loop I make sure that all tested points will be within limits.
Then I guess what you suggest would not work? However that kind of matlab-one line is what I'm looking for.
Isn't there a simple way to get the 3D coordinates of a cube of size N? (i.e 3x3x3 = [1 1 1; 1 1 2; ... ;2 1 3...;3 3 3])
Thanks again!

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Graphics Object Programming 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by