;+
; NAME:
;	MEANCENT
; PURPOSE:
;	To locate the center and determine the size of the image of the
;       Portable Low Brightness Source in a Polar Camera image.
; CATEGORY:
;	Calibration utility
; CALLING SEQUENCE:
;	MEANCENT, IMAGE, C_COL, C_ROW, TOTPIX $
;               , DEBUG = debug, RADIUS = radius, NOTV = notv
; INPUTS:
;	IMAGE:  The Polar Camera image containing the PLBS image.
; OPTIONAL INPUTS:
;	None.
; KEYWORD PARAMETERS:
;	/DEBUG: If set, additional information is printed.
;       RADIUS: If a named variable is supplied, the mean radius
;               (pixels) of the PLBS screen image is returned.
;       /NOTV:  If set, defeats the display of the image as it is
;               analyzed.
; OUTPUTS:
;	C_COL:  The image column of the center of the screen image.
;       C_ROW:  The image row of the center of the screen image.
;       TOTPIX: The number of pixels in the screen image.
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	A window is created (unless /NOTV is set) and the supplied image
;       is displayed and manipulated.
; RESTRICTIONS:
;	Unpredictable results are returned if the image doesn't show
;       the PLBS screen.
; PROCEDURE:
;	Only pixels > 25% of the 99.9th percentile pixel value are
;       considered.  Outlying pixels are removed using the DILATE and
;       ERODE operators.  The geometric centroid of the remaining pixels
;       is returned as the center of the PLBS screen image.  Pixels more
;       than 50 pixels away from the center are eliminated.  If the
;       number of pixels so eliminated is more than 1% of the total an
;       error message is displayed.  The SOBEL operator is applied to
;       find the edge of the PLBS screen image.  The radius of the
;       screen image is calculated from the mean distance of the edge
;       pixels from the center.  If this distance is inconsistent with
;       the number of pixels in the image an error is reported.  If
;       either error condition occurs null results are returned.
; EXAMPLE:
;
; SEE ALSO:
;	REDUCALS.PRO
; MODIFICATION HISTORY:
; 	Written by:	DPSteele, ISR, spring 1994
;       Documented by:  DPSteele, ISAS, May 1996.
;-
PRO meancent,img,c,r,tm,debug=debug,radius=rad,notv=notv

IF N_PARAMS() EQ 0 THEN BEGIN
    doc_library,'meancent'
    RETURN
    ENDIF

; Decide whether to display images
dispim=NOT KEYWORD_SET(notv)

; save !ORDER
oldorder=!ORDER
!ORDER=0

; work with 256x256 images
s=SIZE(img)
IF s(1) NE 256 THEN tmp=REBIN(img,256,256) ELSE tmp=img

; display the image
IF dispim THEN BEGIN
    WINDOW,/FREE,XSIZE=256,YSIZE=256
    wind=!D.WINDOW
    TVSCL,tmp
    ENDIF

; select pixels viewing LBS
;maxpix=MAX(tmp)
maxpix=pcentile(tmp,0.999)
sf=0.25

IF Max(tmp) LE sf*maxpix THEN BEGIN
    Stop
    c=-1
    r=-2
    tm=-1
    rad=-1
    Return
    ENDIF
    
mk,omask			; original mask
omask(WHERE(tmp GT sf*maxpix))=1

; remove outlying pixels by opening the mask
cross=[[0,1,0],[1,1,1],[0,1,0]]
mask=DILATE(omask,cross)	; altered mask
box=REPLICATE(1,3,3)
mask=ERODE(mask,box)
mask=ERODE(mask,box)

; generate column and row arrays
x=BYTE(INDGEN(256)#REPLICATE(1,1,256))
y=TRANSPOSE(x)

; form sums for calculating mean column and mean row
tmx=TOTAL(mask*x)		; altered mask
tmy=TOTAL(mask*y)
tm=TOTAL(mask)
IF KEYWORD_SET(debug) THEN PRINT,tmx,tmy,tm
IF tm GE 10 THEN BEGIN	; image actually shows an extended bright object
    c=tmx/tm
    r=tmy/tm

; set outlying bright pixels to zero
    cols=omask*x			; original mask
    rows=omask*y
    outliers=WHERE((SQRT((cols-c)^2+(rows-r)^2) GT 50) AND $
                   omask,noutlier)
    IF noutlier GT 0.01*tm THEN BEGIN
        outstr=STRING(noutlier,tm $
                     ,FORMAT='(I0," of ",I0," pixels are outliers!")')
        MESSAGE,outstr,/INFORMATIONAL
        ENDIF
    IF noutlier GT 0 THEN tmp(cols(outliers),rows(outliers))=0
    IF dispim THEN TVSCL,tmp

; estimate radius of LBS image using edge of LBS image
    stmp=SOBEL(tmp)
    edge=WHERE(stmp GT MAX(stmp)/2.)
    one2two,edge,stmp,ce,re

; calculate mean distance from center to edge
    rad=mean(SQRT((ce-c)^2+(re-r)^2))
    IF dispim THEN arcs,rad,0,360,c,r,/dev

; Check validity of radius
    rad_total=!PI*rad*rad
    error=(ABS((rad_total-tm)/tm) GT 0.2)

    IF error THEN BEGIN
        errstr=STRING(tm,rad_total $
                     ,FORMAT='("Bright pixels: ",I0,"; area of image: ",I0)')
        MESSAGE,errstr,/INFORMATIONAL
        ENDIF

; If the image model is invalid, return null results.

;    IF (N_ELEMENTS(outstr) GT 0) OR error THEN BEGIN
    IF error THEN BEGIN
        c=-1
        r=-1
        tm=-1
        rad=-1
        ENDIF

    ENDIF ELSE BEGIN
    MESSAGE,'"Bright feature" contains <10 pixels!  Returned center meaningless',/INFORMATIONAL
    c=-999
    r=-999
    rad=-999
    ENDELSE

; restore !ORDER
!ORDER=oldorder
IF dispim THEN BEGIN
    WAIT,5
    WDELETE,wind
    ENDIF

RETURN
END