Coverage for python / lsst / pipe / tasks / tractBackground.py: 0%
135 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-15 00:08 +0000
1# This file is part of pipe_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/>.
22__all__ = [
23 "TractBackground",
24 "TractBackgroundConfig",
25]
27import itertools
29import numpy
31import lsst.afw.geom as afwGeom
32import lsst.afw.image as afwImage
33import lsst.afw.math as afwMath
34import lsst.geom as geom
35from lsst.pex.config import ChoiceField, Config, Field, ListField
36from lsst.pipe.tasks.background import interpolateBadPixels, smoothArray
39class TractBackgroundConfig(Config):
40 """Configuration for TractBackground
42 Parameters `xBin` and `yBin` are in pixels, as translation from warps to
43 tract and back only requires geometric transformations in the warped pixel
44 plane.
45 """
47 xBin = Field(dtype=float, default=128, doc="Bin size in x")
48 yBin = Field(dtype=float, default=128, doc="Bin size in y")
49 minFrac = Field(dtype=float, default=0.1, doc="Minimum fraction of bin size for good measurement")
50 # TODO: not masking detections could impact reference visit selection.
51 # Either could make the field configurable, or swap to pre-made references
52 mask = ListField(
53 dtype=str,
54 doc="Mask planes to treat as bad",
55 default=["BAD", "SAT", "EDGE", "NO_DATA"],
56 )
57 interpolation = ChoiceField(
58 doc="How to interpolate the background values. This maps to an enum; see afw::math::Background",
59 dtype=str,
60 default="AKIMA_SPLINE",
61 optional=True,
62 allowed={
63 "CONSTANT": "Use a single constant value",
64 "LINEAR": "Use linear interpolation",
65 "NATURAL_SPLINE": "Cubic spline with zero second derivative at endpoints",
66 "AKIMA_SPLINE": "Higher-level nonlinear spline that is more robust to outliers",
67 "NONE": "No background estimation is to be attempted",
68 },
69 )
70 doSmooth = Field(dtype=bool, default=False, doc="Do smoothing?")
71 smoothScale = Field(dtype=float, default=2.0, doc="Smoothing scale, as a multiple of the bin size")
72 binning = Field(dtype=int, default=64, doc="Binning to use for warp background model (pixels)")
75class TractBackground:
76 """Background model for a tract, comprised of warped exposures.
78 Similar to lsst.pipe.background.FocalPlaneBackground class, in that we
79 model the background using the "superpixel" method, measuring the
80 background in each superpixel and interpolating between them.
82 There is one use pattern for building a background model: create a
83 `TractBackground`, then `addWarp` for each of the patches in a tract.
85 Once you've built the background model, you can apply it to individual
86 warps with the `toWarpBackground` method.
88 Parameters
89 ----------
90 config : `TractBackgroundConfig`
91 Configuration for measuring tract backgrounds.
92 skymap : `lsst.skymap.ringsSkyMap.RingsSkyMap`
93 Skymap object
94 tract : `int`
95 Working Tract ID.
96 transform : `lsst.afw.geom.TransformPoint2ToPoint2`
97 Transformation from tract coordinates to warp coordinates.
98 values : `lsst.afw.image.ImageF`
99 Measured background values.
100 numbers : `lsst.afw.image.ImageF`
101 Number of pixels in each background measurement.
102 variances : `lsst.afw.image.ImageF`
103 Measured background variances.
104 """
106 def __init__(self, config, skymap, tract, values=None, numbers=None, variances=None):
107 # Developers should note that changes to the signature of this method
108 # require coordinated changes to the `__reduce__` and `clone` methods.
109 self.config = config
110 self.skymap = skymap
111 self.tract = tract
112 self.tractInfo = self.skymap[tract]
113 tractDimX, tractDimY = self.tractInfo.getBBox().getDimensions()
114 self.dims = geom.Extent2I(tractDimX / self.config.xBin, tractDimY / self.config.yBin)
116 if values is None:
117 values = afwImage.ImageF(self.dims)
118 values.set(0.0)
119 else:
120 values = values.clone()
121 assert values.getDimensions() == self.dims
122 self._values = values
123 if numbers is None:
124 numbers = afwImage.ImageF(self.dims) # float for dynamic range and convenience
125 numbers.set(0.0)
126 else:
127 numbers = numbers.clone()
128 assert numbers.getDimensions() == self.dims
129 self._numbers = numbers
130 if variances is None:
131 variances = afwImage.ImageF(self.dims)
132 variances.set(0.0)
133 else:
134 variances = variances.clone()
135 assert variances.getDimensions() == self.dims
136 self._variances = variances
138 def __reduce__(self):
139 return self.__class__, (
140 self.config,
141 self.skymap,
142 self.tract,
143 self._values,
144 self._numbers,
145 self._variances,
146 )
148 def clone(self):
149 return self.__class__(
150 self.config, self.skymap, self.tract, self._values, self._numbers, self._variances
151 )
153 def addWarp(self, warp):
154 """Bins masked images of warps and adds these values into a blank image
155 with the binned tract dimensions at the location of the warp in the
156 tract.
158 Parameters
159 ----------
160 warp : `lsst.afw.image.ExposureF`
161 Warped image corresponding to a single patch in a single visit.
162 """
163 image = warp.getMaskedImage()
164 maskVal = image.getMask().getPlaneBitMask(self.config.mask)
166 warped = afwImage.ImageF(self._values.getBBox())
167 warpedCounts = afwImage.ImageF(self._numbers.getBBox())
168 warpedVariances = afwImage.ImageF(self._variances.getBBox())
170 stats = afwMath.StatisticsControl()
171 stats.setAndMask(maskVal)
172 stats.setNanSafe(True)
174 # Pixel locations in binned tract-scale image
175 pixels = itertools.product(
176 warped.getBBox().x.arange(),
177 warped.getBBox().y.arange(),
178 )
179 for xx, yy in pixels:
180 llc = geom.Point2D((xx - 0.5) * self.config.xBin, (yy - 0.5) * self.config.yBin)
181 urc = geom.Point2D(
182 (xx + 0.5) * self.config.xBin + self.config.xBin - 1,
183 (yy + 0.5) * self.config.yBin + self.config.yBin - 1,
184 )
185 bbox = geom.Box2I(geom.Point2I(llc), geom.Point2I(urc))
186 bbox.clip(image.getBBox()) # Works in tract coordinates
187 if bbox.isEmpty():
188 continue
189 subImage = image.Factory(image, bbox)
190 result = afwMath.makeStatistics(
191 subImage, afwMath.MEANCLIP | afwMath.NPOINT | afwMath.VARIANCECLIP, stats
192 )
193 mean = result.getValue(afwMath.MEANCLIP)
194 num = result.getValue(afwMath.NPOINT)
195 var = result.getValue(afwMath.VARIANCECLIP)
196 if not numpy.isfinite(mean) or not numpy.isfinite(num):
197 continue
198 warped[xx, yy, afwImage.LOCAL] = mean * num
199 warpedCounts[xx, yy, afwImage.LOCAL] = num
200 warpedVariances[xx, yy, afwImage.LOCAL] = var * num
202 self._values += warped
203 self._numbers += warpedCounts
204 self._variances += warpedVariances
206 def __iadd__(self, other):
207 """Merges with another TractBackground
209 This allows multiple background models to be constructed from
210 different warps, and then merged to form a single consistent
211 background model for the entire tract.
213 Parameters
214 ----------
215 other : `TractBackground`
216 Another background model to merge.
218 Returns
219 -------
220 self : `TractBackground`
221 The merged background model.
222 """
223 if (self.config.xBin, self.config.yBin) != (other.config.xBin, other.config.yBin):
224 raise RuntimeError(
225 "Size mismatch: %s vs %s"
226 % ((self.config.xBin, self.config.yBin), (other.config.xBin, other.config.yBin))
227 )
228 if self.dims != other.dims:
229 raise RuntimeError("Dimensions mismatch: %s vs %s" % (self.dims, other.dims))
230 self._values += other._values
231 self._numbers += other._numbers
232 self._variances += other._variances
233 return self
235 def toWarpBackground(self, warp):
236 """Produce a background model for a warp
238 The superpixel model is transformed back to the native pixel
239 resolution, for application to the background of an individual warp.
241 Parameters
242 ----------
243 warp : `lsst.afw.image.ExposureF`
244 Warped image corresponding to a single patch in a single visit.
246 Returns
247 -------
248 bg : `lsst.afw.math.BackgroundList`
249 Background model for warp.
250 """
251 # Transform to binned warp plane
252 binTransform = geom.AffineTransform.makeScaling(self.config.binning)
254 # Transform from binned tract plane to tract plane
255 # Start at the patch corner, not the warp corner overlap region
256 wcs = self.tractInfo.getWcs()
257 coo = wcs.pixelToSky(1, 1)
258 ptch = self.tractInfo.findPatch(coo)
259 ptchDimX, ptchDimY = ptch.getInnerBBox().getDimensions()
260 if ptchDimX != ptchDimY:
261 raise ValueError(
262 "Patch dimensions %d,%d are unequal: cannot proceed as written." % (ptchDimX, ptchDimY)
263 )
264 overlap = (warp.getBBox().getDimensions().getX() - ptchDimX) // 2
265 corner = warp.getBBox().getMin()
266 if corner[0] % ptchDimX != 0:
267 corner[0] += overlap
268 corner[1] += overlap
269 offset = geom.Extent2D(corner[0], corner[1])
270 tractTransform = (
271 geom.AffineTransform.makeTranslation(geom.Extent2D(-0.5, -0.5))
272 * geom.AffineTransform.makeScaling(1.0 / self.config.xBin, 1.0 / self.config.yBin)
273 * geom.AffineTransform.makeTranslation(offset)
274 )
275 transform = afwGeom.makeTransform(tractTransform)
277 # Full transform
278 toSample = afwGeom.makeTransform(binTransform).then(transform)
280 # Full tract sky model and normalization array
281 tractPlane = self.getStatsImage()
282 tpNorm = afwImage.ImageF(tractPlane.getBBox())
283 tpNorm.set(1.0)
285 # Binned warp image and normalization array
286 bbox = warp.getBBox()
287 image = afwImage.ImageF(bbox.getDimensions() // self.config.binning)
288 norm = afwImage.ImageF(image.getBBox())
290 # Warps full tract model to warp image scale
291 ctrl = afwMath.WarpingControl("bilinear")
292 afwMath.warpImage(image, tractPlane, toSample.inverted(), ctrl)
293 afwMath.warpImage(norm, tpNorm, toSample.inverted(), ctrl)
294 image /= norm
296 # Only sky background model, so include only null values in mask plane
297 mask = afwImage.Mask(image.getBBox())
298 isBad = numpy.isnan(image.getArray())
299 mask.getArray()[isBad] = mask.getPlaneBitMask("BAD")
300 image.getArray()[isBad] = image.getArray()[~isBad].mean()
302 return afwMath.BackgroundList(
303 (
304 afwMath.BackgroundMI(warp.getBBox(), afwImage.makeMaskedImage(image, mask)),
305 afwMath.stringToInterpStyle(self.config.interpolation),
306 afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"),
307 afwMath.ApproximateControl.UNKNOWN,
308 0,
309 0,
310 False,
311 )
312 )
314 def getStatsImage(self):
315 """Return the background model data
317 This is the measurement of the background for each of the superpixels.
318 """
319 values = self._values.clone()
320 values /= self._numbers
322 thresh = self.config.minFrac * (self.config.xBin) * (self.config.yBin)
323 isBad = self._numbers.getArray() < thresh
324 if self.config.doSmooth:
325 array = values.getArray()
326 array[:] = smoothArray(array, isBad, self.config.smoothScale)
327 isBad = numpy.isnan(values.array)
328 # Note: this extrapolates outside the focal plane to the tract edges
329 if numpy.any(isBad):
330 interpolateBadPixels(values.getArray(), isBad, self.config.interpolation)
332 return values
334 def getVarianceImage(self):
335 """Return a variance image with the background model dimensions
337 Notes
338 -----
339 Does not interpolate or extrapolate over bad values.
340 """
341 variances = self._variances.clone()
342 variances /= self._numbers
344 return variances