;+
; NAME:
;	smei_star_lsqbins
; PURPOSE:
;	Refine PSF and background area definition based on
;	properties of the original star.
; CATEGORY:
; CALLING SEQUENCE:
	PRO smei_star_lsqbins, original_star, rr, rr0, in_psf, in_bkgnd, bkgnd_count=bkgnd_count
; INPUTS:
;	original_star	array[nx,ny]; type: float
;						area from original skymap around the star to be subtracted
;						(current a 50x50 box i.e 5x5 degrees)
;	rr				array[2,nx*ny]; type: integer??
;						column and row indices into original skymap for all bins
;						in original_star
;	rr0 			array[2,nx*ny]; type: float
;						x,y locations relative to the centroid of the star
;						(in the same units as rr, i.e. bin spacings)
;	in_psf			array[nx,ny]; type: byte
;						1: bin is inside the PSF
;						0: bin is outside the PSF
;	in_bkgnd		array[nx,ny]; type: byte
;						1: bin is included in background area
;						0: bin is not included in background area
; OPTIONAL INPUT PARAMETERS:
;	bkgnd_count=bkgnd_count
;					scalar; type: float (<=1.0) or integer (n > 1)
;						sets the maximum number of pixels to be retained for
;						the background area (in the in_bkgnd array).
;						bkgnd_count >  1  : maximum number of bins to be retained
;						bkgnd_count <= 1.0: maximum fraction of current number of
;											background bins to be retained
; OUTPUTS:
; OPTIONAL OUTPUT PARAMETERS:
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
; PROCEDURE:
;	The first filter on in_psf and in_bkgnd is that only bins with finite
;	values in original star are left non-zero.
;
;	If bkgnd_count is defined than the background array in_bkgnd is modified
;	so that at most bkgnd_count bins are left non-zero. These will be the bins
;	in the background with the lowest value in original_star.
;	Note that in_bkgnd is NOT modified if there are already less than
;	bkgnd_count bins in it.
; MODIFICATION HISTORY:
;	AUG-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

finite_star = finite(original_star)

in_psf   = finite_star AND in_psf  
in_bkgnd = finite_star AND in_bkgnd

IF IsType(bkgnd_count,/defined) THEN BEGIN

	i_bkgnd = where(in_bkgnd, n_bkgnd)

	; bkgnd_count is either a number of bins or a fraction <= 1.0
	; Translate fraction to number of bins

	n_cutoff = bkgnd_count
	IF n_cutoff LE 1 THEN n_cutoff = floor( (n_cutoff < 1)*total(in_bkgnd) )

	IF n_bkgnd GT n_cutoff THEN BEGIN	; More than enough bins in background

		a = original_star[i_bkgnd]		; Values in background bins
		a = (sort(a))[n_cutoff:*]		; Sort and find indices of pixel with high values
		in_bkgnd[ i_bkgnd[a] ] = 0		; Remove bins with high values from background

	ENDIF

ENDIF

RETURN  &  END