LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
secondMomentStarSelector.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 import collections
23 import math
24 
25 import numpy
26 
27 import lsst.pex.config as pexConfig
28 import lsst.afw.detection as afwDetection
29 import lsst.afw.display.ds9 as ds9
30 import lsst.afw.image as afwImage
31 import lsst.afw.math as afwMath
32 import lsst.afw.table as afwTable
33 import lsst.afw.geom as afwGeom
34 import lsst.afw.geom.ellipses as geomEllip
35 import lsst.afw.cameraGeom as cameraGeom
36 from . import algorithmsLib
37 from lsst.meas.base import SingleFrameMeasurementTask, SingleFrameMeasurementConfig
38 
39 class SecondMomentStarSelectorConfig(pexConfig.Config):
40  fluxLim = pexConfig.Field(
41  doc = "specify the minimum psfFlux for good Psf Candidates",
42  dtype = float,
43  default = 12500.0,
44  check = lambda x: x >= 0.0,
45  )
46  fluxMax = pexConfig.Field(
47  doc = "specify the maximum psfFlux for good Psf Candidates (ignored if == 0)",
48  dtype = float,
49  default = 0.0,
50  check = lambda x: x >= 0.0,
51  )
52  clumpNSigma = pexConfig.Field(
53  doc = "candidate PSF's shapes must lie within this many sigma of the average shape",
54  dtype = float,
55  default = 2.0,
56  check = lambda x: x >= 0.0,
57  )
58  kernelSize = pexConfig.Field(
59  doc = "size of the kernel to create",
60  dtype = int,
61  default = 21,
62  )
63  borderWidth = pexConfig.Field(
64  doc = "number of pixels to ignore around the edge of PSF candidate postage stamps",
65  dtype = int,
66  default = 0,
67  )
68  badFlags = pexConfig.ListField(
69  doc = "List of flags which cause a source to be rejected as bad",
70  dtype = str,
71  default = ["base_PixelFlags_flag_edge",
72  "base_PixelFlags_flag_interpolatedCenter",
73  "base_PixelFlags_flag_saturatedCenter",
74  "base_PixelFlags_flag_crCenter",
75  ]
76  )
77  histSize = pexConfig.Field(
78  doc = "Number of bins in moment histogram",
79  dtype = int,
80  default = 64,
81  check = lambda x: x > 0,
82  )
83  histMomentMax = pexConfig.Field(
84  doc = "Maximum moment to consider",
85  dtype = float,
86  default = 100.0,
87  check = lambda x: x > 0,
88  )
89  histMomentMaxMultiplier = pexConfig.Field(
90  doc = "Multiplier of mean for maximum moments histogram range",
91  dtype = float,
92  default = 5.0,
93  check = lambda x: x > 0,
94  )
95  histMomentClip = pexConfig.Field(
96  doc = "Clipping threshold for moments histogram range",
97  dtype = float,
98  default = 5.0,
99  check = lambda x: x > 0,
100  )
101  histMomentMinMultiplier = pexConfig.Field(
102  doc = "Multiplier of mean for minimum moments histogram range",
103  dtype = float,
104  default = 2.0,
105  check = lambda x: x > 0,
106  )
107 
108 Clump = collections.namedtuple('Clump', ['peak', 'x', 'y', 'ixx', 'ixy', 'iyy', 'a', 'b', 'c'])
109 
110 class CheckSource(object):
111  """A functor to check whether a source has any flags set that should cause it to be labeled bad."""
112 
113  def __init__(self, table, badFlags, fluxLim, fluxMax):
114  self.keys = [table.getSchema().find(name).key for name in badFlags]
115  self.keys.append(table.getCentroidFlagKey())
116  self.fluxLim = fluxLim
117  self.fluxMax = fluxMax
118 
119  def __call__(self, source):
120  for k in self.keys:
121  if source.get(k):
122  return False
123  if self.fluxLim != None and source.getPsfFlux() < self.fluxLim: # ignore faint objects
124  return False
125  if self.fluxMax != 0.0 and source.getPsfFlux() > self.fluxMax: # ignore bright objects
126  return False
127  return True
128 
130  ConfigClass = SecondMomentStarSelectorConfig
131 
132  def __init__(self, config):
133  """Construct a star selector that uses second moments
134 
135  This is a naive algorithm and should be used with caution.
136 
137  @param[in] config: An instance of SecondMomentStarSelectorConfig
138  """
139  self.config = config
140 
141  def selectStars(self, exposure, catalog, matches=None):
142  """Return a list of PSF candidates that represent likely stars
143 
144  A list of PSF candidates may be used by a PSF fitter to construct a PSF.
145 
146  @param[in] exposure: the exposure containing the sources
147  @param[in] catalog: a SourceCatalog containing sources that may be stars
148  @param[in] matches: astrometric matches; ignored by this star selector
149 
150  @return psfCandidateList: a list of PSF candidates.
151  """
152  import lsstDebug
153  display = lsstDebug.Info(__name__).display
154 
155  displayExposure = lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
156 
157  isGoodSource = CheckSource(catalog.getTable(), self.config.badFlags, self.config.fluxLim,
158  self.config.fluxMax)
159 
160  detector = exposure.getDetector()
161 
162  mi = exposure.getMaskedImage()
163  #
164  # Create an Image of Ixx v. Iyy, i.e. a 2-D histogram
165  #
166 
167  # Use stats on our Ixx/yy values to determine the xMax/yMax range for clump image
168  iqqList = []
169  for s in catalog:
170  ixx, iyy = s.getIxx(), s.getIyy()
171  # ignore NaN and unrealistically large values
172  if (ixx == ixx and ixx < self.config.histMomentMax and
173  iyy == iyy and iyy < self.config.histMomentMax and
174  isGoodSource(s)):
175  iqqList.append(s.getIxx())
176  iqqList.append(s.getIyy())
177  stat = afwMath.makeStatistics(iqqList, afwMath.MEANCLIP | afwMath.STDEVCLIP | afwMath.MAX)
178  iqqMean = stat.getValue(afwMath.MEANCLIP)
179  iqqStd = stat.getValue(afwMath.STDEVCLIP)
180  iqqMax = stat.getValue(afwMath.MAX)
181 
182  iqqLimit = max(iqqMean + self.config.histMomentClip*iqqStd,
183  self.config.histMomentMaxMultiplier*iqqMean)
184  # if the max value is smaller than our range, use max as the limit, but don't go below N*mean
185  if iqqLimit > iqqMax:
186  iqqLimit = max(self.config.histMomentMinMultiplier*iqqMean, iqqMax)
187 
188  psfHist = _PsfShapeHistogram(detector=detector, xSize=self.config.histSize, ySize=self.config.histSize,
189  ixxMax=iqqLimit, iyyMax=iqqLimit)
190 
191  if display and displayExposure:
192  frame = 0
193  ds9.mtv(mi, frame=frame, title="PSF candidates")
194 
195  with ds9.Buffering():
196  for source in catalog:
197  if isGoodSource(source):
198  if psfHist.insert(source): # n.b. this call has the side effect of inserting
199  ctype = ds9.GREEN # good
200  else:
201  ctype = ds9.MAGENTA # rejected
202  else:
203  ctype = ds9.RED # bad
204 
205  if display and displayExposure:
206  ds9.dot("o", source.getX() - mi.getX0(),
207  source.getY() - mi.getY0(), frame=frame, ctype=ctype)
208 
209  clumps = psfHist.getClumps(display=display)
210 
211  #
212  # Go through and find all the PSF-like objects
213  #
214  # We'll split the image into a number of cells, each of which contributes only
215  # one PSF candidate star
216  #
217  psfCandidateList = []
218 
219  pixToTanXYTransform = None
220  if detector is not None:
221  tanSys = detector.makeCameraSys(cameraGeom.TAN_PIXELS)
222  pixToTanXYTransform = detector.getTransformMap().get(tanSys)
223 
224  # psf candidate shapes must lie within this many RMS of the average shape
225  # N.b. if Ixx == Iyy, Ixy = 0 the criterion is
226  # dx^2 + dy^2 < self.config.clumpNSigma*(Ixx + Iyy) == 2*self.config.clumpNSigma*Ixx
227  for source in catalog:
228  if not isGoodSource(source): continue
229  Ixx, Ixy, Iyy = source.getIxx(), source.getIxy(), source.getIyy()
230  if pixToTanXYTransform:
231  p = afwGeom.Point2D(source.getX(), source.getY())
232  linTransform = pixToTanXYTransform.linearizeForwardTransform(p).getLinear()
233  m = geomEllip.Quadrupole(Ixx, Iyy, Ixy)
234  m.transform(linTransform)
235  Ixx, Iyy, Ixy = m.getIxx(), m.getIyy(), m.getIxy()
236 
237  x, y = psfHist.momentsToPixel(Ixx, Iyy)
238  for clump in clumps:
239  dx, dy = (x - clump.x), (y - clump.y)
240 
241  if math.sqrt(clump.a*dx*dx + 2*clump.b*dx*dy + clump.c*dy*dy) < 2*self.config.clumpNSigma:
242  # A test for > would be confused by NaN
243  if not isGoodSource(source):
244  continue
245  try:
246  psfCandidate = algorithmsLib.makePsfCandidate(source, exposure)
247 
248  # The setXXX methods are class static, but it's convenient to call them on
249  # an instance as we don't know Exposure's pixel type
250  # (and hence psfCandidate's exact type)
251  if psfCandidate.getWidth() == 0:
252  psfCandidate.setBorderWidth(self.config.borderWidth)
253  psfCandidate.setWidth(self.config.kernelSize + 2*self.config.borderWidth)
254  psfCandidate.setHeight(self.config.kernelSize + 2*self.config.borderWidth)
255 
256  im = psfCandidate.getMaskedImage().getImage()
257  if not numpy.isfinite(afwMath.makeStatistics(im, afwMath.MAX).getValue()):
258  continue
259  psfCandidateList.append(psfCandidate)
260 
261  if display and displayExposure:
262  ds9.dot("o", source.getX() - mi.getX0(), source.getY() - mi.getY0(),
263  size=4, frame=frame, ctype=ds9.CYAN)
264  except Exception as err:
265  pass # FIXME: should log this!
266  break
267 
268  return psfCandidateList
269 
270 class _PsfShapeHistogram(object):
271  """A class to represent a histogram of (Ixx, Iyy)
272  """
273  def __init__(self, xSize=32, ySize=32, ixxMax=30, iyyMax=30, detector=None, xy0=afwGeom.Point2D(0,0)):
274  """Construct a _PsfShapeHistogram
275 
276  The maximum seeing FWHM that can be tolerated is [xy]Max/2.35 pixels.
277  The 'resolution' of stars vs galaxies/CRs is provided by [xy]Size/[xy]Max.
278  A larger (better) resolution may thresh the peaks, but a smaller (worse)
279  resolution will allow stars and galaxies/CRs to mix. The disadvantages of
280  a larger (better) resolution can be compensated (some) by using multiple
281  histogram peaks.
282 
283  @input[in] [xy]Size: the size of the psfImage (in pixels)
284  @input[in] ixxMax, iyyMax: the maximum values for I[xy][xy]
285  """
286  self._xSize, self._ySize = xSize, ySize
287  self._xMax, self._yMax = ixxMax, iyyMax
288  self._psfImage = afwImage.ImageF(afwGeom.ExtentI(xSize, ySize), 0)
289  self._num = 0
290  self.detector = detector
291  self.xy0 = xy0
292 
293  def getImage(self):
294  return self._psfImage
295 
296  def insert(self, source):
297  """Insert source into the histogram."""
298 
299  ixx, iyy, ixy = source.getIxx(), source.getIyy(), source.getIxy()
300  if self.detector:
301  tanSys = self.detector.makeCameraSys(cameraGeom.TAN_PIXELS)
302  if tanSys in self.detector.getTransformMap():
303  pixToTanXYTransform = self.detector.getTransformMap()[tanSys]
304  p = afwGeom.Point2D(source.getX(), source.getY())
305  linTransform = pixToTanXYTransform.linearizeForwardTransform(p).getLinear()
306  m = geomEllip.Quadrupole(ixx, iyy, ixy)
307  m.transform(linTransform)
308  ixx, iyy, ixy = m.getIxx(), m.getIyy(), m.getIxy()
309 
310  try:
311  pixel = self.momentsToPixel(ixx, iyy)
312  i = int(pixel[0])
313  j = int(pixel[1])
314  except:
315  return 0
316 
317  if i in range(0, self._xSize) and j in range(0, self._ySize):
318  if i != 0 or j != 0:
319  self._psfImage.set(i, j, self._psfImage.get(i, j) + 1)
320  self._num += 1
321  return 1 # success
322 
323  return 0 # failure
324 
325  def momentsToPixel(self, ixx, iyy):
326  #x = math.sqrt(ixx) * self._xSize / self._xMax
327  #y = math.sqrt(iyy) * self._ySize / self._yMax
328  x = ixx * self._xSize / self._xMax
329  y = iyy * self._ySize / self._yMax
330  return x, y
331 
332  def pixelToMoments(self, x, y):
333  """Given a peak position in self._psfImage, return the corresponding (Ixx, Iyy)"""
334 
335  #ixx = (x*self._xMax/self._xSize)**2
336  #iyy = (y*self._yMax/self._ySize)**2
337  ixx = x*self._xMax/self._xSize
338  iyy = y*self._yMax/self._ySize
339  return ixx, iyy
340 
341  def getClumps(self, sigma=1.0, display=False):
342  if self._num <= 0:
343  raise RuntimeError("No candidate PSF sources")
344 
345  psfImage = self.getImage()
346  #
347  # Embed psfImage into a larger image so we can smooth when measuring it
348  #
349  width, height = psfImage.getWidth(), psfImage.getHeight()
350  largeImg = psfImage.Factory(afwGeom.ExtentI(2*width, 2*height))
351  largeImg.set(0)
352 
353  bbox = afwGeom.BoxI(afwGeom.PointI(width, height), afwGeom.ExtentI(width, height))
354  subLargeImg = psfImage.Factory(largeImg, bbox, afwImage.LOCAL)
355  subLargeImg <<= psfImage
356  del subLargeImg
357  #
358  # Now measure that image, looking for the highest peak. Start by building an Exposure
359  #
360  msk = afwImage.MaskU(largeImg.getDimensions())
361  msk.set(0)
362  var = afwImage.ImageF(largeImg.getDimensions())
363  var.set(1)
364  mpsfImage = afwImage.MaskedImageF(largeImg, msk, var)
365  mpsfImage.setXY0(afwGeom.PointI(-width, -height))
366  del msk
367  del var
368  exposure = afwImage.makeExposure(mpsfImage)
369 
370  #
371  # Next run an object detector
372  #
373  maxVal = afwMath.makeStatistics(psfImage, afwMath.MAX).getValue()
374  threshold = maxVal - sigma*math.sqrt(maxVal)
375  if threshold <= 0.0:
376  threshold = maxVal
377 
378  threshold = afwDetection.Threshold(threshold)
379 
380  ds = afwDetection.FootprintSet(mpsfImage, threshold, "DETECTED")
381  #
382  # And measure it. This policy isn't the one we use to measure
383  # Sources, it's only used to characterize this PSF histogram
384  #
385  schema = afwTable.SourceTable.makeMinimalSchema()
386  psfImageConfig = SingleFrameMeasurementConfig()
387  psfImageConfig.doApplyApCorr = "no"
388  psfImageConfig.slots.centroid = "base_SdssCentroid"
389  psfImageConfig.slots.psfFlux = None #"base_PsfFlux"
390  psfImageConfig.slots.apFlux = "base_CircularApertureFlux_3_0"
391  psfImageConfig.slots.modelFlux = None
392  psfImageConfig.slots.instFlux = None
393  psfImageConfig.slots.calibFlux = None
394  psfImageConfig.slots.shape = "base_SdssShape"
395  # Formerly, this code had centroid.sdss, flux.psf, flux.naive,
396  # flags.pixel, and shape.sdss
397  psfImageConfig.algorithms.names = ["base_SdssCentroid", "base_CircularApertureFlux", "base_SdssShape"]
398  psfImageConfig.algorithms["base_CircularApertureFlux"].radii = [3.0]
399  psfImageConfig.validate()
400  task = SingleFrameMeasurementTask(schema, config=psfImageConfig)
401 
402  catalog = afwTable.SourceCatalog(schema)
403 
404  gaussianWidth = 1.5 # Gaussian sigma for detection convolution
405  exposure.setPsf(algorithmsLib.DoubleGaussianPsf(11, 11, gaussianWidth))
406 
407  ds.makeSources(catalog)
408  #
409  # Show us the Histogram
410  #
411  if display:
412  frame = 1
413  dispImage = mpsfImage.Factory(mpsfImage, afwGeom.BoxI(afwGeom.PointI(width, height),
414  afwGeom.ExtentI(width, height)),
415  afwImage.LOCAL)
416  ds9.mtv(dispImage,title="PSF Selection Image", frame=frame)
417 
418 
419  clumps = list() # List of clumps, to return
420  e = None # thrown exception
421  IzzMin = 1.0 # Minimum value for second moments
422  IzzMax = (self._xSize/8.0)**2 # Max value ... clump r < clumpImgSize/8
423  # diameter should be < 1/4 clumpImgSize
424  apFluxes = []
425  task.run(exposure, catalog) # notes that this is backwards for the new framework
426  for i, source in enumerate(catalog):
427  if source.getCentroidFlag():
428  continue
429  x, y = source.getX(), source.getY()
430 
431  apFluxes.append(source.getApFlux())
432 
433  val = mpsfImage.getImage().get(int(x) + width, int(y) + height)
434 
435  psfClumpIxx = source.getIxx()
436  psfClumpIxy = source.getIxy()
437  psfClumpIyy = source.getIyy()
438 
439  if display:
440  if i == 0:
441  ds9.pan(x, y, frame=frame)
442 
443  ds9.dot("+", x, y, ctype=ds9.YELLOW, frame=frame)
444  ds9.dot("@:%g,%g,%g" % (psfClumpIxx, psfClumpIxy, psfClumpIyy), x, y,
445  ctype=ds9.YELLOW, frame=frame)
446 
447  if psfClumpIxx < IzzMin or psfClumpIyy < IzzMin:
448  psfClumpIxx = max(psfClumpIxx, IzzMin)
449  psfClumpIyy = max(psfClumpIyy, IzzMin)
450  if display:
451  ds9.dot("@:%g,%g,%g" % (psfClumpIxx, psfClumpIxy, psfClumpIyy), x, y,
452  ctype=ds9.RED, frame=frame)
453 
454  det = psfClumpIxx*psfClumpIyy - psfClumpIxy*psfClumpIxy
455  try:
456  a, b, c = psfClumpIyy/det, -psfClumpIxy/det, psfClumpIxx/det
457  except ZeroDivisionError:
458  a, b, c = 1e4, 0, 1e4
459 
460  clumps.append(Clump(peak=val, x=x, y=y, a=a, b=b, c=c,
461  ixx=psfClumpIxx, ixy=psfClumpIxy, iyy=psfClumpIyy))
462 
463  if len(clumps) == 0:
464  msg = "Failed to determine center of PSF clump"
465  if e:
466  msg += ": %s" % e
467  raise RuntimeError(msg)
468 
469  # if it's all we got return it
470  if len(clumps) == 1:
471  return clumps
472 
473  # which clump is the best?
474  # if we've undistorted the moments, stars should only have 1 clump
475  # use the apFlux from the clump measurement, and take the highest
476  # ... this clump has more psf star candidate neighbours than the others.
477 
478  # get rid of any that are huge, and thus poorly defined
479  goodClumps = []
480  for clump in clumps:
481  if clump.ixx < IzzMax and clump.iyy < IzzMax:
482  goodClumps.append(clump)
483 
484  # if culling > IzzMax cost us all clumps, we'll have to take what we have
485  if len(goodClumps) == 0:
486  goodClumps = clumps
487 
488  # use the 'brightest' clump
489  iBestClump = numpy.argsort(apFluxes)[0]
490  clumps = [clumps[iBestClump]]
491  return clumps
A Threshold is used to pass a threshold value to detection algorithms.
Definition: Threshold.h:44
An integer coordinate rectangle.
Definition: Box.h:53
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
Exposure< ImagePixelT, MaskPixelT, VariancePixelT >::Ptr makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, boost::shared_ptr< Wcs const > wcs=boost::shared_ptr< Wcs const >())
Definition: Exposure.h:308