Coverage for python/lsst/ip/diffim/detectAndMeasure.py: 91%
697 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-11-08 01:43 -0800
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/>.
22import dataclasses
23import numpy as np
24import requests
25import os
28import lsst.afw.detection as afwDetection
29import lsst.afw.image as afwImage
30import lsst.afw.math as afwMath
31import lsst.afw.table as afwTable
32import lsst.daf.base as dafBase
33import lsst.geom
34from lsst.images import Mask, get_legacy_difference_image_mask_planes
35from lsst.images.fields import field_from_legacy_background
36from lsst.ip.diffim.utils import (evaluateMaskFraction, computeDifferenceImageMetrics,
37 populate_sattle_visit_cache,
38 get_difference_image_provenance)
39from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask
40from lsst.meas.algorithms import FindGlintTrailsTask, FindCosmicRaysConfig, findCosmicRays
41from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig
42import lsst.meas.deblender
43import lsst.meas.extensions.trailedSources # noqa: F401
44import lsst.meas.extensions.shapeHSM
45import lsst.pex.config as pexConfig
46from lsst.pex.exceptions import InvalidParameterError, LengthError
47import lsst.pipe.base as pipeBase
48import lsst.utils
49from lsst.utils.timer import timeMethod
51from . import DipoleFitTask
53__all__ = ["DetectAndMeasureConfig", "DetectAndMeasureTask",
54 "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"]
56_TEMPLATE_EPOCH_NAMES = {
57 "TEMPLATE MJD MIN": "templateMjdMin",
58 "TEMPLATE MJD MAX": "templateMjdMax",
59 "TEMPLATE MJD MEAN": "templateMjdMean",
60}
61"""Image metadata key of each summary of the template epochs of the
62diaSources, mapped to the name of the diaSource field and task metadata
63entry it summarizes.
64"""
67class BadSubtractionError(pipeBase.AlgorithmError):
68 """Raised when the residuals in footprints of stars used to compute the
69 psf-matching kernel exceeds the configured maximum.
70 """
71 def __init__(self, *, ratio, threshold):
72 msg = ("The ratio of residual power in source footprints on the"
73 " difference image to the power in the footprints on the"
74 f" science image was {ratio}, which exceeds the maximum"
75 f" threshold of {threshold}")
76 super().__init__(msg)
77 self.ratio = ratio
78 self.threshold = threshold
80 @property
81 def metadata(self):
82 return {"ratio": self.ratio,
83 "threshold": self.threshold
84 }
87class NoDiaSourcesError(pipeBase.AlgorithmError):
88 """Raised when there are no diaSources detected on an image difference.
89 """
90 def __init__(self):
91 msg = ("No diaSources detected!")
92 super().__init__(msg)
94 @property
95 def metadata(self):
96 return {}
99class TooManyCosmicRays(pipeBase.AlgorithmError):
100 """Raised if the cosmic ray task fails with too many cosmics.
102 Parameters
103 ----------
104 maxCosmicRays : `int`
105 Maximum number of cosmic rays allowed.
106 """
107 def __init__(self, maxCosmicRays, **kwargs):
108 msg = f"Cosmic ray task found more than {maxCosmicRays} cosmics."
109 self.msg = msg
110 self._metadata = kwargs
111 super().__init__(msg, **kwargs)
112 self._metadata["maxCosmicRays"] = maxCosmicRays
114 def __str__(self):
115 # Exception doesn't handle **kwargs, so we need a custom str.
116 return f"{self.msg}: {self.metadata}"
118 @property
119 def metadata(self):
120 for key, value in self._metadata.items():
121 if not (isinstance(value, int) or isinstance(value, float) or isinstance(value, str)):
122 raise TypeError(f"{key} is of type {type(value)}, but only (int, float, str) are allowed.")
123 return self._metadata
126class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections,
127 dimensions=("instrument", "visit", "detector"),
128 defaultTemplates={"coaddName": "deep",
129 "warpTypeSuffix": "",
130 "fakesType": ""}):
131 science = pipeBase.connectionTypes.Input(
132 doc="Input science exposure.",
133 dimensions=("instrument", "visit", "detector"),
134 storageClass="ExposureF",
135 name="{fakesType}calexp"
136 )
137 matchedTemplate = pipeBase.connectionTypes.Input(
138 doc="Warped and PSF-matched template used to create the difference image.",
139 dimensions=("instrument", "visit", "detector"),
140 storageClass="ExposureF",
141 name="{fakesType}{coaddName}Diff_matchedExp",
142 )
143 difference = pipeBase.connectionTypes.Input(
144 doc="Result of subtracting template from science.",
145 dimensions=("instrument", "visit", "detector"),
146 storageClass="ExposureF",
147 name="{fakesType}{coaddName}Diff_differenceTempExp",
148 )
149 kernelSources = pipeBase.connectionTypes.Input(
150 doc="Final selection of sources used for psf matching.",
151 dimensions=("instrument", "visit", "detector"),
152 storageClass="SourceCatalog",
153 name="{fakesType}{coaddName}Diff_psfMatchSources"
154 )
155 outputSchema = pipeBase.connectionTypes.InitOutput(
156 doc="Schema (as an example catalog) for output DIASource catalog.",
157 storageClass="SourceCatalog",
158 name="{fakesType}{coaddName}Diff_diaSrc_schema",
159 )
160 diaSources = pipeBase.connectionTypes.Output(
161 doc="Detected diaSources on the difference image.",
162 dimensions=("instrument", "visit", "detector"),
163 storageClass="SourceCatalog",
164 name="{fakesType}{coaddName}Diff_diaSrc",
165 )
166 subtractedMeasuredExposure = pipeBase.connectionTypes.Output(
167 doc="Difference image with detection mask plane filled in.",
168 dimensions=("instrument", "visit", "detector"),
169 storageClass="ExposureF",
170 name="{fakesType}{coaddName}Diff_differenceExp",
171 )
172 differenceBackground = pipeBase.connectionTypes.Output(
173 doc="Background model that was subtracted from the difference image.",
174 dimensions=("instrument", "visit", "detector"),
175 storageClass="Background",
176 name="difference_background",
177 )
178 maskedStreaks = pipeBase.connectionTypes.Output(
179 doc='Catalog of streak fit parameters for the difference image.',
180 storageClass="ArrowNumpyDict",
181 dimensions=("instrument", "visit", "detector"),
182 name="{fakesType}{coaddName}Diff_streaks",
183 )
184 glintTrailInfo = pipeBase.connectionTypes.Output(
185 doc='Dict of fit parameters for glint trails in the catalog.',
186 storageClass="ArrowNumpyDict",
187 dimensions=("instrument", "visit", "detector"),
188 name="trailed_glints",
189 )
191 def __init__(self, *, config):
192 super().__init__(config=config)
193 if not (self.config.doCalculateResidualMetics): 193 ↛ 194line 193 didn't jump to line 194 because the condition on line 193 was never true
194 self.inputs.remove("kernelSources")
195 if not (self.config.writeStreakInfo and self.config.doMaskStreaks): 195 ↛ 198line 195 didn't jump to line 198 because the condition on line 195 was always true
196 self.outputs.remove("maskedStreaks")
197 # In future mode the background is recorded on the output images.
198 if (self.config.image_type == "future"
199 or not (self.config.doSubtractBackground and self.config.doWriteBackground)):
200 self.outputs.remove("differenceBackground")
201 if not (self.config.writeGlintInfo): 201 ↛ 202line 201 didn't jump to line 202 because the condition on line 201 was never true
202 self.outputs.remove("glintTrailInfo")
203 if self.config.image_type == "future":
204 self.difference = dataclasses.replace(self.difference, storageClass="DifferenceImage")
205 self.subtractedMeasuredExposure = dataclasses.replace(self.subtractedMeasuredExposure,
206 storageClass="DifferenceImage")
209class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig,
210 pipelineConnections=DetectAndMeasureConnections):
211 """Config for DetectAndMeasureTask
212 """
213 doMerge = pexConfig.Field(
214 dtype=bool,
215 default=True,
216 doc="Merge positive and negative diaSources with grow radius "
217 "set by growFootprint"
218 )
219 doForcedMeasurement = pexConfig.Field(
220 dtype=bool,
221 default=True,
222 doc="Force photometer diaSource locations on PVI?"
223 )
224 doAddMetrics = pexConfig.Field(
225 dtype=bool,
226 default=False,
227 doc="Add columns to the source table to hold analysis metrics?"
228 )
229 doSubtractBackground = pexConfig.Field(
230 dtype=bool,
231 doc="Subtract a background model from the image before detection?",
232 default=True,
233 )
234 doWriteBackground = pexConfig.Field(
235 dtype=bool,
236 doc="Persist the fitted background model? Ignored when"
237 " image_type='future', which records it on the output images.",
238 default=False,
239 deprecated="Use image_type='future', which records the background on the output images."
240 " Will be removed after v31.",
241 )
242 doCalculateResidualMetics = pexConfig.Field(
243 dtype=bool,
244 doc="Calculate metrics to assess image subtraction quality for the task"
245 "metadata?",
246 default=True,
247 )
248 subtractInitialBackground = pexConfig.ConfigurableField(
249 target=lsst.meas.algorithms.SubtractBackgroundTask,
250 doc="Task to perform intial background subtraction, before first detection pass.",
251 )
252 subtractFinalBackground = pexConfig.ConfigurableField(
253 target=lsst.meas.algorithms.SubtractBackgroundTask,
254 doc="Task to perform final background subtraction, after first detection pass.",
255 )
256 doFindCosmicRays = pexConfig.Field(
257 dtype=bool,
258 doc="Detect and mask cosmic rays on the difference image?"
259 "CRs can be interpolated over by setting cosmicray.keepCRs=False",
260 default=True,
261 )
262 cosmicray = pexConfig.ConfigField(
263 dtype=FindCosmicRaysConfig,
264 doc="Options for finding and masking cosmic rays",
265 )
266 detection = pexConfig.ConfigurableField(
267 target=SourceDetectionTask,
268 doc="Final source detection for diaSource measurement",
269 )
270 streakDetection = pexConfig.ConfigurableField(
271 target=SourceDetectionTask,
272 doc="Separate source detection used only for streak masking",
273 )
274 doDeblend = pexConfig.Field(
275 dtype=bool,
276 default=False,
277 doc="Deblend DIASources after detection?"
278 )
279 deblend = pexConfig.ConfigurableField(
280 target=lsst.meas.deblender.SourceDeblendTask,
281 doc="Task to split blended sources into their components."
282 )
283 measurement = pexConfig.ConfigurableField(
284 target=DipoleFitTask,
285 doc="Task to measure sources on the difference image.",
286 )
287 doApCorr = lsst.pex.config.Field(
288 dtype=bool,
289 default=True,
290 doc="Run subtask to apply aperture corrections"
291 )
292 applyApCorr = lsst.pex.config.ConfigurableField(
293 target=ApplyApCorrTask,
294 doc="Task to apply aperture corrections"
295 )
296 forcedMeasurement = pexConfig.ConfigurableField(
297 target=ForcedMeasurementTask,
298 doc="Task to force photometer science image at diaSource locations.",
299 )
300 growFootprint = pexConfig.Field(
301 dtype=int,
302 default=2,
303 doc="Grow positive and negative footprints by this many pixels before merging"
304 )
305 diaSourceMatchRadius = pexConfig.Field(
306 dtype=float,
307 default=0.5,
308 doc="Match radius (in arcseconds) for DiaSource to Source association"
309 )
310 doSkySources = pexConfig.Field(
311 dtype=bool,
312 default=False,
313 doc="Generate sky sources?",
314 )
315 skySources = pexConfig.ConfigurableField(
316 target=SkyObjectsTask,
317 doc="Generate sky sources",
318 )
319 doMaskStreaks = pexConfig.Field(
320 dtype=bool,
321 default=True,
322 doc="Turn on streak masking",
323 )
324 maskStreaks = pexConfig.ConfigurableField(
325 target=MaskStreaksTask,
326 doc="Subtask for masking streaks. Only used if doMaskStreaks is True. "
327 "Adds a mask plane to an exposure, with the mask plane name set by streakMaskName.",
328 )
329 streakBinFactor = pexConfig.Field(
330 dtype=int,
331 default=4,
332 doc="Bin scale factor to use when rerunning detection for masking streaks. "
333 "Only used if doMaskStreaks is True.",
334 )
335 writeStreakInfo = pexConfig.Field(
336 dtype=bool,
337 default=False,
338 doc="Record the parameters of any detected streaks. For LSST, this should be turned off except for "
339 "development work."
340 )
341 findGlints = pexConfig.ConfigurableField(
342 target=FindGlintTrailsTask,
343 doc="Subtask for finding glint trails, usually caused by satellites or debris."
344 )
345 writeGlintInfo = pexConfig.Field(
346 dtype=bool,
347 default=True,
348 doc="Record the parameters of any detected glint trails."
349 )
350 setPrimaryFlags = pexConfig.ConfigurableField(
351 target=SetPrimaryFlagsTask,
352 doc="Task to add isPrimary and deblending-related flags to the catalog."
353 )
354 badSourceFlags = lsst.pex.config.ListField(
355 dtype=str,
356 doc="Sources with any of these flags set are removed before writing the output catalog.",
357 default=("base_PixelFlags_flag_offimage",
358 "base_PixelFlags_flag_interpolatedCenterAll",
359 "base_PixelFlags_flag_badCenter",
360 "base_PixelFlags_flag_edgeCenter",
361 "base_PixelFlags_flag_nodataCenter",
362 "base_PixelFlags_flag_saturatedCenter",
363 "base_PixelFlags_flag_sat_templateCenter",
364 "base_PixelFlags_flag_spikeCenter",
365 ),
366 )
367 clearMaskPlanes = lsst.pex.config.ListField(
368 dtype=str,
369 doc="Mask planes to clear before running detection.",
370 default=("DETECTED", "DETECTED_NEGATIVE", "NOT_DEBLENDED", "STREAK"),
371 )
372 doRejectBadMaskPlaneDetections = pexConfig.Field(
373 dtype=bool,
374 default=True,
375 doc="Reject any peaks detected on ``badMaskPlanes`` before measurement."
376 "These should all be rejected downstream by ``badSourceFlags``, "
377 "but filtering earlier can save time."
378 )
379 badMaskPlanes = lsst.pex.config.ListField(
380 dtype=str,
381 doc="Detections whose footprint peak lies on a pixel with any of these"
382 " mask planes set will be rejected before measurement."
383 " Any missing mask planes will be silently ignored.",
384 default=("NO_DATA", "BAD", "SAT", "EDGE"),
385 )
386 raiseOnBadSubtractionRatio = pexConfig.Field(
387 dtype=bool,
388 default=True,
389 doc="Raise an error if the ratio of power in detected footprints"
390 " on the difference image to the power in footprints on the science"
391 " image exceeds ``badSubtractionRatioThreshold``",
392 )
393 badSubtractionRatioThreshold = pexConfig.Field(
394 dtype=float,
395 default=0.2,
396 doc="Maximum ratio of power in footprints on the difference image to"
397 " the same footprints on the science image."
398 "Only used if ``raiseOnBadSubtractionRatio`` is set",
399 )
400 badSubtractionVariationThreshold = pexConfig.Field(
401 dtype=float,
402 default=0.4,
403 doc="Maximum standard deviation of the ratio of power in footprints on"
404 " the difference image to the same footprints on the science image."
405 "Only used if ``raiseOnBadSubtractionRatio`` is set",
406 )
407 raiseOnNoDiaSources = pexConfig.Field(
408 dtype=bool,
409 default=True,
410 doc="Raise an algorithm error if no diaSources are detected.",
411 )
412 run_sattle = pexConfig.Field(
413 dtype=bool,
414 default=False,
415 doc="If true, dia source bounding boxes will be sent for verification"
416 "to the sattle service."
417 )
418 sattle_historical = pexConfig.Field(
419 dtype=bool,
420 default=False,
421 doc="If re-running a pipeline that requires sattle, this should be set "
422 "to True. This will populate sattle's cache with the historic data "
423 "closest in time to the exposure."
424 )
425 idGenerator = DetectorVisitIdGeneratorConfig.make_field()
426 image_type = pexConfig.ChoiceField[str](
427 "Image type of the input and measured difference images (and, for"
428 " DetectAndMeasureScoreTask, the input and measured score images).",
429 allowed={
430 "legacy": "Read and write lsst.afw.image.ExposureF.",
431 "future": "Read and write lsst.images.DifferenceImage.",
432 },
433 optional=False,
434 default="legacy",
435 )
437 def setDefaults(self):
438 # Background subtraction
439 # Use a small binsize for the first pass to reduce detections on glints
440 # and extended structures. Should not affect the detectability of
441 # faint diaSources
442 self.subtractInitialBackground.binSize = 8
443 self.subtractInitialBackground.useApprox = False
444 self.subtractInitialBackground.statisticsProperty = "MEDIAN"
445 self.subtractInitialBackground.doFilterSuperPixels = True
446 self.subtractInitialBackground.ignoredPixelMask = ["BAD",
447 "EDGE",
448 "DETECTED",
449 "DETECTED_NEGATIVE",
450 "NO_DATA",
451 ]
452 # Use a larger binsize for the final background subtraction, to reduce
453 # over-subtraction of bright objects.
454 self.subtractFinalBackground.binSize = 40
455 self.subtractFinalBackground.useApprox = False
456 self.subtractFinalBackground.statisticsProperty = "MEDIAN"
457 self.subtractFinalBackground.doFilterSuperPixels = True
458 self.subtractFinalBackground.ignoredPixelMask = ["BAD",
459 "EDGE",
460 "DETECTED",
461 "DETECTED_NEGATIVE",
462 "NO_DATA",
463 ]
464 # Cosmic ray detection
465 self.cosmicray.keepCRs = True # do not interpolate over detected CRs
466 # DiaSource Detection
467 self.detection.thresholdPolarity = "both"
468 self.detection.thresholdValue = 5.0
469 self.detection.reEstimateBackground = False
470 self.detection.thresholdType = "pixel_stdev"
471 self.detection.excludeMaskPlanes = []
473 # Copy configs for binned streak detection from the base detection task
474 self.streakDetection.thresholdType = self.detection.thresholdType
475 self.streakDetection.reEstimateBackground = False
476 self.streakDetection.excludeMaskPlanes = self.detection.excludeMaskPlanes
477 self.streakDetection.thresholdValue = self.detection.thresholdValue
478 # Only detect positive streaks
479 self.streakDetection.thresholdPolarity = "positive"
480 # Do not grow detected mask for streaks
481 self.streakDetection.nSigmaToGrow = 0
482 # Set the streak mask along the entire fit line, not only where the
483 # detected mask is set.
484 self.maskStreaks.onlyMaskDetected = False
485 # Restrict streak masking from growing too large
486 self.maskStreaks.maxStreakWidth = 100
487 # Restrict the number of iterations allowed for fitting streaks
488 # When the fit is good it should solve quickly, and exit a bad fit quickly
489 self.maskStreaks.maxFitIter = 10
490 # Only mask to 2 sigma in width
491 self.maskStreaks.nSigmaMask = 2
492 # Threshold for including streaks after the Hough Transform.
493 # A lower value will detect more features that are less linear.
494 self.maskStreaks.absMinimumKernelHeight = 2
496 self.measurement.plugins.names |= ["ext_trailedSources_Naive",
497 "base_LocalPhotoCalib",
498 "base_LocalWcs",
499 "ext_shapeHSM_HsmSourceMoments",
500 "ext_shapeHSM_HsmPsfMoments",
501 "base_ClassificationSizeExtendedness",
502 ]
503 self.measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments"
504 self.measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
505 self.measurement.plugins["base_SdssCentroid"].maxDistToPeak = 5.0
506 self.forcedMeasurement.plugins = ["base_TransformedCentroid", "base_PsfFlux"]
507 self.forcedMeasurement.copyColumns = {
508 "id": "objectId", "parent": "parentObjectId", "coord_ra": "coord_ra", "coord_dec": "coord_dec"}
509 self.forcedMeasurement.slots.centroid = "base_TransformedCentroid"
510 self.forcedMeasurement.slots.shape = None
512 # Keep track of which footprints contain streaks
513 self.measurement.plugins["base_PixelFlags"].masksFpAnywhere = [
514 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SAT_TEMPLATE", "SPIKE"]
515 self.measurement.plugins["base_PixelFlags"].masksFpCenter = [
516 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SAT_TEMPLATE", "SPIKE"]
517 self.skySources.avoidMask = ["DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA", "EDGE"]
519 def validate(self):
520 super().validate()
522 if self.run_sattle:
523 if not os.getenv("SATTLE_URI_BASE"):
524 raise pexConfig.FieldValidationError(DetectAndMeasureConfig.run_sattle, self,
525 "Sattle requested but SATTLE_URI_BASE "
526 "environment variable not set.")
529class DetectAndMeasureTask(lsst.pipe.base.PipelineTask):
530 """Detect and measure sources on a difference image.
531 """
532 ConfigClass = DetectAndMeasureConfig
533 _DefaultName = "detectAndMeasure"
535 def __init__(self, **kwargs):
536 super().__init__(**kwargs)
537 self.schema = afwTable.SourceTable.makeMinimalSchema()
539 self.algMetadata = dafBase.PropertyList()
540 if self.config.doSubtractBackground:
541 self.makeSubtask("subtractInitialBackground")
542 self.makeSubtask("subtractFinalBackground")
543 self.makeSubtask("detection", schema=self.schema)
544 if self.config.doDeblend:
545 self.makeSubtask("deblend", schema=self.schema)
546 self.makeSubtask("setPrimaryFlags", schema=self.schema, isSingleFrame=True)
547 self.makeSubtask("measurement", schema=self.schema,
548 algMetadata=self.algMetadata)
549 if self.config.doApCorr: 549 ↛ 551line 549 didn't jump to line 551 because the condition on line 549 was always true
550 self.makeSubtask("applyApCorr", schema=self.measurement.schema)
551 if self.config.doForcedMeasurement:
552 self.schema.addField(
553 "ip_diffim_forced_PsfFlux_instFlux", "D",
554 "Forced PSF flux measured on the direct image.",
555 units="count")
556 self.schema.addField(
557 "ip_diffim_forced_PsfFlux_instFluxErr", "D",
558 "Forced PSF flux error measured on the direct image.",
559 units="count")
560 self.schema.addField(
561 "ip_diffim_forced_PsfFlux_area", "F",
562 "Forced PSF flux effective area of PSF.",
563 units="pixel")
564 self.schema.addField(
565 "ip_diffim_forced_PsfFlux_flag", "Flag",
566 "Forced PSF flux general failure flag.")
567 self.schema.addField(
568 "ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag",
569 "Forced PSF flux not enough non-rejected pixels in data to attempt the fit.")
570 self.schema.addField(
571 "ip_diffim_forced_PsfFlux_flag_edge", "Flag",
572 "Forced PSF flux object was too close to the edge of the image to use the full PSF model.")
573 self.schema.addField(
574 "ip_diffim_forced_template_PsfFlux_instFlux", "D",
575 "Forced PSF flux measured on the template image.",
576 units="count")
577 self.schema.addField(
578 "ip_diffim_forced_template_PsfFlux_instFluxErr", "D",
579 "Forced PSF flux error measured on the template image.",
580 units="count")
581 self.schema.addField(
582 "ip_diffim_forced_template_PsfFlux_area", "F",
583 "Forced template PSF flux effective area of PSF.",
584 units="pixel")
585 self.schema.addField(
586 "ip_diffim_forced_template_PsfFlux_flag", "Flag",
587 "Forced template PSF flux general failure flag.")
588 self.schema.addField(
589 "ip_diffim_forced_template_PsfFlux_flag_noGoodPixels", "Flag",
590 "Forced template PSF flux not enough non-rejected pixels in data to attempt the fit.")
591 self.schema.addField(
592 "ip_diffim_forced_template_PsfFlux_flag_edge", "Flag",
593 """Forced template PSF flux object was too close to the edge of the image """
594 """to use the full PSF model.""")
595 self.makeSubtask("forcedMeasurement", refSchema=self.schema)
597 self.schema.addField("refMatchId", "L", "unique id of reference catalog match")
598 self.schema.addField("srcMatchId", "L", "unique id of source match")
599 # Create the sky source task for use by metrics,
600 # even if sky sources are not added to the diaSource catalog
601 self.makeSubtask("skySources", schema=self.schema)
602 if self.config.doMaskStreaks: 602 ↛ 605line 602 didn't jump to line 605 because the condition on line 602 was always true
603 self.makeSubtask("maskStreaks")
604 self.makeSubtask("streakDetection")
605 self.makeSubtask("findGlints")
606 self.schema.addField("glint_trail", "Flag", "DiaSource is part of a glint trail.")
607 self.schema.addField("reliability", type="F", doc="Reliability score of the DiaSource")
608 self.schema.addField("reliabilityVersion", type=str, size=7, doc="Version of the reliability model")
609 self.schema.addField("templateMjdMean", type="D", units="d",
610 doc="Weighted mean time (MJD, TAI) of the template inputs that"
611 " cover the diaSource centroid.")
612 self.schema.addField("templateMjdMin", type="D", units="d",
613 doc="Earliest time (MJD, TAI) of the template inputs that"
614 " cover the diaSource centroid.")
615 self.schema.addField("templateMjdMax", type="D", units="d",
616 doc="Latest time (MJD, TAI) of the template inputs that"
617 " cover the diaSource centroid.")
619 # To get the "merge_*" fields in the schema; have to re-initialize
620 # this later, once we have a peak schema post-detection.
621 lsst.afw.detection.FootprintMergeList(self.schema, ["positive", "negative"])
623 # Check that the schema and config are consistent
624 for flag in self.config.badSourceFlags:
625 if flag not in self.schema:
626 raise pipeBase.InvalidQuantumError("Field %s not in schema" % flag)
628 # initialize InitOutputs
629 self.outputSchema = afwTable.SourceCatalog(self.schema)
630 self.outputSchema.getTable().setMetadata(self.algMetadata)
632 def runQuantum(self, butlerQC: pipeBase.QuantumContext,
633 inputRefs: pipeBase.InputQuantizedConnection,
634 outputRefs: pipeBase.OutputQuantizedConnection):
635 inputs = butlerQC.get(inputRefs)
636 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
637 idFactory = idGenerator.make_table_id_factory()
639 # Input image that each output is measured from, keyed by output name.
640 inputImages = {}
641 templateInfo = None
642 if self.config.image_type == "future":
643 dataId = butlerQC.quantum.dataId
644 for output, connection in (("subtractedMeasuredExposure", "difference"),
645 ("scoreMeasuredExposure", "scoreExposure")):
646 if connection not in inputs:
647 continue
648 image = inputs[connection]
649 # The output is written back into this image, so it must
650 # already carry the kernel and the template records.
651 get_difference_image_provenance(image, "kernel", connection, dataId)
652 imageTemplates = get_difference_image_provenance(image, "templates", connection, dataId)
653 if connection == "difference":
654 templateInfo = imageTemplates
655 inputImages[output] = image
656 inputs[connection] = image.to_legacy()
658 # Specify the fields that `annotate` and `convert_outputs_to_future`
659 # need below, to ensure they exist, even as None.
660 measurementResults = pipeBase.Struct(
661 subtractedMeasuredExposure=None,
662 scoreMeasuredExposure=None,
663 diaSources=None,
664 maskedStreaks=None,
665 differenceBackground=None,
666 templateEpochs=None,
667 )
668 try:
669 self.run(**inputs, idFactory=idFactory, measurementResults=measurementResults,
670 templateInfo=templateInfo)
671 except pipeBase.AlgorithmError as e:
672 # Annotate before converting
673 error = pipeBase.AnnotatedPartialOutputsError.annotate(
674 e,
675 self,
676 measurementResults.subtractedMeasuredExposure,
677 measurementResults.scoreMeasuredExposure,
678 measurementResults.diaSources,
679 measurementResults.maskedStreaks,
680 log=self.log
681 )
682 if self.config.image_type == "future":
683 self.convert_outputs_to_future(measurementResults, inputImages)
684 butlerQC.put(measurementResults, outputRefs)
685 raise error from e
686 if self.config.image_type == "future":
687 self.convert_outputs_to_future(measurementResults, inputImages)
688 butlerQC.put(measurementResults, outputRefs)
690 def convert_outputs_to_future(self, results, inputImages):
691 """Write the output images in a result struct back into the
692 `lsst.images` images they were measured from.
694 The image, variance, and mask planes of each input image are replaced
695 by those of the corresponding legacy output, and the background
696 subtracted by this task is attached. The kernel, template records,
697 and observation metadata of the input are kept. Outputs that are
698 absent or `None` are skipped.
700 Parameters
701 ----------
702 results : `lsst.pipe.base.Struct`
703 Output struct to read and modify in place.
704 inputImages : `dict` [`str`, `lsst.images.DifferenceImage`]
705 Input image each output was measured from, keyed by the name of
706 the output. Each is modified in place.
708 Notes
709 -----
710 This uses the default legacy mask plane map of
711 `lsst.images.get_legacy_difference_image_mask_planes` and adds any
712 planes from `~lsst.images.get_legacy_optional_mask_planes` that are set
713 for any pixels in the image.
714 """
715 background = getattr(results, "differenceBackground", None)
716 templateEpochs = getattr(results, "templateEpochs", None) or {}
717 for name, image in inputImages.items():
718 exposure = getattr(results, name, None)
719 if exposure is None: 719 ↛ 720line 719 didn't jump to line 720 because the condition on line 719 was never true
720 continue
721 # These are no-ops when the legacy output is a view of the input,
722 # as returned by `~lsst.images.DifferenceImage.to_legacy`.
723 image.image.array[...] = exposure.image.array
724 image.variance.array[...] = exposure.variance.array
725 image.mask = Mask.from_legacy(exposure.mask, get_legacy_difference_image_mask_planes(),
726 sky_projection=image.sky_projection)
727 image.metadata.update(templateEpochs)
728 if background is not None and len(background) > 0:
729 # The same background is subtracted from the difference image
730 # and the score image, so it applies to both.
731 image.backgrounds.add(
732 "subtracted",
733 field_from_legacy_background(background, bounds=image.bbox, unit=image.unit),
734 description="Background subtracted from this image before detection.",
735 is_subtracted=True,
736 )
737 setattr(results, name, image)
738 if background is not None and len(background) == 0:
739 # An empty BackgroundList cannot be converted to a field, and means
740 # that no background was subtracted.
741 self.log.debug("No background model to attach to the output images.")
743 @timeMethod
744 def run(self, science, matchedTemplate, difference, kernelSources=None,
745 idFactory=None, measurementResults=None, templateInfo=None):
746 """Detect and measure sources on a difference image.
748 The difference image will be convolved with a gaussian approximation of
749 the PSF to form a maximum likelihood image for detection.
750 Close positive and negative detections will optionally be merged into
751 dipole diaSources.
752 Sky sources, or forced detections in background regions, will optionally
753 be added, and the configured measurement algorithm will be run on all
754 detections.
756 Parameters
757 ----------
758 science : `lsst.afw.image.ExposureF`
759 Science exposure that the template was subtracted from.
760 matchedTemplate : `lsst.afw.image.ExposureF`
761 Warped and PSF-matched template that was used produce the
762 difference image.
763 difference : `lsst.afw.image.ExposureF`
764 Result of subtracting template from the science image.
765 kernelSources : `lsst.afw.table.SourceCatalog`, optional
766 Final selection of sources that was used for psf matching.
767 idFactory : `lsst.afw.table.IdFactory`, optional
768 Generator object used to assign ids to detected sources in the
769 difference image. Ids from this generator are not set until after
770 deblending and merging positive/negative peaks.
771 measurementResults : `lsst.pipe.base.Struct`, optional
772 Result struct that is modified to allow saving of partial outputs
773 for some failure conditions. If the task completes successfully,
774 this is also returned.
775 templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \
776 optional
777 Information about the coadds and their input images that the
778 template was built from. Without it the template epochs of every
779 diaSource are NaN.
781 Returns
782 -------
783 measurementResults : `lsst.pipe.base.Struct`
785 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
786 Subtracted exposure with detection mask applied.
787 ``diaSources`` : `lsst.afw.table.SourceCatalog`
788 The catalog of detected sources.
789 ``differenceBackground`` : `lsst.afw.math.BackgroundList`
790 Background that was subtracted from the difference image.
791 ``templateEpochs`` : `dict` [`str`, `float`]
792 Summary of the template epochs of the diaSources, keyed by
793 image metadata key; empty if they are not known.
794 """
795 if measurementResults is None:
796 measurementResults = pipeBase.Struct()
797 if idFactory is None:
798 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
800 # Check image properties and clear detection mask planes.
801 self._prepareInputs(difference)
802 if self.config.doSubtractBackground:
803 background = self._fitAndSubtractBackground(differenceExposure=difference)
804 else:
805 background = afwMath.BackgroundList()
807 if self.config.doFindCosmicRays: 807 ↛ 813line 807 didn't jump to line 813 because the condition on line 807 was always true
808 self.findAndMaskCosmicRays(difference)
810 # Don't use the idFactory until after deblend+merge, so that we aren't
811 # generating ids that just get thrown away (footprint merge doesn't
812 # know about past ids).
813 table = afwTable.SourceTable.make(self.schema)
814 results = self.detection.run(
815 table=table,
816 exposure=difference,
817 doSmooth=True,
818 background=background,
819 clearMask=True,
820 )
821 measurementResults.differenceBackground = background
823 if self.config.doDeblend:
824 sources, positives, negatives = self._deblend(difference,
825 results.positive,
826 results.negative)
828 else:
829 positives = afwTable.SourceCatalog(self.schema)
830 results.positive.makeSources(positives)
831 negatives = afwTable.SourceCatalog(self.schema)
832 results.negative.makeSources(negatives)
833 sources = results.sources
835 self.processResults(science, matchedTemplate, difference,
836 sources, idFactory, kernelSources,
837 positives=positives,
838 negatives=negatives,
839 measurementResults=measurementResults,
840 templateInfo=templateInfo)
841 return measurementResults
843 def _prepareInputs(self, difference):
844 """Ensure that we start with an empty detection and deblended mask.
846 Parameters
847 ----------
848 difference : `lsst.afw.image.ExposureF`
849 The difference image that will be used for detecting diaSources.
850 The mask plane will be modified in place.
852 Raises
853 ------
854 lsst.pipe.base.UpstreamFailureNoWorkFound
855 If the PSF is not usable for measurement.
856 """
857 # Check that we have a valid PSF now before we do more work
858 sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius()
859 if np.isnan(sigma):
860 raise pipeBase.UpstreamFailureNoWorkFound("Invalid PSF detected! PSF width evaluates to NaN.")
861 # Ensure that we start with an empty detection and deblended mask.
862 mask = difference.mask
863 for mp in self.config.clearMaskPlanes:
864 if mp not in mask.getMaskPlaneDict():
865 mask.addMaskPlane(mp)
866 mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes)
868 def _fitAndSubtractBackground(self, differenceExposure, scoreExposure=None):
869 """Fit the final background to ``differenceExposure``.
871 A rough background is first fit on the difference and an internal
872 first-pass detection is run on the (background-subtracted) clone; the
873 resulting detection mask is then transferred to ``differenceExposure``
874 so the final background fit can mask out real sources. This two-pass
875 scheme is what suppresses spurious detections from glints and
876 extended structures: the small-binsize initial fit aggressively
877 over-subtracts those features, and the resulting peak mask lets the
878 large-binsize final fit produce a clean background. When
879 ``scoreExposure`` is supplied (score-image variant), the final
880 background is also subtracted from its image so the caller's final
881 detection runs on a properly background-subtracted score.
883 The first-pass detection's `Struct` is discarded; only the mask it
884 leaves on the clone is used downstream.
886 Parameters
887 ----------
888 differenceExposure : `lsst.afw.image.ExposureF`
889 Difference image. The final background is subtracted from its
890 image in place but its mask is not modified.
891 scoreExposure : `lsst.afw.image.ExposureF`, optional
892 Maximum-likelihood score image. When supplied, the final
893 background is also subtracted from its image.
895 Returns
896 -------
897 background : `lsst.afw.math.BackgroundList`
898 The fit final background.
899 """
900 if scoreExposure is None:
901 detectionExposure = differenceExposure.clone()
902 background = self.subtractInitialBackground.run(detectionExposure).background
903 doSmooth = True
904 else:
905 # Use a clone of differenceExposure because the background is
906 # subtracted in place.
907 background = self.subtractInitialBackground.run(differenceExposure.clone()).background
908 detectionExposure = scoreExposure.clone()
909 detectionExposure.image -= background.getImage()
910 doSmooth = False
911 table = afwTable.SourceTable.make(self.schema)
912 self.detection.run(
913 table=table,
914 exposure=detectionExposure,
915 doSmooth=doSmooth,
916 background=background,
917 clearMask=True,
918 )
919 # Use the temporary detection mask for the final background subtraction,
920 # then restore the original mask. The score image variant never
921 # clears the detection planes of the difference image.
922 mask = differenceExposure.mask
923 originalMask = mask.array.copy()
924 detectedBit = mask.getPlaneBitMask(["DETECTED"])
925 detectedPix = detectionExposure.mask.array & detectedBit > 0
926 detectedNegativeBit = mask.getPlaneBitMask(["DETECTED_NEGATIVE"])
927 detectedNegativePix = detectionExposure.mask.array & detectedNegativeBit > 0
928 mask.array[detectedPix] |= detectedBit
929 mask.array[detectedNegativePix] |= detectedNegativeBit
930 background = self.subtractFinalBackground.run(differenceExposure).background
931 mask.array[...] = originalMask
932 if scoreExposure is not None:
933 # The preconvolution kernel is normalized to 1, so the same
934 # background level applies to the difference and score images.
935 scoreExposure.image -= background.getImage()
936 return background
938 def processResults(self, science, matchedTemplate, difference, sources, idFactory,
939 kernelSources=None, positives=None, negatives=None, measurementResults=None,
940 templateInfo=None):
941 """Measure and process the results of source detection.
943 Parameters
944 ----------
945 science : `lsst.afw.image.ExposureF`
946 Science exposure that the template was subtracted from.
947 matchedTemplate : `lsst.afw.image.ExposureF`
948 Warped and PSF-matched template that was used produce the
949 difference image.
950 difference : `lsst.afw.image.ExposureF`
951 Result of subtracting template from the science image.
952 sources : `lsst.afw.table.SourceCatalog`
953 Detected sources on the difference exposure.
954 idFactory : `lsst.afw.table.IdFactory`
955 Generator object used to assign ids to detected sources in the
956 difference image.
957 kernelSources : `lsst.afw.table.SourceCatalog`, optional
958 Final selection of sources that was used for psf matching.
959 positives : `lsst.afw.table.SourceCatalog`, optional
960 Positive polarity footprints.
961 negatives : `lsst.afw.table.SourceCatalog`, optional
962 Negative polarity footprints.
963 measurementResults : `lsst.pipe.base.Struct`, optional
964 Result struct that is modified to allow saving of partial outputs
965 for some failure conditions. If the task completes successfully,
966 this is also returned.
967 templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \
968 optional
969 Information about the coadds and their input images that the
970 template was built from.
972 Returns
973 -------
974 measurementResults : `lsst.pipe.base.Struct`
976 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
977 Subtracted exposure with detection mask applied.
978 ``diaSources`` : `lsst.afw.table.SourceCatalog`
979 The catalog of detected sources.
980 ``templateEpochs`` : `dict` [`str`, `float`]
981 Summary of the template epochs of the diaSources, keyed by
982 image metadata key; empty if they are not known.
983 """
984 if measurementResults is None: 984 ↛ 985line 984 didn't jump to line 985 because the condition on line 984 was never true
985 measurementResults = pipeBase.Struct()
986 self.metadata["nUnmergedDiaSources"] = len(sources)
987 if self.config.doMerge:
988 # preserve peak schema, if there are any footprints
989 if len(positives) > 0:
990 peakSchema = positives[0].getFootprint().peaks.schema
991 elif len(negatives) > 0: 991 ↛ 994line 991 didn't jump to line 994 because the condition on line 991 was always true
992 peakSchema = negatives[0].getFootprint().peaks.schema
993 else:
994 peakSchema = afwDetection.PeakTable.makeMinimalSchema()
995 mergeList = afwDetection.FootprintMergeList(self.schema,
996 ["positive", "negative"], peakSchema)
997 initialDiaSources = afwTable.SourceCatalog(self.schema)
998 # Start with positive, as FootprintMergeList will self-merge the
999 # subsequent added catalogs, and we want to try to preserve
1000 # deblended positive sources.
1001 mergeList.addCatalog(initialDiaSources.table, positives, "positive", minNewPeakDist=0)
1002 mergeList.addCatalog(initialDiaSources.table, negatives, "negative", minNewPeakDist=0)
1003 mergeList.getFinalSources(initialDiaSources)
1004 # Flag as negative those sources that *only* came from the negative
1005 # footprint set.
1006 initialDiaSources["is_negative"] = initialDiaSources["merge_footprint_negative"] & \
1007 ~initialDiaSources["merge_footprint_positive"]
1008 self.log.info("Merging detections into %d sources", len(initialDiaSources))
1009 # All positive peaks were added before any negative peaks, so
1010 # re-order the peaks in each footprint by significance.
1011 self._reorderPeaksBySignificance(initialDiaSources)
1012 else:
1013 initialDiaSources = sources
1015 # Assign source ids at the end: deblend/merge mean that we don't keep
1016 # track of parents and children, we only care about the final ids.
1017 for source in initialDiaSources:
1018 source.setId(idFactory())
1019 # Ensure sources added after this get correct ids.
1020 initialDiaSources.getTable().setIdFactory(idFactory)
1021 initialDiaSources.setMetadata(self.algMetadata)
1023 self.metadata["nMergedDiaSources"] = len(initialDiaSources)
1025 if self.config.doMaskStreaks: 1025 ↛ 1028line 1025 didn't jump to line 1028 because the condition on line 1025 was always true
1026 streakInfo = self._runStreakMasking(difference)
1028 if self.config.doSkySources:
1029 self.addSkySources(initialDiaSources, difference.mask, difference.info.id)
1031 if self.config.doRejectBadMaskPlaneDetections: 1031 ↛ 1036line 1031 didn't jump to line 1036 because the condition on line 1031 was always true
1032 # Save time by rejecting peaks on bad mask planes prior to
1033 # measurement.
1034 initialDiaSources = self._rejectBadMaskedDetections(initialDiaSources, difference.mask)
1036 if not initialDiaSources.isContiguous(): 1036 ↛ 1037line 1036 didn't jump to line 1037 because the condition on line 1036 was never true
1037 initialDiaSources = initialDiaSources.copy(deep=True)
1039 self.measureDiaSources(initialDiaSources, science, difference, matchedTemplate)
1041 # Remove unphysical diaSources per config.badSourceFlags
1042 diaSources = self._removeBadSources(initialDiaSources)
1044 if self.config.run_sattle:
1045 diaSources = self.filterSatellites(diaSources, science)
1047 # Flag diaSources in glint trails, but do not remove them
1048 diaSources, trail_parameters = self._find_glint_trails(diaSources)
1049 if self.config.writeGlintInfo: 1049 ↛ 1052line 1049 didn't jump to line 1052 because the condition on line 1049 was always true
1050 measurementResults.mergeItems(trail_parameters, 'glintTrailInfo')
1052 self._measureTemplateEpochs(diaSources, templateInfo)
1053 measurementResults.templateEpochs = self._summarizeTemplateEpochs(diaSources)
1055 if self.config.doForcedMeasurement:
1056 self.measureForcedSources(diaSources, science, difference.getWcs())
1057 self.measureForcedSources(diaSources, matchedTemplate, difference.getWcs(),
1058 template=True)
1060 # Clear the image plane for regions with NO_DATA.
1061 # These regions are most often caused by insufficient template coverage.
1062 # Do this for the final difference image after detection and measurement
1063 # since the subtasks should all be configured to handle NO_DATA properly
1064 difference.image.array[difference.mask.array & difference.mask.getPlaneBitMask('NO_DATA') > 0] = 0
1066 measurementResults.subtractedMeasuredExposure = difference
1068 if self.config.doMaskStreaks and self.config.writeStreakInfo: 1068 ↛ 1069line 1068 didn't jump to line 1069 because the condition on line 1068 was never true
1069 measurementResults.maskedStreaks = streakInfo.maskedStreaks
1071 if kernelSources is not None:
1072 self.calculateMetrics(science, difference, diaSources, kernelSources)
1074 if np.count_nonzero(~diaSources["sky_source"]) > 0:
1075 measurementResults.diaSources = diaSources
1076 elif self.config.raiseOnNoDiaSources:
1077 raise NoDiaSourcesError()
1078 elif len(diaSources) > 0: 1078 ↛ 1082line 1078 didn't jump to line 1082 because the condition on line 1078 was always true
1079 # This option allows returning sky sources,
1080 # even if there are no diaSources
1081 measurementResults.diaSources = diaSources
1082 self.log.info("Measured %d diaSources and %d sky sources",
1083 np.count_nonzero(~diaSources["sky_source"]),
1084 np.count_nonzero(diaSources["sky_source"])
1085 )
1086 return measurementResults
1088 def _reorderPeaksBySignificance(self, diaSources):
1089 """Sort each merged source's footprint peaks by significance.
1091 Parameters
1092 ----------
1093 diaSources : `lsst.afw.table.SourceCatalog`
1094 Merged sources whose footprint peaks are re-ordered in place.
1095 """
1096 for source in diaSources:
1097 peaks = source.getFootprint().getPeaks()
1098 if len(peaks) < 2 or "significance" not in peaks.schema.getNames():
1099 continue
1100 ordered = sorted(peaks, key=lambda peak: peak["significance"], reverse=True)
1101 newPeaks = afwDetection.PeakCatalog(peaks.table)
1102 newPeaks.extend(ordered)
1103 peaks.clear()
1104 source.getFootprint().setPeakCatalog(newPeaks)
1106 def _deblend(self, difference, positiveFootprints, negativeFootprints):
1107 """Deblend the positive and negative footprints and return a catalog
1108 containing just the children, and the deblended footprints.
1110 Parameters
1111 ----------
1112 difference : `lsst.afw.image.Exposure`
1113 Result of subtracting template from the science image.
1114 positiveFootprints, negativeFootprints : `lsst.afw.detection.FootprintSet`
1115 Positive and negative polarity footprints measured on
1116 ``difference`` to be deblended separately.
1118 Returns
1119 -------
1120 sources : `lsst.afw.table.SourceCatalog`
1121 Positive and negative deblended children.
1122 positives, negatives : `lsst.afw.table.SourceCatalog`
1123 Deblended positive and negative polarity sources with footprints
1124 detected on ``difference``.
1125 """
1127 def deblend(footprints, negative=False):
1128 """Deblend a positive or negative footprint set,
1129 and return the deblended children.
1131 Parameters
1132 ----------
1133 footprints : `lsst.afw.detection.FootprintSet`
1134 negative : `bool`
1135 Set True if the footprints contain negative fluxes
1137 Returns
1138 -------
1139 sources : `lsst.afw.table.SourceCatalog`
1140 """
1141 sources = afwTable.SourceCatalog(self.schema)
1142 footprints.makeSources(sources)
1143 if negative:
1144 # Invert the image so the deblender can run on positive peaks
1145 difference_inverted = difference.clone()
1146 difference_inverted.image *= -1
1147 self.deblend.run(exposure=difference_inverted, sources=sources)
1148 children = sources[sources["parent"] != 0]
1149 # Set the heavy footprint pixel values back to reality
1150 for child in children:
1151 footprint = child.getFootprint()
1152 array = footprint.getImageArray()
1153 array *= -1
1154 else:
1155 self.deblend.run(exposure=difference, sources=sources)
1156 self.setPrimaryFlags.run(sources)
1157 children = sources["detect_isDeblendedSource"] == 1
1158 sources = sources[children].copy(deep=True)
1159 # Clear parents, so that measurement plugins behave correctly.
1160 sources['parent'] = 0
1161 return sources.copy(deep=True)
1163 positives = deblend(positiveFootprints)
1164 negatives = deblend(negativeFootprints, negative=True)
1166 sources = afwTable.SourceCatalog(self.schema)
1167 sources.reserve(len(positives) + len(negatives))
1168 sources.extend(positives, deep=True)
1169 sources.extend(negatives, deep=True)
1170 if len(negatives) > 0:
1171 sources[-len(negatives):]["is_negative"] = True
1172 return sources, positives, negatives
1174 def _rejectBadMaskedDetections(self, initialDiaSources, mask):
1175 """Remove detections whose footprint peak lies on a bad-masked pixel.
1177 Detections are rejected if any peak of the source footprint falls on a
1178 pixel with any of the ``config.badMaskPlanes`` bits set.
1180 Parameters
1181 ----------
1182 initialDiaSources : `lsst.afw.table.SourceCatalog`
1183 The catalog of detected peaks.
1184 mask : `lsst.afw.image.Mask`
1185 Mask of the image used for detection.
1187 Returns
1188 -------
1189 initialDiaSources : `lsst.afw.table.SourceCatalog`
1190 The catalog with bad-masked detections removed.
1191 """
1192 if not self.config.badMaskPlanes or len(initialDiaSources) == 0: 1192 ↛ 1193line 1192 didn't jump to line 1193 because the condition on line 1192 was never true
1193 self.metadata["nRejectedBadMaskPlanes"] = 0
1194 return initialDiaSources
1196 presentPlanes = [mp for mp in self.config.badMaskPlanes
1197 if mp in mask.getMaskPlaneDict()]
1198 if not presentPlanes: 1198 ↛ 1199line 1198 didn't jump to line 1199 because the condition on line 1198 was never true
1199 self.metadata["nRejectedBadMaskPlanes"] = 0
1200 return initialDiaSources
1202 badBitMask = mask.getPlaneBitMask(presentPlanes)
1203 x0, y0 = mask.getXY0()
1205 selector = np.ones(len(initialDiaSources), dtype=bool)
1206 for i, source in enumerate(initialDiaSources):
1207 peaks = source.getFootprint().getPeaks()
1208 for peak in peaks:
1209 ix = peak.getIx() - x0
1210 iy = peak.getIy() - y0
1211 if mask.array[iy, ix] & badBitMask:
1212 selector[i] = False
1213 break
1214 nRejected = np.count_nonzero(~selector)
1215 self.metadata["nRejectedBadMaskPlanes"] = nRejected
1216 if nRejected > 0:
1217 self.log.info("Rejected %d detections on pixels with bad mask "
1218 "planes %s before measurement.",
1219 nRejected, presentPlanes)
1220 return initialDiaSources[selector].copy(deep=True)
1222 def _removeBadSources(self, diaSources):
1223 """Remove unphysical diaSources from the catalog.
1225 Parameters
1226 ----------
1227 diaSources : `lsst.afw.table.SourceCatalog`
1228 The catalog of detected sources.
1230 Returns
1231 -------
1232 diaSources : `lsst.afw.table.SourceCatalog`
1233 The updated catalog of detected sources, with any source that has a
1234 flag in ``config.badSourceFlags`` set removed.
1235 """
1236 selector = np.ones(len(diaSources), dtype=bool)
1237 for flag in self.config.badSourceFlags:
1238 flags = diaSources[flag]
1239 nBad = np.count_nonzero(flags)
1240 if nBad > 0:
1241 self.log.debug("Found %d unphysical sources with flag %s.", nBad, flag)
1242 selector &= ~flags
1243 nBadTotal = np.count_nonzero(~selector)
1244 self.metadata["nRemovedBadFlaggedSources"] = nBadTotal
1245 self.log.info("Removed %d unphysical sources.", nBadTotal)
1246 return diaSources[selector].copy(deep=True)
1248 def _find_glint_trails(self, diaSources):
1249 """Define a new flag column for diaSources that are in a glint trail.
1251 Parameters
1252 ----------
1253 diaSources : `lsst.afw.table.SourceCatalog`
1254 The catalog of detected sources.
1256 Returns
1257 -------
1258 diaSources : `lsst.afw.table.SourceCatalog`
1259 The updated catalog of detected sources, with a new bool column
1260 called 'glint_trail' added.
1262 trail_parameters : `dict`
1263 Parameters of all the trails that were found.
1264 """
1265 if self.config.doSkySources:
1266 # Do not include sky sources in glint detection
1267 candidateDiaSources = diaSources[~diaSources["sky_source"]].copy(deep=True)
1268 else:
1269 candidateDiaSources = diaSources
1270 trailed_glints = self.findGlints.run(candidateDiaSources)
1271 glint_mask = [True if id in trailed_glints.trailed_ids else False for id in diaSources['id']]
1272 if np.any(glint_mask):
1273 diaSources['glint_trail'] = np.array(glint_mask)
1275 slopes = np.array([trail.slope for trail in trailed_glints.parameters])
1276 intercepts = np.array([trail.intercept for trail in trailed_glints.parameters])
1277 stderrs = np.array([trail.stderr for trail in trailed_glints.parameters])
1278 lengths = np.array([trail.length for trail in trailed_glints.parameters])
1279 angles = np.array([trail.angle for trail in trailed_glints.parameters])
1280 center_xs = np.array([trail.center_x for trail in trailed_glints.parameters])
1281 center_ys = np.array([trail.center_y for trail in trailed_glints.parameters])
1282 parameters = {'slopes': slopes, 'intercepts': intercepts, 'stderrs': stderrs, 'lengths': lengths,
1283 'angles': angles, 'center_xs': center_xs, 'center_ys': center_ys}
1285 trail_parameters = pipeBase.Struct(glintTrailInfo=parameters)
1287 return diaSources, trail_parameters
1289 def _measureTemplateEpochs(self, diaSources, templateInfo):
1290 """Determine the observed epochs of the template for each diaSource.
1292 Parameters
1293 ----------
1294 diaSources : `lsst.afw.table.SourceCatalog`
1295 Contiguous catalog to fill in place.
1296 templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`] or \
1297 `None`
1298 Information about the coadds and their input images that the
1299 template was built from.
1301 Notes
1302 -----
1303 An input covers a diaSource if both its coadd's bounds and its own
1304 bounds contain the centroid. Each coadd's inputs give a mean weighted
1305 by their coadd weights, and where coadds overlap their means are
1306 averaged with weights ``sqrt(W)``, ``W`` being the total weight of a
1307 coadd's inputs there. `GetTemplateTask` merges the patches of a tract
1308 with weights ``variance**-0.5``, and the variance of a coadd is about
1309 ``1/W``.
1310 """
1311 x = diaSources.getX()
1312 y = diaSources.getY()
1313 numerator = np.zeros(len(diaSources))
1314 denominator = np.zeros(len(diaSources))
1315 epochMin = np.full(len(diaSources), np.inf)
1316 epochMax = np.full(len(diaSources), -np.inf)
1317 for info in templateInfo or []:
1318 inCoadd = info.bounds.contains(x=x, y=y)
1319 if not inCoadd.any():
1320 continue
1321 weightSum = np.zeros(len(diaSources))
1322 weightedSum = np.zeros(len(diaSources))
1323 for record in info.inputs:
1324 covers = inCoadd & record.bounds.contains(x=x, y=y)
1325 weightSum[covers] += record.weight
1326 weightedSum[covers] += record.weight*record.mid_mjd_tai
1327 epochMin[covers] = np.minimum(epochMin[covers], record.mid_mjd_tai)
1328 epochMax[covers] = np.maximum(epochMax[covers], record.mid_mjd_tai)
1329 inputs = weightSum > 0
1330 numerator[inputs] += weightedSum[inputs]/np.sqrt(weightSum[inputs])
1331 denominator[inputs] += np.sqrt(weightSum[inputs])
1332 known = denominator > 0
1333 epoch = np.full(len(diaSources), np.nan)
1334 epoch[known] = numerator[known]/denominator[known]
1335 epochMin[~known] = np.nan
1336 epochMax[~known] = np.nan
1337 diaSources["templateMjdMean"] = epoch
1338 diaSources["templateMjdMin"] = epochMin
1339 diaSources["templateMjdMax"] = epochMax
1341 def _summarizeTemplateEpochs(self, diaSources):
1342 """Summarize the template epochs of a catalog in the task metadata.
1344 The earliest, latest and mean values over the rows with a template
1345 epoch, sky sources included, are set in the task metadata under the
1346 names of the diaSource fields they summarize. They are NaN if no row
1347 has a template epoch.
1349 Parameters
1350 ----------
1351 diaSources : `lsst.afw.table.SourceCatalog`
1352 Contiguous catalog with the template epoch fields set.
1354 Returns
1355 -------
1356 templateEpochs : `dict` [`str`, `float`]
1357 The same values keyed by their image metadata key in
1358 ``_TEMPLATE_EPOCH_NAMES``, or empty if no row has a template epoch.
1359 """
1360 known = np.isfinite(diaSources["templateMjdMean"])
1361 if not known.any():
1362 for name in _TEMPLATE_EPOCH_NAMES.values():
1363 self.metadata[name] = np.nan
1364 return {}
1365 values = {"templateMjdMin": float(np.min(diaSources["templateMjdMin"][known])),
1366 "templateMjdMax": float(np.max(diaSources["templateMjdMax"][known])),
1367 "templateMjdMean": float(np.mean(diaSources["templateMjdMean"][known]))}
1368 for name, value in values.items():
1369 self.metadata[name] = value
1370 return {key: values[name] for key, name in _TEMPLATE_EPOCH_NAMES.items()}
1372 def addSkySources(self, diaSources, mask, seed,
1373 subtask=None):
1374 """Add sources in empty regions of the difference image
1375 for measuring the background.
1377 Parameters
1378 ----------
1379 diaSources : `lsst.afw.table.SourceCatalog`
1380 The catalog of detected sources.
1381 mask : `lsst.afw.image.Mask`
1382 Mask plane for determining regions where Sky sources can be added.
1383 seed : `int`
1384 Seed value to initialize the random number generator.
1385 """
1386 if subtask is None: 1386 ↛ 1388line 1386 didn't jump to line 1388 because the condition on line 1386 was always true
1387 subtask = self.skySources
1388 if subtask.config.nSources <= 0:
1389 self.metadata[f"n_{subtask.getName()}"] = 0
1390 return
1391 skySourceFootprints = subtask.run(mask=mask, seed=seed, catalog=diaSources)
1392 self.metadata[f"n_{subtask.getName()}"] = len(skySourceFootprints)
1394 def measureDiaSources(self, diaSources, science, difference, matchedTemplate):
1395 """Use (matched) template and science image to constrain dipole fitting.
1397 Parameters
1398 ----------
1399 diaSources : `lsst.afw.table.SourceCatalog`
1400 The catalog of detected sources.
1401 science : `lsst.afw.image.ExposureF`
1402 Science exposure that the template was subtracted from.
1403 difference : `lsst.afw.image.ExposureF`
1404 Result of subtracting template from the science image.
1405 matchedTemplate : `lsst.afw.image.ExposureF`
1406 Warped and PSF-matched template that was used produce the
1407 difference image.
1408 """
1409 # Ensure that the required mask planes are present
1410 for mp in self.config.measurement.plugins["base_PixelFlags"].masksFpAnywhere:
1411 difference.mask.addMaskPlane(mp)
1412 # Note that this may not be correct if we convolved the science image.
1413 # In the future we may wish to persist the matchedScience image.
1414 self.measurement.run(diaSources, difference, science, matchedTemplate)
1415 if self.config.doApCorr: 1415 ↛ exitline 1415 didn't return from function 'measureDiaSources' because the condition on line 1415 was always true
1416 apCorrMap = difference.getInfo().getApCorrMap()
1417 if apCorrMap is None: 1417 ↛ 1420line 1417 didn't jump to line 1420 because the condition on line 1417 was always true
1418 self.log.warning("Difference image does not have valid aperture correction; skipping.")
1419 else:
1420 self.applyApCorr.run(
1421 catalog=diaSources,
1422 apCorrMap=apCorrMap,
1423 )
1425 def measureForcedSources(self, diaSources, image, wcs, template=False):
1426 """Perform forced measurement of the diaSources on a direct image.
1428 Parameters
1429 ----------
1430 diaSources : `lsst.afw.table.SourceCatalog`
1431 The catalog of detected sources.
1432 image: `lsst.afw.image.ExposureF`
1433 Exposure that the forced measurement is being performed on
1434 wcs : `lsst.afw.geom.SkyWcs`
1435 Coordinate system definition (wcs) for the exposure.
1436 template : `bool`
1437 Is the forced measurement being performed on the template?
1438 """
1439 # Run forced psf photometry on the image at the diaSource locations.
1440 # Copy the measured flux and error into the diaSource.
1441 forcedSources = self.forcedMeasurement.generateMeasCat(image, diaSources, wcs)
1442 self.forcedMeasurement.run(forcedSources, image, diaSources, wcs)
1444 if template:
1445 base_key = 'ip_diffim_forced_template_PsfFlux'
1446 else:
1447 base_key = 'ip_diffim_forced_PsfFlux'
1448 mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema)
1449 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFlux")[0],
1450 f"{base_key}_instFlux", True)
1451 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFluxErr")[0],
1452 f"{base_key}_instFluxErr", True)
1453 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_area")[0],
1454 f"{base_key}_area", True)
1455 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag")[0],
1456 f"{base_key}_flag", True)
1457 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_noGoodPixels")[0],
1458 f"{base_key}_flag_noGoodPixels", True)
1459 mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_edge")[0],
1460 f"{base_key}_flag_edge", True)
1461 for diaSource, forcedSource in zip(diaSources, forcedSources):
1462 diaSource.assign(forcedSource, mapper)
1464 def findAndMaskCosmicRays(self, difference):
1465 """Detect and mask cosmic rays on the difference image.
1467 Parameters
1468 ----------
1469 difference : `lsst.afw.image.Exposure`
1470 The background-subtracted difference image.
1471 The mask plane will be modified in place.
1472 """
1473 # This is run on a background-subtracted difference image, so the
1474 # remaining background should be ~0.
1475 medianBg = 0.
1476 try:
1477 crs = findCosmicRays(
1478 difference.maskedImage,
1479 difference.psf,
1480 medianBg,
1481 pexConfig.makePropertySet(self.config.cosmicray),
1482 self.config.cosmicray.keepCRs,
1483 )
1484 except LengthError:
1485 raise TooManyCosmicRays(self.config.cosmicray.nCrPixelMax) from None
1486 num = 0
1487 if crs is not None: 1487 ↛ 1493line 1487 didn't jump to line 1493 because the condition on line 1487 was always true
1488 mask = difference.maskedImage.mask
1489 crBit = mask.getPlaneBitMask("CR")
1490 afwDetection.setMaskFromFootprintList(mask, crs, crBit)
1491 num = len(crs)
1493 text = "masked" if self.config.cosmicray.keepCRs else "interpolated over"
1494 self.log.info("Identified and %s %s cosmic rays.", text, num)
1495 self.metadata["cosmic_ray_count"] = num
1497 def calculateMetrics(self, science, difference, diaSources, kernelSources):
1498 """Add difference image QA metrics to the Task metadata.
1500 This may be used to produce corresponding metrics (see
1501 lsst.analysis.tools.tasks.diffimTaskDetectorVisitMetricAnalysis).
1503 Parameters
1504 ----------
1505 science : `lsst.afw.image.ExposureF`
1506 Science exposure that was subtracted.
1507 difference : `lsst.afw.image.Exposure`
1508 The target difference image to calculate metrics for.
1509 diaSources : `lsst.afw.table.SourceCatalog`
1510 The catalog of detected sources.
1511 kernelSources : `lsst.afw.table.SourceCatalog`
1512 Final selection of sources that was used for psf matching.
1513 """
1514 mask = difference.mask
1515 badPix = (mask.array & mask.getPlaneBitMask(self.config.detection.excludeMaskPlanes)) > 0
1516 self.metadata["nGoodPixels"] = np.sum(~badPix)
1517 self.metadata["nBadPixels"] = np.sum(badPix)
1518 detPosPix = (mask.array & mask.getPlaneBitMask("DETECTED")) > 0
1519 detNegPix = (mask.array & mask.getPlaneBitMask("DETECTED_NEGATIVE")) > 0
1520 self.metadata["nPixelsDetectedPositive"] = np.sum(detPosPix)
1521 self.metadata["nPixelsDetectedNegative"] = np.sum(detNegPix)
1522 detPosPix &= badPix
1523 detNegPix &= badPix
1524 self.metadata["nBadPixelsDetectedPositive"] = np.sum(detPosPix)
1525 self.metadata["nBadPixelsDetectedNegative"] = np.sum(detNegPix)
1527 metricsMaskPlanes = list(mask.getMaskPlaneDict().keys())
1528 for maskPlane in metricsMaskPlanes:
1529 try:
1530 self.metadata["%s_mask_fraction"%maskPlane.lower()] = evaluateMaskFraction(mask, maskPlane)
1531 except InvalidParameterError:
1532 self.metadata["%s_mask_fraction"%maskPlane.lower()] = -1
1533 self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane)
1535 if self.config.doSkySources:
1536 skySources = diaSources[diaSources["sky_source"]]
1537 else:
1538 skySources = None
1539 metrics = computeDifferenceImageMetrics(science, difference, kernelSources, sky_sources=skySources)
1541 self.metadata["residualFootprintRatioMean"] = metrics.differenceFootprintRatioMean
1542 self.metadata["residualFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev
1543 self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean
1544 self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev
1545 self.log.info("Mean, stdev of ratio of difference to science "
1546 "pixels in star footprints: %5.4f, %5.4f",
1547 self.metadata["residualFootprintRatioMean"],
1548 self.metadata["residualFootprintRatioStdev"])
1549 if self.config.raiseOnBadSubtractionRatio:
1550 if metrics.differenceFootprintRatioMean > self.config.badSubtractionRatioThreshold: 1550 ↛ 1551line 1550 didn't jump to line 1551 because the condition on line 1550 was never true
1551 raise BadSubtractionError(ratio=metrics.differenceFootprintRatioMean,
1552 threshold=self.config.badSubtractionRatioThreshold)
1553 if metrics.differenceFootprintRatioStdev > self.config.badSubtractionVariationThreshold: 1553 ↛ 1554line 1553 didn't jump to line 1554 because the condition on line 1553 was never true
1554 raise BadSubtractionError(ratio=metrics.differenceFootprintRatioStdev,
1555 threshold=self.config.badSubtractionVariationThreshold)
1557 def getSattleDiaSourceAllowlist(self, diaSources, science):
1558 """Query the sattle service and determine which diaSources are allowed.
1560 Parameters
1561 ----------
1562 diaSources : `lsst.afw.table.SourceCatalog`
1563 The catalog of detected sources.
1564 science : `lsst.afw.image.ExposureF`
1565 Science exposure that was subtracted.
1567 Returns
1568 ----------
1569 allow_list : `list` of `int`
1570 diaSourceIds of diaSources that can be made public.
1572 Raises
1573 ------
1574 requests.HTTPError
1575 Raised if sattle call does not return success.
1576 """
1577 wcs = science.getWcs()
1578 visit_info = science.getInfo().getVisitInfo()
1579 visit_id = visit_info.getId()
1580 sattle_uri_base = os.getenv('SATTLE_URI_BASE')
1582 dia_sources_json = []
1583 for source in diaSources:
1584 source_bbox = source.getFootprint().getBBox()
1585 corners = wcs.pixelToSky([lsst.geom.Point2D(c) for c in source_bbox.getCorners()])
1586 bbox_radec = [[pt.getRa().asDegrees(), pt.getDec().asDegrees()] for pt in corners]
1587 dia_sources_json.append({"diasource_id": source["id"], "bbox": bbox_radec})
1589 payload = {"visit_id": visit_id, "detector_id": science.getDetector().getId(),
1590 "diasources": dia_sources_json, "historical": self.config.sattle_historical}
1592 sattle_output = requests.put(f'{sattle_uri_base}/diasource_allow_list',
1593 json=payload)
1595 # retry once if visit cache is not populated
1596 if sattle_output.status_code == 404:
1597 self.log.warning(f'Visit {visit_id} not found in sattle cache, re-sending')
1598 populate_sattle_visit_cache(visit_info, historical=self.config.sattle_historical)
1599 sattle_output = requests.put(f'{sattle_uri_base}/diasource_allow_list', json=payload)
1601 sattle_output.raise_for_status()
1603 return sattle_output.json()['allow_list']
1605 @timeMethod
1606 def filterSatellites(self, diaSources, science):
1607 """Remove diaSources overlapping predicted satellite positions.
1609 Parameters
1610 ----------
1611 diaSources : `lsst.afw.table.SourceCatalog`
1612 The catalog of detected sources.
1613 science : `lsst.afw.image.ExposureF`
1614 Science exposure that was subtracted.
1616 Returns
1617 ----------
1618 filterdDiaSources : `lsst.afw.table.SourceCatalog`
1619 Filtered catalog of diaSources
1620 """
1622 allow_list = self.getSattleDiaSourceAllowlist(diaSources, science)
1624 if allow_list:
1625 allow_set = set(allow_list)
1626 allowed_ids = [source['id'] in allow_set for source in diaSources]
1627 diaSources = diaSources[np.array(allowed_ids)].copy(deep=True)
1628 else:
1629 self.log.warning('Sattle allowlist is empty, all diaSources removed')
1630 diaSources = diaSources[0:0].copy(deep=True)
1631 return diaSources
1633 def _runStreakMasking(self, difference):
1634 """Do streak masking and optionally save the resulting streak
1635 fit parameters in a catalog.
1637 Only returns non-empty streakInfo if self.config.writeStreakInfo
1638 is set. The difference image is binned by self.config.streakBinFactor
1639 (and detection is run a second time) so that regions with lower
1640 surface brightness streaks are more likely to fall above the
1641 detection threshold.
1643 Parameters
1644 ----------
1645 difference: `lsst.afw.image.Exposure`
1646 The exposure in which to search for streaks. Must have a detection
1647 mask.
1649 Returns
1650 -------
1651 streakInfo: `lsst.pipe.base.Struct`
1652 ``rho`` : `np.ndarray`
1653 Angle of detected streak.
1654 ``theta`` : `np.ndarray`
1655 Distance from center of detected streak.
1656 ``sigma`` : `np.ndarray`
1657 Width of streak profile.
1658 ``reducedChi2`` : `np.ndarray`
1659 Reduced chi2 of the best-fit streak profile.
1660 ``modelMaximum`` : `np.ndarray`
1661 Peak value of the fit line profile.
1662 """
1663 maskedImage = difference.maskedImage
1664 # Bin the diffim to enhance low surface brightness streaks
1665 binnedMaskedImage = afwMath.binImage(maskedImage,
1666 self.config.streakBinFactor,
1667 self.config.streakBinFactor)
1668 binnedExposure = afwImage.ExposureF(binnedMaskedImage.getBBox())
1669 binnedExposure.setMaskedImage(binnedMaskedImage)
1670 # Clear the DETECTED mask plane before streak detection
1671 binnedExposure.mask &= ~binnedExposure.mask.getPlaneBitMask('DETECTED')
1672 # Rerun detection to set the DETECTED mask plane on binnedExposure
1673 sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius()
1674 _table = afwTable.SourceTable.make(afwTable.SourceTable.makeMinimalSchema())
1675 self.streakDetection.run(table=_table, exposure=binnedExposure, doSmooth=True,
1676 sigma=sigma/self.config.streakBinFactor)
1677 binnedDetectedMaskPlane = binnedExposure.mask.array & binnedExposure.mask.getPlaneBitMask('DETECTED')
1678 rescaledDetectedMaskPlane = binnedDetectedMaskPlane.repeat(self.config.streakBinFactor,
1679 axis=0).repeat(self.config.streakBinFactor,
1680 axis=1)
1681 # Create new version of a diffim with DETECTED based on binnedExposure
1682 streakMaskedImage = maskedImage.clone()
1683 ysize, xsize = rescaledDetectedMaskPlane.shape
1684 streakMaskedImage.mask.array[:ysize, :xsize] |= rescaledDetectedMaskPlane
1685 # Detect streaks on this new version of the diffim
1686 streaks = self.maskStreaks.run(streakMaskedImage)
1687 streakMaskPlane = streakMaskedImage.mask.array & streakMaskedImage.mask.getPlaneBitMask('STREAK')
1688 # Apply the new STREAK mask to the original diffim
1689 maskedImage.mask.array |= streakMaskPlane
1691 if self.config.writeStreakInfo: 1691 ↛ 1692line 1691 didn't jump to line 1692 because the condition on line 1691 was never true
1692 rhos = np.array([line.rho for line in streaks.lines])
1693 thetas = np.array([line.theta for line in streaks.lines])
1694 sigmas = np.array([line.sigma for line in streaks.lines])
1695 chi2s = np.array([line.reducedChi2 for line in streaks.lines])
1696 modelMaximums = np.array([line.modelMaximum for line in streaks.lines])
1697 streakInfo = {'rho': rhos, 'theta': thetas, 'sigma': sigmas, 'reducedChi2': chi2s,
1698 'modelMaximum': modelMaximums}
1699 else:
1700 streakInfo = {'rho': np.array([]), 'theta': np.array([]), 'sigma': np.array([]),
1701 'reducedChi2': np.array([]), 'modelMaximum': np.array([])}
1702 return pipeBase.Struct(maskedStreaks=streakInfo)
1705class DetectAndMeasureScoreConnections(DetectAndMeasureConnections):
1706 scoreExposure = pipeBase.connectionTypes.Input(
1707 doc="Maximum likelihood image for detection.",
1708 dimensions=("instrument", "visit", "detector"),
1709 storageClass="ExposureF",
1710 name="{fakesType}{coaddName}Diff_scoreTempExp",
1711 )
1712 scoreMeasuredExposure = pipeBase.connectionTypes.Output(
1713 doc="Score image after backgrond subtraction and detection.",
1714 dimensions=("instrument", "visit", "detector"),
1715 storageClass="ExposureF",
1716 name="{fakesType}{coaddName}Diff_scoreExp",
1717 )
1719 def __init__(self, *, config):
1720 super().__init__(config=config)
1721 if self.config.image_type == "future":
1722 self.scoreExposure = dataclasses.replace(self.scoreExposure, storageClass="DifferenceImage")
1723 self.scoreMeasuredExposure = dataclasses.replace(self.scoreMeasuredExposure,
1724 storageClass="DifferenceImage")
1727class DetectAndMeasureScoreConfig(DetectAndMeasureConfig,
1728 pipelineConnections=DetectAndMeasureScoreConnections):
1729 pass
1732class DetectAndMeasureScoreTask(DetectAndMeasureTask):
1733 """Detect DIA sources using a score image,
1734 and measure the detections on the difference image.
1736 Source detection is run on the supplied score, or maximum likelihood,
1737 image. Note that no additional convolution will be done in this case.
1738 Close positive and negative detections will optionally be merged into
1739 dipole diaSources.
1740 Sky sources, or forced detections in background regions, will optionally
1741 be added, and the configured measurement algorithm will be run on all
1742 detections.
1743 """
1744 ConfigClass = DetectAndMeasureScoreConfig
1745 _DefaultName = "detectAndMeasureScore"
1747 @timeMethod
1748 def run(self, science, matchedTemplate, difference, scoreExposure, kernelSources,
1749 idFactory=None, measurementResults=None, templateInfo=None):
1750 """Detect and measure sources on a score image.
1752 Parameters
1753 ----------
1754 science : `lsst.afw.image.ExposureF`
1755 Science exposure that the template was subtracted from.
1756 matchedTemplate : `lsst.afw.image.ExposureF`
1757 Warped and PSF-matched template that was used produce the
1758 difference image.
1759 difference : `lsst.afw.image.ExposureF`
1760 Result of subtracting template from the science image.
1761 scoreExposure : `lsst.afw.image.ExposureF`
1762 Score or maximum likelihood difference image
1763 kernelSources : `lsst.afw.table.SourceCatalog`
1764 Final selection of sources that was used for psf matching.
1765 idFactory : `lsst.afw.table.IdFactory`, optional
1766 Generator object used to assign ids to detected sources in the
1767 difference image. Ids from this generator are not set until after
1768 deblending and merging positive/negative peaks.
1769 measurementResults : `lsst.pipe.base.Struct`, optional
1770 Result struct that is modified to allow saving of partial outputs
1771 for some failure conditions. If the task completes successfully,
1772 this is also returned.
1773 templateInfo : `list` [`lsst.images.DifferenceImageTemplateInfo`], \
1774 optional
1775 Information about the coadds and their input images that the
1776 template was built from. Without it the template epochs of every
1777 diaSource are NaN.
1779 Returns
1780 -------
1781 measurementResults : `lsst.pipe.base.Struct`
1783 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF`
1784 Subtracted exposure with detection mask applied.
1785 ``scoreMeasuredExposure`` : `lsst.afw.image.ExposureF`
1786 Score exposure with detection mask applied.
1787 ``diaSources`` : `lsst.afw.table.SourceCatalog`
1788 The catalog of detected sources.
1789 ``differenceBackground`` : `lsst.afw.math.BackgroundList`
1790 Background that was subtracted from the difference image.
1791 ``templateEpochs`` : `dict` [`str`, `float`]
1792 Summary of the template epochs of the diaSources, keyed by
1793 image metadata key; empty if they are not known.
1794 """
1795 if measurementResults is None:
1796 measurementResults = pipeBase.Struct()
1797 if idFactory is None:
1798 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
1800 # Check image properties and clear detection mask planes.
1801 self._prepareInputs(difference)
1802 self._prepareInputs(scoreExposure)
1803 if self.config.doSubtractBackground:
1804 background = self._fitAndSubtractBackground(
1805 differenceExposure=difference, scoreExposure=scoreExposure,
1806 )
1807 else:
1808 background = afwMath.BackgroundList()
1810 if self.config.doFindCosmicRays: 1810 ↛ 1819line 1810 didn't jump to line 1819 because the condition on line 1810 was always true
1811 # Cosmic rays are detected on the difference and the resulting
1812 # mask is propagated to the score image used for detection.
1813 self.findAndMaskCosmicRays(difference)
1814 scoreExposure.mask.array |= difference.mask.array
1816 # Don't use the idFactory until after deblend+merge, so that we aren't
1817 # generating ids that just get thrown away (footprint merge doesn't
1818 # know about past ids).
1819 table = afwTable.SourceTable.make(self.schema)
1820 detectResults = self.detection.run(
1821 table=table,
1822 exposure=scoreExposure,
1823 doSmooth=False,
1824 background=background,
1825 clearMask=True,
1826 )
1827 measurementResults.differenceBackground = background
1828 measurementResults.scoreMeasuredExposure = scoreExposure
1830 # Copy the final detection masks from the Score image to the difference.
1831 # Do not copy all mask planes, since the convolution will grow them by
1832 # the size of the convolution kernel.
1833 detectedBit = difference.mask.getPlaneBitMask(["DETECTED"])
1834 detectedPix = scoreExposure.mask.array & detectedBit > 0
1835 detectedNegativeBit = difference.mask.getPlaneBitMask(["DETECTED_NEGATIVE"])
1836 detectedNegativePix = scoreExposure.mask.array & detectedNegativeBit > 0
1837 difference.mask.array[detectedPix] |= detectedBit
1838 difference.mask.array[detectedNegativePix] |= detectedNegativeBit
1839 # Also set the EDGE mask, since the science image was convolved
1840 edgeBit = difference.mask.getPlaneBitMask(["EDGE"])
1841 edgePix = scoreExposure.mask.array & edgeBit > 0
1842 difference.mask.array[edgePix] |= edgeBit
1844 if self.config.doDeblend:
1845 sources, positives, negatives = self._deblend(difference,
1846 detectResults.positive,
1847 detectResults.negative)
1849 self.processResults(science, matchedTemplate, difference,
1850 sources, idFactory, kernelSources,
1851 positives=positives,
1852 negatives=negatives,
1853 measurementResults=measurementResults,
1854 templateInfo=templateInfo)
1856 else:
1857 positives = afwTable.SourceCatalog(self.schema)
1858 detectResults.positive.makeSources(positives)
1859 negatives = afwTable.SourceCatalog(self.schema)
1860 detectResults.negative.makeSources(negatives)
1861 self.processResults(science, matchedTemplate, difference,
1862 detectResults.sources, idFactory, kernelSources,
1863 positives=positives,
1864 negatives=negatives,
1865 measurementResults=measurementResults,
1866 templateInfo=templateInfo)
1867 return measurementResults