How to vectorize this code

The goal here is to find peaks in the xy plane. This is a general example but my particular utilization uses a much bigger dataset and is slower than I would like. Would vectorizing it help? Another method? I want there to be a thresholding.

[x,y,z] = peaks;

prom = 1;

% Find dimensions to set up loop
xdim = size(x,1);
ydim = size(x,2);

% Loop through x dimension to find peaks of each row
xpeaks = zeros(size(z));

for i = 1:xdim
    [~,locs] = findpeaks(z(i,:), 'MinPeakProminence', prom);
    xpeaks(i,locs) = 1;
    
end

% Loop through y dimension to find peaks of each row
ypeaks = zeros(size(z));

for i = 1:ydim
    [~,locs] = findpeaks(z(:,i), 'MinPeakProminence', prom);
    ypeaks(locs,i) = 1;
    
end

% Find indices that were peaks in both x and y
peak_inds = xpeaks+ypeaks == 2;

% Plot
figure
peaks
hold on
plot3(x(peak_inds),y(peak_inds),z(peak_inds),'r*','MarkerSize',24)

Answer 1

[Mara]

[x,y,z] = peaks;

prom = 1;

% Find dimensions to set up loop
xdim = size(x,1);
ydim = size(x,2);

% vectorize (the vectorization always concatenates the columns but you can
% just transpose the matrix to "concatenate your rows").
zv1 = z(:);
z_transposed = z';
zv2 = z_transposed(:);

% find the peaks in each of the two vectors and mark indices with 1 in xpeaks, ypeaks
[xpeaks, ypeaks] = deal(zeros(size(zv1)));

[~,pks] = findpeaks(zv1, 'MinPeakProminence', prom);
xpeaks(pks) = 1;

[~,pks] = findpeaks(zv2, 'MinPeakProminence', prom);
ypeaks(pks) = 1;

% reshape the vectors back into the original matrix size
xpeaks = reshape(xpeaks, xdim, []);
ypeaks = reshape(ypeaks, ydim, [])'; % here it is transposed back (')

% Find indices that were peaks in both x and y
peak_inds = xpeaks+ypeaks == 2;

% Plot
figure
peaks
hold on
plot3(x(peak_inds),y(peak_inds),z(peak_inds),'r*','MarkerSize',24)

[Stephen Thompson]

Excellent, this runs much faster too.

Answer 2

[Utkarsh]

Hi Stephen Thompson,

From your question, it seems like you want to find peaks in a 2D matrix without using a for loop

For this you may look at imregionalmax function which finds peak and returns a logical matrix for the input.

For example,

img = randn(3,3) 
imregionalmax(img) 

[Stephen Thompson]

I did looks at that - however I need the thresholding, not merely all peaks.