pro fit_xlimb, image, x0, y0, r0, r_err, $ index=index, $ summed_pixels=summed_pixels, $ pixel_size=pixel_size, $ verbose=verbose, limb_points=limb_points, $ xlimb=xlimb,$ ylimb=ylimb, 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, $ resolution=resolution,n_points=n_points,no_ask=no_ask ;+ ;NAME: ; fit_xlimb ;PURPOSE: ; Find the center and radius of an image. NOTE: For SXT images, the ; default is to return the result in full resolution pixels. ;CALLING SEQUENCE: ; fit_xlimb, image, x, y, r, pixel_size=2.5, date='05-Oct-93' ; fit_xlimb, image, x, y, r, pixel_size=2.5, intial=[434,512,403*2.5] ; fit_xlimb, data(*,*,nimage), x, y, r, index=index(nimage),/norfit ;OUTPUT ; x0, y0, image center in FR pixels. ;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. ; n_points Number of points (+/-) along radius to check for inflection. ; Default is 15 for HR and (if index input) n_points/2 for QR. ; 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 ; 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 ; Suppress final print to screen of fit result. ; /limb_points = if set, show limb points, even if /quiet ; /no_ask = go ahead and fit, don't ask first. ; initial_values = 3-element vector giving initial guess: [xc,yc,rsun] ; xc and yc are in pixels, rsun 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. ; NOTE: If xfit ERRORs off then x=y=-1.0 are returned. ; 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 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 FIT_LIMB.PRO ; 17-May-05 (LWA) - If /quiet don't print final result to screen. ; 2-Aug-05 (LWA) - Increased n_points from 10 to 25. ; Added keyword limb_points. ; 3-Aug-05 (LWA) - Stripped unused variables and changed name to fit_xlimb.pro ; 7-Aug-05 (LWA) - Return x=y=-1.0 on fitting error. ; 15-Aug-05 (LWA) - Do test on (txfit GE 3) in RLOOP. ; 2-Nov-05 (LWA) - Increased n_points from 25 to 35 ; 8-Nov-05 (LWA) - If /interactive then n_points=5. ; 12-Dec-05 (LWA) - n_points=5, experimental. ; 17-Jan-05 (LWA) - Clarified header. ; Fixed xlimb output bug. ; 6-Feb-06 (LWA) - Set n_points=15. ; 8-Oct-07 (LWA) - Made n_points a keyword. ; If qindex and QR then n_points=npoints/2. ; Use sxt_shift_res to obtain FR, pix defined at LL corner. ; 28-Jul-08 (LWA) - Clarified header and internal comments. ; 8-May-13 (LWA) - Added /quiet in calls to get_rb0p. ;- ;on_error,2 CIRCLE_ITER = 50 ; iteration limit for circle fitting x0 = 0. y0 = 0. r0 = 0. r_err = 0. verbose = keyword_set(verbose) no_ask=keyword_set(no_ask) if keyword_set(noiterate) then noiterate=1 else noiterate=0 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,/quiet) resolution = 2^gt_res(index) corner = gt_corner(index, /from_sc) xcen = -corner(0) ycen = -corner(1) ; ---- Output parameters 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,/quiet) $ 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)) ;print,'!ERR GE 4' ;stop if not no_ask then begin yesnox,' Do you want to fit_limb?',in,'yes' endif else begin in=1 endelse if not in then begin r0=radius x0=xcen y0=ycen goto,noxfit endif 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. ; if keyword_set(interactive) then n_points=5 else $ ; n_points = 35 ;LWA 11/8/05 if NOT keyword_set(n_points) then n_points=15 else n_points=n_points ;LWA 10/8/07 if qindex then begin if gt_res(index) eq 2 then n_points=n_points/2 ;LWA 10/8/07 endif ; 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_xlimb: Could not get enough limb points" x0 = -1. y0 = -1. r0 = 0. r_err = 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(txfit) GE 3 then 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) 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_xlimb: Could not get enough limb points" x0 = -1. y0 = -1. r0 = 0. r_err = 0. endif else begin if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) $ OR (keyword_set(limb_points)) 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 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(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) if NOT check then begin print,"ERROR: fit_xlimb: No convergence in circle fit" x0 = -1. y0 = -1. r0 = 0. r_err = 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) ; 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) if (not keyword_set(summed_pixels) and qindex) then begin if qindex then begin case 1 of gt_res(index) eq 2 : xy=sxt_shift_res([x0,y0],/xy,/qf) else : xy=sxt_shift_res([x0,y0],/xy,/hf) endcase x0=xy(0)-0.5 & y0=xy(1)-0.5 ;LWA 10/9/07, shift for pix_center coord. endif endif ; 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 if NOT keyword_set(quiet) then 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 !P.MULTI = save_P_multi xlimb=xfit ylimb=yfit noxfit: end