Coverage for tests / test_addToCoadd.py: 19%
170 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-28 08:45 +0000
1#
2# LSST Data Management System
3# Copyright 2008, 2009, 2010 LSST Corporation.
4#
5# This product includes software developed by the
6# LSST Project (http://www.lsst.org/).
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <http://www.lsstcorp.org/LegalNotices/>.
21#
23"""Test lsst.coadd.utils.addToCoadd
24"""
25import os
26import unittest
28import numpy as np
30import lsst.utils
31import lsst.utils.tests
32import lsst.geom as geom
33import lsst.afw.image as afwImage
34import lsst.afw.math as afwMath
35import lsst.afw.display.ds9 as ds9
36import lsst.pex.exceptions as pexExcept
37import lsst.coadd.utils as coaddUtils
39try:
40 display
41except NameError:
42 display = False
44try:
45 AfwdataDir = lsst.utils.getPackageDir('afwdata')
46except Exception:
47 AfwdataDir = None
48# path to a medium-sized MaskedImage, relative to afwdata package root
49MedMiSubpath = os.path.join("data", "med.fits")
52def slicesFromBox(box, image):
53 """Computes the numpy slice in x and y associated with a parent bounding box
54 given an image/maskedImage/exposure
55 """
56 startInd = (box.getMinX() - image.getX0(), box.getMinY() - image.getY0())
57 stopInd = (startInd[0] + box.getWidth(), startInd[1] + box.getHeight())
58# print "slicesFromBox: box=(min=%s, dim=%s), imxy0=%s, imdim=%s, startInd=%s, stopInd=%s" %\
59# (box.getMin(), box.getDimensions(), image.getXY0(), image.getDimensions(), startInd, stopInd)
60 return (
61 slice(startInd[0], stopInd[0]),
62 slice(startInd[1], stopInd[1]),
63 )
66def referenceAddToCoadd(coadd, weightMap, maskedImage, badPixelMask, weight):
67 """Reference implementation of lsst.coadd.utils.addToCoadd
69 Unlike lsst.coadd.utils.addToCoadd this one does not update the input coadd and weightMap,
70 but instead returns the new versions (as numpy arrays).
72 Inputs:
73 - coadd: coadd before adding maskedImage (a MaskedImage)
74 - weightMap: weight map before adding maskedImage (an Image)
75 - maskedImage: masked image to add to coadd (a MaskedImage)
76 - badPixelMask: mask of bad pixels to ignore (an int)
77 - weight: relative weight of this maskedImage (a float)
79 Returns three items:
80 - overlapBBox: the overlap region relative to the parent (geom.Box2I)
81 - coaddArrayList: new coadd as a list of image, mask, variance numpy arrays
82 - weightMapArray: new weight map, as a numpy array
83 """
84 overlapBBox = coadd.getBBox()
85 overlapBBox.clip(maskedImage.getBBox())
87 coaddArrayList = coadd.image.array.copy(), coadd.mask.array.copy(), coadd.variance.array.copy()
88 weightMapArray = weightMap.array.copy()
90 if overlapBBox.isEmpty():
91 return (overlapBBox, coaddArrayList, weightMapArray)
93 maskedImageArrayList = maskedImage.image.array, maskedImage.mask.array, maskedImage.variance.array
94 badMaskArr = (maskedImageArrayList[1] & badPixelMask) != 0
96 coaddSlices = slicesFromBox(overlapBBox, coadd)
97 imageSlices = slicesFromBox(overlapBBox, maskedImage)
99 badMaskView = badMaskArr[imageSlices[1], imageSlices[0]]
100 for ind in range(3):
101 coaddView = coaddArrayList[ind][coaddSlices[1], coaddSlices[0]]
102 maskedImageView = maskedImageArrayList[ind][imageSlices[1], imageSlices[0]]
103 if ind == 1: # mask plane
104 coaddView |= np.where(badMaskView, 0, maskedImageView)
105 else: # image or variance plane
106 if ind == 0: # image
107 weightFac = weight
108 else: # variance
109 weightFac = weight**2
110 coaddView += np.where(badMaskView, 0, maskedImageView)*weightFac
111 weightMapView = weightMapArray[coaddSlices[1], coaddSlices[0]]
112 weightMapView += np.where(badMaskView, 0, 1)*weight
113 return overlapBBox, coaddArrayList, weightMapArray
116class AddToCoaddTestCase(unittest.TestCase):
117 """A test case for addToCoadd
118 """
120 def _testAddToCoaddImpl(self, useMask, uniformWeight=True):
121 """Test coadd"""
123 trueImageValue = 10.0
124 imBBox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(10, 20))
125 if useMask:
126 coadd = afwImage.MaskedImageF(imBBox)
127 weightMap = coadd.getImage().Factory(coadd.getBBox())
129 badBits = 0x1
130 badPixel = (float("NaN"), badBits, 0)
131 truth = (trueImageValue, 0x0, 0)
132 else:
133 coadd = afwImage.ImageF(imBBox)
134 weightMap = coadd.Factory(coadd.getBBox())
136 badPixel = float("NaN")
137 truth = trueImageValue
139 for i in range(0, 20, 3):
140 image = coadd.Factory(coadd.getDimensions())
141 image.set(badPixel)
143 subBBox = geom.Box2I(geom.Point2I(0, i),
144 image.getDimensions() - geom.Extent2I(0, i))
145 subImage = image.Factory(image, subBBox, afwImage.LOCAL)
146 subImage.set(truth)
147 del subImage
149 weight = 1.0 if uniformWeight else 1.0 + 0.1*i
150 if useMask:
151 coaddUtils.addToCoadd(coadd, weightMap, image, badBits, weight)
152 else:
153 coaddUtils.addToCoadd(coadd, weightMap, image, weight)
155 self.assertEqual(image[-1, -1, afwImage.LOCAL], truth)
157 coadd /= weightMap
159 if display:
160 ds9.mtv(image, title="image", frame=1)
161 ds9.mtv(coadd, title="coadd", frame=2)
162 ds9.mtv(weightMap, title="weightMap", frame=3)
164 stats = afwMath.makeStatistics(coadd, afwMath.MEAN | afwMath.STDEV)
166 return [trueImageValue, stats.getValue(afwMath.MEAN), 0.0, stats.getValue(afwMath.STDEV)]
168 def testAddToCoaddMask(self):
169 """Test coadded MaskedImages"""
171 for uniformWeight in (False, True):
172 truth_mean, mean, truth_stdev, stdev = self._testAddToCoaddImpl(True, uniformWeight)
174 self.assertEqual(truth_mean, mean)
175 self.assertEqual(truth_stdev, stdev)
177 def testAddToCoaddNaN(self):
178 """Test coadded Images with NaN"""
180 for uniformWeight in (False, True):
181 truth_mean, mean, truth_stdev, stdev = self._testAddToCoaddImpl(False, uniformWeight)
183 self.assertEqual(truth_mean, mean)
184 self.assertEqual(truth_stdev, stdev)
187class AddToCoaddAfwdataTestCase(unittest.TestCase):
188 """A test case for addToCoadd using afwdata
189 """
191 def referenceTest(self, coadd, weightMap, image, badPixelMask, weight):
192 """Compare lsst implemenation of addToCoadd to a reference implementation.
194 Returns the overlap bounding box
195 """
196 # this call leaves coadd and weightMap alone:
197 overlapBBox, refCoaddArrayList, refweightMapArray = \
198 referenceAddToCoadd(coadd, weightMap, image, badPixelMask, weight)
199 # this updates coadd and weightMap:
200 afwOverlapBox = coaddUtils.addToCoadd(coadd, weightMap, image, badPixelMask, weight)
201 self.assertEqual(overlapBBox, afwOverlapBox)
203 coaddArrayList = coadd.image.array, coadd.mask.array, coadd.variance.array
204 weightMapArray = weightMap.array
206 for name, ind in (("image", 0), ("mask", 1), ("variance", 2)):
207 if not np.allclose(coaddArrayList[ind], refCoaddArrayList[ind]):
208 errMsgList = (
209 "Computed %s does not match reference for badPixelMask=%s:" % (name, badPixelMask),
210 "computed= %s" % (coaddArrayList[ind],),
211 "reference= %s" % (refCoaddArrayList[ind],),
212 )
213 errMsg = "\n".join(errMsgList)
214 self.fail(errMsg)
215 if not np.allclose(weightMapArray, refweightMapArray):
216 errMsgList = (
217 "Computed weight map does not match reference for badPixelMask=%s:" % (badPixelMask,),
218 "computed= %s" % (weightMapArray,),
219 "reference= %s" % (refweightMapArray,),
220 )
221 errMsg = "\n".join(errMsgList)
222 self.fail(errMsg)
223 return overlapBBox
225 @unittest.skipUnless(AfwdataDir, "afwdata not available")
226 def testMed(self):
227 """Test addToCoadd by adding an image with known bad pixels using varying masks
228 """
229 medBBox = geom.Box2I(geom.Point2I(130, 315), geom.Extent2I(20, 21))
230 medMIPath = os.path.join(AfwdataDir, MedMiSubpath)
231 maskedImage = afwImage.MaskedImageF(afwImage.MaskedImageF(medMIPath), medBBox)
232 coadd = afwImage.MaskedImageF(medBBox)
233 weightMap = afwImage.ImageF(medBBox)
234 weight = 0.9
235 for badPixelMask in (0x00, 0xFF):
236 self.referenceTest(coadd, weightMap, maskedImage, badPixelMask, weight)
238 @unittest.skipUnless(AfwdataDir, "afwdata not available")
239 def testMultSizes(self):
240 """Test addToCoadd by adding various subregions of the med image
241 to a coadd that's a slightly different shape
242 """
243 bbox = geom.Box2I(geom.Point2I(130, 315), geom.Extent2I(30, 31))
244 medMIPath = os.path.join(AfwdataDir, MedMiSubpath)
245 fullMaskedImage = afwImage.MaskedImageF(medMIPath)
246 maskedImage = afwImage.MaskedImageF(fullMaskedImage, bbox)
247 coaddBBox = geom.Box2I(
248 maskedImage.getXY0() + geom.Extent2I(-6, +4),
249 maskedImage.getDimensions() + geom.Extent2I(10, -10))
250 coadd = afwImage.MaskedImageF(coaddBBox)
251 weightMap = afwImage.ImageF(coaddBBox)
252 badPixelMask = 0x0
254 # add masked image that extends beyond coadd in y
255 overlapBBox = self.referenceTest(coadd, weightMap, maskedImage, badPixelMask, 0.5)
256 self.assertFalse(overlapBBox.isEmpty())
258 # add masked image that extends beyond coadd in x
259 bbox = geom.Box2I(geom.Point2I(120, 320), geom.Extent2I(50, 10))
260 maskedImage = afwImage.MaskedImageF(fullMaskedImage, bbox)
261 overlapBBox = self.referenceTest(coadd, weightMap, maskedImage, badPixelMask, 0.5)
262 self.assertFalse(overlapBBox.isEmpty())
264 # add masked image that is fully within the coadd
265 bbox = geom.Box2I(geom.Point2I(130, 320), geom.Extent2I(10, 10))
266 maskedImage = afwImage.MaskedImageF(fullMaskedImage, bbox)
267 overlapBBox = self.referenceTest(coadd, weightMap, maskedImage, badPixelMask, 0.5)
268 self.assertFalse(overlapBBox.isEmpty())
270 # add masked image that does not overlap coadd
271 bbox = geom.Box2I(geom.Point2I(0, 0), geom.Extent2I(10, 10))
272 maskedImage = afwImage.MaskedImageF(fullMaskedImage, bbox)
273 overlapBBox = self.referenceTest(coadd, weightMap, maskedImage, badPixelMask, 0.5)
274 self.assertTrue(overlapBBox.isEmpty())
276 def testAssertions(self):
277 """Test that addToCoadd requires coadd and weightMap to have the same dimensions and xy0"""
278 maskedImage = afwImage.MaskedImageF(geom.Extent2I(10, 10))
279 coadd = afwImage.MaskedImageF(geom.Extent2I(11, 11))
280 coadd.setXY0(5, 6)
281 for dw, dh in (1, 0), (0, 1), (-1, 0), (0, -1):
282 weightMapBBox = geom.Box2I(coadd.getXY0(), coadd.getDimensions() + geom.Extent2I(dw, dh))
283 weightMap = afwImage.ImageF(weightMapBBox)
284 weightMap.setXY0(coadd.getXY0())
285 try:
286 coaddUtils.addToCoadd(coadd, weightMap, maskedImage, 0x0, 0.1)
287 self.fail("should have raised exception")
288 except pexExcept.Exception:
289 pass
290 for dx0, dy0 in (1, 0), (0, 1), (-1, 0), (0, -1):
291 weightMapBBox = geom.Box2I(coadd.getXY0() + geom.Extent2I(dx0, dy0), coadd.getDimensions())
292 weightMap = afwImage.ImageF(weightMapBBox)
293 try:
294 coaddUtils.addToCoadd(coadd, weightMap, maskedImage, 0x0, 0.1)
295 self.fail("should have raised exception")
296 except pexExcept.Exception:
297 pass
300class MemoryTester(lsst.utils.tests.MemoryTestCase):
301 pass
304def setup_module(module):
305 lsst.utils.tests.init()
308if __name__ == "__main__": 308 ↛ 309line 308 didn't jump to line 309 because the condition on line 308 was never true
309 lsst.utils.tests.init()
310 unittest.main()