% function  pks = findpeaks2d(d, thres, filt)
% %FINDPEAKS2D Find local peaks in 2D data
% %   PKS = FINDPEAKS2D(I) finds local peaks in intensity matrix I. A local
% %   peak is defined as a data sample which is either larger than the
% %   neighboring samples or is equal to Inf.  PKS is as an Nx2 matrix [R C],
% %   with row vector R and column vector C indicating peak locations.
% %
% %   FINDPEAKS2D(Y,TH) specifies a threshold for removing noise in the data.
% %
% %   FINDPEAKS2D(Y,TH,F) specifies a 2D filter kernel F to smooth data prior
% %   to finding the peaks.
% %
% %   Applies an optional 2-D median filter if Image Processing Toolbox
% %   is available.
% 
% % Algorithm derived from FastPeakFind on MATLAB Central by Adi Natan
% 
% validateattributes(d,{'numeric'}, ...
%     {'2d','real'},'d','findpeaks2d');
% 
% if isfloat(d)
%     max_d = max(d(:));
%     min_d = min(d(:));
%     d = uint16((d-min_d)./(max_d-min_d) .* 65535);
% end
% if nargin < 2
%     thres = max([min(max(d,[],1)), min(max(d,[],2))]);
% end
% if nargin < 3
%     filt = fgaussian(7,1);
% end
% % Use 2D median filter if available (Image Processing Toolbox)
% if license('test','image_toolbox')
%     d = medfilt2(d,[3,3]);
% end
% if isa(d,'uint8')
%     d = d.*uint8(d > thres);
% else
%     d = d.*uint16(d > thres);
% end
% d = conv2(single(d),filt,'same') ;
% % Re-apply threshold:
% d = d .* (d > 0.9*thres);
% 
% [x,y] = find(d);
% pks = [];
% for j = 1:numel(y)
%     dxy = d(x(j),y(j));
%     if ...
%             (dxy >= d(x(j)-1,y(j)-1)) &&...
%             (dxy >  d(x(j)-1,y(j)  )) &&...
%             (dxy >= d(x(j)-1,y(j)+1)) &&...
%             (dxy >  d(x(j)  ,y(j)-1)) && ...
%             (dxy >  d(x(j)  ,y(j)+1)) && ...
%             (dxy >= d(x(j)+1,y(j)-1)) && ...
%             (dxy >  d(x(j)+1,y(j)  )) && ...
%             (dxy >= d(x(j)+1,y(j)+1))
%         pks = [pks; x(j) y(j)]; %#ok<AGROW>
%     end
% end
% 
% function h = fgaussian(sz,std)
% % 2D Gaussian filter kernel, square
% % std is gaussian variance, e.g., 0.5
% % sz must be scalar, odd, full-kernel width
% %
% % hsz is kernel half-width:
% hsz = (sz-1)/2;
% [x,y] = meshgrid(-hsz:hsz,-hsz:hsz);
% h = exp(-(x.*x+y.*y)/(2*std*std));
% h(h<eps*max(h(:))) = 0;
% sumh = sum(h(:));
% if sumh ~= 0
%     h  = h/sumh;
% end
