6 views (last 30 days)

Simple as that, in my code I am first determining the neighbourhood of a voxel and then determining the most frequent value in the neighbourhood. I am doing this for millions of voxels, so the computation time adds up. I am expecting this to take quite some time, however my entire code has been running for about 2.5 days now with no end in sight and I've determined that the bottleneck is this section of my code. Using profiler it appears that the function "mode" takes up almost 50% of the computation time of the this function, and about 25% of my entire code.

My method of determining the neighbourhood is the next slowest part of the code using: neighbourhood = im(i-1:i+1,j-1:j+1,k-1:k+1);

Does anyone have a good alternative to "mode" i can use in this situation? If anyone has any more efficient suggestions for finding the neighbourhood of the voxels that would also be appreciated.

Thanks!

for i = 1:size(lcoords,1)

%Get neighbourhood of voxel i from list of image voxel subscripts

neighbourhood = im(lcoords(i,1)-1:lcoords(i,1)+1,lcoords(i,2)-1:lcoords(i,2)+1, lcoords(i,3)-1:lcoords(i,3)+1);

neighbourhood(neighbourhood == imcoords(i,2)) = []; %Ignore parts of the neighbourhood that have same value as current voxel

if ~isempty(neighbourhood) %If the voxel is surrounded by any values not including its own region, find the mode of the neighbourhood

imcoords(i,3) = mode(neighbourhood(:)); %record the mode

end

end

Matt J
on 7 Apr 2020

Edited: Matt J
on 7 Apr 2020

Yes, actually there are no zeros in this matrix. -1s could be changed to zeros as they represent void space that I am uninsterested in. However the point of the code is to identify voxels that touch void space (at least one -1) and determine which value that voxel is touching the most.

Starting another answer for this line of solution. I have a class on the File Exchange that can store 3D arrays in sparse form.

Not only could this save you a lot of memory, but it sounds like at least some of what you're trying to do can be done via sparse 3D convolution with a 3x3x3 kernel of ones. This can be done using the convn method of the class. For example, below I create a 3D binary image A the same dimensions as yours. It contains about 1.7 million non-zeros and consumes 26 MB. Only one of the 1's is not touching any zeros. Using convolution, the code calculates a separate binary matrix B with all elements not touching any zeros removed. Since there is only one such element, B contains only 1 fewer element than A:

>> A=spfun(@ceil, ndSparse.sprand([2854,2906,2013],1e-4) );

>> A(1:3,1:3,1:3)=1;

>>

>> whos A

Name Size Kilobytes Class Attributes

A 2854x2906x2013 26102 ndSparse

>>

>> tic; B=A-(convn(A,ones(3,3,3),'same')>26.999); toc;

Elapsed time is 12.614820 seconds.

>>

>> nnz(A)

ans =

1669474

>> nnz(B)

ans =

1669473

Matt J
on 15 Apr 2020

Ah well. Sadly, it's not nearly as sparse as I'd hoped and not enough for you to get an advantage out of ndSparse representation. The mode calculation is going to allocate at least 130 GB with that many non-zeros.

I would say your best bet is to do the processing in chunks like I recommended in my first answer. What I envision would be to process im(:,:,1:100), then im(:,:,100:199), then im(199:298) and so on. The overlap is deliberate, so that every 3x3x3 neighborhood is captured in the partitioning.

Sign in to comment.

Matt J
on 6 Apr 2020

Edited: Matt J
on 6 Apr 2020

It would be better not to implement the mode calculation repeatedly in a loop over the voxels. If you have enough RAM to hold 27 copies of your image volume, then this is one way to replace the loop with a vectorized calculation:

[di,dj,dk]=ndgrid(-1:+1);

Ishifts=repmat(im,1,1,1,27);

for n=1:27

Ishifts(:,:,:,n)=circshift(im,[di(n),dj(n),dk(n)]);

end

result=mode( Ishifts , 4);

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 5 Comments

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/515739-faster-alternative-to-mode#comment_822277

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/515739-faster-alternative-to-mode#comment_822277

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/515739-faster-alternative-to-mode#comment_822911

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/515739-faster-alternative-to-mode#comment_822911

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/515739-faster-alternative-to-mode#comment_822940

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/515739-faster-alternative-to-mode#comment_822940

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/515739-faster-alternative-to-mode#comment_822945

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/515739-faster-alternative-to-mode#comment_822945

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/515739-faster-alternative-to-mode#comment_823057

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/515739-faster-alternative-to-mode#comment_823057

Sign in to comment.