Coverage for tests/test_build_camera.py: 98%
156 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-11-10 01:57 -0800
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/>.
21import os.path
22import re
23import unittest
25import astropy.units as u
26import numpy as np
27import yaml
28from scipy.linalg import lstsq
29from scipy.optimize import minimize
31import lsst.drp.tasks
32import lsst.drp.tasks.gbdesAstrometricFit
33import lsst.utils
34from lsst.afw.cameraGeom.testUtils import FIELD_ANGLE, PIXELS
35from lsst.afw.table import ExposureCatalog
36from lsst.drp.tasks.build_camera import (
37 BuildCameraFromAstrometryConfig,
38 BuildCameraFromAstrometryTask,
39 _z_function,
40)
41from lsst.obs.subaru import HyperSuprimeCam
43TESTDIR = os.path.abspath(os.path.dirname(__file__))
46class TestBuildCameraFromAstrometry(lsst.utils.tests.TestCase):
48 def setUp(self):
50 rng = np.random.default_rng(seed=54321)
52 # The tests will be done with data simulated to look like HSC data.
53 self.camera = HyperSuprimeCam().getCamera()
54 self.detectorList = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
55 self.visitList = [1176, 17900, 17930, 17934]
56 self.rotationAngle = (270 * u.degree).to(u.radian).value
57 self.task = BuildCameraFromAstrometryTask()
59 datadir = os.path.join(TESTDIR, "data")
60 with open(os.path.join(datadir, "sample_wcs.yaml"), "r") as f:
61 modelParams = yaml.load(f, Loader=yaml.Loader)
62 self.mapParams = {}
63 for key, value in modelParams.items():
64 if isinstance(key, str) and re.fullmatch(".+/poly", key):
65 if value["Type"] == "Identity":
66 continue
67 polyCoefficients = np.concatenate(
68 [value["XPoly"]["Coefficients"], value["YPoly"]["Coefficients"]]
69 )
70 self.mapParams[key] = polyCoefficients
72 # Make semi-random visit maps:
73 for visit in self.visitList:
74 positionOffset = rng.random(2)
75 angleOffset = rng.random() * (2 * np.pi)
76 params = np.array(
77 [
78 positionOffset[0],
79 np.cos(angleOffset),
80 -np.sin(angleOffset),
81 positionOffset[1],
82 np.sin(angleOffset),
83 np.cos(angleOffset),
84 ]
85 )
86 self.mapParams[f"{visit}/poly"] = params
88 visitSummary = ExposureCatalog.readFits(os.path.join(datadir, "visitSummary_1176.fits"))
89 gbdesTask = lsst.drp.tasks.gbdesAstrometricFit.GbdesAstrometricFitTask()
90 _, self.mapTemplate = gbdesTask.make_yaml(visitSummary)
91 self.mapTemplate["BAND/DEVICE/poly"]["XMax"] = self.camera[10].getBBox().endX
92 self.mapTemplate["BAND/DEVICE/poly"]["YMax"] = self.camera[10].getBBox().endY
94 recenteredModel = self.task._recenter_model(
95 self.mapParams, self.mapTemplate, self.detectorList, self.camera
96 )
97 self.originalFieldAngleX = (recenteredModel.newTPx * u.degree).to(u.radian).value
98 self.originalFieldAngleY = (recenteredModel.newTPy * u.degree).to(u.radian).value
100 bbox = self.camera[self.detectorList[0]].getBBox()
101 self.xPix = (self.task.x + 1) * bbox.endX / 2
102 self.yPix = (self.task.y + 1) * bbox.endY / 2
104 def test_recenter_model(self):
105 deviceParams = []
106 for element in self.mapTemplate["BAND/DEVICE"]["Elements"]:
107 for detector in self.detectorList:
108 detectorTemplate = element.replace("DEVICE", str(detector))
109 detectorTemplate = detectorTemplate.replace("BAND", ".+")
110 for k, params in self.mapParams.items():
111 if re.fullmatch(detectorTemplate, k):
112 deviceParams.append(params)
113 deviceArray = np.vstack(deviceParams)
115 visitParams = []
116 for element in self.mapTemplate["EXPOSURE"]["Elements"]:
117 for visit in self.visitList:
118 visitTemplate = element.replace("EXPOSURE", str(visit))
119 for k, params in self.mapParams.items():
120 if re.fullmatch(visitTemplate, k):
121 visitParams.append(params)
122 identityVisitParams = np.array([0, 1, 0, 0, 0, 1])
123 visitParams.append(identityVisitParams)
124 expoArray = np.vstack(visitParams)
126 # Get the tangent plane positions from the original device maps:
127 origTpX = []
128 origTpY = []
129 for deviceMap in deviceArray:
130 nCoeffsDev = len(deviceMap) // 2
131 deviceDegree = lsst.drp.tasks.gbdesAstrometricFit._degreeFromNCoeffs(nCoeffsDev)
132 intX = _z_function(deviceMap[:nCoeffsDev], self.task.x, self.task.y, order=deviceDegree)
133 intY = _z_function(deviceMap[nCoeffsDev:], self.task.x, self.task.y, order=deviceDegree)
134 origTpX.append(intX)
135 origTpY.append(intY)
137 origTpX = np.array(origTpX).ravel()
138 origTpY = np.array(origTpY).ravel()
140 for method in ["affine", "rigid"]:
141 # Get the interim positions with the new device maps:
142 recenteredModel = self.task._recenter_model(
143 self.mapParams,
144 self.mapTemplate,
145 self.detectorList,
146 self.camera,
147 transform=method,
148 )
149 origTpX, origTpY = recenteredModel.inputCoords
150 newIntX = recenteredModel.newTPx
151 newIntY = recenteredModel.newTPy
153 # Now fit the per-visit parameters with the constraint that the new
154 # tangent plane positions match the old tangent plane positions,
155 # and then check they are sufficiently close to the old values.
156 newExpoArrayEmp = np.zeros(expoArray.shape)
157 for e, expo in enumerate(expoArray):
158 origMapX = expo[0] + origTpX * expo[1] + origTpY * expo[2]
159 origMapY = expo[3] + origTpX * expo[4] + origTpY * expo[5]
161 def min_function(params):
162 tpX = params[0] + params[1] * newIntX + params[2] * newIntY
163 tpY = params[3] + params[4] * newIntX + params[5] * newIntY
165 diff = ((origMapX - tpX)) ** 2 + ((origMapY - tpY)) ** 2
166 return 1000 * diff.sum()
168 def jac(params):
169 tpX = params[0] + params[1] * newIntX + params[2] * newIntY
170 tpY = params[3] + params[4] * newIntX + params[5] * newIntY
171 dX = 2 * 1000 * (origMapX - tpX)
172 dY = 2 * 1000 * (origMapY - tpY)
173 jacobian = [
174 -dX.sum(),
175 -(dX * newIntX).sum(),
176 -(dX * newIntY).sum(),
177 -dY.sum(),
178 -(dY * newIntX).sum(),
179 -(dY * newIntY).sum(),
180 ]
181 return jacobian
183 res = minimize(min_function, expo, method="Newton-CG", jac=jac, options={"xtol": 1e-14})
184 newExpoArrayEmp[e] = res.x
186 newMapX = (
187 newExpoArrayEmp[e][0] + newIntX * newExpoArrayEmp[e][1] + newIntY * newExpoArrayEmp[e][2]
188 )
189 newMapY = (
190 newExpoArrayEmp[e][3] + newIntX * newExpoArrayEmp[e][4] + newIntY * newExpoArrayEmp[e][5]
191 )
192 self.assertFloatsAlmostEqual(origMapX, newMapX, atol=1e-10)
193 self.assertFloatsAlmostEqual(origMapY, newMapY, atol=1e-10)
195 baselineCameraTPx, baselineCameraTPy = recenteredModel.targetCoords
197 # Check that the average offset between the input camera and the
198 # new camera is close to zero.
199 dX = newIntX - baselineCameraTPx
200 dY = newIntY - baselineCameraTPy
202 self.assertAlmostEqual(dX.mean(), 0)
203 self.assertAlmostEqual(dY.mean(), 0)
205 # Check that the average rotation the input camera and the new
206 # camera is close to zero.
207 norm = (baselineCameraTPx**2 + baselineCameraTPy**2) ** 0.5 * (newIntX**2 + newIntY**2) ** 0.5
208 cosTheta = (baselineCameraTPx * newIntX + baselineCameraTPy * newIntY) / norm
209 self.assertAlmostEqual(cosTheta.mean(), 1, places=4)
211 if method == "affine":
212 # Check that there is no slope on average between input camera
213 # and new camera.
214 A = np.column_stack([np.ones_like(baselineCameraTPx), baselineCameraTPx, baselineCameraTPy])
215 xSolution, *_ = lstsq(A, newIntX)
216 ySolution, *_ = lstsq(A, newIntY)
217 self.assertFloatsAlmostEqual(xSolution, np.array([0, 1, 0]), atol=1e-6)
218 self.assertFloatsAlmostEqual(ySolution, np.array([0, 0, 1]), atol=1e-6)
220 def test_run_with_basic_model(self):
222 config = BuildCameraFromAstrometryConfig()
223 config.modelSplitting = "basic"
224 task = BuildCameraFromAstrometryTask(config=config)
226 camera = task.run(
227 self.mapParams,
228 self.mapTemplate,
229 self.detectorList,
230 self.visitList,
231 self.camera,
232 )
233 testX, testY = [], []
234 for dev in camera:
235 faX, faY = (
236 dev.getTransform(PIXELS, FIELD_ANGLE)
237 .getMapping()
238 .applyForward(np.array([self.xPix, self.yPix]))
239 )
240 testX.append(faX)
241 testY.append(faY)
242 testX = np.concatenate(testX)
243 testY = np.concatenate(testY)
245 self.assertFloatsAlmostEqual(self.originalFieldAngleX, testX, atol=1e-12)
246 self.assertFloatsAlmostEqual(self.originalFieldAngleY, testY, atol=1e-12)
248 def test_run_with_splitModel(self):
250 config = BuildCameraFromAstrometryConfig()
251 config.modelSplitting = "physical"
252 config.modelSplittingTolerance = 1e-6
253 task = BuildCameraFromAstrometryTask(config=config)
254 camera = task.run(
255 self.mapParams,
256 self.mapTemplate,
257 self.detectorList,
258 self.visitList,
259 self.camera,
260 )
262 testX, testY = [], []
263 for dev in camera:
264 faX, faY = (
265 dev.getTransform(PIXELS, FIELD_ANGLE)
266 .getMapping()
267 .applyForward(np.array([self.xPix, self.yPix]))
268 )
269 testX.append(faX)
270 testY.append(faY)
271 testX = np.concatenate(testX)
272 testY = np.concatenate(testY)
274 # The "physical" model splitting is not expected to
275 # reconstruct the input field angles perfectly.
276 self.assertFloatsAlmostEqual(self.originalFieldAngleX, testX, atol=1e-4)
277 self.assertFloatsAlmostEqual(self.originalFieldAngleY, testY, atol=1e-4)
280def setup_module(module):
281 lsst.utils.tests.init()
284if __name__ == "__main__": 284 ↛ 285line 284 didn't jump to line 285 because the condition on line 284 was never true
285 lsst.utils.tests.init()
286 unittest.main()