function lm = locmax(input_file,threshold);

%LOCMAX
%
% LM = LOCMAX( INOUT_FILE, THRESHOLD)
%
% Finds local maxima of INPUT_FILE above THRESHOLD. A local maximum 
% is a voxel which is greater than or equal to all its 6 neighbours,
% and strictly greater than at least one of them.
%
% LM is a matrix with 4 columns. 
% Col 1: values of local maxima, sorted in descending order.
% Cols 2-4: x,y,z voxels of local maxima in C notation i.e. starting at 0.
% If x or y step sizes are negative, then their voxel indices are reversed 
% so that they agree with 'register'.

in_handle=openimage(deblank(input_file));

numslices=getimageinfo(in_handle, 'NumSlices')
numys=getimageinfo(in_handle, 'ImageHeight')
numxs=getimageinfo(in_handle, 'ImageWidth')
steps=getimageinfo(in_handle, 'Steps')

xm1=[1 1:(numxs-1)];
xp1=[2:numxs numxs];
ym1=[1 1:(numys-1)];
yp1=[2:numys numys];
zp1=[2:numslices numslices];
X=reshape(getimages(in_handle,1),numxs,numys);
X2=X;
imat=(1:numxs)'*ones(1,numys)-1;
if steps(1)<0
   imat=numxs-1-imat;
end
jmat=ones(numxs,1)*(1:numys)-1;
if steps(2)<0
   jmat=numys-1-jmat;
end
lm=[];
for slice=1:numslices
   X1=X;
   X=X2;
   X2=reshape(getimages(in_handle,zp1(slice)),numxs,numys);
   X3=X(xm1,:);
   X4=X(xp1,:);
   X5=X(:,ym1);
   X6=X(:,yp1);
   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)  ...
      &(X>threshold);
   i=imat(islm);
   j=jmat(islm);
   k=(slice-1)*ones(length(i),1);
   Xlm=X(islm);
   lm=[lm; [Xlm i j k]];
end
lm=flipud(sortrows(lm,1));

% end