;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009 ;function find_centroid_sat_xrg ; ;PURPOSE: ; To find a centroid of a small rectangular region that is known to ; contain a star. This program is to be used in conjunction with the ; script examine_persistent_pixels.pro. find_centroid.pro should be used ; to simply return a centroid. ; ;INPUTS: ; Image - The 2d image containing the star ; ;OUTPUTS: A 4 element array ; 1 - X Centroid Coord ; 2 - Y Centroid Coord ; 3 - Second Moment in X ; 4 - Second Moment in Y ; ; SatRad = Radius at which the star becomes unsaturated ;KEYWORDS: ; TVFLAG ; 0 - Don't plot or Tv ; 1 - Tv the star and print out the centroids as iterated ; ; SATURATIONLIM: ; Set this keyword to the saturation value for the star. ; When the mean number of counts per pixel falls below this ; for a given aperture radius, this radius will be a good estimate ; for the radius that contains the saturated pixels. ; IsStar: ; 0 if the object is not a star. 1 if it is a star ; ElecStr: ; 'ASIC' or 'LEACH' ; RefPix: ; 4 for HxRG multiplexers ; function find_centroid_sat_xRG, Image, XStart, XStop, YStart, YStop, IsStar, $ SatRad, TvFlag = TvFlag, $ SaturationLim = SaturationLim, Verbose = Verbose, $ ElecStr = ElecStr, RefPix=RefPix If N_Elements(TvFlag) eq 0 then TvFlag = 0 If N_Elements(SaturationLim) eq 0 then SaturationLim = 1 If N_Elements(Verbose) eq 0 then Verbose = 0 If N_Elements(ElecStr) eq 0 then ElecStr = 'ASIC' If N_Elements(RefPix) eq 0 then RefPix = 4 SatRad = 0 IsStar = 1 Naxis1 = (Size(Image))[1] & Naxis2 = (Size(Image))[2] X = 0D & Y = 0D XStartL = XStart & XStopL = XStop YStartL = YStart & YStopL = YStop LoopNum = 0 Repeat Begin LoopNum+=1 SubImage = Image(XStartL:XStopL, YStartL: YStopL) ;Make the arrays with negative and positive indices for easy weighting SubImSize = Size(SubImage) XCoords = Indgen(SubImSize(1),SubImSize(2)) mod SubImSize(1) XCoords = XCoords - Mean(XCoords) YCoords = reverse(Indgen(SubImSize(1),SubImSize(2)),2) / SubImSize(1) YCoords = YCoords - Mean(YCoords) ;Calculate Flux, First and Second Moments Flux = Total(SubImage) X = Total(XCoords*SubImage)/Flux Y = Total(YCoords*SubImage)/Flux X2 = Total(XCoords^2*SubImage)/Flux Y2 = Total(YCoords^2*SubImage)/Flux if ElecStr eq 'ASIC' then begin ;Shift the image over and try to recentroid XStartL = XStartL + round(X) > 0 XStopL = XStopL + round(X) < (Naxis1-2*RefPix) YStartL = YStartL - round(Y) > 0 YStopL = YStopL - round(Y) < (Naxis2-2*RefPix) Endif Else if ElecStr eq 'LEACH' then begin XStartL = XStartL + round(X) > 0 XStopL = XStopL + round(X) < (Naxis1-2*RefPix-1) YStartL = YStartL - round(Y) > 0 YStopL = YStopL - round(Y) < (Naxis2-2*RefPix) EndiF ;Set TvFlag = 1 to check if centroiding is being done correctly If TvFlag eq 1 then begin If TvFlag eq 1 then window,xsize=768,ysize=768 tvimage, congrid(bytscl(SubImage),512,512), $ background= fsc_color('white'),/erase print, 'X Centroid: '+Strtrim(X,2) print, 'Y Centroid: '+Strtrim(Y,2) print, 'X Sec Mom : '+Strtrim(sqrt(X2),2) print, 'Y Sec Mom : '+Strtrim(sqrt(Y2),2) EndIf EndRep Until Abs(X) lt 1 and Abs(Y) lt 1 or (LoopNum gt 40) ;Now find out where the outer edge of the star/saturated mask is ;by checking flux in different apertures RCoords = XCoords^2 +YCoords^2 MaxRad = fix(Max(Sqrt(RCoords))) FluxIn = Fltarr(MaxRad+1) & FluxOut = Fltarr(MaxRad+1) ;Keep track of mean count number In, Out, and in 1 pixel Annulus MeanCountOut = Fltarr(MaxRad+1) MeanCountAnn = Fltarr(MaxRad+1) MeanCountIn = Fltarr(MaxRad+1) For Rad = 1, fix(MaxRad) do begin InRadPix = where(sqrt(RCoords) le Rad, complement = OutRadPix) AnnRadPix = where(sqrt(RCoords) ge Rad $ and sqrt(RCoords) lt Rad+1) FluxIn(Rad) = total(SubImage(InRadPix)) FluxOut(Rad) = total(SubImage(OutRadPix)) MeanCountIn(Rad) = Mean(SubImage(InRadPix)) MeanCountOut(Rad) = Mean(SubImage(OutRadPix)) MeanCountAnn(Rad) = Mean(SubImage(AnnRadPix)) If SatRad eq 0 and MeanCountAnn(Rad) lt SaturationLim then $ SatRad = Rad EndFor ;Determine whether this is really a star If SatRad eq 0 or Flux lt 4*Max(SubImage) or $ X2 lt 1 or Y2 lt 1 then IsStar = 0 ;The centroid might be different if on the edge XCen = XStartL + Round(Float(XStopL-XStartL)/2) YCen = YStartL + Round(Float(YStopL-YStartL)/2) If Verbose then begin If IsStar then print, 'This is a star' If N_Elements(SatRad) gt 0 then print, 'Saturation Radius :'+Strtrim(SatRad,2) If IsStar then print, 'XCen : ', XCen, 'YCen : ', YCen EndIf If TvFlag then stop return, [ XCen, YCen, X2, Y2] end