function [x, y, v] = peak(inp,threshold,blocksize,max_peaks)
s = size(inp);
if nargin<4
max_peaks = 20;
end
if nargin<3
blocksize = 5;
end
if nargin<2
threshold = 0.05;
end
np = 0;
fprintf('image size: %d x %d\n',s(1),s(2));
fprintf('block size: %d\n', blocksize);
fprintf('threshold: %g\n',threshold);
fprintf('max # peaks: %d\n',max_peaks);
while np < max_peaks
[col_val, row_idx] = max(inp);
[val, col] = max(col_val);
row = row_idx(col);
if val<threshold
break;
end
ra = max(row-blocksize,1);
rb = min(row+blocksize,s(1));
ca = max(col-blocksize,1);
cb = min(col+blocksize,s(2));
inp(ra:rb,ca:cb) = 0;
np = np + 1;
x(np) = row;
y(np) = col;
v(np) = val;
if (np==1)
fprintf('%4s\t%8s\t%4s\t%4s\n','peak','value','row','col');
end
fprintf('%4d\t%8.6f\t%4d\t%4d\n',np,val,row,col);
end
if (np==0)
disp('no peaks found exceeding threshold')
end