LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
detectAndMeasure.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
22import numpy as np
23
24import lsst.afw.detection as afwDetection
25import lsst.afw.table as afwTable
26import lsst.daf.base as dafBase
27import lsst.geom
28from lsst.ip.diffim.utils import evaluateMaskFraction
29from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask
30from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig
33import lsst.meas.extensions.shapeHSM
34import lsst.pex.config as pexConfig
35from lsst.pex.exceptions import InvalidParameterError
36import lsst.pipe.base as pipeBase
37import lsst.utils
38from lsst.utils.timer import timeMethod
39
40from . import DipoleFitTask
41
42__all__ = ["DetectAndMeasureConfig", "DetectAndMeasureTask",
43 "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"]
44
45
46class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections,
47 dimensions=("instrument", "visit", "detector"),
48 defaultTemplates={"coaddName": "deep",
49 "warpTypeSuffix": "",
50 "fakesType": ""}):
51 science = pipeBase.connectionTypes.Input(
52 doc="Input science exposure.",
53 dimensions=("instrument", "visit", "detector"),
54 storageClass="ExposureF",
55 name="{fakesType}calexp"
56 )
57 matchedTemplate = pipeBase.connectionTypes.Input(
58 doc="Warped and PSF-matched template used to create the difference image.",
59 dimensions=("instrument", "visit", "detector"),
60 storageClass="ExposureF",
61 name="{fakesType}{coaddName}Diff_matchedExp",
62 )
63 difference = pipeBase.connectionTypes.Input(
64 doc="Result of subtracting template from science.",
65 dimensions=("instrument", "visit", "detector"),
66 storageClass="ExposureF",
67 name="{fakesType}{coaddName}Diff_differenceTempExp",
68 )
69 outputSchema = pipeBase.connectionTypes.InitOutput(
70 doc="Schema (as an example catalog) for output DIASource catalog.",
71 storageClass="SourceCatalog",
72 name="{fakesType}{coaddName}Diff_diaSrc_schema",
73 )
74 diaSources = pipeBase.connectionTypes.Output(
75 doc="Detected diaSources on the difference image.",
76 dimensions=("instrument", "visit", "detector"),
77 storageClass="SourceCatalog",
78 name="{fakesType}{coaddName}Diff_diaSrc",
79 )
80 subtractedMeasuredExposure = pipeBase.connectionTypes.Output(
81 doc="Difference image with detection mask plane filled in.",
82 dimensions=("instrument", "visit", "detector"),
83 storageClass="ExposureF",
84 name="{fakesType}{coaddName}Diff_differenceExp",
85 )
86 maskedStreaks = pipeBase.connectionTypes.Output(
87 doc='Streak profile information.',
88 storageClass="ArrowNumpyDict",
89 dimensions=("instrument", "visit", "detector"),
90 name="{fakesType}{coaddName}Diff_streaks",
91 )
92
93 def __init__(self, *, config):
94 super().__init__(config=config)
95 if not (self.config.writeStreakInfo and self.config.doMaskStreaks):
96 self.outputs.remove("maskedStreaks")
97
98
99class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
100 pipelineConnections=DetectAndMeasureConnections):
101 """Config for DetectAndMeasureTask
102 """
103 doMerge = pexConfig.Field(
104 dtype=bool,
105 default=True,
106 doc="Merge positive and negative diaSources with grow radius "
107 "set by growFootprint"
108 )
109 doForcedMeasurement = pexConfig.Field(
110 dtype=bool,
111 default=True,
112 doc="Force photometer diaSource locations on PVI?")
113 doAddMetrics = pexConfig.Field(
114 dtype=bool,
115 default=False,
116 doc="Add columns to the source table to hold analysis metrics?"
117 )
118 detection = pexConfig.ConfigurableField(
119 target=SourceDetectionTask,
120 doc="Final source detection for diaSource measurement",
121 )
122 deblend = pexConfig.ConfigurableField(
124 doc="Task to split blended sources into their components."
125 )
126 measurement = pexConfig.ConfigurableField(
127 target=DipoleFitTask,
128 doc="Task to measure sources on the difference image.",
129 )
130 doApCorr = lsst.pex.config.Field(
131 dtype=bool,
132 default=True,
133 doc="Run subtask to apply aperture corrections"
134 )
136 target=ApplyApCorrTask,
137 doc="Task to apply aperture corrections"
138 )
139 forcedMeasurement = pexConfig.ConfigurableField(
140 target=ForcedMeasurementTask,
141 doc="Task to force photometer science image at diaSource locations.",
142 )
143 growFootprint = pexConfig.Field(
144 dtype=int,
145 default=2,
146 doc="Grow positive and negative footprints by this many pixels before merging"
147 )
148 diaSourceMatchRadius = pexConfig.Field(
149 dtype=float,
150 default=0.5,
151 doc="Match radius (in arcseconds) for DiaSource to Source association"
152 )
153 doSkySources = pexConfig.Field(
154 dtype=bool,
155 default=False,
156 doc="Generate sky sources?",
157 )
158 skySources = pexConfig.ConfigurableField(
159 target=SkyObjectsTask,
160 doc="Generate sky sources",
161 )
162 doMaskStreaks = pexConfig.Field(
163 dtype=bool,
164 default=False,
165 doc="Turn on streak masking",
166 )
167 maskStreaks = pexConfig.ConfigurableField(
168 target=MaskStreaksTask,
169 doc="Subtask for masking streaks. Only used if doMaskStreaks is True. "
170 "Adds a mask plane to an exposure, with the mask plane name set by streakMaskName.",
171 )
172 writeStreakInfo = pexConfig.Field(
173 dtype=bool,
174 default=False,
175 doc="Record the parameters of any detected streaks. For LSST, this should be turned off except for "
176 "development work."
177 )
178 setPrimaryFlags = pexConfig.ConfigurableField(
179 target=SetPrimaryFlagsTask,
180 doc="Task to add isPrimary and deblending-related flags to the catalog."
181 )
182 badSourceFlags = lsst.pex.config.ListField(
183 dtype=str,
184 doc="Sources with any of these flags set are removed before writing the output catalog.",
185 default=("base_PixelFlags_flag_offimage",
186 "base_PixelFlags_flag_interpolatedCenterAll",
187 "base_PixelFlags_flag_badCenterAll",
188 "base_PixelFlags_flag_edgeCenterAll",
189 "base_PixelFlags_flag_saturatedCenterAll",
190 ),
191 )
192 idGenerator = DetectorVisitIdGeneratorConfig.make_field()
193
194 def setDefaults(self):
195 # DiaSource Detection
196 self.detection.thresholdPolarity = "both"
197 self.detection.thresholdValue = 5.0
198 self.detection.reEstimateBackground = False
199 self.detection.thresholdType = "pixel_stdev"
200 self.detection.excludeMaskPlanes = ["EDGE"]
201
202 # Add filtered flux measurement, the correct measurement for pre-convolved images.
203 self.measurement.algorithms.names.add("base_PeakLikelihoodFlux")
204 self.measurement.plugins.names |= ["ext_trailedSources_Naive",
205 "base_LocalPhotoCalib",
206 "base_LocalWcs",
207 "ext_shapeHSM_HsmSourceMoments",
208 "ext_shapeHSM_HsmPsfMoments",
209 ]
210 self.measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments"
211 self.measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
212 self.measurement.plugins["base_SdssCentroid"].maxDistToPeak = 5.0
213 self.forcedMeasurement.plugins = ["base_TransformedCentroid", "base_PsfFlux"]
214 self.forcedMeasurement.copyColumns = {
215 "id": "objectId", "parent": "parentObjectId", "coord_ra": "coord_ra", "coord_dec": "coord_dec"}
216 self.forcedMeasurement.slots.centroid = "base_TransformedCentroid"
217 self.forcedMeasurement.slots.shape = None
218
219 # Keep track of which footprints contain streaks
220 self.measurement.plugins["base_PixelFlags"].masksFpAnywhere = [
221 "STREAK", "INJECTED", "INJECTED_TEMPLATE"]
222 self.measurement.plugins["base_PixelFlags"].masksFpCenter = [
223 "STREAK", "INJECTED", "INJECTED_TEMPLATE"]
224 self.skySources.avoidMask = ["DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA", "EDGE"]
225
226
227class DetectAndMeasureTask(lsst.pipe.base.PipelineTask):
228 """Detect and measure sources on a difference image.
229 """
230 ConfigClass = DetectAndMeasureConfig
231 _DefaultName = "detectAndMeasure"
232
233 def __init__(self, **kwargs):
234 super().__init__(**kwargs)
235 self.schema = afwTable.SourceTable.makeMinimalSchema()
236 # Add coordinate error fields:
237 afwTable.CoordKey.addErrorFields(self.schema)
238
239 self.algMetadata = dafBase.PropertyList()
240 self.makeSubtask("detection", schema=self.schema)
241 self.makeSubtask("deblend", schema=self.schema)
242 self.makeSubtask("setPrimaryFlags", schema=self.schema, isSingleFrame=True)
243 self.makeSubtask("measurement", schema=self.schema,
244 algMetadata=self.algMetadata)
245 if self.config.doApCorr:
246 self.makeSubtask("applyApCorr", schema=self.measurement.schema)
247 if self.config.doForcedMeasurement:
248 self.schema.addField(
249 "ip_diffim_forced_PsfFlux_instFlux", "D",
250 "Forced PSF flux measured on the direct image.",
251 units="count")
252 self.schema.addField(
253 "ip_diffim_forced_PsfFlux_instFluxErr", "D",
254 "Forced PSF flux error measured on the direct image.",
255 units="count")
256 self.schema.addField(
257 "ip_diffim_forced_PsfFlux_area", "F",
258 "Forced PSF flux effective area of PSF.",
259 units="pixel")
260 self.schema.addField(
261 "ip_diffim_forced_PsfFlux_flag", "Flag",
262 "Forced PSF flux general failure flag.")
263 self.schema.addField(
264 "ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag",
265 "Forced PSF flux not enough non-rejected pixels in data to attempt the fit.")
266 self.schema.addField(
267 "ip_diffim_forced_PsfFlux_flag_edge", "Flag",
268 "Forced PSF flux object was too close to the edge of the image to use the full PSF model.")
269 self.makeSubtask("forcedMeasurement", refSchema=self.schema)
270
271 self.schema.addField("refMatchId", "L", "unique id of reference catalog match")
272 self.schema.addField("srcMatchId", "L", "unique id of source match")
273 if self.config.doSkySources:
274 self.makeSubtask("skySources", schema=self.schema)
275 if self.config.doMaskStreaks:
276 self.makeSubtask("maskStreaks")
277
278 # Check that the schema and config are consistent
279 for flag in self.config.badSourceFlags:
280 if flag not in self.schema:
281 raise pipeBase.InvalidQuantumError("Field %s not in schema" % flag)
282
283 # initialize InitOutputs
284 self.outputSchema = afwTable.SourceCatalog(self.schema)
285 self.outputSchema.getTable().setMetadata(self.algMetadata)
286
287 def runQuantum(self, butlerQC: pipeBase.QuantumContext,
288 inputRefs: pipeBase.InputQuantizedConnection,
289 outputRefs: pipeBase.OutputQuantizedConnection):
290 inputs = butlerQC.get(inputRefs)
291 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
292 idFactory = idGenerator.make_table_id_factory()
293 outputs = self.run(**inputs, idFactory=idFactory)
294 butlerQC.put(outputs, outputRefs)
295
296 @timeMethod
297 def run(self, science, matchedTemplate, difference,
298 idFactory=None):
299 """Detect and measure sources on a difference image.
300
301 The difference image will be convolved with a gaussian approximation of
302 the PSF to form a maximum likelihood image for detection.
303 Close positive and negative detections will optionally be merged into
304 dipole diaSources.
305 Sky sources, or forced detections in background regions, will optionally
306 be added, and the configured measurement algorithm will be run on all
307 detections.
308
309 Parameters
310 ----------
311 science : `lsst.afw.image.ExposureF`
312 Science exposure that the template was subtracted from.
313 matchedTemplate : `lsst.afw.image.ExposureF`
314 Warped and PSF-matched template that was used produce the
315 difference image.
316 difference : `lsst.afw.image.ExposureF`
317 Result of subtracting template from the science image.
318 idFactory : `lsst.afw.table.IdFactory`, optional
319 Generator object used to assign ids to detected sources in the
320 difference image. Ids from this generator are not set until after
321 deblending and merging positive/negative peaks.
322
323 Returns
324 -------
325 measurementResults : `lsst.pipe.base.Struct`
326
327 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
328 Subtracted exposure with detection mask applied.
329 ``diaSources`` : `lsst.afw.table.SourceCatalog`
330 The catalog of detected sources.
331 """
332 if idFactory is None:
333 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
334
335 # Ensure that we start with an empty detection mask.
336 mask = difference.mask
337 mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
338
339 # Don't use the idFactory until after deblend+merge, so that we aren't
340 # generating ids that just get thrown away (footprint merge doesn't
341 # know about past ids).
342 table = afwTable.SourceTable.make(self.schema)
343 results = self.detection.run(
344 table=table,
345 exposure=difference,
346 doSmooth=True,
347 )
348
349 sources, positives, negatives = self._deblend(difference,
350 results.positive,
351 results.negative)
352
353 return self.processResults(science, matchedTemplate, difference, sources, idFactory,
354 positiveFootprints=positives,
355 negativeFootprints=negatives)
356
357 def processResults(self, science, matchedTemplate, difference, sources, idFactory,
358 positiveFootprints=None, negativeFootprints=None,):
359 """Measure and process the results of source detection.
360
361 Parameters
362 ----------
363 science : `lsst.afw.image.ExposureF`
364 Science exposure that the template was subtracted from.
365 matchedTemplate : `lsst.afw.image.ExposureF`
366 Warped and PSF-matched template that was used produce the
367 difference image.
368 difference : `lsst.afw.image.ExposureF`
369 Result of subtracting template from the science image.
370 sources : `lsst.afw.table.SourceCatalog`
371 Detected sources on the difference exposure.
372 idFactory : `lsst.afw.table.IdFactory`
373 Generator object used to assign ids to detected sources in the
374 difference image.
375 positiveFootprints : `lsst.afw.detection.FootprintSet`, optional
376 Positive polarity footprints.
377 negativeFootprints : `lsst.afw.detection.FootprintSet`, optional
378 Negative polarity footprints.
379
380 Returns
381 -------
382 measurementResults : `lsst.pipe.base.Struct`
383
384 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
385 Subtracted exposure with detection mask applied.
386 ``diaSources`` : `lsst.afw.table.SourceCatalog`
387 The catalog of detected sources.
388 """
389 self.metadata.add("nUnmergedDiaSources", len(sources))
390 if self.config.doMerge:
391 fpSet = positiveFootprints
392 fpSet.merge(negativeFootprints, self.config.growFootprint,
393 self.config.growFootprint, False)
394 initialDiaSources = afwTable.SourceCatalog(self.schema)
395 fpSet.makeSources(initialDiaSources)
396 self.log.info("Merging detections into %d sources", len(initialDiaSources))
397 else:
398 initialDiaSources = sources
399
400 # Assign source ids at the end: deblend/merge mean that we don't keep
401 # track of parents and children, we only care about the final ids.
402 for source in initialDiaSources:
403 source.setId(idFactory())
404 # Ensure sources added after this get correct ids.
405 initialDiaSources.getTable().setIdFactory(idFactory)
406 initialDiaSources.setMetadata(self.algMetadata)
407
408 self.metadata.add("nMergedDiaSources", len(initialDiaSources))
409
410 if self.config.doMaskStreaks:
411 streakInfo = self._runStreakMasking(difference.maskedImage)
412
413 if self.config.doSkySources:
414 self.addSkySources(initialDiaSources, difference.mask, difference.info.id)
415
416 if not initialDiaSources.isContiguous():
417 initialDiaSources = initialDiaSources.copy(deep=True)
418
419 self.measureDiaSources(initialDiaSources, science, difference, matchedTemplate)
420 diaSources = self._removeBadSources(initialDiaSources)
421
422 if self.config.doForcedMeasurement:
423 self.measureForcedSources(diaSources, science, difference.getWcs())
424
425 self.calculateMetrics(difference)
426
427 measurementResults = pipeBase.Struct(
428 subtractedMeasuredExposure=difference,
429 diaSources=diaSources,
430 )
431 if self.config.doMaskStreaks and self.config.writeStreakInfo:
432 measurementResults.mergeItems(streakInfo, 'maskedStreaks')
433
434 return measurementResults
435
436 def _deblend(self, difference, positiveFootprints, negativeFootprints):
437 """Deblend the positive and negative footprints and return a catalog
438 containing just the children, and the deblended footprints.
439
440 Parameters
441 ----------
442 difference : `lsst.afw.image.Exposure`
443 Result of subtracting template from the science image.
444 positiveFootprints, negativeFootprints : `lsst.afw.detection.FootprintSet`
445 Positive and negative polarity footprints measured on
446 ``difference`` to be deblended separately.
447
448 Returns
449 -------
450 sources : `lsst.afw.table.SourceCatalog`
451 Positive and negative deblended children.
452 positives, negatives : `lsst.afw.detection.FootprintSet`
453 Deblended positive and negative polarity footprints measured on
454 ``difference``.
455 """
456 def makeFootprints(sources):
457 footprints = afwDetection.FootprintSet(difference.getBBox())
458 footprints.setFootprints([src.getFootprint() for src in sources])
459 return footprints
460
461 def deblend(footprints):
462 """Deblend a positive or negative footprint set,
463 and return the deblended children.
464 """
465 sources = afwTable.SourceCatalog(self.schema)
466 footprints.makeSources(sources)
467 self.deblend.run(exposure=difference, sources=sources)
468 self.setPrimaryFlags.run(sources)
469 children = sources["detect_isDeblendedSource"] == 1
470 sources = sources[children].copy(deep=True)
471 # Clear parents, so that measurement plugins behave correctly.
472 sources['parent'] = 0
473 return sources.copy(deep=True)
474
475 positives = deblend(positiveFootprints)
476 negatives = deblend(negativeFootprints)
477
478 sources = afwTable.SourceCatalog(self.schema)
479 sources.reserve(len(positives) + len(negatives))
480 sources.extend(positives, deep=True)
481 sources.extend(negatives, deep=True)
482 return sources, makeFootprints(positives), makeFootprints(negatives)
483
484 def _removeBadSources(self, diaSources):
485 """Remove unphysical diaSources from the catalog.
486
487 Parameters
488 ----------
489 diaSources : `lsst.afw.table.SourceCatalog`
490 The catalog of detected sources.
491
492 Returns
493 -------
494 diaSources : `lsst.afw.table.SourceCatalog`
495 The updated catalog of detected sources, with any source that has a
496 flag in ``config.badSourceFlags`` set removed.
497 """
498 selector = np.ones(len(diaSources), dtype=bool)
499 for flag in self.config.badSourceFlags:
500 flags = diaSources[flag]
501 nBad = np.count_nonzero(flags)
502 if nBad > 0:
503 self.log.debug("Found %d unphysical sources with flag %s.", nBad, flag)
504 selector &= ~flags
505 nBadTotal = np.count_nonzero(~selector)
506 self.metadata.add("nRemovedBadFlaggedSources", nBadTotal)
507 self.log.info("Removed %d unphysical sources.", nBadTotal)
508 return diaSources[selector].copy(deep=True)
509
510 def addSkySources(self, diaSources, mask, seed,
511 subtask=None):
512 """Add sources in empty regions of the difference image
513 for measuring the background.
514
515 Parameters
516 ----------
517 diaSources : `lsst.afw.table.SourceCatalog`
518 The catalog of detected sources.
519 mask : `lsst.afw.image.Mask`
520 Mask plane for determining regions where Sky sources can be added.
521 seed : `int`
522 Seed value to initialize the random number generator.
523 """
524 if subtask is None:
525 subtask = self.skySources
526 skySourceFootprints = subtask.run(mask=mask, seed=seed, catalog=diaSources)
527 self.metadata.add(f"n_{subtask.getName()}", len(skySourceFootprints))
528
529 def measureDiaSources(self, diaSources, science, difference, matchedTemplate):
530 """Use (matched) template and science image to constrain dipole fitting.
531
532 Parameters
533 ----------
534 diaSources : `lsst.afw.table.SourceCatalog`
535 The catalog of detected sources.
536 science : `lsst.afw.image.ExposureF`
537 Science exposure that the template was subtracted from.
538 difference : `lsst.afw.image.ExposureF`
539 Result of subtracting template from the science image.
540 matchedTemplate : `lsst.afw.image.ExposureF`
541 Warped and PSF-matched template that was used produce the
542 difference image.
543 """
544 # Ensure that the required mask planes are present
545 for mp in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere:
546 difference.mask.addMaskPlane(mp)
547 # Note that this may not be correct if we convolved the science image.
548 # In the future we may wish to persist the matchedScience image.
549 self.measurement.run(diaSources, difference, science, matchedTemplate)
550 if self.config.doApCorr:
551 apCorrMap = difference.getInfo().getApCorrMap()
552 if apCorrMap is None:
553 self.log.warning("Difference image does not have valid aperture correction; skipping.")
554 else:
555 self.applyApCorr.run(
556 catalog=diaSources,
557 apCorrMap=apCorrMap,
558 )
559
560 def measureForcedSources(self, diaSources, science, wcs):
561 """Perform forced measurement of the diaSources on the science image.
562
563 Parameters
564 ----------
565 diaSources : `lsst.afw.table.SourceCatalog`
566 The catalog of detected sources.
567 science : `lsst.afw.image.ExposureF`
568 Science exposure that the template was subtracted from.
569 wcs : `lsst.afw.geom.SkyWcs`
570 Coordinate system definition (wcs) for the exposure.
571 """
572 # Run forced psf photometry on the PVI at the diaSource locations.
573 # Copy the measured flux and error into the diaSource.
574 forcedSources = self.forcedMeasurement.generateMeasCat(science, diaSources, wcs)
575 self.forcedMeasurement.run(forcedSources, science, diaSources, wcs)
576 mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema)
577 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFlux")[0],
578 "ip_diffim_forced_PsfFlux_instFlux", True)
579 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFluxErr")[0],
580 "ip_diffim_forced_PsfFlux_instFluxErr", True)
581 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_area")[0],
582 "ip_diffim_forced_PsfFlux_area", True)
583 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag")[0],
584 "ip_diffim_forced_PsfFlux_flag", True)
585 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_noGoodPixels")[0],
586 "ip_diffim_forced_PsfFlux_flag_noGoodPixels", True)
587 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_edge")[0],
588 "ip_diffim_forced_PsfFlux_flag_edge", True)
589 for diaSource, forcedSource in zip(diaSources, forcedSources):
590 diaSource.assign(forcedSource, mapper)
591
592 def calculateMetrics(self, difference):
593 """Add image QA metrics to the Task metadata.
594
595 Parameters
596 ----------
597 difference : `lsst.afw.image.Exposure`
598 The target image to calculate metrics for.
599
600 """
601 mask = difference.mask
602 badPix = (mask.array & mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0
603 self.metadata.add("nGoodPixels", np.sum(~badPix))
604 self.metadata.add("nBadPixels", np.sum(badPix))
605 detPosPix = (mask.array & mask.getPlaneBitMask("DETECTED")) > 0
606 detNegPix = (mask.array & mask.getPlaneBitMask("DETECTED_NEGATIVE")) > 0
607 self.metadata.add("nPixelsDetectedPositive", np.sum(detPosPix))
608 self.metadata.add("nPixelsDetectedNegative", np.sum(detNegPix))
609 detPosPix &= badPix
610 detNegPix &= badPix
611 self.metadata.add("nBadPixelsDetectedPositive", np.sum(detPosPix))
612 self.metadata.add("nBadPixelsDetectedNegative", np.sum(detNegPix))
613
614 metricsMaskPlanes = list(mask.getMaskPlaneDict().keys())
615 for maskPlane in metricsMaskPlanes:
616 try:
617 self.metadata.add("%s_mask_fraction"%maskPlane.lower(), evaluateMaskFraction(mask, maskPlane))
618 except InvalidParameterError:
619 self.metadata.add("%s_mask_fraction"%maskPlane.lower(), -1)
620 self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane)
621
622 def _runStreakMasking(self, maskedImage):
623 """Do streak masking at put results into catalog.
624
625 Parameters
626 ----------
627 maskedImage: `lsst.afw.image.maskedImage`
628 The image in which to search for streaks. Must have a detection
629 mask.
630
631 Returns
632 -------
633 streakInfo: `lsst.pipe.base.Struct`
634 ``rho`` : `np.ndarray`
635 Angle of detected streak.
636 ``theta`` : `np.ndarray`
637 Distance from center of detected streak.
638 ``sigma`` : `np.ndarray`
639 Width of streak profile.
640 """
641 streaks = self.maskStreaks.run(maskedImage)
642 if self.config.writeStreakInfo:
643 rhos = np.array([line.rho for line in streaks.lines])
644 thetas = np.array([line.theta for line in streaks.lines])
645 sigmas = np.array([line.sigma for line in streaks.lines])
646 streakInfo = {'rho': rhos, 'theta': thetas, 'sigma': sigmas}
647 else:
648 streakInfo = {'rho': np.array([]), 'theta': np.array([]), 'sigma': np.array([])}
649 return pipeBase.Struct(maskedStreaks=streakInfo)
650
651
652class DetectAndMeasureScoreConnections(DetectAndMeasureConnections):
653 scoreExposure = pipeBase.connectionTypes.Input(
654 doc="Maximum likelihood image for detection.",
655 dimensions=("instrument", "visit", "detector"),
656 storageClass="ExposureF",
657 name="{fakesType}{coaddName}Diff_scoreExp",
658 )
659
660
661class DetectAndMeasureScoreConfig(DetectAndMeasureConfig,
662 pipelineConnections=DetectAndMeasureScoreConnections):
663 pass
664
665
666class DetectAndMeasureScoreTask(DetectAndMeasureTask):
667 """Detect DIA sources using a score image,
668 and measure the detections on the difference image.
669
670 Source detection is run on the supplied score, or maximum likelihood,
671 image. Note that no additional convolution will be done in this case.
672 Close positive and negative detections will optionally be merged into
673 dipole diaSources.
674 Sky sources, or forced detections in background regions, will optionally
675 be added, and the configured measurement algorithm will be run on all
676 detections.
677 """
678 ConfigClass = DetectAndMeasureScoreConfig
679 _DefaultName = "detectAndMeasureScore"
680
681 @timeMethod
682 def run(self, science, matchedTemplate, difference, scoreExposure,
683 idFactory=None):
684 """Detect and measure sources on a score image.
685
686 Parameters
687 ----------
688 science : `lsst.afw.image.ExposureF`
689 Science exposure that the template was subtracted from.
690 matchedTemplate : `lsst.afw.image.ExposureF`
691 Warped and PSF-matched template that was used produce the
692 difference image.
693 difference : `lsst.afw.image.ExposureF`
694 Result of subtracting template from the science image.
695 scoreExposure : `lsst.afw.image.ExposureF`
696 Score or maximum likelihood difference image
697 idFactory : `lsst.afw.table.IdFactory`, optional
698 Generator object used to assign ids to detected sources in the
699 difference image. Ids from this generator are not set until after
700 deblending and merging positive/negative peaks.
701
702 Returns
703 -------
704 measurementResults : `lsst.pipe.base.Struct`
705
706 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
707 Subtracted exposure with detection mask applied.
708 ``diaSources`` : `lsst.afw.table.SourceCatalog`
709 The catalog of detected sources.
710 """
711 if idFactory is None:
712 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
713
714 # Ensure that we start with an empty detection mask.
715 mask = scoreExposure.mask
716 mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
717
718 # Don't use the idFactory until after deblend+merge, so that we aren't
719 # generating ids that just get thrown away (footprint merge doesn't
720 # know about past ids).
721 table = afwTable.SourceTable.make(self.schema)
722 results = self.detection.run(
723 table=table,
724 exposure=scoreExposure,
725 doSmooth=False,
726 )
727 # Copy the detection mask from the Score image to the difference image
728 difference.mask.assign(scoreExposure.mask, scoreExposure.getBBox())
729
730 sources, positives, negatives = self._deblend(difference,
731 results.positive,
732 results.negative)
733
734 return self.processResults(science, matchedTemplate, difference, sources, idFactory,
735 positiveFootprints=positives, negativeFootprints=negatives)
A mapping between the keys of two Schemas, used to copy data between them.
Class for storing ordered metadata with comments.
run(self, coaddExposures, bbox, wcs, dataIds, physical_filter=None, **kwargs)