Suppose we have an array like this :

A = peaks
A(A<0)=0   % Remove negative values to make things a bit more simple.

I would like to write a script that is able to select the whole volume of the peaks contained in A.

For example, include all the values that surround a local maximum, until it hits values which are either 0 or a local minimum.

Below is a link to an answered question on a similar but more complicated topic.

In these results, 6 peaks are identified, but what I want is 3 separated peaks, as described in the picture below:

Hope this is clear enough.

CodePudding user response：

I am posting another separate answer to this question because it is a completely different approach, and it is complementary to the first one depending on the specific problem you are trying to solve

This approach is based on another definition of "Being inside a peak".

We will say that a point (Xi,Yi) belongs to the peak #k if, when doing a gradient ascent, the position (Xi_,Yi_) (Iterating during the ascent), converges to the position of the maximum of the peak #k.

This process can be visualized by using

We will use the approach presented

Note that the matrix filtMtx allows you to index each peak, X(filtMtx == ii),Y(filtMtx == ii),im(filtMtx == ii) are the (X,Y,Z) coordinates of all the points related to peak ii.

Please also note that this way doesn't guarantee you that all points will be linked to a maxima (for example, local minima have nil gradient and will not be linked to anything), but it guarantees that most of your points will.

CodePudding user response：

There might be better (faster) approaches, but the following one seems to be working fine.

So the problem with what you want to do is that although it is possible to find the positions of the local maximas

If you want to extract the values in your original array corresponding to each region:

[sz1,sz2] = size(im);
Output_Matrices = zeros(sz1,sz2,numel(s));

for i=1:numel(s)
    
    tmp=im*0;
    tmp(s(i).PixelIdxList)=1;
    tmp2=tmp.*im;
    
    % The maximum amplitude and location
    
    [refV,b]=max(tmp2(:));
    [x2,y2]=ind2sub(size(im),b);
    
    % select the region around local max amplitude
    tmp=bwselect(im>refV*Best_tresholds(i),y2,x2,4);
    
    [xi,yi]=find(tmp);
    
    hold on, plot(yi,xi,'r.')
    hold on, text(y2 10,x2,num2str(i),'Color','white','FontSize',16)
    
    % Your output matrix i-th slice holds the values for the i_th peak
    Output_Matrices(xi,yi,i) = im(xi,yi); 


end