MATLAB Answers

Faster alternative to mode?

6 views (last 30 days)
Eric Chadwick
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

  5 Comments

Show 2 older comments
darova
darova on 7 Apr 2020
The only idea i have is to change these lines
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
On
tmp = neighbourhood == imcoords(i,2)
if any(~tmp)
I think it would be faster
What about mode? Is there something to simplify there?
darova
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
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.

Sign in to comment.

Accepted Answer

Matt J
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

  21 Comments

Eric Chadwick
Eric Chadwick on 15 Apr 2020
I'm not sure what to tell you but when I create A its 30 GB, when i reduce the pre-allocation as you said it becomes 15 GB.
On the computer I was running the code on, the allocation of A has not completed because another process is actually taking longer to finish which is stranged because in the smaller image that process does not usually take a signifcant amount of time. That something I will have to look into.
Back to this problem though, I was actually getting an out of memory error when trying those lines but I was accidentally using the non-sparse version of the image. I will let you know when and if it finishes.
Eric Chadwick
Eric Chadwick on 15 Apr 2020
A finished and had 205,097,786 non-zeros
Matt J
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.

More Answers (1)

Matt J
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);

  6 Comments

Show 3 older comments
Eric Chadwick
Eric Chadwick on 7 Apr 2020
That's inconvenient but it is still no reason to go voxel-by-voxel. You can break the image into (overlapping) sub-chunks and apply my proposal to that.
Okay, I haven't had time to try to implement your code in that way. I will try that.
What are the dimensions of this image?
2854x2906x2013
That's not possible. If the original data is 8 bit, there is no way the neighborhood modes require more than 8 bits to represent.
The 32-bit image has replaced the formerly binary values of the image with numbers ranging from -1 to greater than 255.
Matt - I'd bet a decent sum of money the problem is mainly pre-allocation, or the lack thereof.
John - even if a missing pre-allocation step is remedied, looping element by element through a 16GB image will still make things impractically slow.
There is preallocation for every variable shown, I just only included the parts that I felt were relevant to the problem.
Matt J
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
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?

Sign in to comment.

Sign in to answer this question.