Coverage for tests / test_coaddBoundedField.py: 16%
137 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-14 23:55 +0000
1#
2# LSST Data Management System
3# Copyright 2008-2017 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
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#
22import os
23import unittest
25import numpy as np
27import lsst.geom
28import lsst.afw.math
29import lsst.afw.geom
30import lsst.afw.image
31import lsst.meas.algorithms
32import lsst.pex.exceptions
33import lsst.utils.tests
36class CoaddBoundedFieldTestCase(lsst.utils.tests.TestCase):
38 def makeRandomWcs(self, crval, maxOffset=10.0):
39 """Make a random TAN Wcs that's complex enough to create an interesting test, by combining
40 random offsets and rotations. We don't include any random scaling, because we don't want
41 the Jacobian to be significantly different from one in these tests, as we want to compare
42 warped images (which implicitly include the Jacobian) against the behavior of CoaddBoundedField
43 (which intentionally does not).
44 """
45 crpix = lsst.geom.Point2D(*self.rng.uniform(low=-maxOffset, high=maxOffset, size=2))
46 scale = 0.01*lsst.geom.arcseconds
47 orientation = np.pi*self.rng.random()*lsst.geom.radians
48 cdMatrix = lsst.afw.geom.makeCdMatrix(scale=scale, orientation=orientation)
49 return lsst.afw.geom.makeSkyWcs(crpix=crpix, crval=crval, cdMatrix=cdMatrix)
51 def makeRandomField(self, bbox):
52 """Create a random ChebyshevBoundedField"""
53 coefficients = self.rng.standard_normal((3, 3))
54 # We add a positive DC offset to make sure our fields more likely to combine constructively;
55 # this makes the test more stringent, because we get less accidental small numbers.
56 coefficients[0, 0] = self.rng.uniform(low=4, high=6)
57 return lsst.afw.math.ChebyshevBoundedField(bbox, coefficients)
59 def constructElements(self, validBox):
60 """Construct the elements of a CoaddBoundedField."""
61 rng = np.random.Generator(np.random.MT19937(50))
62 validPolygon = lsst.afw.geom.Polygon(lsst.geom.Box2D(validBox)) if validBox else None
63 elements = []
64 validBoxes = []
66 for i in range(10):
67 elements.append(
68 lsst.meas.algorithms.CoaddBoundedField.Element(
69 self.makeRandomField(self.elementBBox),
70 self.makeRandomWcs(self.crval),
71 validPolygon,
72 rng.uniform(low=0.5, high=1.5)
73 )
74 )
75 validBoxes.append(validBox)
76 return elements, validBoxes
78 def setUp(self):
79 self.rng = np.random.Generator(np.random.MT19937(50))
80 self.crval = lsst.geom.SpherePoint(45.0, 45.0, lsst.geom.degrees)
81 self.elementBBox = lsst.geom.Box2I(lsst.geom.Point2I(-50, -50), lsst.geom.Point2I(50, 50))
82 self.coaddWcs = self.makeRandomWcs(self.crval, maxOffset=0.0)
83 self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(-75, -75), lsst.geom.Point2I(75, 75))
84 self.possibleValidBoxes = (None, lsst.geom.Box2I(lsst.geom.Point2I(-25, -25),
85 lsst.geom.Point2I(25, 25)))
87 def testEvaluate(self):
88 """Test the main implementation of CoaddBoundedField::evaluate() by creating images of
89 each constituent field, warping them, and coadding them to form a check image.
91 We run this test twice: once with a regular ValidPolygon, and once
92 with the polygon undefined (ie, set to None); see DM-11008.
93 """
94 def runTest(elements, validBoxes):
95 field = lsst.meas.algorithms.CoaddBoundedField(self.bbox, self.coaddWcs, elements, 0.0)
96 coaddImage = lsst.afw.image.ImageF(self.bbox)
97 warpCtrl = lsst.afw.math.WarpingControl("bilinear")
98 weightMap = lsst.afw.image.ImageF(self.bbox)
99 for element, validBox in zip(elements, validBoxes):
100 elementImage = (lsst.afw.image.ImageF(validBox) if validBox
101 else lsst.afw.image.ImageF(element.field.getBBox()))
102 element.field.fillImage(elementImage, True)
103 warp = lsst.afw.image.ImageF(self.bbox)
104 lsst.afw.math.warpImage(warp, self.coaddWcs, elementImage, element.wcs, warpCtrl, 0.0)
105 coaddImage.scaledPlus(element.weight, warp)
106 warp.getArray()[warp.getArray() != 0.0] = element.weight
107 weightMap += warp
108 coaddImage /= weightMap
109 coaddImage.getArray()[np.isnan(coaddImage.getArray())] = 0.0
110 fieldImage = lsst.afw.image.ImageF(self.bbox)
111 field.fillImage(fieldImage)
113 # The coaddImage we're comparing to has artifacts at the edges of all the input exposures,
114 # due to the interpolation, so we just check that the number of unequal pixels is small (<10%)
115 # This can depend on the random number seed, so if a failure is seen, uncomment the line below
116 # to enable a plot of the differences, and verify by eye that they look like edge artifacts. If
117 # they do, just modify the seed (at the top of this file) or change number-of-pixels threshold
118 # until the test passes.
120 diff = np.abs(fieldImage.getArray() - coaddImage.getArray())
121 relTo = np.abs(fieldImage.getArray())
122 rtol = 1E-2
123 atol = 1E-7
124 bad = np.logical_and(diff > rtol*relTo, diff > atol)
126 # enable this to see a plot of the comparison (but it will always
127 # fail, since it doesn't take into account the artifacts in
128 # coaddImage)
129 if False:
130 self.assertFloatsAlmostEqual(fieldImage.getArray(), coaddImage.getArray(), plotOnFailure=True,
131 rtol=rtol, atol=atol, relTo=relTo)
133 self.assertLess(bad.sum(), 0.10*self.bbox.getArea())
135 for validBox in self.possibleValidBoxes:
136 elements, validBoxes = self.constructElements(validBox)
137 self.assertNotEqual(elements[0], elements[1])
138 runTest(elements, validBoxes)
140 def testEquality(self):
141 def runTest(elements, validBoxes):
142 field1 = lsst.meas.algorithms.CoaddBoundedField(self.bbox, self.coaddWcs, elements, 0.0)
143 field2 = lsst.meas.algorithms.CoaddBoundedField(self.bbox, self.coaddWcs, elements, 0.0)
144 self.assertEqual(field1, field2)
146 bbox = lsst.geom.Box2I(lsst.geom.Point2I(-75, -75), lsst.geom.Point2I(75, 75))
147 field3 = lsst.meas.algorithms.CoaddBoundedField(bbox, self.coaddWcs, elements, 0.0)
148 self.assertEqual(field1, field3)
150 field4 = lsst.meas.algorithms.CoaddBoundedField(self.bbox, self.coaddWcs, elements, 0.0)
151 self.assertEqual(field1, field4)
153 # NOTE: make a copy of the list; [:] to copy segfaults,
154 # copy.copy() doesn't behave nicely on py2 w/our pybind11 objects,
155 # and .elements.copy() doesn't exist on py2.
156 elementsCopy = list(elements)
157 field5 = lsst.meas.algorithms.CoaddBoundedField(self.bbox, self.coaddWcs, elementsCopy, 0.0)
158 self.assertEqual(field1, field5)
160 # inequality tests below here
161 field6 = lsst.meas.algorithms.CoaddBoundedField(self.bbox, self.coaddWcs, elements, 3.0)
162 self.assertNotEqual(field1, field6)
163 field7 = lsst.meas.algorithms.CoaddBoundedField(self.bbox, self.coaddWcs, [], 0.0)
164 self.assertNotEqual(field1, field7)
166 bbox = lsst.geom.Box2I(lsst.geom.Point2I(-74, -75), lsst.geom.Point2I(75, 75))
167 field8 = lsst.meas.algorithms.CoaddBoundedField(bbox, self.coaddWcs, elements, 0.0)
168 self.assertNotEqual(field1, field8)
170 crval = lsst.geom.SpherePoint(45.0, 45.0, lsst.geom.degrees)
171 coaddWcs = self.makeRandomWcs(crval, maxOffset=2.0)
172 field9 = lsst.meas.algorithms.CoaddBoundedField(self.bbox, coaddWcs, elements, 0.0)
173 self.assertNotEqual(field1, field9)
175 elementsCopy = list(elements)
176 elementsCopy[2].weight = 1000000
177 field10 = lsst.meas.algorithms.CoaddBoundedField(self.bbox, self.coaddWcs, elementsCopy, 0.0)
178 self.assertNotEqual(field1, field10)
179 elementsCopy = list(elements)
180 elementsCopy.pop(0)
181 field11 = lsst.meas.algorithms.CoaddBoundedField(self.bbox, self.coaddWcs, elementsCopy, 0.0)
182 self.assertNotEqual(field1, field11)
184 for validBox in self.possibleValidBoxes:
185 runTest(*self.constructElements(validBox))
187 def testPersistence(self):
188 """Test that we can round-trip a CoaddBoundedField through FITS."""
189 def runTest(elements, validBoxes):
190 filename = "testCoaddBoundedField.fits"
191 field1 = lsst.meas.algorithms.CoaddBoundedField(self.bbox, self.coaddWcs, elements, 0.0)
192 field1.writeFits(filename)
193 field2 = lsst.meas.algorithms.CoaddBoundedField.readFits(filename)
194 self.assertEqual(field1, field2)
195 elements2 = field2.getElements()
196 self.assertEqual(elements, elements2)
197 for el1, el2 in zip(elements, elements2):
198 self.assertEqual(el1, el2)
199 self.assertEqual(el1.field, el2.field)
200 self.assertEqual(el1.wcs, el2.wcs)
201 self.assertWcsAlmostEqualOverBBox(el1.wcs, el2.wcs, self.bbox,
202 maxDiffPix=0, maxDiffSky=0*lsst.geom.arcseconds)
203 self.assertEqual(el1.validPolygon, el2.validPolygon)
204 self.assertEqual(el1.weight, el2.weight)
206 self.assertEqual(field1.getCoaddWcs(), field2.getCoaddWcs())
207 self.assertWcsAlmostEqualOverBBox(self.coaddWcs, field2.getCoaddWcs(), self.bbox,
208 maxDiffPix=0, maxDiffSky=0*lsst.geom.arcseconds)
209 self.assertEqual(field1.getDefault(), field2.getDefault())
210 image1 = lsst.afw.image.ImageD(self.bbox)
211 image2 = lsst.afw.image.ImageD(self.bbox)
212 field1.fillImage(image1)
213 field2.fillImage(image2)
214 self.assertImagesEqual(image1, image2)
215 os.remove(filename)
217 for validBox in self.possibleValidBoxes:
218 runTest(*self.constructElements(validBox))
220 def tearDown(self):
221 del self.coaddWcs
222 del self.bbox
225class TestMemory(lsst.utils.tests.MemoryTestCase):
226 pass
229def setup_module(module):
230 lsst.utils.tests.init()
233if __name__ == "__main__": 233 ↛ 234line 233 didn't jump to line 234 because the condition on line 233 was never true
234 lsst.utils.tests.init()
235 unittest.main()