Faster alternative to mode?

5 次查看（过去 30 天）
Eric Chadwick
Eric Chadwick 2020-4-6
评论： Matt J 2020-4-15
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 个评论
显示 3更早的评论
darova
darova 2020-4-7
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 2020-4-7
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
Matt J 2020-4-7
编辑：Matt J 2020-4-7
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.
https://www.mathworks.com/matlabcentral/fileexchange/29832-n-dimensional-sparse-arrays
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 个评论
显示 19更早的评论
Eric Chadwick
Eric Chadwick 2020-4-15
A finished and had 205,097,786 non-zeros
Matt J
Matt J 2020-4-15
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.

请先登录，再进行评论。

更多回答（1 个）

Matt J
Matt J 2020-4-6
编辑：Matt J 2020-4-6
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 个评论
显示 4更早的评论
Matt J
Matt J 2020-4-7
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 2020-4-7
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?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by