LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
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 and deblended mask.
336 mask = difference.mask
337 clearMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE", "NOT_DEBLENDED"]
338 for mp in clearMaskPlanes:
339 if mp not in mask.getMaskPlaneDict():
340 mask.addMaskPlane(mp)
341 mask &= ~mask.getPlaneBitMask(clearMaskPlanes)
342
343 # Don't use the idFactory until after deblend+merge, so that we aren't
344 # generating ids that just get thrown away (footprint merge doesn't
345 # know about past ids).
346 table = afwTable.SourceTable.make(self.schema)
347 results = self.detection.run(
348 table=table,
349 exposure=difference,
350 doSmooth=True,
351 )
352
353 sources, positives, negatives = self._deblend(difference,
354 results.positive,
355 results.negative)
356
357 return self.processResults(science, matchedTemplate, difference, sources, idFactory,
358 positiveFootprints=positives,
359 negativeFootprints=negatives)
360
361 def processResults(self, science, matchedTemplate, difference, sources, idFactory,
362 positiveFootprints=None, negativeFootprints=None,):
363 """Measure and process the results of source detection.
364
365 Parameters
366 ----------
367 science : `lsst.afw.image.ExposureF`
368 Science exposure that the template was subtracted from.
369 matchedTemplate : `lsst.afw.image.ExposureF`
370 Warped and PSF-matched template that was used produce the
371 difference image.
372 difference : `lsst.afw.image.ExposureF`
373 Result of subtracting template from the science image.
374 sources : `lsst.afw.table.SourceCatalog`
375 Detected sources on the difference exposure.
376 idFactory : `lsst.afw.table.IdFactory`
377 Generator object used to assign ids to detected sources in the
378 difference image.
379 positiveFootprints : `lsst.afw.detection.FootprintSet`, optional
380 Positive polarity footprints.
381 negativeFootprints : `lsst.afw.detection.FootprintSet`, optional
382 Negative polarity footprints.
383
384 Returns
385 -------
386 measurementResults : `lsst.pipe.base.Struct`
387
388 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
389 Subtracted exposure with detection mask applied.
390 ``diaSources`` : `lsst.afw.table.SourceCatalog`
391 The catalog of detected sources.
392 """
393 self.metadata.add("nUnmergedDiaSources", len(sources))
394 if self.config.doMerge:
395 fpSet = positiveFootprints
396 fpSet.merge(negativeFootprints, self.config.growFootprint,
397 self.config.growFootprint, False)
398 initialDiaSources = afwTable.SourceCatalog(self.schema)
399 fpSet.makeSources(initialDiaSources)
400 self.log.info("Merging detections into %d sources", len(initialDiaSources))
401 else:
402 initialDiaSources = sources
403
404 # Assign source ids at the end: deblend/merge mean that we don't keep
405 # track of parents and children, we only care about the final ids.
406 for source in initialDiaSources:
407 source.setId(idFactory())
408 # Ensure sources added after this get correct ids.
409 initialDiaSources.getTable().setIdFactory(idFactory)
410 initialDiaSources.setMetadata(self.algMetadata)
411
412 self.metadata.add("nMergedDiaSources", len(initialDiaSources))
413
414 if self.config.doMaskStreaks:
415 streakInfo = self._runStreakMasking(difference.maskedImage)
416
417 if self.config.doSkySources:
418 self.addSkySources(initialDiaSources, difference.mask, difference.info.id)
419
420 if not initialDiaSources.isContiguous():
421 initialDiaSources = initialDiaSources.copy(deep=True)
422
423 self.measureDiaSources(initialDiaSources, science, difference, matchedTemplate)
424 diaSources = self._removeBadSources(initialDiaSources)
425
426 if self.config.doForcedMeasurement:
427 self.measureForcedSources(diaSources, science, difference.getWcs())
428
429 self.calculateMetrics(difference)
430
431 measurementResults = pipeBase.Struct(
432 subtractedMeasuredExposure=difference,
433 diaSources=diaSources,
434 )
435 if self.config.doMaskStreaks and self.config.writeStreakInfo:
436 measurementResults.mergeItems(streakInfo, 'maskedStreaks')
437
438 return measurementResults
439
440 def _deblend(self, difference, positiveFootprints, negativeFootprints):
441 """Deblend the positive and negative footprints and return a catalog
442 containing just the children, and the deblended footprints.
443
444 Parameters
445 ----------
446 difference : `lsst.afw.image.Exposure`
447 Result of subtracting template from the science image.
448 positiveFootprints, negativeFootprints : `lsst.afw.detection.FootprintSet`
449 Positive and negative polarity footprints measured on
450 ``difference`` to be deblended separately.
451
452 Returns
453 -------
454 sources : `lsst.afw.table.SourceCatalog`
455 Positive and negative deblended children.
456 positives, negatives : `lsst.afw.detection.FootprintSet`
457 Deblended positive and negative polarity footprints measured on
458 ``difference``.
459 """
460 def makeFootprints(sources):
461 footprints = afwDetection.FootprintSet(difference.getBBox())
462 footprints.setFootprints([src.getFootprint() for src in sources])
463 return footprints
464
465 def deblend(footprints):
466 """Deblend a positive or negative footprint set,
467 and return the deblended children.
468 """
469 sources = afwTable.SourceCatalog(self.schema)
470 footprints.makeSources(sources)
471 self.deblend.run(exposure=difference, sources=sources)
472 self.setPrimaryFlags.run(sources)
473 children = sources["detect_isDeblendedSource"] == 1
474 sources = sources[children].copy(deep=True)
475 # Clear parents, so that measurement plugins behave correctly.
476 sources['parent'] = 0
477 return sources.copy(deep=True)
478
479 positives = deblend(positiveFootprints)
480 negatives = deblend(negativeFootprints)
481
482 sources = afwTable.SourceCatalog(self.schema)
483 sources.reserve(len(positives) + len(negatives))
484 sources.extend(positives, deep=True)
485 sources.extend(negatives, deep=True)
486 return sources, makeFootprints(positives), makeFootprints(negatives)
487
488 def _removeBadSources(self, diaSources):
489 """Remove unphysical diaSources from the catalog.
490
491 Parameters
492 ----------
493 diaSources : `lsst.afw.table.SourceCatalog`
494 The catalog of detected sources.
495
496 Returns
497 -------
498 diaSources : `lsst.afw.table.SourceCatalog`
499 The updated catalog of detected sources, with any source that has a
500 flag in ``config.badSourceFlags`` set removed.
501 """
502 selector = np.ones(len(diaSources), dtype=bool)
503 for flag in self.config.badSourceFlags:
504 flags = diaSources[flag]
505 nBad = np.count_nonzero(flags)
506 if nBad > 0:
507 self.log.debug("Found %d unphysical sources with flag %s.", nBad, flag)
508 selector &= ~flags
509 nBadTotal = np.count_nonzero(~selector)
510 self.metadata.add("nRemovedBadFlaggedSources", nBadTotal)
511 self.log.info("Removed %d unphysical sources.", nBadTotal)
512 return diaSources[selector].copy(deep=True)
513
514 def addSkySources(self, diaSources, mask, seed,
515 subtask=None):
516 """Add sources in empty regions of the difference image
517 for measuring the background.
518
519 Parameters
520 ----------
521 diaSources : `lsst.afw.table.SourceCatalog`
522 The catalog of detected sources.
523 mask : `lsst.afw.image.Mask`
524 Mask plane for determining regions where Sky sources can be added.
525 seed : `int`
526 Seed value to initialize the random number generator.
527 """
528 if subtask is None:
529 subtask = self.skySources
530 skySourceFootprints = subtask.run(mask=mask, seed=seed, catalog=diaSources)
531 self.metadata.add(f"n_{subtask.getName()}", len(skySourceFootprints))
532
533 def measureDiaSources(self, diaSources, science, difference, matchedTemplate):
534 """Use (matched) template and science image to constrain dipole fitting.
535
536 Parameters
537 ----------
538 diaSources : `lsst.afw.table.SourceCatalog`
539 The catalog of detected sources.
540 science : `lsst.afw.image.ExposureF`
541 Science exposure that the template was subtracted from.
542 difference : `lsst.afw.image.ExposureF`
543 Result of subtracting template from the science image.
544 matchedTemplate : `lsst.afw.image.ExposureF`
545 Warped and PSF-matched template that was used produce the
546 difference image.
547 """
548 # Ensure that the required mask planes are present
549 for mp in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere:
550 difference.mask.addMaskPlane(mp)
551 # Note that this may not be correct if we convolved the science image.
552 # In the future we may wish to persist the matchedScience image.
553 self.measurement.run(diaSources, difference, science, matchedTemplate)
554 if self.config.doApCorr:
555 apCorrMap = difference.getInfo().getApCorrMap()
556 if apCorrMap is None:
557 self.log.warning("Difference image does not have valid aperture correction; skipping.")
558 else:
559 self.applyApCorr.run(
560 catalog=diaSources,
561 apCorrMap=apCorrMap,
562 )
563
564 def measureForcedSources(self, diaSources, science, wcs):
565 """Perform forced measurement of the diaSources on the science image.
566
567 Parameters
568 ----------
569 diaSources : `lsst.afw.table.SourceCatalog`
570 The catalog of detected sources.
571 science : `lsst.afw.image.ExposureF`
572 Science exposure that the template was subtracted from.
573 wcs : `lsst.afw.geom.SkyWcs`
574 Coordinate system definition (wcs) for the exposure.
575 """
576 # Run forced psf photometry on the PVI at the diaSource locations.
577 # Copy the measured flux and error into the diaSource.
578 forcedSources = self.forcedMeasurement.generateMeasCat(science, diaSources, wcs)
579 self.forcedMeasurement.run(forcedSources, science, diaSources, wcs)
580 mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema)
581 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFlux")[0],
582 "ip_diffim_forced_PsfFlux_instFlux", True)
583 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFluxErr")[0],
584 "ip_diffim_forced_PsfFlux_instFluxErr", True)
585 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_area")[0],
586 "ip_diffim_forced_PsfFlux_area", True)
587 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag")[0],
588 "ip_diffim_forced_PsfFlux_flag", True)
589 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_noGoodPixels")[0],
590 "ip_diffim_forced_PsfFlux_flag_noGoodPixels", True)
591 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_edge")[0],
592 "ip_diffim_forced_PsfFlux_flag_edge", True)
593 for diaSource, forcedSource in zip(diaSources, forcedSources):
594 diaSource.assign(forcedSource, mapper)
595
596 def calculateMetrics(self, difference):
597 """Add image QA metrics to the Task metadata.
598
599 Parameters
600 ----------
601 difference : `lsst.afw.image.Exposure`
602 The target image to calculate metrics for.
603
604 """
605 mask = difference.mask
606 badPix = (mask.array & mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0
607 self.metadata.add("nGoodPixels", np.sum(~badPix))
608 self.metadata.add("nBadPixels", np.sum(badPix))
609 detPosPix = (mask.array & mask.getPlaneBitMask("DETECTED")) > 0
610 detNegPix = (mask.array & mask.getPlaneBitMask("DETECTED_NEGATIVE")) > 0
611 self.metadata.add("nPixelsDetectedPositive", np.sum(detPosPix))
612 self.metadata.add("nPixelsDetectedNegative", np.sum(detNegPix))
613 detPosPix &= badPix
614 detNegPix &= badPix
615 self.metadata.add("nBadPixelsDetectedPositive", np.sum(detPosPix))
616 self.metadata.add("nBadPixelsDetectedNegative", np.sum(detNegPix))
617
618 metricsMaskPlanes = list(mask.getMaskPlaneDict().keys())
619 for maskPlane in metricsMaskPlanes:
620 try:
621 self.metadata.add("%s_mask_fraction"%maskPlane.lower(), evaluateMaskFraction(mask, maskPlane))
622 except InvalidParameterError:
623 self.metadata.add("%s_mask_fraction"%maskPlane.lower(), -1)
624 self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane)
625
626 def _runStreakMasking(self, maskedImage):
627 """Do streak masking at put results into catalog.
628
629 Parameters
630 ----------
631 maskedImage: `lsst.afw.image.maskedImage`
632 The image in which to search for streaks. Must have a detection
633 mask.
634
635 Returns
636 -------
637 streakInfo: `lsst.pipe.base.Struct`
638 ``rho`` : `np.ndarray`
639 Angle of detected streak.
640 ``theta`` : `np.ndarray`
641 Distance from center of detected streak.
642 ``sigma`` : `np.ndarray`
643 Width of streak profile.
644 """
645 streaks = self.maskStreaks.run(maskedImage)
646 if self.config.writeStreakInfo:
647 rhos = np.array([line.rho for line in streaks.lines])
648 thetas = np.array([line.theta for line in streaks.lines])
649 sigmas = np.array([line.sigma for line in streaks.lines])
650 streakInfo = {'rho': rhos, 'theta': thetas, 'sigma': sigmas}
651 else:
652 streakInfo = {'rho': np.array([]), 'theta': np.array([]), 'sigma': np.array([])}
653 return pipeBase.Struct(maskedStreaks=streakInfo)
654
655
656class DetectAndMeasureScoreConnections(DetectAndMeasureConnections):
657 scoreExposure = pipeBase.connectionTypes.Input(
658 doc="Maximum likelihood image for detection.",
659 dimensions=("instrument", "visit", "detector"),
660 storageClass="ExposureF",
661 name="{fakesType}{coaddName}Diff_scoreExp",
662 )
663
664
665class DetectAndMeasureScoreConfig(DetectAndMeasureConfig,
666 pipelineConnections=DetectAndMeasureScoreConnections):
667 pass
668
669
670class DetectAndMeasureScoreTask(DetectAndMeasureTask):
671 """Detect DIA sources using a score image,
672 and measure the detections on the difference image.
673
674 Source detection is run on the supplied score, or maximum likelihood,
675 image. Note that no additional convolution will be done in this case.
676 Close positive and negative detections will optionally be merged into
677 dipole diaSources.
678 Sky sources, or forced detections in background regions, will optionally
679 be added, and the configured measurement algorithm will be run on all
680 detections.
681 """
682 ConfigClass = DetectAndMeasureScoreConfig
683 _DefaultName = "detectAndMeasureScore"
684
685 @timeMethod
686 def run(self, science, matchedTemplate, difference, scoreExposure,
687 idFactory=None):
688 """Detect and measure sources on a score image.
689
690 Parameters
691 ----------
692 science : `lsst.afw.image.ExposureF`
693 Science exposure that the template was subtracted from.
694 matchedTemplate : `lsst.afw.image.ExposureF`
695 Warped and PSF-matched template that was used produce the
696 difference image.
697 difference : `lsst.afw.image.ExposureF`
698 Result of subtracting template from the science image.
699 scoreExposure : `lsst.afw.image.ExposureF`
700 Score or maximum likelihood difference image
701 idFactory : `lsst.afw.table.IdFactory`, optional
702 Generator object used to assign ids to detected sources in the
703 difference image. Ids from this generator are not set until after
704 deblending and merging positive/negative peaks.
705
706 Returns
707 -------
708 measurementResults : `lsst.pipe.base.Struct`
709
710 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
711 Subtracted exposure with detection mask applied.
712 ``diaSources`` : `lsst.afw.table.SourceCatalog`
713 The catalog of detected sources.
714 """
715 if idFactory is None:
716 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
717
718 # Ensure that we start with an empty detection mask.
719 mask = scoreExposure.mask
720 mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
721
722 # Don't use the idFactory until after deblend+merge, so that we aren't
723 # generating ids that just get thrown away (footprint merge doesn't
724 # know about past ids).
725 table = afwTable.SourceTable.make(self.schema)
726 results = self.detection.run(
727 table=table,
728 exposure=scoreExposure,
729 doSmooth=False,
730 )
731 # Copy the detection mask from the Score image to the difference image
732 difference.mask.assign(scoreExposure.mask, scoreExposure.getBBox())
733
734 sources, positives, negatives = self._deblend(difference,
735 results.positive,
736 results.negative)
737
738 return self.processResults(science, matchedTemplate, difference, sources, idFactory,
739 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)