FUNCTION lextrem_box,ima,mm,sb
;
; Returns coordinates and values for local extrema in the image
;
; INPUTS: IMA = input image (2D array)
;	  MM  = number; if positive, maxima are searched; 
;		if negative or zero, minima are searched
;	  SB  = side-length of a box (pixels)
; OUTPUT: 2D long-integer array (table) with 3 columns
;	  0: x-coord., 1: y-coord, 2: value
;
; EXAMPLE: extrema_table = LEXTREM_BOX(ima,-1,15) 
;		... searching for minima with a box of 15x15 pix.
;
; Method: Image is divided into boxes sb x sb pixels and
;	  an extremum is found in each box. Extrema found
;	  at the borders of boxes are discarded because
;	  they do not correspond to the definition.
;
on_error,1

si=size(ima)
if si(0) ne 2 then message,'Input image is not a 2D array !!'

; Division into boxes

sb=fix(sb)
nx=si(1)/sb	;integer division, no. of boxes in x
ny=si(2)/sb	;integer division, no. of boxes in y
if (si(1) mod sb) lt 2 then nx=nx-1
if (si(2) mod sb) lt 2 then ny=ny-1
; right and upper borders remain unexplored...

out=lonarr(3,nx*ny)
cnt=0

; Loop over boxes

for j=0,ny-1 do begin
  for i=0,nx-1 do begin
	xcor=i*sb
	ycor=j*sb
	box=ima(xcor:xcor+sb+1,ycor:ycor+sb+1)  ; overlap by 2 pixels
	if mm gt 0 then val=max(box,pos) else val=min(box,pos)
	x=pos mod (sb+2)
	y=pos/(sb+2)
; Discard extrema at the edges of boxes
	if (x ne 0 and x ne sb+1) and (y ne 0 and y ne sb+1) then begin
	  out(0,cnt)=x+xcor
	  out(1,cnt)=y+ycor
	  out(2,cnt)=val
	  cnt=cnt+1
	endif
  endfor
endfor

w=where(out(2,*))	;remove empty rest of out
RETURN,out(*,w)

END