function FindStars, Image, threshold=Threshold, maxstars=MaxStars, nan=NaN,		$
		starpix=StarPix, peakpix=PeakPix, starpos=StarPos, peakpos=PeakPos, btest=bTest

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	FindStars
; PURPOSE:
;	Locate stars in 2D image based on a 2D template defining the shape of a star
;	in the image
; CATEGORY:
; CALLING SEQUENCE:
; INPUTS:
;	Image		2D array, any type
;					Image containing stars
; OPTIONAL INPUT PARAMETERS:
;	Threshold		scalar, any type
;						only maxima above Threshold are considered
;	MaxStars		scalar, integer
;						maximum number of stars to be removed
;	/nan			if set pixels in stars will be replaced by ValuesNaN
; OUTPUTS:
;	Image		2D array, same type as input
;						image with stars removed
; OPTIONAL OUTPUT PARAMETERS:
;	StarPix		1D array, type long integer
;					indices of pixels in stars
;	BadPix		1D array, type long integer
;					indices of areas around local maxima not matching the star criterium
; CALLS:
;	InitVar, IsType, ArrayLocation, echo, FitStar, twin
; PROCEDURE:
; >	The image is scanned from the maximum down. A box around the maximum is
;	matched to a star template to decide whether the maximum is a star.
; >	If a star is located an area covering the star is set to the median of the
;	box used to test for a star (or to NaN if the /nan switch is set).
;	The star location are accumulated in StarPix
; >	Local maxima not matching the star template are collected in BadPix.
;	The intensities of these maxima are not modified.
;	Ideally these would be cosmic rays. Some of them may be stars which for some
;	reason did not fit well enough to qualify as a star.
; >	The search for stars can be limited in two ways: the 'maxstars' and the 'threshold'
;	keywords. The 'threshold' keyword is probably preferable and should be set to
;	a value sufficiently far above the noise in the image.
; MODIFICATION HISTORY:
;	OCT-1998, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

sz = size(Image)
Frame = Image


InitVar, bTest, /key
InitVar, NaN, /key
InitVar, MaxStars , n_elements(Frame)

if IsType(Threshold, /undefined) then Threshold = median(Image)

BadPix  = -1
StarPix = -1
StarPos = -1
PeakPix = -1
PeakPos = -1

nP = 0
nQ = 0
nR = 0
nRMax = 25

bScan = finite(Frame)
iScan = where (bScan)

while nR lt nRMax and iScan[0] ne -1 do begin
	Maximum = max(Frame[iScan],/nan,Star)
	Star = iScan[Star]

	NoFit = nP ge MaxStars or Maximum lt Threshold
	bStar = FitStar(Frame,bScan,Star,Box,Mask,BoxBase,BoxDev,BoxMed,Residue,devratio=DevRatio,nofit=NoFit,nearbystar=NearbyStar)
	if not NoFit and not bStar then	begin
		if DevRatio lt .7 then	$
			bStar = FitStar(Frame,bScan,[Star,Star],Box,Mask,BoxBase,BoxDev,BoxMed,Residue,devratio=DevRatio)
	endif

	S  = ArrayLocation(Star, size=sz)				; X,Y position of maximum

	W  = (size(Box))[1:2]
	XS = lindgen(W[0])# replicate(1L,W[1])			; X-coordinate in Box
	YS = replicate(1L,W[0])# lindgen(W[1])			; Y-coordinate

	W	= (W-1)/2
	XS	= S[0]-W[0]+XS								; X-coordinate in Frame around maximum in S
	YS	= S[1]-W[1]+YS								; Y-coordinate

	if bTest then begin
		twin,/show,/bl
		surface, Frame[XS,YS]
	endif

;	Frame[XS[Mask],YS[Mask]] = BoxMed
	Frame[XS[Mask],YS[Mask]] = !VALUES.F_NAN
	bScan[XS[Mask],YS[Mask]] = 0B

	MaskPix = XS[Mask]+YS[Mask]*sz[1]

	if bStar then begin								; Identified a star
		nP = nP+1
		print, 'Star'+strcompress(nP)+' at ['+strcompress(S[0])+','+ $
			strcompress(S[1])+'], peak intensity '+strcompress(Maximum)+', '+strcompress(DevRatio)

		if NearbyStar then begin
			D = min(((StarPos mod sz[1])-S[0])^2+(StarPos/sz[1]-S[1])^2,i)
			D = ArrayLocation(StarPos[i],size=sz)
			message, /info, 'Nearby star'+strcompress(i+1)+' at ['+strcompress(D[0])+','+strcompress(D[1])+']
		endif

		if BadPix[0] ne -1 then begin
			if PeakPix[0] eq -1 then begin
				PeakPix = BadPix
				PeakPos = BadPos
			endif else begin
				PeakPix = [PeakPix,BadPix]
				PeakPos = [PeakPos,BadPos]
			endelse
			BadPix = -1
		endif

		if StarPix[0] eq -1 then begin
			StarPix = MaskPix
			StarPos = Star
		endif else begin
			StarPix = [StarPix,MaskPix]
			StarPos = [StarPos,Star]
		endelse

	    nR = 0
	endif else begin

		nQ = nQ+1
		print, 'Peak'+strcompress(nQ)+' at ['+strcompress(S[0])+','+ $
			strcompress(S[1])+'], peak intensity '+strcompress(Maximum)+', '+strcompress(DevRatio)

		if BadPix[0] eq -1 then begin
			BadPix = MaskPix
			BadPos = Star
		endif else begin
			BadPix = [BadPix,MaskPix]
			BadPos = [BadPos,Star]
		endelse

	    nR = nR+1
	endelse

	if bTest then begin
		twin,/show,/tl
		surface, Frame[XS,YS]
		i = 0
		echo, 'pause', i
		if i eq 1 then message, 'stop
	endif

	iScan = where(bScan)
endwhile

if PeakPix[0] ne -1 then begin
	PeakPix = PeakPix[uniq(PeakPix,sort(PeakPix))]	; Remove double entries

	BadPix  = replicate(1B,n_elements(PeakPix))		; Remove entries that are also on StarPix list
	for i=0,n_elements(PeakPix)-1 do if (where(PeakPix[i] eq StarPix))[0] ne -1 then BadPix[i] = 0B
	PeakPix = PeakPix[where(BadPix)]

	Frame[PeakPix] = Image[PeakPix]
endif

; At this point 'StarPix' contains all the indices assigned to stars.
; Frame[StarPix] is set to the frame value after removal of the star. All other
; pixels are the same as in the input array Image

if StarPix[0] ne -1 then begin
	StarPix = StarPix[uniq(StarPix,sort(StarPix))]
	if NaN then Image[StarPix] = ValuesNan(sz) else Image = Frame
endif

return, nP  &  end