Trying to avoid a simple triple for loop.

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

Answer 1

[Daniel Shub]

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;

[fred sampedro]

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!