


pro center_sobel, img, p,para
;+
; 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 verion 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-07-20: Michael Steinegger, BBSO/NJIT
;        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.
;-

;----- find the limb by using a non-linear edge enhancement operator
mask=img
;mask = smooth(img,4)
;lim = MEDIAN ( mask ) / 4.
;lim = MEDIAN( mask ) / 2.	;----- this works better for YNAO images
;lim = MEAN( mask ) / 3.	;----- now this seems to work better...
;lim=mean(mask)/6.
lim=total(mask)/n_elements(mask)/6.
vmin=5
mask = ( SOBEL( TEMPORARY( mask ) ) GT lim )


;----- elminate problems with borders and a few bad colums

b = 4
;mask( *, p.ny-10 : * ) = 0
;mask( *, 0 : b - 1 ) = 0
;mask( 0 : b - 1 + xo, * ) = 0
;mask( p.nx - 8 : *, * ) = 0
;----- get x-coordinates and diameter
;;xx = TOTAL( mask, 2 ) GT 0    ;----- horizontal sum
xx = TOTAL( mask, 2 ) GT vmin	; this should give a better radius
;----- solar image is centered in x
IF ( xx( b ) EQ 0 ) AND ( xx( p.nx - b - 1 ) 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 ) EQ 1 ) AND ( xx( p.nx - b - 1 ) EQ 0 ) THEN BEGIN
    xx( 0 : b - 1 ) = 1    ;----- set border to 1
    index = WHERE( 1 - xx )
    x2 = index( 0 ) - 1
    xd = -1
    xc = x2 - p.xd / 2
ENDIF

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

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

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

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

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

IF ( yy( b ) EQ 0 ) AND ( yy( p.ny - b - 1 ) 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 - 1 ) EQ 0 ) THEN BEGIN
  yy( 0 : b - 1 ) = 1    ;----- set border to 1
  index = WHERE( 1 - yy )
  y2 = index( 0 ) - 1
  yd = -1
  yc = y2 - p.yd / 2
  IF yc LT 0 THEN yc = y2 + yc		;---- solves some occacional problems

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 = -1
  yc = y1 + p.yd / 2
ENDIF

;------ here starts the new part (M. Steinegger, 2000/07/20)

;; if the xc and yc we've calculated are obviously bogus, put in a guess
;if (xc < 900 OR xc > 1100) then xc = 1000
;if (yc < 900 OR yc > 1100) then yc = 1000
;----- check initial guesses for diameter
if (n_elements(xd) eq 0 or n_elements(xc) eq 0 or n_elements(yd) eq 0 $
	or n_elements(yc) eq 0) then hr_pos_flatfield,img,p,xd,yd,xc,yc
IF xd EQ -1 AND yd EQ -1 THEN xd = p.xd
IF xd EQ -1 AND yd NE -1 THEN xd = yd
IF yd EQ -1 AND xd NE -1 THEN yd = xd

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

r  = FIX( ROUND( ( p.xd + p.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
mask( WHERE( d LT r * 0.9 OR d GT r * 1.1 ) ) = 0	; for KSO images

;-----  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.
      r  = ( p.xd + p.yd ) / 4.

;;      print,sigma,bad(0),N_ELEMENTS(bad), xc, yc, xd

   ENDIF
ENDWHILE	
yd = xd

para=[ xc, yc, xd,yd]

END





