LSSTApplications  18.1.0
LSSTDataManagementBasePackage
dipoleMeasurement.py
Go to the documentation of this file.
1 # This file is part of ip_diffim.
2 #
3 # Developed for the LSST Data Management System.
4 # This product includes software developed by the LSST Project
5 # (https://www.lsst.org).
6 # See the COPYRIGHT file at the top-level directory of this distribution
7 # for details of code ownership.
8 #
9 # This program is free software: you can redistribute it and/or modify
10 # it under the terms of the GNU General Public License as published by
11 # the Free Software Foundation, either version 3 of the License, or
12 # (at your option) any later version.
13 #
14 # This program is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU General Public License for more details.
18 #
19 # You should have received a copy of the GNU General Public License
20 # along with this program. If not, see <https://www.gnu.org/licenses/>.
21 
22 import numpy as np
23 
24 import lsst.afw.geom as afwGeom
25 import lsst.afw.image as afwImage
26 import lsst.pex.config as pexConfig
27 from lsst.log import Log
28 import lsst.meas.deblender.baseline as deblendBaseline
29 from lsst.meas.base.pluginRegistry import register
30 from lsst.meas.base import SingleFrameMeasurementTask, SingleFrameMeasurementConfig, \
31  SingleFramePluginConfig, SingleFramePlugin
32 import lsst.afw.display as afwDisplay
33 
34 __all__ = ("DipoleMeasurementConfig", "DipoleMeasurementTask", "DipoleAnalysis", "DipoleDeblender",
35  "SourceFlagChecker", "ClassificationDipoleConfig", "ClassificationDipolePlugin")
36 
37 
38 class ClassificationDipoleConfig(SingleFramePluginConfig):
39  """Configuration for classification of detected diaSources as dipole or not"""
40  minSn = pexConfig.Field(
41  doc="Minimum quadrature sum of positive+negative lobe S/N to be considered a dipole",
42  dtype=float, default=np.sqrt(2) * 5.0,
43  )
44  maxFluxRatio = pexConfig.Field(
45  doc="Maximum flux ratio in either lobe to be considered a dipole",
46  dtype=float, default=0.65
47  )
48 
49 
50 @register("ip_diffim_ClassificationDipole")
51 class ClassificationDipolePlugin(SingleFramePlugin):
52  """A plugin to classify whether a diaSource is a dipole.
53  """
54 
55  ConfigClass = ClassificationDipoleConfig
56 
57  @classmethod
59  """
60  Returns
61  -------
62  result : `callable`
63  """
64  return cls.APCORR_ORDER
65 
66  def __init__(self, config, name, schema, metadata):
67  SingleFramePlugin.__init__(self, config, name, schema, metadata)
69  self.keyProbability = schema.addField(name + "_value", type="D",
70  doc="Set to 1 for dipoles, else 0.")
71  self.keyFlag = schema.addField(name + "_flag", type="Flag", doc="Set to 1 for any fatal failure.")
72 
73  def measure(self, measRecord, exposure):
74  passesSn = self.dipoleAnalysis.getSn(measRecord) > self.config.minSn
75  negFlux = np.abs(measRecord.get("ip_diffim_PsfDipoleFlux_neg_instFlux"))
76  negFluxFlag = measRecord.get("ip_diffim_PsfDipoleFlux_neg_flag")
77  posFlux = np.abs(measRecord.get("ip_diffim_PsfDipoleFlux_pos_instFlux"))
78  posFluxFlag = measRecord.get("ip_diffim_PsfDipoleFlux_pos_flag")
79 
80  if negFluxFlag or posFluxFlag:
81  self.fail(measRecord)
82  # continue on to classify
83 
84  totalFlux = negFlux + posFlux
85 
86  # If negFlux or posFlux are NaN, these evaluate to False
87  passesFluxNeg = (negFlux / totalFlux) < self.config.maxFluxRatio
88  passesFluxPos = (posFlux / totalFlux) < self.config.maxFluxRatio
89  if (passesSn and passesFluxPos and passesFluxNeg):
90  val = 1.0
91  else:
92  val = 0.0
93 
94  measRecord.set(self.keyProbability, val)
95 
96  def fail(self, measRecord, error=None):
97  measRecord.set(self.keyFlag, True)
98 
99 
101  """Measurement of detected diaSources as dipoles"""
102 
103  def setDefaults(self):
104  SingleFrameMeasurementConfig.setDefaults(self)
105  self.plugins = ["base_CircularApertureFlux",
106  "base_PixelFlags",
107  "base_SkyCoord",
108  "base_PsfFlux",
109  "ip_diffim_NaiveDipoleCentroid",
110  "ip_diffim_NaiveDipoleFlux",
111  "ip_diffim_PsfDipoleFlux",
112  "ip_diffim_ClassificationDipole",
113  ]
114 
115  self.slots.calibFlux = None
116  self.slots.modelFlux = None
117  self.slots.gaussianFlux = None
118  self.slots.shape = None
119  self.slots.centroid = "ip_diffim_NaiveDipoleCentroid"
120  self.doReplaceWithNoise = False
121 
122 
124  """Measurement of Sources, specifically ones from difference images, for characterization as dipoles
125 
126  Parameters
127  ----------
128  sources : 'lsst.afw.table.SourceCatalog'
129  Sources that will be measured
130  badFlags : `list` of `dict`
131  A list of flags that will be used to determine if there was a measurement problem
132 
133  Notes
134  -----
135  The list of badFlags will be used to make a list of keys to check for measurement flags on. By
136  default the centroid keys are added to this list
137 
138  Description
139 
140  This class provides a default configuration for running Source measurement on image differences.
141 
142  .. code-block:: py
143 
144  class DipoleMeasurementConfig(SingleFrameMeasurementConfig):
145  "Measurement of detected diaSources as dipoles"
146  def setDefaults(self):
147  SingleFrameMeasurementConfig.setDefaults(self)
148  self.plugins = ["base_CircularApertureFlux",
149  "base_PixelFlags",
150  "base_SkyCoord",
151  "base_PsfFlux",
152  "ip_diffim_NaiveDipoleCentroid",
153  "ip_diffim_NaiveDipoleFlux",
154  "ip_diffim_PsfDipoleFlux",
155  "ip_diffim_ClassificationDipole",
156  ]
157  self.slots.calibFlux = None
158  self.slots.modelFlux = None
159  self.slots.instFlux = None
160  self.slots.shape = None
161  self.slots.centroid = "ip_diffim_NaiveDipoleCentroid"
162  self.doReplaceWithNoise = False
163 
164  These plugins enabled by default allow the user to test the hypothesis that the Source is a dipole.
165  This includes a set of measurements derived from intermediate base classes
166  DipoleCentroidAlgorithm and DipoleFluxAlgorithm.
167  Their respective algorithm control classes are defined in
168  DipoleCentroidControl and DipoleFluxControl.
169  Each centroid and flux measurement will have _neg (negative)
170  and _pos (positive lobe) fields.
171 
172  The first set of measurements uses a "naive" alrogithm
173  for centroid and flux measurements, implemented in
174  NaiveDipoleCentroidControl and NaiveDipoleFluxControl.
175  The algorithm uses a naive 3x3 weighted moment around
176  the nominal centroids of each peak in the Source Footprint. These algorithms fill the table fields
177  ip_diffim_NaiveDipoleCentroid* and ip_diffim_NaiveDipoleFlux*
178 
179  The second set of measurements undertakes a joint-Psf model on the negative
180  and positive lobe simultaneously. This fit simultaneously solves for the negative and positive
181  lobe centroids and fluxes using non-linear least squares minimization.
182  The fields are stored in table elements ip_diffim_PsfDipoleFlux*.
183 
184  Because this Task is just a config for SingleFrameMeasurementTask, the same result may be acheived by
185  manually editing the config and running SingleFrameMeasurementTask. For example:
186 
187  .. code-block:: py
188 
189  config = SingleFrameMeasurementConfig()
190  config.plugins.names = ["base_PsfFlux",
191  "ip_diffim_PsfDipoleFlux",
192  "ip_diffim_NaiveDipoleFlux",
193  "ip_diffim_NaiveDipoleCentroid",
194  "ip_diffim_ClassificationDipole",
195  "base_CircularApertureFlux",
196  "base_SkyCoord"]
197 
198  config.slots.calibFlux = None
199  config.slots.modelFlux = None
200  config.slots.instFlux = None
201  config.slots.shape = None
202  config.slots.centroid = "ip_diffim_NaiveDipoleCentroid"
203  config.doReplaceWithNoise = False
204 
205  schema = afwTable.SourceTable.makeMinimalSchema()
206  task = SingleFrameMeasurementTask(schema, config=config)-
207 
208  Debug variables
209 
210  The ``lsst.pipe.base.cmdLineTask.CmdLineTask`` command line task interface supports a
211  flag-d/--debug to import debug.py from your PYTHONPATH. The relevant contents of debug.py
212  for this Task include:
213 
214  .. code-block:: py
215 
216  import sys
217  import lsstDebug
218  def DebugInfo(name):
219  di = lsstDebug.getInfo(name)
220  if name == "lsst.ip.diffim.dipoleMeasurement":
221  di.display = True # enable debug output
222  di.maskTransparency = 90 # display mask transparency
223  di.displayDiaSources = True # show exposure with dipole results
224  return di
225  lsstDebug.Info = DebugInfo
226  lsstDebug.frame = 1
227 
228  config.slots.calibFlux = None
229  config.slots.modelFlux = None
230  config.slots.gaussianFlux = None
231  config.slots.shape = None
232  config.slots.centroid = "ip_diffim_NaiveDipoleCentroid"
233  config.doReplaceWithNoise = False
234 
235  This code is dipoleMeasTask.py in the examples directory, and can be run as e.g.
236 
237  .. code-block:: none
238 
239  examples/dipoleMeasTask.py
240  examples/dipoleMeasTask.py --debug
241  examples/dipoleMeasTask.py --debug --image /path/to/image.fits
242 
243 
244 
245  Start the processing by parsing the command line, where the user has the option of
246  enabling debugging output and/or sending their own image for demonstration
247  (in case they have not downloaded the afwdata package).
248 
249  .. code-block:: py
250 
251  if __name__ == "__main__":
252  import argparse
253  parser = argparse.ArgumentParser(
254  description="Demonstrate the use of SourceDetectionTask and DipoleMeasurementTask")
255  parser.add_argument('--debug', '-d', action="store_true", help="Load debug.py?", default=False)
256  parser.add_argument("--image", "-i", help="User defined image", default=None)
257  args = parser.parse_args()
258  if args.debug:
259  try:
260  import debug
261  debug.lsstDebug.frame = 2
262  except ImportError as e:
263  print(e, file=sys.stderr)
264  run(args)
265 
266  The processing occurs in the run function. We first extract an exposure from disk or afwdata, displaying
267  it if requested:
268 
269  .. code-block:: py
270 
271  def run(args):
272  exposure = loadData(args.image)
273  if args.debug:
274  afwDisplay.Display(frame=1).mtv(exposure)
275 
276  Create a default source schema that we will append fields to as we add more algorithms:
277 
278  .. code-block:: py
279 
280  schema = afwTable.SourceTable.makeMinimalSchema()
281 
282  Create the detection and measurement Tasks, with some minor tweaking of their configs:
283 
284  .. code-block:: py
285 
286  # Create the detection task
287  config = SourceDetectionTask.ConfigClass()
288  config.thresholdPolarity = "both"
289  config.background.isNanSafe = True
290  config.thresholdValue = 3
291  detectionTask = SourceDetectionTask(config=config, schema=schema)
292  # And the measurement Task
293  config = DipoleMeasurementTask.ConfigClass()
294  config.plugins.names.remove('base_SkyCoord')
295  algMetadata = dafBase.PropertyList()
296  measureTask = DipoleMeasurementTask(schema, algMetadata, config=config)
297 
298  Having fully initialied the schema, we create a Source table from it:
299 
300  .. code-block:: py
301 
302  # Create the output table
303  tab = afwTable.SourceTable.make(schema)
304 
305  Run detection:
306 
307  .. code-block:: py
308 
309  # Process the data
310  results = detectionTask.run(tab, exposure)
311 
312  Because we are looking for dipoles, we need to merge the positive and negative detections:
313 
314  .. code-block:: py
315 
316  # Merge the positve and negative sources
317  fpSet = results.fpSets.positive
318  growFootprint = 2
319  fpSet.merge(results.fpSets.negative, growFootprint, growFootprint, False)
320  diaSources = afwTable.SourceCatalog(tab)
321  fpSet.makeSources(diaSources)
322  print("Merged %s Sources into %d diaSources (from %d +ve, %d -ve)" % (len(results.sources),
323  len(diaSources),
324  results.fpSets.numPos,
325  results.fpSets.numNeg))
326 
327  Finally, perform measurement (both standard and dipole-specialized) on the merged sources:
328 
329  .. code-block:: py
330 
331  measureTask.run(diaSources, exposure)
332 
333  Optionally display debugging information:
334 
335  .. code-block:: py
336 
337  # Display dipoles if debug enabled
338  if args.debug:
339  dpa = DipoleAnalysis()
340  dpa.displayDipoles(exposure, diaSources)
341 
342  """
343  ConfigClass = DipoleMeasurementConfig
344  _DefaultName = "dipoleMeasurement"
345 
346 
347 
350 
352  """Functor class to check whether a diaSource has flags set that should cause it to be labeled bad."""
353 
354  def __init__(self, sources, badFlags=None):
355  self.badFlags = ['base_PixelFlags_flag_edge', 'base_PixelFlags_flag_interpolatedCenter',
356  'base_PixelFlags_flag_saturatedCenter']
357  if badFlags is not None:
358  for flag in badFlags:
359  self.badFlags.append(flag)
360  self.keys = [sources.getSchema().find(name).key for name in self.badFlags]
361  self.keys.append(sources.table.getCentroidFlagKey())
362 
363  def __call__(self, source):
364  """Call the source flag checker on a single Source
365 
366  Parameters
367  ----------
368  source :
369  Source that will be examined
370  """
371  for k in self.keys:
372  if source.get(k):
373  return False
374  return True
375 
376 
378  """Functor class that provides (S/N, position, orientation) of measured dipoles"""
379 
380  def __init__(self):
381  pass
382 
383  def __call__(self, source):
384  """Parse information returned from dipole measurement
385 
386  Parameters
387  ----------
388  source : `lsst.afw.table.SourceRecord`
389  The source that will be examined"""
390  return self.getSn(source), self.getCentroid(source), self.getOrientation(source)
391 
392  def getSn(self, source):
393  """Get the total signal-to-noise of the dipole; total S/N is from positive and negative lobe
394 
395  Parameters
396  ----------
397  source : `lsst.afw.table.SourceRecord`
398  The source that will be examined"""
399 
400  posflux = source.get("ip_diffim_PsfDipoleFlux_pos_instFlux")
401  posfluxErr = source.get("ip_diffim_PsfDipoleFlux_pos_instFluxErr")
402  negflux = source.get("ip_diffim_PsfDipoleFlux_neg_instFlux")
403  negfluxErr = source.get("ip_diffim_PsfDipoleFlux_neg_instFluxErr")
404 
405  # Not a dipole!
406  if (posflux < 0) is (negflux < 0):
407  return 0
408 
409  return np.sqrt((posflux/posfluxErr)**2 + (negflux/negfluxErr)**2)
410 
411  def getCentroid(self, source):
412  """Get the centroid of the dipole; average of positive and negative lobe
413 
414  Parameters
415  ----------
416  source : `lsst.afw.table.SourceRecord`
417  The source that will be examined"""
418 
419  negCenX = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_x")
420  negCenY = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_y")
421  posCenX = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_x")
422  posCenY = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_y")
423  if (np.isinf(negCenX) or np.isinf(negCenY) or np.isinf(posCenX) or np.isinf(posCenY)):
424  return None
425 
426  center = afwGeom.Point2D(0.5*(negCenX+posCenX),
427  0.5*(negCenY+posCenY))
428  return center
429 
430  def getOrientation(self, source):
431  """Calculate the orientation of dipole; vector from negative to positive lobe
432 
433  Parameters
434  ----------
435  source : `lsst.afw.table.SourceRecord`
436  The source that will be examined"""
437 
438  negCenX = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_x")
439  negCenY = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_y")
440  posCenX = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_x")
441  posCenY = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_y")
442  if (np.isinf(negCenX) or np.isinf(negCenY) or np.isinf(posCenX) or np.isinf(posCenY)):
443  return None
444 
445  dx, dy = posCenX-negCenX, posCenY-negCenY
446  angle = afwGeom.Angle(np.arctan2(dx, dy), afwGeom.radians)
447  return angle
448 
449  def displayDipoles(self, exposure, sources):
450  """Display debugging information on the detected dipoles
451 
452  Parameters
453  ----------
454  exposure : `lsst.afw.image.Exposure`
455  Image the dipoles were measured on
456  sources : `lsst.afw.table.SourceCatalog`
457  The set of diaSources that were measured"""
458 
459  import lsstDebug
460  display = lsstDebug.Info(__name__).display
461  displayDiaSources = lsstDebug.Info(__name__).displayDiaSources
462  maskTransparency = lsstDebug.Info(__name__).maskTransparency
463  if not maskTransparency:
464  maskTransparency = 90
465  disp = afwDisplay.Display(frame=lsstDebug.frame)
466  disp.setMaskTransparency(maskTransparency)
467  disp.mtv(exposure)
468 
469  if display and displayDiaSources:
470  with disp.Buffering():
471  for source in sources:
472  cenX, cenY = source.get("ipdiffim_DipolePsfFlux_centroid")
473  if np.isinf(cenX) or np.isinf(cenY):
474  cenX, cenY = source.getCentroid()
475 
476  isdipole = source.get("ip_diffim_ClassificationDipole_value")
477  if isdipole and np.isfinite(isdipole):
478  # Dipole
479  ctype = afwDisplay.GREEN
480  else:
481  # Not dipole
482  ctype = afwDisplay.RED
483 
484  disp.dot("o", cenX, cenY, size=2, ctype=ctype)
485 
486  negCenX = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_x")
487  negCenY = source.get("ip_diffim_PsfDipoleFlux_neg_centroid_y")
488  posCenX = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_x")
489  posCenY = source.get("ip_diffim_PsfDipoleFlux_pos_centroid_y")
490  if (np.isinf(negCenX) or np.isinf(negCenY) or np.isinf(posCenX) or np.isinf(posCenY)):
491  continue
492 
493  disp.line([(negCenX, negCenY), (posCenX, posCenY)], ctype=afwDisplay.YELLOW)
494 
495  lsstDebug.frame += 1
496 
497 
499  """Functor to deblend a source as a dipole, and return a new source with deblended footprints.
500 
501  This necessarily overrides some of the functionality from
502  meas_algorithms/python/lsst/meas/algorithms/deblend.py since we
503  need a single source that contains the blended peaks, not
504  multiple children sources. This directly calls the core
505  deblending code deblendBaseline.deblend (optionally _fitPsf for
506  debugging).
507 
508  Not actively being used, but there is a unit test for it in
509  dipoleAlgorithm.py.
510  """
511 
512  def __init__(self):
513  # Set up defaults to send to deblender
514 
515  # Always deblend as Psf
516  self.psfChisqCut1 = self.psfChisqCut2 = self.psfChisqCut2b = np.inf
517  self.log = Log.getLogger('ip.diffim.DipoleDeblender')
518  self.sigma2fwhm = 2. * np.sqrt(2. * np.log(2.))
519 
520  def __call__(self, source, exposure):
521  fp = source.getFootprint()
522  peaks = fp.getPeaks()
523  peaksF = [pk.getF() for pk in peaks]
524  fbb = fp.getBBox()
525  fmask = afwImage.Mask(fbb)
526  fmask.setXY0(fbb.getMinX(), fbb.getMinY())
527  fp.spans.setMask(fmask, 1)
528 
529  psf = exposure.getPsf()
530  psfSigPix = psf.computeShape().getDeterminantRadius()
531  psfFwhmPix = psfSigPix * self.sigma2fwhm
532  subimage = afwImage.ExposureF(exposure, bbox=fbb, deep=True)
533  cpsf = deblendBaseline.CachingPsf(psf)
534 
535  # if fewer than 2 peaks, just return a copy of the source
536  if len(peaks) < 2:
537  return source.getTable().copyRecord(source)
538 
539  # make sure you only deblend 2 peaks; take the brighest and faintest
540  speaks = [(p.getPeakValue(), p) for p in peaks]
541  speaks.sort()
542  dpeaks = [speaks[0][1], speaks[-1][1]]
543 
544  # and only set these peaks in the footprint (peaks is mutable)
545  peaks.clear()
546  for peak in dpeaks:
547  peaks.append(peak)
548 
549  if True:
550  # Call top-level deblend task
551  fpres = deblendBaseline.deblend(fp, exposure.getMaskedImage(), psf, psfFwhmPix,
552  log=self.log,
553  psfChisqCut1=self.psfChisqCut1,
554  psfChisqCut2=self.psfChisqCut2,
555  psfChisqCut2b=self.psfChisqCut2b)
556  else:
557  # Call lower-level _fit_psf task
558 
559  # Prepare results structure
560  fpres = deblendBaseline.DeblenderResult(fp, exposure.getMaskedImage(), psf, psfFwhmPix, self.log)
561 
562  for pki, (pk, pkres, pkF) in enumerate(zip(dpeaks, fpres.deblendedParents[0].peaks, peaksF)):
563  self.log.debug('Peak %i', pki)
564  deblendBaseline._fitPsf(fp, fmask, pk, pkF, pkres, fbb, dpeaks, peaksF, self.log,
565  cpsf, psfFwhmPix,
566  subimage.getMaskedImage().getImage(),
567  subimage.getMaskedImage().getVariance(),
568  self.psfChisqCut1, self.psfChisqCut2, self.psfChisqCut2b)
569 
570  deblendedSource = source.getTable().copyRecord(source)
571  deblendedSource.setParent(source.getId())
572  peakList = deblendedSource.getFootprint().getPeaks()
573  peakList.clear()
574 
575  for i, peak in enumerate(fpres.deblendedParents[0].peaks):
576  if peak.psfFitFlux > 0:
577  suffix = "pos"
578  else:
579  suffix = "neg"
580  c = peak.psfFitCenter
581  self.log.info("deblended.centroid.dipole.psf.%s %f %f",
582  suffix, c[0], c[1])
583  self.log.info("deblended.chi2dof.dipole.%s %f",
584  suffix, peak.psfFitChisq / peak.psfFitDof)
585  self.log.info("deblended.flux.dipole.psf.%s %f",
586  suffix, peak.psfFitFlux * np.sum(peak.templateImage.getArray()))
587  peakList.append(peak.peak)
588  return deblendedSource
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
A class representing an angle.
Definition: Angle.h:127
Definition: Log.h:691
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:78
def __init__(self, config, name, schema, metadata)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...