FUNCTION lextrem_seg,ima,sm,mm,sh,mask
;
; Returns coordinates and values for local extrema in the image
;
; INPUTS: IMA = input image (2D array), might be smoothed by a
;         FWHM = SM Gaussian
;	  MM  = number; if positive, maxima are searched;
;		if negative or zero, minima are searched
;	  SH  = step (pix) parameter for STROAL.PRO routine
;	  MASK = 2D 0/1 mask of the region of interest
; OUTPUT: 2D long-integer array (table) with 4 columns
;	  0: x-coord., 1: y-coord, 2: value, 3: area (pix) of the
;				surrounding peak/dip region
;
; EXAMPLE: extrema_table = LEXTREM_SEG(ima,3,1,2,mask)
;	... smoothing by FWHM=3 pix. Gaussian, then searching for maxima
;	    with a step of 2 pixels inside the region given by 'mask'.
;
; Called routines: SCONVOL, STROAL, OBJEKT2, (TVW)
;
on_error,1

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

if n_params() lt 5 then mask=1	;default - no masking
if n_params() lt 4 then sh=1	;default sh=1
if n_params() lt 3 then mm=1	;default - maxima
if n_params() lt 2 then sm=1	;default - no smoothing

; Smoothing, segmentation, and object labelling

if sm eq 1 then seg=ima else seg=SCONVOL(ima,fwhm=sm)
seg=STROAL(seg,sh)
if mm gt 0 then seg=seg gt 0 else seg=seg lt 0
seg=ima*seg*mask		;segmented image

;  TVW,seg

OBJEKT2,seg,tab		; objects = regions with local extrema

; Finding local extrema

tab=tab(*,where(tab(0,*)))		;exclude 0 (1-pix) objects
nobj=max(tab(0,*))			;number of objects
out=lonarr(4,nobj)

for j=1,nobj do begin			;loop over objects
  pcn=tab(1,where(tab(0,*) eq j,c))	;set of j-th object coords
  if mm gt 0 then out(2,j-1)=max(seg(pcn),pos) else  $
	out(2,j-1)=min(seg(pcn),pos)	;value and subscript of max/min
  extpos=pcn(pos)			;1D position of max/min
  out(0,j-1)=extpos mod si(1)		;x-coordinate of max/min
  out(1,j-1)=extpos/si(1)		;y-coordinate of max/min
  out(3,j-1)=c				;area of the object
endfor

RETURN,out
END