5454- The code uses a new asymmetry weighting function developed with help from Connie Rockosi
5555- This code is adapted from the SDSS centroiding code, which was written by Jim Gunn
5656 and cleaned up by Connie Rockosi.
57-
57+
5858History:
59592004-03-22 ROwen First release.
60602004-04-07 ROwen Packaged as part of PyGuide and moved test code elsewhere.
@@ -133,7 +133,7 @@ def _fmtList(alist):
133133
134134class CentroidData :
135135 """Centroid data, including the following fields:
136-
136+
137137 flags; check before paying attention to the remaining data:
138138 - isOK if False then centroiding failed; see msgStr for more info
139139 - msgStr warning or error message, if any
@@ -142,7 +142,7 @@ class CentroidData:
142142 basic info:
143143 - rad radius for centroid search (pix)
144144 - imStats image statistics such as median and std. dev. (if known); an ImUtil.ImStats object.
145-
145+
146146 star data:
147147 - xyCtr the x,y centroid (pixels); None if unspecified
148148 - xyErr the predicted 1-sigma uncertainty in xyCtr (pixels); [nan, nan] if unspecified
@@ -157,12 +157,12 @@ class CentroidData:
157157 pixNoise(rad) = sqrt((readNoise/ccdGain)^2 + (meanVal(rad)-bias)/ccdGain)
158158 - pix the total number of unmasked pixels (ADU)
159159 - counts the total number of counts (ADU)
160-
160+
161161 Warning: asymm is supposed to be normalized, but it gets large
162162 for bright objects with lots of masked pixels. This may be
163163 simply because the value is only computed at the nearest integer pixel
164164 or because the noise is assumed gaussian, or some error.
165-
165+
166166 Suggested use:
167167 - check isOK; if False do not use the data
168168 - check nSat(); if not None and more than a few then be cautious in using the data
@@ -183,7 +183,7 @@ def __init__(self,
183183 self .isOK = isOK
184184 self .msgStr = msgStr
185185 self .nSat = nSat
186-
186+
187187 self .rad = rad
188188 if imStats is None :
189189 imStats = ImUtil .ImStats ()
@@ -193,11 +193,11 @@ def __init__(self,
193193 if xyErr is None :
194194 xyErr = numpy .array ([numpy .nan , numpy .nan ])
195195 self .xyErr = xyErr
196-
196+
197197 self .asymm = asymm
198198 self .pix = pix
199199 self .counts = counts
200-
200+
201201 def __repr__ (self ):
202202 dataList = []
203203 for arg in ("isOK" , "msgStr" , "nSat" , "xyCtr" , "xyErr" , "asymm" , "pix" , "counts" , "rad" , "imStats" ):
@@ -231,10 +231,10 @@ def basicCentroid(
231231 3: print basic iteration info, 4: print detailed iteration info.
232232 Note: there are no warnings at this time because the relevant info is returned.
233233 - doDS9 if True, diagnostic images are displayed in ds9
234-
234+
235235 Masks are optional. If specified, they must be the same shape as "data"
236236 and should be of type Bool. None means no mask (all data is OK).
237-
237+
238238 Returns a CentroidData object (which see for more info),
239239 but with no imStats info.
240240 """
@@ -272,21 +272,21 @@ def basicCentroid(
272272 # display circle showing the centroider input in frame 1
273273 args = list (xyGuess ) + [rad ]
274274 ds9Win .xpaset ("regions" , "image; circle %s # group=ctrcirc" % _fmtList (args ))
275-
275+
276276 try :
277277 # OK, use this as first guess at maximum. Extract radial profiles in
278278 # a 3x3 gridlet about this, and walk to find minimum fitting error
279279 maxi , maxj = ijIndGuess
280280 asymmArr = numpy .zeros ([3 ,3 ], float )
281281 totPtsArr = numpy .zeros ([3 ,3 ], int )
282282 totCountsArr = numpy .zeros ([3 ,3 ], float )
283-
283+
284284 niter = 0
285285 while True :
286286 niter += 1
287287 if niter > _MaxIter :
288288 raise RuntimeError ("could not find a star in %s iterations" % (niter ,))
289-
289+
290290 for i in range (3 ):
291291 ii = maxi + i - 1
292292 for j in range (3 ):
@@ -299,42 +299,42 @@ def basicCentroid(
299299# (warning: the error estimate will be invalid and chiSq will not be normalized)
300300# asymmArr[i, j], totCountsArr[i, j], totPtsArr[i, j] = radProf.radAsymm(
301301# data, mask, (ii, jj), rad)
302-
302+
303303 if verbosity > 3 :
304304 print ("basicCentroid: ind=[%s, %s] ctr=(%s, %s) asymm=%10.1f, totPts=%s, totCounts=%s" % \
305305 (i , j , ii , jj , asymmArr [i , j ], totPtsArr [i , j ], totCountsArr [i , j ]))
306-
306+
307307 # have error matrix. Find minimum
308308 ii , jj = scipy .ndimage .minimum_position (asymmArr )
309309 ii -= 1
310310 jj -= 1
311-
311+
312312 if verbosity > 2 :
313313 print ("basicCentroid: error matrix min ii=%d, jj=%d, errmin=%5.1f" % (ii , jj , asymmArr [ii ,jj ]))
314314 if verbosity > 3 :
315315 print ("basicCentroid: asymm matrix =\n " , asymmArr )
316-
316+
317317 if (ii != 0 or jj != 0 ):
318318 # minimum error not in center; walk and try again
319319 maxi += ii
320320 maxj += jj
321321 if verbosity > 2 :
322322 print ("shift by" , - ii , - jj , "to" , maxi , maxj )
323-
323+
324324 if ((maxi - ijIndGuess [0 ])** 2 + (maxj - ijIndGuess [1 ])** 2 ) >= rad ** 2 :
325325 raise RuntimeError ("could not find star within %r pixels" % (rad ,))
326-
326+
327327 # shift asymmArr and totPtsArr to minimum is in center again
328328 asymmArr = scipy .ndimage .shift (asymmArr , (- ii , - jj ))
329329 totCountsArr = scipy .ndimage .shift (totCountsArr , (- ii , - jj ))
330330 totPtsArr = scipy .ndimage .shift (totPtsArr , (- ii , - jj ))
331331 else :
332332 # Have minimum. Get out and go home.
333333 break
334-
334+
335335 if verbosity > 2 :
336336 print ("basicCentroid: found ijMax=%s after %r iterations" % ((maxi , maxj ), niter ,))
337-
337+
338338 # perform a parabolic fit to find true centroid
339339 # and compute the error estimate
340340 # y(x) = ymin + a(x-xmin)^2
@@ -346,32 +346,32 @@ def basicCentroid(
346346 bi = 0.5 * (asymmArr [2 , 1 ] - asymmArr [0 , 1 ])
347347 aj = 0.5 * (asymmArr [1 , 2 ] - 2.0 * asymmArr [1 , 1 ] + asymmArr [1 , 0 ])
348348 bj = 0.5 * (asymmArr [1 , 2 ] - asymmArr [1 , 0 ])
349-
349+
350350 di = - 0.5 * bi / ai
351351 dj = - 0.5 * bj / aj
352352 ijCtr = (
353353 maxi + di ,
354354 maxj + dj ,
355355 )
356356 xyCtr = ImUtil .xyPosFromIJPos (ijCtr )
357-
357+
358358 if verbosity > 2 :
359359 print ("basicCentroid: asymmArr[0:3, 1]=%s, ai=%s, bi=%s, di=%s, iCtr=%s" % (asymmArr [0 :3 , 1 ], ai , bi , di , ijCtr [0 ]))
360360 print ("basicCentroid: asymmArr[1, 0:3]=%s, aj=%s, bj=%s, dj=%s, jCtr=%s" % (asymmArr [1 , 0 :3 ], aj , bj , dj , ijCtr [1 ]))
361-
361+
362362 # crude error estimate, based on measured asymmetry
363363 # note: I also tried using the minimum along i,j but that sometimes is negative
364364 # and this is already so crude that it's not likely to help
365365 radAsymmSigma = asymmArr [1 ,1 ]
366366 iErr = math .sqrt (radAsymmSigma / ai )
367367 jErr = math .sqrt (radAsymmSigma / aj )
368368 xyErr = (jErr , iErr )
369-
369+
370370 if ds9Win :
371371 # display x at centroid
372372 ds9Win .xpaset ("regions" , "image; x point %s # group=centroid" % \
373373 _fmtList (xyCtr ))
374-
374+
375375
376376 # count # saturated pixels, if satMask available
377377 if satMask is None :
@@ -400,10 +400,10 @@ def makeDisk(i, j):
400400 )
401401 subMask = subMaskObj .getSubFrame ()
402402 numpy .logical_and (maybeSatPixel , numpy .logical_not (subMask ), maybeSatPixel )
403-
403+
404404 numpy .logical_and (subSatMask , maybeSatPixel , maybeSatPixel )
405405 nSat = scipy .ndimage .sum (maybeSatPixel )
406-
406+
407407 ctrData = CentroidData (
408408 isOK = True ,
409409 rad = rad ,
@@ -444,7 +444,7 @@ def centroid(
444444 checkSig = (True , True ),
445445):
446446 """Centroid and then confirm that there is usable signal at the location.
447-
447+
448448 Inputs:
449449 - data image data [i,j]
450450 - mask a mask of invalid data (1 if invalid, 0 if valid); None if no mask.
@@ -462,14 +462,14 @@ def centroid(
462462 - doDS9 if True, display diagnostic images in ds9
463463 - checkSig Verify usable signal for circle at (xyGuess, xy centroid)?
464464 If both are false then imStats is not computed.
465-
465+
466466 Returns a CentroidData object (which see for more info).
467467 """
468468 if verbosity > 2 :
469469 print ("data =" , data )
470470 print ("mask =" , mask )
471471 print ("centroid(xyGuess=%s, rad=%s, ccdInfo=%s, thresh=%s)" % (xyGuess , rad , ccdInfo , thresh ))
472-
472+
473473 if checkSig [0 ]:
474474 signalOK , imStats = checkSignal (
475475 data = data ,
@@ -524,7 +524,7 @@ def checkSignal(
524524 verbosity = 0 ,
525525):
526526 """Check that there is usable signal in a given circle.
527-
527+
528528 Inputs:
529529 - data image data [i,j]
530530 - mask a mask [i,j] of 0's (valid data) or 1's (invalid); None if no mask.
@@ -537,11 +537,11 @@ def checkSignal(
537537 values less than PyGuide.Constants.MinThresh are silently increased
538538 - doSmooth if True apply a 3x3 median filter to smooth the data
539539 - verbosity 0: no output, 1: print warnings, 2: print information, 3: print iteration info.
540-
540+
541541 Return:
542542 - signalOK True if usable signal is present
543543 - imStats background statistics: a PyGuide.ImStats object
544-
544+
545545 Details of usable signal:
546546 - Computes median and stdDev in a region extending from a circle of radius "rad"
547547 to a box of size (rad+_OuterRadAdd)*2 on a side.
@@ -572,7 +572,7 @@ def checkSignal(
572572 nPts = subData .size ,
573573 )
574574 subCtrIJ = subDataObj .subIJFromFullIJ (ImUtil .ijPosFromXYPos (xyCtr ))
575-
575+
576576 if mask is not None :
577577 subMaskObj = ImUtil .subFrameCtr (
578578 mask ,
@@ -588,7 +588,7 @@ def checkSignal(
588588 def makeCircle (i , j ):
589589 return ((i - subCtrIJ [0 ])** 2 + (j - subCtrIJ [1 ])** 2 ) > rad ** 2
590590 circleMask = numpy .fromfunction (makeCircle , subData .shape )
591-
591+
592592 # make a copy of the data outside a circle of radius "rad";
593593 # use this to compute background stats
594594 bkgndPixels = numpy .extract (numpy .logical_and (circleMask , numpy .logical_not (subMask )), subData )
@@ -604,7 +604,7 @@ def makeCircle(i, j):
604604 return False , ImUtil .ImStats (
605605 nPts = bkgndPixels .size ,
606606 )
607-
607+
608608 imStats = ImUtil .skyStats (bkgndPixels , thresh )
609609 del (bkgndPixels )
610610
@@ -618,7 +618,7 @@ def makeCircle(i, j):
618618 if doSmooth :
619619 scipy .ndimage .median_filter (smoothedData , 3 , output = smoothedData )
620620 del (dataPixels )
621-
621+
622622 # look for a blob of at least 2x2 adjacent pixels with smoothed value >= dataCut
623623 # note: it'd be much simpler but less safe to simply test:
624624 # if max(smoothedData) < dataCut: # have signal
@@ -641,15 +641,15 @@ def makeCircle(i, j):
641641def conditionData (data ):
642642 """Convert dataArr to the correct type
643643 such that basicCentroid can operate most efficiently on it.
644-
644+
645645 Warning: does not copy the data unless necessary.
646646 """
647647 return conditionArr (data , desType = numpy .float32 )
648648
649649def conditionMask (mask ):
650650 """Convert mask to the correct type
651651 such that basicCentroid can operate most efficiently on it.
652-
652+
653653 Mask is optional, so a value of None returns None.
654654
655655 Warning: does not copy the data unless necessary.
0 commit comments