FUNCTION HALPHA_CENTER, img, p

;FUNCTION CENTER_SOBEL, img, p, ccd 
;+
; NAME:
;    CENTER_SOBEL
;
; PURPOSE:
;    Calculate diameter and center coordinates of synoptic full disk
;    images. If the solar image is not centered, CENTER_SOBEL returns
;    the disk center coordinates and the diameter values are set to -1.
;
; CATEGORY:
;    BBSO Archiving System.
;    Image Processing.
;
; CALLING SEQUENCE:
;    co = CENTER_SOBEL( img, p, ccd )
;
; INPUTS:
;    IMG: Image frame.
;    P:   STRUCTURE image specs.
;
; OUTPUTS:
;    CO:  LONG array [ center coordinate X, center coordinate Y, diameter X,
;         diameter Y ].
;
; KEYWORDS:
;    None. 
;
; MODIFICATION HISTORY:
;    1995-12-04: Anders Johanneson, Big Bear Solar Observatory
;        Original version of AJEDCENT.PRO.
;    1999-12-13: Carsten Denker, Big Bear Solar Observatory
;        Intergration in the new BBSO archiving system, variables are passed 
;        as structures, off-center images via CO keyword, information via
;        MESSAGE procedure.
;    2000-05-18: Michael Steinegger, Big Bear Solar Observatory
;        The Sobel function is only used to obtain the first guess for the
;        center and radius, the exact values are obtained by fitting a
;        circle to the Sobel function using FIT_CIRCLE.PRO.
;    2000-05-23: Michael Steinegger, Big Bear Solar Observatory
;        Version HALPHA_CENTER.PRO to be used by other programs like e.g.
;        HALPHA_FLAT.PRO. Does not need structure CCD.
;-
        
;-----  find the limb by using a non-linear edge enhancement operator
 
;lim = MEDIAN ( img ) / 4. This might be a potential problem!!!
lim = MEDIAN ( img ) / 5.
mask = ( SOBEL( img ) GT lim ) 
   
;-----  elminate problems with borders and a few bad columns 

;type = STRMID( ccd.fdr, 0, 1 )
type = 'h'	; to be used for H-alpha images
CASE type OF
   'h':  xo = 10
   ELSE: xo = 50
ENDCASE
b = 4
mask( *, p.ny - b -1 : * ) = 0 
mask( *, 0 : b - 1 ) = 0 
mask( 0 : b - 1 + xo, * ) = 0 
mask( p.nx - b - 1 : *, * ) = 0
                       
;-----  get x-coordinates and diameter

;xx = TOTAL( mask, 2 ) GT 0	;-----  horizontal sum
xx = TOTAL( mask, 2 ) GT 3	; this should give a better radius

;-----  solar image is centered in x   
                  
IF ( xx( b + xo ) EQ 0 ) AND ( xx( p.nx - b - 2 ) EQ 0 ) THEN BEGIN 
    index = WHERE( xx, n )
    x1 = index( 0 )
    x2 = index( n - 1 )
    xd = x2 - x1 + 1
    xc = xd / 2 + x1
ENDIF

;-----  solar image is on the left

IF ( xx( b + xo ) EQ 1 ) AND ( xx( p.nx - b - 2 ) EQ 0 ) THEN BEGIN
    xx( 0 : b + xo - 1 ) = 1	;-----  set border to 1 
    index = WHERE( 1 - xx )
    x2 = index( 0 ) - 1 
    xd = p.xd
    xc = x2 - p.xd / 2 
ENDIF 

;-----  solar image is on the right

IF ( xx( b + xo ) EQ 0 ) AND ( xx( p.nx - b - 2 ) EQ 1 ) THEN BEGIN
  index = WHERE( xx )
  x1 = index( 0 )
  xd = p.xd
  xc = x1 + p.xd / 2 
ENDIF

;-----  get y-coordinates and diameter

;yy = TOTAL( mask, 1 ) GT 0	;-----  vertical sum
yy = TOTAL( mask, 1 ) GT 3	; this should give a better radius

;-----  solar image is centered in y

IF ( yy( b ) EQ 0 ) AND ( yy( p.ny - b - 2 ) EQ 0 ) THEN BEGIN 
    index = WHERE( yy, n )
    y1 = index( 0 )
    y2 = index( n - 1 )
    yd = y2 - y1 + 1
    yc = yd / 2 + y1 
ENDIF 

;-----  solar image is on the bottom

IF ( yy( b ) EQ 1 ) AND ( yy ( p.ny - b - 2 ) EQ 0 ) THEN BEGIN 
  yy( 0 : b - 1 ) = 1		;-----  set border to 1 
  index = WHERE( 1 - yy )
  y2 = index( 0 ) - 1 
  yd = p.yd
  yc = y2 - p.yd / 2
ENDIF

;-----  solar image is on the top

IF ( yy( b ) EQ 0 ) AND ( yy( p.ny - b - 2 ) EQ 1 ) THEN BEGIN 
  index = WHERE( yy )
  y1 = index( 0 ) 
  yd = p.yd
  yc = y1 + p.yd / 2 
ENDIF


;-----  eliminate center of the disk (active regions!)

r  = FIX( ROUND( ( xd + yd ) / 4. ) )
DIST_CIRCLE, d, [ p.nx, p.ny ], xc, yc
mask( WHERE( d LT r * 0.95 OR d GT r * 1.1 ) ) = 0

;-----  fit a circle to obtain radius and center coordinates,
;       iterate until all points are within 2 * sigma of radius
            
bad = [0L]
WHILE bad( 0 ) NE -1 DO BEGIN 
   wo = WHERE( mask GT 0 )
   xfit = wo MOD p.nx
   yfit = wo / FLOAT( p.nx )
   result = FIT_CIRCLE( xfit, yfit, [ xc, yc, r ] )
   radius = SQRT( ( xfit - result( 0 ) ) ^ 2. + ( yfit - result( 1 ) ) ^ 2. )
   sigma  = SQRT ( TOTAL( ( radius - result( 2 ) ) ^ 2. ) / $
                 ( N_ELEMENTS ( wo ) - 1 ) )
   bad = WHERE( radius GT ( result( 2 ) + 2. * sigma ) OR $
                radius LT ( result( 2 ) - 2. * sigma ) ) 
   IF bad( 0 ) NE -1 THEN BEGIN
      mask( xfit( bad ), yfit( bad) ) = 0
      xc = FIX( ROUND ( result( 0 ) ) )
      yc = FIX( ROUND ( result( 1 ) ) )
      xd = FIX( ROUND ( result( 2 ) * 2. ) )
      yd = xd
      r  = ( xd + yd ) / 4.
;;      print,sigma,bad(0),N_ELEMENTS(bad), xc, yc, xd
   ENDIF
ENDWHILE

yd = xd 

RETURN, [ xc, yc, xd, yd ]
 
END