How to use "imregionalmax" without using for loop

Hello, I have an 51×20×20 Array and I wish to find the second biggest peak in each 20×20 array (I find peaks by performing "imregionalmax" over each 20×20 arrays). As the dimension of my array represents, I have 51 of these 20×20 arrays. I want to do this without using for loop to iterate over all these 51 arrays and also, increase my speed. Can anyone help me with this?

UPDATE

This is my original time-consuming part of code:

z = 0;
for i=1:size(AF_Nuni,1),
    peaks = sort(AF_Nuni(i,imregionalmax(squeeze(AF_Nuni(i,:,:)))));
    peaksi = peaks(end - 1);
    z = z + peaksi;
end

"AF_Nuni" is a 1156*34*34 gpuArray. it takes about 8.2 seconds which is not acceptable for me.

[Matt J]

It is doubtful that a for-loop would slow you down.

[Image Analyst]

I agree with @Matt J. A paltry 20 iterations is nothing. It will take virtually no time at all.

tic;
for k = 1 : 20
    ;
end
toc

Elapsed time is 0.003824 seconds.

It's the call to imregionalmax that takes the majority of time because it has to loop over all rows and columns in the image to find the max values. That will take much more time than the 20 iterations.

[moh mor]

Thank you @Image Analyst and @Matt J , the number of for iteration is more than 1000 in my real program and I have to implement this thousands of times. So, I think for loop is still an important factor as well. Can you mention efficient replacements for "imregionalmax" function to increase speed?

[Matt J]

Can you mention efficient replacements for "imregionalmax" function to increase speed?

What about my Answer below?

[DGM]

A 1000x1000x1000 uint8 array takes about 80 seconds to process on my computer with a 12 year old motherboard that literally came out of a dumpster. Any modern system would take even less time.

tic
sz = [1000 1000 1000];
inpict = uint8(randi([0 255],sz));
mk = false(sz);
toc

Elapsed time is 12.071484 seconds.

tic
for c = 1:sz(3)
    mk(:,:,c) = imregionalmax(inpict(:,:,c));
end
toc

Elapsed time is 35.938297 seconds.

In all likelihood, a large portion (if not the majority) of the time spent will be reading the possibly multi-gigabyte files from disk, array reshaping, and other preprocessing and display tasks, not just the actual imregionalmax() call alone. You still haven't posted any code or notes about actual image sizes, filetypes, or numeric classes, so nobody can do anything but guess how to make your code faster.

[moh mor]

thank you guys, i thought that the problem is easy to understand completely. so sorry that I have not posted my code yet.

This is my original time-consuming part of code:

z = 0;
for i=1:size(AF_Nuni,1),
    peaks = sort(AF_Nuni(i,imregionalmax(squeeze(AF_Nuni(i,:,:)))));
    peaksi = peaks(end - 1);
    z = z + peaksi;
end

"AF_Nuni" is a 1156*34*34 gpuArray. it takes about 8.2 seconds which is not acceptable for me.

[Matt J]

Why peaks(end-1)? You want the 2nd largest peak? WHat if there is only 1 peak?

[moh mor]

@Matt J

no, I am sure that there are a lot of peaks in each one of these 2D array.

[Matt J]

So you are confirming that you want the second largest peak? Does it matter then if the maximum peak occurs twice? In such a case, peaks(end-1) will be the same as the max peak.

[moh mor]

@Matt J

this is not going to happen. anyway, if there are two maximum points, we'll consider one of them as our second biggest peak. no probem!

[Catalytic]

Are the values of AF_Nuni non-negative?

[moh mor]

yes it is nonnegative and something between 0 and 1.

Answer 1

[Matt J]

Perhaps as follows,

AF_Nuni = permute(AF_Nuni,[2,3,1]); %Avoid this permutation by forming AF_Nuni 
                                    %in an appropriate order

AF_Nuni( ~imregionalmaxSices(AF_Nuni)  ) = nan;

peaks=maxk(   reshape(  AF_nuni,  [],size(AF_Nuni,3)  )  ,2,1);

z=sum(peaks(2,:));


function D=imregionalmaxSices(A)

    B=padarray(A,[1,1],-inf);  %A is the input array
    
    C=reshape( imregionalmax(B(:,:))  ,size(B));
    
    D=C(2:end-1,2:end-1,:);  %the result


end

[moh mor]

brilliant answer. running speed decreased from 8.2s to 6.8s. But, this is not enough. if there is any other solution, please let me know. thank you

[Matt J]

You could also try an alternative definition of imregionalmaxSlices,

function D=imregionalmaxSices(A)
    
   D=false(size(A));
   for i=1:size(A,3)
      D(:,:,i)=imregionalmax(A(:,:,i)); 
   end
   
end

Also, try to avoid the permute() operation, as I recommended in the comments.

[moh mor]

"imregionalmax" is time-consuming and I think it is not avoidable.

[Matt J]

"imregionalmax" is time-consuming and I think it is not avoidable.

