Improving the efficiency of triple nested loop

The code tests a 3D binary array for the largest sphere that can fit into the porous regions. The code works however for a 40x40x40 array for a triple nested loop, the code runs slowly. I was wondering if there were a way that it could be rewritten as to not take so long. An initial idea was not using interp3 every single iteration but I was not sure how to write that.

function poresize=poresizeAlgorithm(tiff,scaffoldLength)

tiffInv=~tiff;

j=size(tiff,1);
k=size(tiff,2);
l=size(tiff,3);

[x,y,z] = meshgrid([1:j],[1:k],[1:l]);
[x1,y1,z1] = meshgrid([1:j],[1:k],[1:l]);

r=0.5;
position=[];
for m=1+ceil(r):size(tiffInv,1)-ceil(r)
    for n=1+ceil(r):size(tiffInv,2)-ceil(r)
        for o=1+ceil(r):size(tiffInv,3)-ceil(r)
            
            centreSphere=tiffInv(m,n,o);
            if centreSphere==1
                %add sphere to array
                equ = ((x-m).^2 + (y-n).^2 + (z-o).^2)./(r.^2);
                z11 = interp3(x,y,z,equ,x1,y1,z1,'nearest'); %was spline
                ix = 1 > z11;
                volSphere=sum(ix,'all');
                
                
            testSphere=ix+tiffInv;
                numTwo=sum(testSphere(:) == 2);
                if volSphere==numTwo
                    output=[m,n,o,r];
                    position=[position output];
                    r=r+0.25;
                    
                    position2=[reshape(position,4,[])]';
                    continue
                end
            end     
        end
    end  
end

poresize=(scaffoldLength/j) * (position2(end,4)*2);
end

[Dyuman Joshi]

x1, y1, z1 are same as x, y, z.

So, using interp3 does not make sense to me, as the output is going to be the same as the input "equ".

Have you tried Profiling your code to see what the bottleneck(s) is(are)?

[Stephen23]

[x,y,z] = meshgrid([1:j],[1:k],[1:l]);
%                  ^   ^ ^   ^ ^   ^      superfluous square brackets
[x1,y1,z1] = meshgrid([1:j],[1:k],[1:l]);
%                     ^   ^ ^   ^ ^   ^   superfluous square brackets
position2=[reshape(position,4,[])]';
%         ^                      ^        superfluous square brackets

Superfluous square brackets do nothing except slow down your code and make it harder to read. Get rid of them.

[Matthew Bedding]

Thank you both for the help!

[Dyuman Joshi]

What is the objective/idea behind the interp3 line?

If you can provide additional details regarding what you are trying to do, we might be able to offer more suggestions.

[Matthew Bedding]

The interp3 line is to generate a binary 3D array where 1 defines voxels that are inside the sphere and 0 defines voxels that are outside the sphere.

Answer 1

[Matt J]

It seems like you could just use bwdist,

poresize = 2*max(bwdist(tiff),[],'all')

[Matthew Bedding]

Since it will be used for a fluid dynamics equation, the algoirthm needs to determine the diameter of the largest sphere rather than just the largest difference between two voxels; but thank you for sharing the function - will be useful for something else!

[Matt J]

I don't see the difiference. The distance of a black voxel to the nearest white voxel is the radius of the largest sphere, centered at the black voxel, that will fit inside the black region. Taking the max of this over all the voxels should be what you say you want.

Note though that you have not posted an image illustrating the problem. It's always recommended to do so with image processing problems.

[Matthew Bedding]

Apologies, thanks for your help anyway.

[Matt J]

So, you are now persuaded that my answer does solve the problem you've posted?

[Matthew Bedding]

It doesn't but thanks for your help