# Faster alternative to mode?

7 views (last 30 days)
Eric Chadwick on 6 Apr 2020
Commented: Matt J on 15 Apr 2020
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
darova on 7 Apr 2020
Example:
tic
for i = 1:1e6
a = ones(5);
a(3) = [];
end
toc
tic
for i = 1:1e6
a = ones(5);
a(3) = 0;
end
toc
Because ofre-sizing the arrays on each iteration
Eric Chadwick on 7 Apr 2020
Thanks, I will change those lines, but the real bottleneck is the mode function. Thats what I would like to simplify as it takes up 50% of the time in my function.

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
Eric Chadwick on 15 Apr 2020
A finished and had 205,097,786 non-zeros
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.

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);
Matt J on 7 Apr 2020
The 32-bit image has replaced the formerly binary values of the image with numbers ranging from -1 to greater than 255.
Is there any sparsity that you can take advantage of? Does the image contain mostly zeros?
Eric Chadwick 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.
For example, if a voxel of value 1 touches at least one -1, I am interested in it and will then determine whether it is touching mostly other voxels of another positive value or mostly -1s.
I have not worked with sparse matrices very much. What is your idea?