pro fit_limb, image, x0, y0, r0, r_err, $ index=index, $ summed_pixels=summed_pixels, $ pixel_size=pixel_size, $ verbose=verbose, $ oblateness=oblateness,ob_angle=ob_angle, bias=bias,xlimb=xfit,$ ylimb=yfit, interactive=interactive, quiet=quiet, $ initial_values=initial,fast=fast, xinitial=hxa_x, $ yinitial=hxa_y, $ decompress=decompress, noiterate=noiterate,iterate=iterate, $ kill=RLOOP,norfit=norfit,date=date, $ noobcorrect=noobcorrect,ellipse=ellipse,resolution=resolution ;+ ;NAME: ; fit_limb ;PURPOSE: ; Find the center and radius of an image. NOTE: For SXT images, the ; default is to return the result in full resolution pixels. ;CATEGORY: ;CALLING SEQUENCE: ; fit_limb, image, x, y, r, pixel_size=2.5, date='05-Oct-93' ; fit_limb, image, x, y, r, pixel_size=2.5, intial=[434,512,403*2.5] ; fit_limb, data(*,*,nimage), x, y, r, index=index(nimage),/norfit ;INPUTS: ; image Full sun image for which the center and radius will ; be computed. It should be a 2-D array (not a 3-D data cube) ; It should be decompressed, and ideally background subtracted. ;OPTIONAL KEYWORD INPUTS: ; index SXT index structure. If passed, it will use this value to ; derive the initial guess of the location. ; summed_pixels - If set, return the pixel coordinate of the sun center in ; the summed pixels, not the full resolution pixel coordinates ; which are the default (only relevant if index is set). ; pixel_size - The pixel size in arcseconds. If the SXT index is passed ; it will be derived from that. ; /verbose = if set, print a bunch of diagnostics ; /ellipse = if set, fit an ellipse rather than a circle. May work ; better for oblate images, but should not be used for circular ; images since it allows another degree of freedom in the fit. ; oblateness = return second harmonic of radius ; ob_angle = return phase of 2nd harmonic ; /noobcorrect = do not correct final sun center value for oblateness ; bias = return distortion of limb due higher harmonics ; xlimb,ylimb = returns the coordinates of the selected limb pixels ; /interactive = allows interactive selection of initial center ; guess and deselection of limb points. ; /quiet = No images are displayed if set and interactive is not set ; initial_values = 3-element vector giving initial guess: [x0,y0,r0] ; x0 and y0 are in pixels, r0 is in arcseconds !! ; If initial values are passed, HXA is NOT used and index ; need not be passed as a parameter. ; /fast = If set, do the fast but slightly less accurate computation of ; the limb position (limb derivative is not smoothed). ; xinitial = returns x coordinate of image center for the initial guess ; yinitial = returns y coordinate of image center for the initial guess ; /decompress = If set, use sxt_decomp to decompress the image ; /iterate = iterate in the circle fitting routine. In principle this ; could improve the accuracy. (default) ; /noiterate = do not iterate. ; kill = number of iterations eliminating bad limb points. The ; default is 0 for non x-ray images and 2 for x-ray images. ; X-ray images often need kill of 1 or 2. ; Smaller numbers (0 or 1) make the algorithm ; correspondingly less sensitive to the initial guess. ; Higher numbers will eliminate more limb points with the potential ; for a better fit in some cases. ; /norfit = if set, do not fit the radius of the sun: use the initial guess ; If you know the pixel size exactly and pass the date or index ; so that the sun's radius can be computed, this keyword should ; be used since it removes one degree of freedom from the circle ; fit. ; date = date for computation of solar radius. Index is used if it is ; present. If date is omitted, a crude value is computed from the ; image. ;OUTPUTS: ; x = x position of the center of the image. ; For SXT images, the units are FULL res pixels unless /SUMMED_PIXELS ; is set. There is also (0, 1 or 3) pixel offset for FR, HR, and QR. ; y = y position of the center of the image ; For SXT images, the units are FULL res pixels unless /SUMMED_PIXELS ; is set. There is also (0, 1 or 3) pixel offset for FR, HR, and QR. ; r = radius of the image in units of (original) pixels ; ; RSun/RFit = The ratio of the Sun's radius computed for the time of the ; image to the fit radius. This is the size of the SXT pixels ; in arcseconds. However, if the index in NOT suppled, the ; Sun's radius can not be computed and the result is the ; ratio of the initial guess of the radius to the fitted value. ; In this case, this is NOT the size of the SXT pixels. ; RSun/RFit is printed as the program ends. It is not returned ; in a variable. ;SIDE EFFECTS: ; Makes plots in the current window. ;RESTRICTIONS: ; o The algorithm relies on a VERY good initital guess to remove bad limb ; points. ; o The correction for oblateness is not done for x-ray images. ; o If part of the limb is obscurred, you may need to use a large kill ; value (around 5). Such data must be treated on a case by case basis ; by a smart user! ;PROCEDURE: ; The initial ; guess is used to find a rough limb position which is used to compute ; the points with the maximum derivative in the intensity along the radii. ; These maximum derivative points are assumed to delineate the limb and a ; circle is fit to these points and the center and the radius are computed. ; The derivative is smoothed using the deriv_lud procedure. ;MODIFICATION HISTORY: ; T. Metcalf 4/1992 ; T. Metcalf 10/1992: Switched from gt_hxa to hxa_suncenter to get the ; initital guess. ; T. Metcalf 2/1993: Switched from hxa_suncenter to get_suncenter ; *** COPIED FROM SXT_CENTER.PRO TO FIT_LIMB.PRO *** ; 21-Dec-93 (MDM) - Modified to work with non-SXT images ; 6-Jan-94 (MDM) - Minor changes ; 14-Jan-94 (MDM) - Removed FID_FLAG keyword (was marking if RD_PNT ; data had fidutial mark) ; 19-Jan-94 (MDM) - Corrected error not allowing it to work on SXT image ; - Modified header. ; 3-Feb-94 (MDM) - Modified to use GT_CORNER with /FROM_SC switch ; to determine the initial sun center guess. ; T. Metcalf 2/1994: o Default kill for non x-ray images is now 0. ; o Added obcorrect keyword. ; o Added ellipse keyword ; 15-Feb-94 (MDM) - Removed on_error,2 statement ; - Corrected error where IF statement was trying to ; operate on a one byte array, not a scalar ; (the problem was in GET_RB0P) ; 21-Feb-94 (MDM) - Modified the header ; 9-Dec-94 (MDM) - MOdified the header ; 22-Mar-95 (MDM) - Modified to add a 0.5 pixel offset for ; the results when not passing an SXT index. ;- ;on_error,2 CIRCLE_ITER = 50 ; iteration limit for circle fitting x0 = 0. y0 = 0. r0 = 0. r_err = 0. ob_angle = 0. oblateness = 0. bias = 0. verbose = keyword_set(verbose) if keyword_set(noiterate) then noiterate=1 else noiterate=0 ;;if keyword_set(iterate) then noiterate=0 else noiterate=1 if keyword_set(nordpnt) then begin if n_elements(initial) NE 3 then interactive=1 endif qindex = keyword_set(index) save_P_multi = !P.MULTI !P.MULTI = [0,3,2,0,0] ;---------- Get the image if (size(image))(0) NE 2 then message,'The image must be a 2-D array' if qindex then begin if n_elements(index) NE 1 then $ message,'Index must have a single element matching the image' endif nx = n_elements(image(*,0)) ny = n_elements(image(0,*)) if keyword_set(decompress) and qindex then image = sxt_decomp(image) isxray = 0 if (qindex) then begin isxray=(gt_filtb(index) GT 1) pixel_size = gt_pix_size(index) ;pixel size in arcsec end else begin if (n_elements(pixel_size) eq 0) then begin print, 'FIT_LIMB: Need the pixel size' return end end max_value = max(image) ;?? Why did tom set to 255 or 4096 ;-------------------- Determine the intial values if (n_elements(initial) EQ 3) then begin ;user passed them in xcen = initial(0) ycen = initial(1) radius = initial(2) endif else begin if (qindex) then begin ; ---- Compute the solar radius from the ephemeris radius = get_rb0p(index,/radius) ;;radius = rb0p(0) ;in arcsec ;;cen = sxt_cen(index) ;MDM removed 3-Feb-94 ;;xcen = cen(0) / (2^gt_res(index)) ;;ycen = cen(1) / (2^gt_res(index)) resolution = 2^gt_res(index) corner = gt_corner(index, /from_sc) xcen = -corner(0) ycen = -corner(1) hxa_x = xcen*resolution hxa_y = ycen*resolution endif else begin pp = where(image GT 0.3*max_value) img = bytarr(nx,ny) img(pp) = 1 centroidw,img,xcen,ycen if keyword_set(date) then radius = get_rb0p(date,/radius) $ else radius = sqrt(total(img(pp))/!pi)*pixel_size end end init_radius = radius radius = radius/pixel_size ;Convert arcsec to pixels ;-------------------- Now get a guess for the center of the circle from the user if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then begin erase tv,bytscl(image) xyouts,.15,.95,"Initial Guess",align=0.5,/norm,charsize=1.5 endif if keyword_set(interactive) then begin print print,"Guess Disk Center as Accurately as Possible:" print print,"Click left mouse button repeatedly to get the center of the image" print,"Click middle mouse button to redraw" print,"Click right mouse button to continue when you get close" endif theta = findgen(361) * !dtor !ERR = 0 REPEAT begin if !ERR LT 4 then begin if !ERR GE 2 then tv,bytscl(image) else if !ERR EQ 1 then begin xcen = x ycen = y endif xcir = xcen + radius * cos(theta) ycir = ycen + radius * sin(theta) if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then plots,xcir,ycir,/device endif if keyword_set(interactive) then cursor,x,y,3,/device endrep UNTIL (!ERR GE 4) OR (NOT keyword_set(interactive)) if keyword_set(interactive) then erase if NOT keyword_set(quiet) then begin print print,"Working ..." print endif ;-------------------- Start the work of fitting the limb ; Number of points on either side of limb to check for the inflection ; point in intensity. Increase if there was a fiducial mark. n_points = 10 ; Compute the positions of the inflection points across the limb for t=0,n_elements(theta)-1 do begin ;Loop through angles, 1 degree per step inten = fltarr(n_points*2+1) ; inten storesthe intensity along a radius off_screen = 0 for r=-n_points,n_points do begin ; Get the intensity along a radius x = fix(xcen + (radius+r)*cos(theta(t))) y = fix(ycen + (radius+r)*sin(theta(t))) if NOT (x GT (nx-1) OR x LT 0 OR y GT (ny-1) or y LT 0) then begin inten(r+n_points) = image(fix(x),fix(y)) endif $ else off_screen=1 endfor if min(inten) EQ max(inten) then off_screen=1 ; Watch for saturated areas ; Now we have the intensity along the radius, so compute the point with ; the maximum derivative and store its position in the xfit and yfit ; vectors. if NOT off_screen then begin rinten = radius-n_points+findgen(n_elements(inten)) ; vector of radii if NOT keyword_set(fast) then begin ; Compute the *smoothed* derivative deriv_lud,rinten,inten,derivative,mid=rmidpoint,reg_param=1.0,/quiet endif $ else begin ; Compute the *unsmoothed* derivative derivative = shift(inten,1)-inten derivative(n_elements(derivative)-1) = 0.0 & derivative(0)=0.0 rmidpoint = 0.5*(shift(rinten,1)+rinten) rmidpoint(n_elements(rmidpoint)-1)=rinten(n_elements(rinten)-1) rmidpoint(0)=rinten(0) endelse test_function = abs(derivative) dmax = max(test_function,inflection_p) ; Find maximum derivative ; Now interpolate to the best value of test_function max ; put inflection_p in range inflection_p = ((inflection_p > 3) < (n_elements(rmidpoint)-4)) contrast = abs(inten(inflection_p+3) - inten(inflection_p-3)) ; line_cent does a parabolic fit to test_function to find the max inflection_p = line_cent(abs(test_function),best=inflection_p,n_p=3) ; put inflection_p in range inflection_p = ((inflection_p > 0) < (n_elements(rmidpoint)-1)) ; Linear interpolation to get the limb radius inflection_p = interpol(rmidpoint,findgen(n_elements(rmidpoint)), $ inflection_p) inflection_p = inflection_p(0) ; convert 1-element array to scalar if keyword_set(verbose) then begin plot,rinten,inten oplot,rmidpoint,test_function,lines=2 oplot,[inflection_p,inflection_p],[-1.e6,1.e6] endif ; If the maximum derivative is at an endpoint, ignore it, ; otherwise add the point to the xfit and yfit vectors. ; Also test the significance of the result, ignore if it's not ; over 3 sigma or if the contrast over the limb is less than a ; threshold value if (inflection_p GT rmidpoint(2)) AND $ (inflection_p LT rmidpoint(n_elements(rmidpoint)-3)) $ AND (dmax GE 2.5*stdev(test_function)) AND $ (contrast GT MAX_VALUE*0.05) AND (contrast LE MAX_VALUE*0.84) $ AND (max(inten) LT MAX_VALUE*0.86) $ then begin xf = xcen + inflection_p*cos(theta(t)) yf = ycen + inflection_p*sin(theta(t)) if n_elements(xfit) LE 0 then xfit = [xf] $ else xfit=[xfit,xf] if n_elements(yfit) LE 0 then yfit = [yf] $ else yfit=[yfit,yf] endif $ else if keyword_set(verbose) then $ xyouts,0.75*(!X.CRANGE(1)-!X.CRANGE(0))+!X.CRANGE(0), $ 0.9*(!Y.CRANGE(1)-!Y.CRANGE(0))+!Y.CRANGE(0),align=0.5, $ 'REJECTED',/data endif endfor if n_elements(xfit) LT 3 then begin print,"ERROR: fit_limb: Could not get enough limb points" x0 = 0. y0 = 0. r0 = 0. r_err = 0. ob_angle = 0. oblateness = 0. bias = 0. endif $ else begin ; Make a cut on the distance from the mean radius. This assumes that ; xcen and ycen are good guesses to the center of the disk! ; fit a circle to the limb points txfit = xfit ; temporary variables tyfit = yfit txcen = xcen tycen = ycen tradius = radius ; Number of iterations killing bad points if n_elements(RLOOP) NE 1 then begin if (isxray) then RLOOP=2 else RLOOP=0 endif $ else RLOOP=fix(abs(RLOOP)) for loop=0,RLOOP-1 do begin if n_elements(xfit) GE 3 then begin if keyword_set(ellipse) then begin circle_fit,txfit,tyfit, txcen, tycen, tradius, xcrude, ycrude, $ rcrude,check,/quiet,/noiterate,norfit=norfit endif else begin ; tolerance of 0.0 means do not iterate fit = fit_circle(txfit,tyfit,[txcen,tycen,tradius], $ tolerance=0.,radius_fix=norfit, $ limit_iter=CIRCLE_ITER,num_iter=num_iter) check = num_iter NE CIRCLE_ITER xcrude=fit(0) & ycrude=fit(1) & rcrude=fit(2) endelse if check then begin rfit = sqrt((xfit-xcrude)^2+(yfit-ycrude)^2) trfit = sqrt((txfit-xcrude)^2+(tyfit-ycrude)^2) endif $ else begin rfit = sqrt((xfit-txcen)^2+(yfit-tycen)^2) trfit = sqrt((txfit-txcen)^2+(tyfit-tycen)^2) rcrude = total(rfit)/n_elements(rfit) endelse r_cutoff = (stdev(trfit)>0.25) good = where( abs(rfit-rcrude) LT r_cutoff,count) if count GT 0 then begin txfit = xfit(good) tyfit = yfit(good) endif txcen = xcrude tycen = ycrude tradius = rcrude endif endfor xfit = txfit yfit = tyfit xcen = txcen ycen = tycen radius = tradius if n_elements(xfit) LT 3 then begin print,"ERROR: fit_limb: Could not get enough limb points" x0 = 0. y0 = 0. r0 = 0. r_err = 0. ob_angle = 0. oblateness = 0. bias = 0. endif $ else begin if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then begin erase tv,bytscl(image,top=!D.N_COLORS-20) xyouts,.15,.95,"Limb Points",align=0.5,/norm,charsize=1.5 plots,xfit,yfit,psym=1,/device ; Show the limb points to the user wait,2. endif if keyword_set(interactive) then begin print print print,"Click left button near obviously bad limb points to erase them" print,"Click right button when finished to continue with the fit" REPEAT begin cursor,x,y,3,/device ; Give the user a chance to edit the points if !ERR LT 4 then begin distance = (xfit-x)^2+(yfit-y)^2 junk = min(distance,closest_p) ; find the closest xfit,yfit point if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then $ plots,xfit(closest_p),yfit(closest_p),/device,psym=1,color=0 if closest_p NE 0 and closest_p NE n_elements(xfit)-1 then begin xfit = $ [xfit(0:closest_p-1),xfit(closest_p+1:n_elements(xfit)-1)] yfit = $ [yfit(0:closest_p-1),yfit(closest_p+1:n_elements(yfit)-1)] endif $ else if closest_p EQ 0 then begin xfit = xfit(1:n_elements(xfit)-1) yfit = yfit(1:n_elements(yfit)-1) endif $ else if closest_p EQ n_elements(xfit)-1 then begin xfit = xfit(0:n_elements(xfit)-2) yfit = yfit(0:n_elements(yfit)-2) endif endif endrep UNTIL !ERR GE 4 endif if keyword_set(interactive) then begin erase if NOT keyword_set(quiet) then begin print print,"Working ..." print endif endif ; fit a circle to the limb points !P.MULTI = save_P_multi if keyword_set(ellipse) then begin circle_fit, xfit, yfit, xcen, ycen, radius, x0, y0, r0, check, $ noiterate=noiterate,norfit=norfit endif else begin if keyword_set(noiterate) then tol = 0.0 else tol = 0.001 fit = fit_circle(xfit,yfit,[xcen,ycen,radius],radius_fix=norfit, $ tolerance=tol,limit_iter=CIRCLE_ITER, $ num_iter=num_iter) check = (num_iter NE CIRCLE_ITER) x0 = fit(0) y0 = fit(1) r0 = fit(2) endelse if NOT check then begin print,"ERROR: fit_limb: No convergence in circle fit" x0 = 0. y0 = 0. r0 = 0. r_err = 0. ob_angle = 0. oblateness = 0. bias = 0. endif $ else begin if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then begin erase tv,bytscl(image) xyouts,.15,.95,"Best Fit",align=0.5,/norm,charsize=1.5 endif x = x0 + r0 * cos(theta) y = y0 + r0 * sin(theta) if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then $ plots,x,y,/device ; display the best fit circle on the image ; Correct for resolution pixel differences (stolen from find_limb) ; Also correct for oblateness of the disk ; ; Analyze limb shape by fitting harmonics (1,2,3,4 * theta). ; ; Stolen from find_limb rad = sqrt((xfit-x0)^2 +(yfit-y0)^2) - r0 theta = atan((yfit-y0),(xfit-x0)) sin_fun = sin(theta) cos_fun = cos(theta) sin_1 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a) cos_1 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a) r_err = sigma/sqrt(n_elements(rad)) ;print, sigma, ' Residuals after fundamental' rad = rad - (sin_1(0) + cos_1(0))/2 - sin_1(1)*sin(theta)- cos_1(1)*cos(theta) sin_fun = sin(2*theta) cos_fun = cos(2*theta) sin_2 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a) cos_2 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a) ;print, sigma, ' Residuals after 1st harmonic' rad = rad-(sin_2(0)+cos_2(0))/2 - sin_2(1)*sin(2*theta) - cos_2(1)*cos(2*theta) sin_fun = sin(3*theta) cos_fun = cos(3*theta) sin_3 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a) cos_3 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a) ;print, sigma, ' Residuals after 2nd harmonic' rad = rad-(sin_3(0)+cos_3(0))/2 - sin_3(1)*sin(3*theta) - cos_3(1)*cos(3*theta) sin_fun = sin(4*theta) cos_fun = cos(4*theta) sin_4 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a) cos_4 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a) ;print, sigma, ' Residuals after 3rd harmonic' rad = rad-(sin_4(0)+cos_4(0))/2 - sin_4(1)*sin(4*theta) - cos_4(1)*cos(4*theta) ob_angle = atan(sin_2(1)/cos_2(1)) oblateness = sqrt(sin_2(1)^2 + cos_2(1)^2) bias = sqrt(sin_3(1)^2 + cos_3(1)^2 + sin_4(1)^2 + cos_4(1)^2) ; No oblateness correction for X-rays if (NOT keyword_set(noobcorrect)) and (NOT isxray) then begin x0 = x0 + sin_1(1) y0 = y0 + cos_1(1) endif if (not keyword_set(summed_pixels) and qindex) then begin off_corr = [-1, 0, 1, -1, 3] ;MDM added 19-Nov-91 ;changing from 1x1 to 1x1,no offset ;changing from 2x2 to 1x1,offset=1 (1x1 pixels 1&2 are 2x2 pixel 0) ;changing from 4x4 to 1x1,offset=3 (1x1 pixels 3,4,5&6 are 4x4 pixel 0) x0 = x0 - 0.5 ; Move the image such that pixel center is defined y0 = y0 - 0.5 ; as the pixel ; -- Output in FULL RESOLUTION pixels x0 = x0 * resolution + off_corr(resolution) y0 = y0 * resolution r0 = r0 * resolution end else begin ; MDM added the ELSE statement 20-Mar-95 x0 = x0 - 0.5 ; Move the image such that pixel center is defined y0 = y0 - 0.5 ; as the pixel end print,x0,y0,r0,init_radius/r0, $ format="('Center: x=',f10.4,' y=',f10.4,' Radius: ',f10.4,' RSun/RFit: ',f10.4)" endelse endelse endelse print !P.MULTI = save_P_multi end