Coverage for python / lsst / drp / tasks / assemble_chi2_coadd.py: 0%
105 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 09:15 +0000
1# This file is part of drp_tasks.
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/>.
23import logging
25import numpy as np
27import lsst.afw.image as afwImage
28import lsst.afw.math as afwMath
29import lsst.afw.table as afwTable
30import lsst.pex.config as pexConfig
31import lsst.pipe.base as pipeBase
32import lsst.pipe.base.connectionTypes as cT
33from lsst.afw.detection import Psf
34from lsst.meas.algorithms import SourceDetectionTask
35from lsst.meas.base import SkyMapIdGeneratorConfig
37log = logging.getLogger(__name__)
40def calculateKernelSize(sigma: float, nSigmaForKernel: float = 7) -> int:
41 """Calculate the size of the smoothing kernel.
43 Parameters
44 ----------
45 sigma:
46 Gaussian sigma of smoothing kernel.
47 nSigmaForKernel:
48 The multiple of `sigma` to use to set the size of the kernel.
49 Note that that is the full width of the kernel bounding box
50 (so a value of 7 means 3.5 sigma on either side of center).
51 The value will be rounded up to the nearest odd integer.
53 Returns
54 -------
55 size:
56 Size of the smoothing kernel.
57 """
58 return (int(sigma * nSigmaForKernel + 0.5) // 2) * 2 + 1 # make sure it is odd
61def convolveImage(image: afwImage.Image, psf: Psf) -> afwImage.Image:
62 """Convolve an image with a psf
64 This methodm and the docstring, is based off the method in
65 `~lsst.meas.algorithms.detection.SourceDetectionTask`.
67 We convolve the image with a Gaussian approximation to the PSF,
68 because this is separable and therefore fast. It's technically a
69 correlation rather than a convolution, but since we use a symmetric
70 Gaussian there's no difference.
72 Parameters
73 ----------
74 image:
75 The image to convovle.
76 psf:
77 The PSF to convolve the `image` with.
79 Returns
80 -------
81 convolved:
82 The result of convolving `image` with the `psf`.
83 """
84 sigma = psf.computeShape(psf.getAveragePosition()).getDeterminantRadius()
85 bbox = image.getBBox()
87 # Smooth using a Gaussian (which is separable, hence fast) of width sigma
88 # Make a SingleGaussian (separable) kernel with the 'sigma'
89 kWidth = calculateKernelSize(sigma)
90 gaussFunc = afwMath.GaussianFunction1D(sigma)
91 gaussKernel = afwMath.SeparableKernel(kWidth, kWidth, gaussFunc, gaussFunc)
93 convolvedImage = image.Factory(bbox)
95 afwMath.convolve(convolvedImage, image, gaussKernel, afwMath.ConvolutionControl())
97 return convolvedImage.Factory(convolvedImage, bbox, afwImage.PARENT, False)
100class AssembleChi2CoaddConnections(
101 pipeBase.PipelineTaskConnections,
102 dimensions=("tract", "patch", "skymap"),
103 defaultTemplates={"inputCoaddName": "deep", "outputCoaddName": "deepChi2"},
104):
105 inputCoadds = cT.Input(
106 doc="Exposure on which to run deblending",
107 name="{inputCoaddName}Coadd_calexp",
108 storageClass="ExposureF",
109 multiple=True,
110 dimensions=("tract", "patch", "band", "skymap"),
111 )
112 chi2Coadd = cT.Output(
113 doc="Chi^2 exposure, produced by merging multiband coadds",
114 name="{outputCoaddName}Coadd_calexp",
115 storageClass="ExposureF",
116 dimensions=("tract", "patch", "skymap"),
117 )
120class AssembleChi2CoaddConfig(pipeBase.PipelineTaskConfig, pipelineConnections=AssembleChi2CoaddConnections):
121 outputPixelatedVariance = pexConfig.Field(
122 dtype=bool,
123 default=False,
124 doc="Whether to output a pixelated variance map for the generated "
125 "chi^2 coadd, or to have a flat variance map defined by combining "
126 "the inverse variance maps of the coadds that were combined.",
127 )
129 useUnionForMask = pexConfig.Field(
130 dtype=bool,
131 default=True,
132 doc="Whether to calculate the union of the mask plane in each band, "
133 "or the intersection of the mask plane in each band.",
134 )
137class AssembleChi2CoaddTask(pipeBase.PipelineTask):
138 """Assemble a chi^2 coadd from a collection of multi-band coadds
140 References
141 ----------
142 .. [1] Szalay, A. S., Connolly, A. J., and Szokoly, G. P., “Simultaneous
143 Multicolor Detection of Faint Galaxies in the Hubble Deep Field”,
144 The Astronomical Journal, vol. 117, no. 1, pp. 68–74,
145 1999. doi:10.1086/300689.
147 .. [2] Kaiser 2001 whitepaper,
148 http://pan-starrs.ifa.hawaii.edu/project/people/kaiser/imageprocessing/im%2B%2B.pdf # noqa: E501, W505
150 .. [3] https://dmtn-015.lsst.io/
152 .. [4] https://project.lsst.org/meetings/law/sites/lsst.org.meetings.law/files/Building%20and%20using%20coadds.pdf # noqa: E501, W505
153 """
155 ConfigClass = AssembleChi2CoaddConfig
156 _DefaultName = "assembleChi2Coadd"
158 def __init__(self, initInputs, **kwargs):
159 super().__init__(initInputs=initInputs, **kwargs)
161 def combinedMasks(self, masks: list[afwImage.MaskX]) -> afwImage.MaskX:
162 """Combine the mask plane in each input coadd
164 Parameters
165 ----------
166 mMask:
167 The MultibandMask in each band.
169 Returns
170 -------
171 result:
172 The resulting single band mask.
173 """
174 refMask = masks[0]
175 bbox = refMask.getBBox()
176 mask = refMask.array
177 for _mask in masks[1:]:
178 if self.config.useUnionForMask:
179 mask = mask | _mask.array
180 else:
181 mask = mask & _mask.array
182 result = refMask.Factory(bbox)
183 result.array[:] = mask
184 return result
186 def runQuantum(self, butlerQC, inputRefs, outputRefs):
187 inputs = butlerQC.get(inputRefs)
188 outputs = self.run(**inputs)
189 butlerQC.put(outputs, outputRefs)
191 def run(self, inputCoadds: list[afwImage.Exposure]) -> pipeBase.Struct:
192 """Assemble the chi2 coadd from the multiband coadds
194 Parameters
195 ----------
196 inputCoadds:
197 The coadds to combine into a single chi2 coadd.
199 Returns
200 -------
201 result:
202 The chi2 coadd created from the input coadds.
203 """
204 convControl = afwMath.ConvolutionControl()
205 convControl.setDoNormalize(False)
206 convControl.setDoCopyEdge(False)
208 # Set a reference exposure to use for creating the new coadd.
209 # It doesn't matter which exposure we use, since we just need the
210 # bounding box information and Factory to create a new expsure with
211 # the same dtype.
212 refExp = inputCoadds[0]
213 bbox = refExp.getBBox()
215 image = refExp.image.Factory(bbox)
216 variance_list = []
217 # Convovle the image in each band and weight by the median variance
218 for calexp in inputCoadds:
219 convolved = convolveImage(calexp.image, calexp.getPsf())
220 _variance = np.median(calexp.variance.array)
221 convolved.array[:] /= _variance
222 image += convolved
223 variance_list.append(_variance)
225 variance = refExp.variance.Factory(bbox)
226 if self.config.outputPixelatedVariance:
227 # Write the per pixel variance to the output coadd
228 variance.array[:] = np.sum([1 / coadd.variance for coadd in inputCoadds], axis=0)
229 else:
230 # Use a flat variance in each band
231 variance.array[:] = np.sum(1 / np.array(variance_list))
232 # Combine the masks planes to calculate the mask plae of the new coadd
233 mask = self.combinedMasks([coadd.mask for coadd in inputCoadds])
234 # Create the exposure
235 maskedImage = refExp.maskedImage.Factory(image, mask=mask, variance=variance)
236 chi2coadd = refExp.Factory(maskedImage, exposureInfo=refExp.getInfo())
237 chi2coadd.info.setFilter(None)
238 return pipeBase.Struct(chi2Coadd=chi2coadd)
241class DetectChi2SourcesConnections(
242 pipeBase.PipelineTaskConnections,
243 dimensions=("tract", "patch", "skymap"),
244 defaultTemplates={"inputCoaddName": "deepChi2", "outputCoaddName": "deepChi2"},
245):
246 detectionSchema = cT.InitOutput(
247 doc="Schema of the detection catalog",
248 name="{outputCoaddName}Coadd_det_schema",
249 storageClass="SourceCatalog",
250 )
251 exposure = cT.Input(
252 doc="Exposure on which detections are to be performed",
253 name="{inputCoaddName}Coadd_calexp",
254 storageClass="ExposureF",
255 dimensions=("tract", "patch", "skymap"),
256 )
257 outputSources = cT.Output(
258 doc="Detected sources catalog",
259 name="{outputCoaddName}Coadd_det",
260 storageClass="SourceCatalog",
261 dimensions=("tract", "patch", "skymap"),
262 )
265class DetectChi2SourcesConfig(pipeBase.PipelineTaskConfig, pipelineConnections=DetectChi2SourcesConnections):
266 detection = pexConfig.ConfigurableField(target=SourceDetectionTask, doc="Detect sources in chi2 coadd")
268 idGenerator = SkyMapIdGeneratorConfig.make_field()
270 def setDefaults(self):
271 super().setDefaults()
272 self.detection.reEstimateBackground = False
273 self.detection.thresholdValue = 3
276class DetectChi2SourcesTask(pipeBase.PipelineTask):
277 _DefaultName = "detectChi2Sources"
278 ConfigClass = DetectChi2SourcesConfig
280 def __init__(self, schema=None, **kwargs):
281 # N.B. Super is used here to handle the multiple inheritance of
282 # PipelineTasks, the init tree call structure has been reviewed
283 # carefully to be sure super will work as intended.
284 super().__init__(**kwargs)
285 if schema is None:
286 schema = afwTable.SourceTable.makeMinimalSchema()
287 self.schema = schema
288 self.makeSubtask("detection", schema=self.schema)
289 self.detectionSchema = afwTable.SourceCatalog(self.schema)
291 def runQuantum(self, butlerQC, inputRefs, outputRefs):
292 inputs = butlerQC.get(inputRefs)
293 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId)
294 inputs["idFactory"] = idGenerator.make_table_id_factory()
295 inputs["expId"] = idGenerator.catalog_id
296 outputs = self.run(**inputs)
297 butlerQC.put(outputs, outputRefs)
299 def run(self, exposure: afwImage.Exposure, idFactory: afwTable.IdFactory, expId: int) -> pipeBase.Struct:
300 """Run detection on a chi2 exposure.
302 Parameters
303 ----------
304 exposure :
305 Exposure on which to detect (maybe backround-subtracted and scaled,
306 depending on configuration).
307 idFactory :
308 IdFactory to set source identifiers.
309 expId :
310 Exposure identifier (integer) for RNG seed.
312 Returns
313 -------
314 result : `lsst.pipe.base.Struct`
315 Results as a struct with attributes:
316 ``outputSources``
317 Catalog of detections (`lsst.afw.table.SourceCatalog`).
318 """
319 table = afwTable.SourceTable.make(self.schema, idFactory)
320 # We override `doSmooth` since the chi2 coadd has already had an
321 # extra PSF convolution applied to decorrelate the images
322 # accross bands.
323 detections = self.detection.run(table, exposure, expId=expId, doSmooth=False)
324 sources = detections.sources
325 return pipeBase.Struct(outputSources=sources)