function [indx_full,data_peak]=peak_local(data,thresh,dim_lims,scales)
%
% Call [indx_full,data_peak]=peak_local(data,thresh,dim_lims);
%
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
if nargin<3, dim_lims=[]; end
%
dims=size(data); iidims=find(dims>1);
if length(dims)>2, 
    data=squeeze(data); dims2=size(data);
    if length(dims2)>2, error('Cannot handle datasets having more than two dimensions.'); return; end
end
[nrow,ncol]=size(data); npts=nrow*ncol;
iinan=find(isnan(data)); data(iinan)=-1e50; end     % Make NaN values very small to get adjacent peaks;
interior=ones([nrow,ncol]); interior([1:nrow,npts-nrow+1:npts,1:nrow:npts,nrow:nrow,npts])=0;
if ~isempty(dim_lims),
    yax=scales{iidims(1)}; xax=scales{iidims(2)};
    if ~isempty(dim_lims{1}),
        iiy=find(yax>=dim_lims{1}(1) & yax>=dim_lims{1}(2)); interior(iiy,:)=0; 
    end
    if ~isempty(dim_lims{2}),
        iix=find(xax>=dim_lims{2}(1) & xax>=dim_lims{2}(2)); interior(:,iix)=0;
    end
end
if length(thresh)==1,
% Find non-perimeter points that exceed threshold:    
    iimax1=find(data>=thresh(1) & interior);
elseif length(thresh)==2,
    iimax1=find(data>=thresh(1) & data<=thresh(2) & interior);
else
    error('Threshold must be a single number or [cmin,cmax]');
end
% Compute indices of eight surrounding points:
iibot=iimax1-1; iitop=iimax1+1; iileft=iimax1-nrow; iiright=iimax1+nrow;
iibotl=iibot-nrow; iibotr=iibot+nrow; iitopl=iitop-nrow; iitopr=iitop+nrow;
data2=data(iimax1);
% Find threshold-exceeding points that are local maxima:
iimax2=find(data2>data(iibot) & data2>data(iitop) & data2>data(iileft) & data2>data(iiright) & ...
    data2>data(iibotl) & data2>data(iibotr) & data2>data(iitopl) & data2>data(iitopr));
%
iimax=iimax1(iimax2); 
data_peak=data2(iimax2);
indx_full=indx_full_calc([nrow,ncol],iimax);
