function lm = locmax(input_file,input_thresh,mask_file,mask_thresh, ...
    fwhm,flip);

%LOCMAX
%
% LM = LOCMAX( INPUT_FILE, INPUT_THRESH [, MASK_FILE, MASK_THRESH ...
%              [,FWHM [, FLIP]]] )
%
% Finds local maxima and clusters of INPUT_FILE above INPUT_THRESH
% and MASK_FILE above MASK_THRESH (ignored if empty; empty by default).
% If MASK_THRESH is a vector [a b], a<=b, then mask is a<MASK_FILE<=b. 
% Clusters are voxels connected in any of the 2*D directions. A local maximum 
% is a voxel which is greater than or equal to all its 2*D neighbours,
% and strictly greater than at least one of them.
%
% FWHM: If a file is provided, cluster volume is measured in resels.
% If FWHM is a scalar, this is taken as the FWHM for the whole volume.
% Default is 1, i.e. cluster volume is measured in mm^D.
%
% FLIP: INPUT_FILE is multiplied by FLIP before processing. For T statistic
% images, FLIP = -1 will look for negative peaks nd clusters. Default is 1.
% 
% LM is a matrix with 7 columns. 
% Col 1:    values of local maxima, sorted in descending order.
% Cols 2-4: x,y,z coords of local maxima in C notation i.e. starting at 0
%           so that they agree with 'register'.
% Col 5:    index of cluster to which they belong, in descending order of 
%           cluster resels (1=largest cluster).
% Col 6:    volume of cluster in mm^D.
% Col 7:    resels of cluster.

%############################################################################
% COPYRIGHT:   Copyright 2002 K.J. Worsley, 
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center, 
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada. 
%              worsley@math.mcgill.ca
%
%              Permission to use, copy, modify, and distribute this
%              software and its documentation for any purpose and without
%              fee is hereby granted, provided that this copyright
%              notice appears in all copies. The author and McGill University
%              make no representations about the suitability of this
%              software for any purpose.  It is provided "as is" without
%              express or implied warranty.
%############################################################################

% Defaults:

if nargin<3
   mask_file=[]
   mask_thresh=[]
end
if nargin < 5
   fwhm=1
end
if nargin < 6
   flip=1
end

d=fmris_read_image(input_file,1,1);
numslices=d.dim(3)
numys=d.dim(2)
numxs=d.dim(1)
numpix=numxs*numys;
D=2+(numslices>1)

i2=2:numxs;
j2=2:numys;
i1=i2-1;
j1=j2-1;
excurset=zeros(numxs,numys);
vox=[];
voxid=[];
edge1=[];
edge2=[];
n=0;

xm1=[1 1:(numxs-1)];
xp1=[2:numxs numxs];
ym1=[1 1:(numys-1)];
yp1=[2:numys numys];
zp1=[2:numslices numslices];
X=d.data*flip;

if ~isempty(mask_file)
    mask_thresh1=mask_thresh(1);
    if length(mask_thresh)>=2
        mask_thresh2=mask_thresh(2);
    else
        mask_thresh2=Inf;
    end
    d=fmris_read_image(mask_file,1,1);
    mask=d.data;
    X(mask<=mask_thresh1 | mask>mask_thresh2)=-Inf;
end

ext=input_file((length(input_file)-2):length(input_file));
isanalyze= all(ext=='img')

X2=X;
imat=(1:numxs)'*ones(1,numys)-1;
if d.vox(1)<0 & ~isanalyze
   imat=numxs-1-imat;
end
jmat=ones(numxs,1)*(1:numys)-1;
if d.vox(2)<0 & ~isanalyze
   jmat=numys-1-jmat;
end
lm=[];
lmvox=[];
resels=[];
for slice=1:numslices
   X1=X;
   X=X2;
   d=fmris_read_image(input_file,zp1(slice),1);
   X2=d.data*flip;
   if ~isempty(mask_file)
      d=fmris_read_image(mask_file,zp1(slice),1);
      mask=d.data;
      X2(mask<=mask_thresh1 | mask>mask_thresh2)=-Inf;
   end
   X3=X(xm1,:);
   X4=X(xp1,:);
   X5=X(:,ym1);
   X6=X(:,yp1);
   
   % Excursion set:
   excursetm=excurset;
   excurset=X>input_thresh;
   
   % Local maxima:
   islm=(X>=X1)&(X>=X2)&(X>=X3)&(X>=X4)&(X>=X5)&(X>=X6) ...
      &((X>X1)+(X>X2)+(X>X3)+(X>X4)+(X>X5)+(X>X6)>0)&excurset;
   lmvox=[lmvox; find(islm)+(slice-1)*numpix];
   i=imat(islm);
   j=jmat(islm);
   k=(slice-1)*ones(length(i),1);
   Xlm=X(islm);
   lm=[lm; [Xlm i j k]];

   if isstr(fwhm)
      d=fmris_read_image(fwhm,slice,2);
      resels=[resels; d.data(excurset)];
   end   
   
   % Clusters:
   voxidm=voxid;
   voxid=reshape(cumsum(excurset(:)),numxs,numys)+n;
   n=n+sum(sum(excurset));
   vox=[vox; find(excurset)+numpix*(slice-1)];
   edge=find([zeros(1,numys); excurset(i1,:).*excurset(i2,:)]);
   edge1=[edge1; voxid(edge-1)];
   edge2=[edge2; voxid(edge)];
   edge=find([zeros(numxs,1)  excurset(:,j1).*excurset(:,j2)]);
   edge1=[edge1; voxid(edge-numxs)];
   edge2=[edge2; voxid(edge)];
   if slice>1
      edge=find(excurset.*excursetm);
      edge1=[edge1; voxidm(edge)];
      edge2=[edge2; voxid(edge)];
   end

end

if n<1
   lm=[];
   return
end

% Find cluster id's in nf (from Numerical Recipes in C, page 346):
nf=1:n;
for l=1:length(edge1)
   j=edge1(l);
   k=edge2(l);
   while nf(j)~=j j=nf(j); end
   while nf(k)~=k k=nf(k); end
   if j~=k nf(j)=k; end
end
for j=1:n
   while nf(j)~=nf(nf(j)) nf(j)=nf(nf(j)); end
end   

% find the unique cluster id's corresponding to the local maxima:
ivox=find(ismember(vox,lmvox));
clmid=nf(ivox);
[uclmid,iclmid,jclmid]=unique(clmid);

% find their volumes:
ucid=unique(nf);
ucvol=hist(nf,ucid)*prod(abs(d.vox));   
if isstr(fwhm)
    ucrsl=zeros(1,length(ucid));
    for i=1:length(ucid)
        ucrsl(i)=sum(resels(nf==ucid(i)));
    end
else
    ucrsl=ucvol/fwhm^D;
end
uclmvol=ucvol(find(ismember(ucid,uclmid)));
uclmrsl=ucrsl(find(ismember(ucid,uclmid)));

% and their ranks (in ascending order):
[sortuclmvol,iuclmvol]=sort(uclmvol);
len=length(iuclmvol);
rankvol=zeros(1,len);
rankvol(iuclmvol)=1:len;

% add these to lm as extra columns:
lm=[lm rankvol(jclmid)' uclmvol(jclmid)' uclmrsl(jclmid)'];

% sort by cluster size, then local maxima, then flip:
lm=flipud(sortrows(lm,[5 1]));
lm(:,5)=len+1-lm(:,5);

return