Do you really need to demand that it be a regional max? Could it just be a local max (i.e., a pixel that is >= its neighbors)? You said earlier that no two maxima will typically have the same value. That being the case, why would you ever even have a regional max that is larger than 1 pixel?

If a local max is enough, my answer can be modified to,

AF_Nuni = permute(AF_Nuni,[2,3,1]); %Avoid this permutation by forming AF_Nuni 
                                    %in an appropriate order

AF_Nuni( AF_Nuni<imdilate(AF_nuni,ones(3)) ) = nan;

peaks=maxk(   reshape(  AF_nuni,  [],size(AF_Nuni,3)  )  ,2,1);

z = sum(peaks(2,:)) ;

[moh mor]

thank @Matt J ,

Actually when I wanted to implement above code, I came across the following error:

Error using morphopInputParser>checkInputImage (line 97)
Underlying class for gpuArray image A must be uint8 or logical.

Error in morphopInputParser (line 12)
A  = checkInputImage(A);

Error in morphop (line 13)
[A,se,padfull,unpad,op_type] = morphopInputParser(varargin{:});

Error in gpuArray/imdilate (line 63)
    B = morphop(A,se,'dilate',mfilename,varargin{:});

Error in CostFunctionf4 (line 27)
AF_Nuni( AF_Nuni<imdilate(AF_Nuni,ones(3)) ) = nan;

I thought that's because of "AF_Nuni" which is a gpuArray. I used "gather(AF_Nuni)" instead which solved the error but the speed did not change. 6.8s.

[Matt J]

Not me. I see a 2x speed-up.

AF_Nuni = rand(34,34,11560);
tic;
AF_Nuni( ~imregionalmaxSices(AF_Nuni)  ) = nan;
peaks=maxk(   reshape(  AF_Nuni,  [],size(AF_Nuni,3)  )  ,2,1);
z=sum(peaks(2,:));
toc

Elapsed time is 0.834718 seconds.

AF_Nuni = rand(34,34,11560);
tic;
AF_Nuni( AF_Nuni<imdilate(AF_Nuni,ones(3)) ) = nan;
peaks=maxk(   reshape(  AF_Nuni,  [],size(AF_Nuni,3)  )  ,2,1);
z = sum(peaks(2,:)) ;
toc

Elapsed time is 0.392757 seconds.

function D=imregionalmaxSices(A)

    B=padarray(A,[1,1],-inf);  %A is the input array
    
    C=reshape( imregionalmax(B(:,:))  ,size(B));
    
    D=C(2:end-1,2:end-1,:);  %the result


end

[Matt J]

Here is a GPU-compatible version. On the GTX 1080 Ti, I see 1.6x speed-up in double float precision, though about a 20% in single floats,

gd=gpuDevice();

Z=double(full(sprand(34,34*11560, 10/34/34)));
Z=reshape(gpuArray(Z),34,34,[]);

AF_Nuni = Z;
tic;
AF_Nuni( ~imregionalmaxSices(AF_Nuni)  ) = nan;
peaks=maxk(   reshape(  AF_Nuni,  [],size(AF_Nuni,3)  )  ,2,1);
z=sum(peaks(2,:));

wait(gd)
toc

Elapsed time is 0.085708 seconds.

AF_Nuni = Z;
tic;
AF_Nuni( AF_Nuni<movmax2(AF_Nuni) ) = nan;
peaks=maxk(   reshape(  AF_Nuni,  [],size(AF_Nuni,3)  )  ,2,1);
z = sum(peaks(2,:)) ;

wait(gd)
toc

Elapsed time is 0.053386 seconds.

function A=movmax2(A)

 A=movmax(A,3,1);
 A=movmax(A,3,2);

end

function D=imregionalmaxSices(A)

    B=padarray(A,[1,1],-inf);  %A is the input array
    
    C=reshape( imregionalmax(B(:,:))  ,size(B));
    
    D=C(2:end-1,2:end-1,:);  %the result

end

[moh mor]

unfortunatley, it gives the following error:

Error using movmax
Invalid data type. First input must be numeric or logical.

[Matt J]

First input must be numeric or logical.

And indeed it should be! Surely if you are giving movmax non-numeric input, it is inadvertent, and something you can easily debug.

Answer 2

[Image Analyst]

peaks = sort(AF_Nuni(i,imregionalmax(squeeze(AF_Nuni(i,:,:)))));

The way you're indexing could be slowing you down. Normally you loop over the last index, not the first one. If you could rotate your volumetric image in advance and then iterate over the last index. You could access memory much more efficiently and you would not need to use squeeze.

[moh mor]

thank you @Image Analyst,

Sorry , but I do not why when I reshape AF_Nuni to a 34*34*1156 and want to do whatever you said, I face with this error in the 1st iteration, i.e. i=1 (I have checked AF_Nuni demention and size already). can you help me with this?

peaks = sort(AF_Nuni(imregionalmax(AF_Nuni(:,:,i)),i));

here the error:

Error using gpuArray/subsref
Index exceeds matrix dimension.

[Image Analyst]

Instead of

for i=1:size(AF_Nuni,1),

you need

for i=1:size(AF_Nuni,3)