Coverage for python / lsst / drp / tasks / build_camera.py: 8%
313 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-22 09:04 +0000
1# This file is part of drp_tasks.
2#
3# LSST Data Management System
4# This product includes software developed by the
5# LSST Project (http://www.lsst.org/).
6# See COPYRIGHT file at the top of the source tree.
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 <https://www.lsstcorp.org/LegalNotices/>.
21#
22import re
23from collections import defaultdict
25import astshim as ast
26import numpy as np
27import scipy.optimize
29import lsst.pex.config
30from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, Camera
31from lsst.afw.geom import TransformPoint2ToPoint2
32from lsst.geom import degrees
33from lsst.pipe.base import Task
35__all__ = [
36 "BuildCameraFromAstrometryConfig",
37 "BuildCameraFromAstrometryTask",
38]
40# Set up global variable to use in minimization callback.
41nFunctionEval = 0
44def _z_function(params, x, y, order=4):
45 """Convenience function for calculating the value of a 2-d polynomial with
46 coefficients given by an array of parameters.
48 Parameters
49 ----------
50 params : `np.ndarray`
51 Coefficientss for polynomial with terms ordered
52 [1, x, y, x^2, xy, y^2, ...].
53 x, y : `np.ndarray`
54 x and y values at which to evaluate the 2-d polynomial.
55 order : `int`, optional
56 Order of the given polynomial.
58 Returns
59 -------
60 z : `np.ndarray`
61 Value of polynomial evaluated at (x, y).
62 """
63 coeffs = np.zeros((order + 1, order + 1))
64 c = 0
65 for i in range(order + 1):
66 for j in range(i + 1):
67 coeffs[j, i - j] = params[c]
68 c += 1
69 z = np.polynomial.polynomial.polyval2d(x, y, coeffs.T)
70 return z
73def _z_function_dx(params, x, y, order=4):
74 """Convenience function for calculating the derivative with respect to x of
75 a 2-d polynomial with coefficients given by an array of parameters.
77 Parameters
78 ----------
79 params : `np.ndarray`
80 Coefficientss for polynomial with terms ordered
81 [1, x, y, x^2, xy, y^2, ...].
82 x, y : `np.ndarray`
83 x and y values at which to evaluate the 2-d polynomial.
84 order : `int`, optional
85 Order of the given polynomial.
87 Returns
88 -------
89 z : `np.ndarray`
90 Derivative of polynomial evaluated at (x, y).
91 """
92 coeffs = np.zeros((order + 1, order + 1))
93 c = 0
94 for i in range(order + 1):
95 for j in range(i + 1):
96 coeffs[j, i - j] = params[c]
97 c += 1
98 derivCoeffs = np.zeros(coeffs.shape)
99 derivCoeffs[:, :-1] = coeffs[:, 1:]
100 derivCoeffs[:, :-1] *= np.arange(1, order + 1)
102 z = np.polynomial.polynomial.polyval2d(x, y, derivCoeffs.T)
103 return z
106def _z_function_dy(params, x, y, order=4):
107 """Convenience function for calculating the derivative with respect to y of
108 a 2-d polynomial with coefficients given by an array of parameters.
110 Parameters
111 ----------
112 params : `np.ndarray`
113 Coefficientss for polynomial with terms ordered
114 [1, x, y, x^2, xy, y^2, ...].
115 x, y : `np.ndarray`
116 x and y values at which to evaluate the 2-d polynomial.
117 order : `int`, optional
118 Order of the given polynomial.
120 Returns
121 -------
122 z : `np.ndarray`
123 Derivative of polynomial evaluated at (x, y).
124 """
125 coeffs = np.zeros((order + 1, order + 1))
126 c = 0
127 for i in range(order + 1):
128 for j in range(i + 1):
129 coeffs[j, i - j] = params[c]
130 c += 1
131 derivCoeffs = np.zeros(coeffs.shape)
132 derivCoeffs[:-1] = coeffs[1:]
133 derivCoeffs[:-1] *= np.arange(1, order + 1).reshape(-1, 1)
135 z = np.polynomial.polynomial.polyval2d(x, y, derivCoeffs.T)
136 return z
139class BuildCameraFromAstrometryConfig(lsst.pex.config.Config):
140 """Configuration for BuildCameraTask."""
142 tangentPlaneDegree = lsst.pex.config.Field(
143 dtype=int,
144 doc=(
145 "Order of polynomial mapping between the focal plane and the tangent plane. Only used if "
146 "modelSplitting='physical'."
147 ),
148 default=9,
149 )
150 focalPlaneDegree = lsst.pex.config.Field(
151 dtype=int,
152 doc=(
153 "Order of polynomial mapping between the pixels and the focal plane. Only used if "
154 "modelSplitting='physical'."
155 ),
156 default=3,
157 )
158 modelSplitting = lsst.pex.config.ChoiceField(
159 dtype=str,
160 doc="How to split the camera model into pixel to focal plane and focal plane to tangent plane parts.",
161 default="basic",
162 allowed={
163 "basic": "Put all the mapping except scaling into the pixel to focal plane part",
164 "physical": (
165 "Fit a more physical mapping where as much as possible of the model variability goes"
166 "into the focal plane to tangent plane part of the model, to imitate the effects of "
167 "the telescope optics"
168 ),
169 },
170 )
171 plateScale = lsst.pex.config.Field(
172 dtype=float,
173 doc=("Scaling between camera coordinates in mm and angle on the sky in" " arcsec."),
174 default=20.005867576692737,
175 )
176 astInversionTolerance = lsst.pex.config.Field(
177 dtype=float,
178 doc="Tolerance for AST map inversion.",
179 default=0.005,
180 )
181 astInversionMaxIter = lsst.pex.config.Field(
182 dtype=int, doc="Maximum iterations fro AST map inversion.", default=10
183 )
184 modelSplittingTolerance = lsst.pex.config.Field(
185 dtype=float,
186 doc="Average error in model splitting minimization acceptable for convergence.",
187 default=1e-8,
188 )
191class BuildCameraFromAstrometryTask(Task):
192 """Build an `lsst.afw.cameraGeom.Camera` object out of the `gbdes`
193 polynomials mapping from pixels to the tangent plane.
195 Parameters
196 ----------
198 """
200 ConfigClass = BuildCameraFromAstrometryConfig
201 _DefaultName = "buildCameraFromAstrometry"
203 def __init__(self, **kwargs):
204 super().__init__(**kwargs)
206 # The gbdes model normalizes the pixel positions to the range -1 - +1.
207 X = np.arange(-1, 1, 0.1)
208 Y = np.arange(-1, 1, 0.1)
209 x, y = np.meshgrid(X, Y)
210 self.x = x.ravel()
211 self.y = y.ravel()
213 def run(self, mapParams, mapTemplate, detectorList, visitList, inputCamera, rotationAngle):
214 """Convert the model parameters into a Camera object.
216 Parameters
217 ----------
218 mapParams : `dict` [`float`]
219 Parameters that describe the per-detector and per-visit parts of
220 the astrometric model.
221 mapTemplate : `dict`
222 Dictionary describing the format of the astrometric model,
223 following the convention in `gbdes`.
224 detectorList : `list` [`int`]
225 List of detector ids.
226 visitList : `list` [`int`]
227 List of ids for visits that were used to train the input model.
228 camera : `lsst.afw.cameraGeom.Camera`
229 Camera object from which to take pupil function, name, and other
230 properties.
231 rotationAngle : `float`
232 Value in radians of the average rotation angle of the input visits.
234 Returns
235 -------
236 camera : `lsst.afw.cameraGeom.Camera`
237 Camera object with transformations set by the input mapParams.
238 """
240 # Normalize the model.
241 newParams, newTPx, newTPy = self._normalize_model(mapParams, mapTemplate, detectorList, visitList)
243 if self.config.modelSplitting == "basic":
244 # Put all of the camera distortion into the pixels->focal plane
245 # part of the distortion model, with the focal plane->tangent plane
246 # part only used for scaling between the focal plane and sky.
247 pixToFocalPlane, focalPlaneToTangentPlane = self._basic_model(newParams, detectorList)
249 else:
250 # Fit two polynomials, such that the first describes the pixel->
251 # focal plane part of the model, and the second describes the focal
252 # plane->tangent plane part of the model, with the goal of
253 # generating a more physically-motivated distortion model.
254 pixToFocalPlane, focalPlaneToTangentPlane = self._split_model(newTPx, newTPy, detectorList)
256 # Turn the mappings into a Camera object.
257 camera = self._translate_to_afw(
258 pixToFocalPlane, focalPlaneToTangentPlane, detectorList, inputCamera, rotationAngle
259 )
261 return camera
263 def _normalize_model(self, mapParams, mapTemplate, detectorList, visitList):
264 """Normalize the camera mappings, such that they correspond to the
265 average visit.
267 In order to break degeneracies, `gbdes` sets one visit to have the
268 identity for the per-visit part of its model. This visit is thus
269 imprinted in the per-detector part of the model. This function
270 normalizes the model such that the parameters correspond instead to the
271 average per-visit behavior. This step assumes that the per-visit part
272 of the model is a first-order polynomial, in which case this conversion
273 can be done without loss of information.
275 Parameters
276 ----------
277 mapParams : `dict` [`float`]
278 Parameters that describe the per-detector and per-visit parts of
279 the astrometric model.
280 mapTemplate : `dict`
281 Dictionary describing the format of the astrometric model,
282 following the convention in `gbdes`.
283 detectorList : `list` [`int`]
284 List of detector ids.
285 visitList : `list` [`int`]
286 List of ids for visits that were used to train the input model.
288 Returns
289 -------
290 newDeviceArray : `np.ndarray`
291 Array of NxM, where N is the number of detectors, and M is the
292 number of coefficients for each per-detector mapping.
293 newTPx, newTPy : `np.ndarray`
294 Projection of `self.x` and `self.y` onto the tangent plane, given
295 the normalized mapping.
296 """
298 # Get the per-device and per-visit parts of the model.
299 deviceParams = []
300 visitParams = []
301 for element in mapTemplate["BAND/DEVICE"]["Elements"]:
302 for detector in detectorList:
303 detectorTemplate = element.replace("DEVICE", str(detector))
304 detectorTemplate = detectorTemplate.replace("BAND", ".+")
305 for k, params in mapParams.items():
306 if re.fullmatch(detectorTemplate, k):
307 deviceParams.append(params)
308 deviceArray = np.vstack(deviceParams)
310 for element in mapTemplate["EXPOSURE"]["Elements"]:
311 for visit in visitList:
312 visitTemplate = element.replace("EXPOSURE", str(visit))
313 for k, params in mapParams.items():
314 if re.fullmatch(visitTemplate, k):
315 visitParams.append(params)
316 identityVisitParams = np.array([0, 1, 0, 0, 0, 1])
317 visitParams.append(identityVisitParams)
318 expoArray = np.vstack(visitParams)
320 expoMean = expoArray.mean(axis=0)
322 # Shift the per-device part of the model to correspond with the mean
323 # per-visit behavior.
324 newDeviceArray = np.zeros(deviceArray.shape)
325 nCoeffsDev = deviceArray.shape[1] // 2
326 newDeviceArray[:, :nCoeffsDev] = (
327 deviceArray[:, :nCoeffsDev] * expoMean[1] + deviceArray[:, nCoeffsDev:] * expoMean[2]
328 )
329 newDeviceArray[:, nCoeffsDev:] = (
330 deviceArray[:, :nCoeffsDev] * expoMean[4] + deviceArray[:, nCoeffsDev:] * expoMean[5]
331 )
332 newDeviceArray[:, 0] += expoMean[0]
333 newDeviceArray[:, nCoeffsDev] += expoMean[3]
335 # Then get the tangent plane positions from the new device model:
336 newTPx = []
337 newTPy = []
338 for deviceMap in newDeviceArray:
339 nCoeffsDev = len(deviceMap) // 2
340 deviceDegree = int(-1.5 + 0.5 * (1 + 8 * nCoeffsDev) ** 0.5)
342 intX = _z_function(deviceMap[:nCoeffsDev], self.x, self.y, order=deviceDegree)
343 intY = _z_function(deviceMap[nCoeffsDev:], self.x, self.y, order=deviceDegree)
344 newTPx.append(intX)
345 newTPy.append(intY)
347 newTPx = np.array(newTPx).ravel()
348 newTPy = np.array(newTPy).ravel()
350 return newDeviceArray, newTPx, newTPy
352 def _basic_model(self, modelParameters, detectorList):
353 """This will just convert the pix->fp parameters into the right format,
354 and return an identity mapping for the fp->tp part.
356 Parameters
357 ----------
358 modelParameters : `np.ndarray`
359 Array of NxM, where N is the number of detectors, and M is the
360 number of coefficients for each per-detector mapping.
361 detectorList : `list` [`int`]
362 List of detector ids.
364 Returns
365 -------
366 pixelsToFocalPlane : `dict` [`int`, [`string`, `float`]]
367 Dictionary giving the per-detector pixel to focal plane mapping,
368 with keys for each detector giving the x and y transformation
369 polynomials.
370 focalPlaneToTangentPlane : `dict` [`string`, `float`]
371 Dictionary giving the focal plane to tangent plane mapping.
372 """
374 nCoeffsFp = modelParameters.shape[1] // 2
375 pixelsToFocalPlane = defaultdict(dict)
376 for d, det in enumerate(detectorList):
377 pixelsToFocalPlane[det]["x"] = modelParameters[d][:nCoeffsFp]
378 pixelsToFocalPlane[det]["y"] = modelParameters[d][nCoeffsFp:]
380 focalPlaneToTangentPlane = {"x": [0, 1, 0], "y": [0, 0, 1]}
382 return pixelsToFocalPlane, focalPlaneToTangentPlane
384 def _split_model(self, targetX, targetY, detectorList, pixToFpGuess=None, fpToTpGuess=None):
385 """Fit a two-step model, with one polynomial per detector fitting the
386 pixels to focal plane part, followed by a polynomial fitting the focal
387 plane to tangent plane part.
389 The linear parameters of the focal plane to tangent plane model are
390 fixed to be an identity map to remove degeneracies.
392 Parameters
393 ----------
394 targetX, targetY : `np.ndarray`
395 Target x and y values in the tangent plane.
396 detectorList : `list` [`int`]
397 List of detector ids.
398 pixToFpGuess : `dict` [`int`, [`string`, `float`]]
399 Initial guess for the pixels to focal plane mapping to use in the
400 model fit, in the form of a dictionary giving the per-detector
401 pixel to focal plane mapping, with keys for each detector giving
402 the x and y transformation polynomials.
403 focalPlaneToTangentPlane : `dict` [`string`, `float`]
404 Initial guess for the focal plane to tangent plane mapping to use
405 in the model fit, in the form of a dictionary giving the focal
406 plane to tangent plane mapping.
408 Returns
409 -------
410 pixelsToFocalPlane : `dict` [`int`, [`string`, `float`]]
411 Dictionary giving the per-detector pixel to focal plane mapping,
412 with keys for each detector giving the x and y transformation
413 polynomials.
414 focalPlaneToTangentPlane : `dict` [`string`, `float`]
415 Dictionary giving the focal plane to tangent plane mapping.
416 """
418 tpDegree = self.config.tangentPlaneDegree
419 fpDegree = self.config.focalPlaneDegree
421 nCoeffsFP = int((fpDegree + 2) * (fpDegree + 1) / 2)
422 nCoeffsFPTot = len(detectorList) * nCoeffsFP * 2
424 nCoeffsTP = int((tpDegree + 2) * (tpDegree + 1) / 2)
425 # The constant and linear terms will be fixed to remove degeneracy with
426 # the pix->fp part of the model
427 nCoeffsFixed = 3
428 nCoeffsFree = nCoeffsTP - nCoeffsFixed
429 fixedParams = np.zeros(nCoeffsTP * 2)
430 fixedParams[1] = 1
431 fixedParams[nCoeffsTP + 2] = 1
433 nX = len(self.x)
434 # We need an array of the form [1, x, y, x^2, xy, y^2, ...] to use when
435 # calculating derivatives.
436 xyArray = np.zeros((nCoeffsFP, nX))
437 c = 0
438 for ii in range(fpDegree + 1):
439 for jj in range(ii + 1):
440 coeffs = np.zeros((fpDegree + 1, fpDegree + 1))
441 coeffs[jj, ii - jj] = 1
442 xyArray[c] = np.polynomial.polynomial.polyval2d(self.x, self.y, coeffs.T)
443 c += 1
445 def two_part_function(params):
446 # The function giving the split model.
447 fpXAll = []
448 fpYAll = []
449 for i in range(len(detectorList)):
450 fpX = _z_function(
451 params[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree
452 )
453 fpY = _z_function(
454 params[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], self.x, self.y, order=fpDegree
455 )
456 fpXAll.append(fpX)
457 fpYAll.append(fpY)
458 fpXAll = np.array(fpXAll).ravel()
459 fpYAll = np.array(fpYAll).ravel()
461 tpParams = fixedParams.copy()
462 tpParams[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree]
463 tpParams[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFPTot + nCoeffsFree :]
465 tpX = _z_function(tpParams[:nCoeffsTP], fpXAll, fpYAll, order=tpDegree)
466 tpY = _z_function(tpParams[nCoeffsTP:], fpXAll, fpYAll, order=tpDegree)
467 return tpX, tpY
469 def min_function(params):
470 # The least-squares function to be minimized.
471 tpX, tpY = two_part_function(params)
472 diff = ((targetX - tpX)) ** 2 + ((targetY - tpY)) ** 2
473 result = diff.sum()
474 return result
476 def jac(params):
477 # The Jacobian of min_function, which is needed for the minimizer.
478 tpX, tpY = two_part_function(params)
479 dX = -2 * (targetX - tpX)
480 dY = -2 * (targetY - tpY)
481 jacobian = np.zeros(len(params))
482 fpParams = params[:nCoeffsFPTot]
483 tpParams = fixedParams.copy()
484 tpParams[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree]
485 tpParams[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFPTot + nCoeffsFree :]
486 fpXAll = []
487 fpYAll = []
488 # Fill in the derivatives for the pix->fp terms.
489 for i in range(len(detectorList)):
490 dX_i = dX[nX * i : nX * (i + 1)]
491 dY_i = dY[nX * i : nX * (i + 1)]
492 fpX = _z_function(
493 fpParams[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree
494 )
495 fpY = _z_function(
496 fpParams[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)],
497 self.x,
498 self.y,
499 order=fpDegree,
500 )
501 fpXAll.append(fpX)
502 fpYAll.append(fpY)
504 dTpX_dFpX = _z_function_dx(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree)
505 dTpX_dFpY = _z_function_dy(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree)
506 dTpY_dFpX = _z_function_dx(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree)
507 dTpY_dFpY = _z_function_dy(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree)
509 dTpX_part = np.concatenate([xyArray * dTpX_dFpX, xyArray * dTpX_dFpY])
510 dTpY_part = np.concatenate([xyArray * dTpY_dFpX, xyArray * dTpY_dFpY])
511 jacobian_i = (dX_i * dTpX_part + dY_i * dTpY_part).sum(axis=1)
512 jacobian[2 * nCoeffsFP * i : 2 * nCoeffsFP * (i + 1)] = jacobian_i
514 fpXAll = np.array(fpXAll).ravel()
515 fpYAll = np.array(fpYAll).ravel()
517 # Fill in the derivatives for the fp->tp terms.
518 for j in range(nCoeffsFree):
519 tParams = np.zeros(nCoeffsTP)
520 tParams[nCoeffsFixed + j] = 1
521 tmpZ = _z_function(tParams, fpXAll, fpYAll, order=tpDegree)
522 jacobian[nCoeffsFPTot + j] = (dX * tmpZ).sum()
523 jacobian[nCoeffsFPTot + nCoeffsFree + j] = (dY * tmpZ).sum()
524 return jacobian
526 def hessian(params):
527 # The Hessian of min_function, which is needed for the minimizer.
528 hessian = np.zeros((len(params), len(params)))
530 fpParams = params[:nCoeffsFPTot]
531 tpParams = fixedParams.copy()
532 tpParams[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree]
533 tpParams[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFPTot + nCoeffsFree :]
534 fpXAll = []
535 fpYAll = []
537 # Loop over the detectors to fill in the d(pix->fp)**2 and
538 # d(pix->fp)d(fp->tp) cross terms.
539 for i in range(len(detectorList)):
540 fpX = _z_function(
541 fpParams[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree
542 )
543 fpY = _z_function(
544 fpParams[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)],
545 self.x,
546 self.y,
547 order=fpDegree,
548 )
549 fpXAll.append(fpX)
550 fpYAll.append(fpY)
551 dTpX_dFpX = _z_function_dx(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree)
552 dTpX_dFpY = _z_function_dy(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree)
553 dTpY_dFpX = _z_function_dx(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree)
554 dTpY_dFpY = _z_function_dy(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree)
556 dTpX_part = np.concatenate([xyArray * dTpX_dFpX, xyArray * dTpX_dFpY])
557 dTpY_part = np.concatenate([xyArray * dTpY_dFpX, xyArray * dTpY_dFpY])
559 for k in range(2 * nCoeffsFP):
560 for m in range(2 * nCoeffsFP):
561 hessian[2 * nCoeffsFP * i + k, 2 * nCoeffsFP * i + m] = (
562 2 * (dTpX_part[k] * dTpX_part[m] + dTpY_part[k] * dTpY_part[m]).sum()
563 )
565 for j in range(nCoeffsFree):
566 tmpParams = np.zeros(nCoeffsTP)
567 tmpParams[nCoeffsFixed + j] = 1
568 tmpZ = _z_function(tmpParams, fpX, fpY, order=tpDegree)
570 hessX = 2 * (dTpX_part[k] * tmpZ).sum()
571 hessY = 2 * (dTpY_part[k] * tmpZ).sum()
572 # dTP_x part
573 hessian[2 * nCoeffsFP * i + k, nCoeffsFPTot + j] = hessX
574 hessian[nCoeffsFPTot + j, 2 * nCoeffsFP * i + k] = hessX
575 # dTP_y part
576 hessian[2 * nCoeffsFP * i + k, nCoeffsFPTot + nCoeffsFree + j] = hessY
577 hessian[nCoeffsFPTot + nCoeffsFree + j, 2 * nCoeffsFP * i + k] = hessY
579 fpXAll = np.array(fpXAll).ravel()
580 fpYAll = np.array(fpYAll).ravel()
581 tmpZArray = np.zeros((nCoeffsFree, nX * len(detectorList)))
582 # Finally, get the d(fp->tp)**2 terms
583 for j in range(nCoeffsFree):
584 tParams = np.zeros(nCoeffsTP)
585 tParams[nCoeffsFixed + j] = 1
586 tmpZArray[j] = _z_function(tParams, fpXAll, fpYAll, order=tpDegree)
587 for j in range(nCoeffsFree):
588 for m in range(nCoeffsFree):
589 # X-Y cross terms are zero
590 hess = 2 * (tmpZArray[j] * tmpZArray[m]).sum()
591 hessian[nCoeffsFPTot + j, nCoeffsFPTot + m] = hess
592 hessian[nCoeffsFPTot + nCoeffsFree + j, nCoeffsFPTot + nCoeffsFree + m] = hess
593 return hessian
595 global nFunctionEval
596 nFunctionEval = 0
598 def callbackMF(params):
599 global nFunctionEval
600 self.log.info(f"Iteration {nFunctionEval}, function value {min_function(params)}")
601 nFunctionEval += 1
603 initGuess = np.zeros(nCoeffsFPTot + 2 * nCoeffsFree)
604 if pixToFpGuess:
605 useVar = min(nCoeffsFP, len(pixToFpGuess[0]["x"]))
606 for i, det in enumerate(detectorList):
607 initGuess[2 * nCoeffsFP * i : 2 * nCoeffsFP * i + useVar] = pixToFpGuess[det]["x"][:useVar]
608 initGuess[(2 * i + 1) * nCoeffsFP : (2 * i + 1) * nCoeffsFP + useVar] = pixToFpGuess[det][
609 "y"
610 ][:useVar]
611 if fpToTpGuess:
612 useVar = min(nCoeffsTP, len(fpToTpGuess["x"]))
613 initGuess[nCoeffsFPTot : nCoeffsFPTot + useVar - nCoeffsFixed] = fpToTpGuess["x"][
614 nCoeffsFixed:useVar
615 ]
616 initGuess[nCoeffsFPTot + nCoeffsFree : nCoeffsFPTot + nCoeffsFree + useVar - nCoeffsFixed] = (
617 fpToTpGuess["y"][nCoeffsFixed:useVar]
618 )
620 self.log.info(f"Initial value of least squares function: {min_function(initGuess)}")
621 res = scipy.optimize.minimize(
622 min_function,
623 initGuess,
624 callback=callbackMF,
625 method="Newton-CG",
626 jac=jac,
627 hess=hessian,
628 options={"xtol": self.config.modelSplittingTolerance},
629 )
631 # Convert parameters to a dictionary.
632 pixToFP = {}
633 for i, det in enumerate(detectorList):
634 pixToFP[det] = {
635 "x": res.x[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP],
636 "y": res.x[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)],
637 }
639 fpToTpAll = fixedParams.copy()
640 fpToTpAll[nCoeffsFixed:nCoeffsTP] = res.x[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree]
641 fpToTpAll[nCoeffsFixed + nCoeffsTP :] = res.x[nCoeffsFPTot + nCoeffsFree :]
642 fpToTp = {"x": fpToTpAll[:nCoeffsTP], "y": fpToTpAll[nCoeffsTP:]}
644 return pixToFP, fpToTp
646 def make_ast_polymap_coeffs(self, order, xCoeffs, yCoeffs):
647 """Convert polynomial coefficients in gbdes format to AST PolyMap input
648 format.
650 Paramaters
651 ----------
652 order: `int`
653 Polynomial order
654 xCoeffs, yCoeffs: `list` of `float`
655 Forward or inverse polynomial coefficients for the x and y axes
656 of output, in this order:
657 x0y0, x1y0, x0y1, x2y0, x1y1, x0y2, ...xNy0, xN-1y1, ...x0yN
658 where N is the polynomial order.
659 Returns
660 -------
661 polyArray : `np.ndarray`
662 Array formatted for AST PolyMap input.
663 """
664 N = len(xCoeffs)
666 polyArray = np.zeros((2 * N, 4))
668 for i in range(order + 1):
669 for j in range(order + 1):
670 if (i + j) > order:
671 continue
672 vectorIndex = int(((i + j) * (i + j + 1)) / 2 + j)
673 polyArray[vectorIndex] = [xCoeffs[vectorIndex], 1, i, j]
674 polyArray[vectorIndex + N] = [yCoeffs[vectorIndex], 2, i, j]
676 return polyArray
678 def _translate_to_afw(
679 self, pixToFocalPlane, focalPlaneToTangentPlane, detectorList, inputCamera, rotationAngle
680 ):
681 """Convert the model parameters to a Camera object.
683 Parameters
684 ----------
685 pixelsToFocalPlane : `dict` [`int`, [`string`, `float`]]
686 Dictionary giving the per-detector pixel to focal plane mapping,
687 with keys for each detector giving the x and y transformation
688 polynomials.
689 focalPlaneToTangentPlane : `dict` [`string`, `float`]
690 Dictionary giving the focal plane to tangent plane mapping.
691 rotationAngle : `float`
692 Value in radians of the average rotation angle of the input visits.
694 Returns
695 -------
696 output_camera : `lsst.afw.cameraGeom.Camera`
697 Camera object containing the pix->fp and fp->tp transforms.
698 """
699 if self.config.modelSplitting == "basic":
700 tpDegree = 1
701 nCoeffsFP = len(pixToFocalPlane[detectorList[0]]["x"])
702 fpDegree = int(-1.5 + 0.5 * (1 + 8 * nCoeffsFP) ** 0.5)
703 else:
704 tpDegree = self.config.tangentPlaneDegree
705 fpDegree = self.config.focalPlaneDegree
707 scaleConvert = (1 * degrees).asArcseconds() / self.config.plateScale
709 cameraBuilder = Camera.Builder(inputCamera.getName())
710 cameraBuilder.setPupilFactoryName(inputCamera.getPupilFactoryName())
712 # Convert fp->tp to AST format:
713 forwardCoeffs = self.make_ast_polymap_coeffs(
714 tpDegree, focalPlaneToTangentPlane["x"], focalPlaneToTangentPlane["y"]
715 )
716 # Reverse rotation from input visits and flip x-axis.
717 cosRot = np.cos(-rotationAngle)
718 sinRot = np.sin(-rotationAngle)
719 rotateAndFlipCoeffs = self.make_ast_polymap_coeffs(1, [0, cosRot, -sinRot], [0, -sinRot, cosRot])
721 ccdZoom = ast.ZoomMap(2, 1 / scaleConvert)
722 ccdToSky = ast.PolyMap(
723 forwardCoeffs,
724 2,
725 "IterInverse=1, TolInverse=%s, NIterInverse=%s"
726 % (self.config.astInversionTolerance / 2.0, self.config.astInversionMaxIter),
727 )
728 rotateAndFlip = ast.PolyMap(
729 rotateAndFlipCoeffs,
730 2,
731 "IterInverse=1, TolInverse=%s, NIterInverse=%s"
732 % (self.config.astInversionTolerance / 2.0, self.config.astInversionMaxIter),
733 )
734 toRadians = ast.ZoomMap(2, (1 * degrees).asRadians())
736 fullMapping = ccdZoom.then(rotateAndFlip).then(ccdToSky).then(rotateAndFlip).then(toRadians)
737 focalPlaneToTPMapping = TransformPoint2ToPoint2(fullMapping)
738 cameraBuilder.setTransformFromFocalPlaneTo(FIELD_ANGLE, focalPlaneToTPMapping)
740 # Convert pix->fp to AST format:
741 for detector in inputCamera:
742 d = detector.getId()
743 if d not in pixToFocalPlane:
744 # This camera will not include detectors that were not used in
745 # astrometric fit.
746 continue
748 detectorBuilder = cameraBuilder.add(detector.getName(), detector.getId())
749 detectorBuilder.setBBox(detector.getBBox())
750 detectorBuilder.setPixelSize(detector.getPixelSize())
751 for amp in detector.getAmplifiers():
752 detectorBuilder.append(amp.rebuild())
754 normCoeffs = self.make_ast_polymap_coeffs(
755 1, [-1.0, 2 / detector.getBBox().getWidth(), 0], [-1.0, 0, 2 / detector.getBBox().getHeight()]
756 )
757 normMap = ast.PolyMap(
758 normCoeffs,
759 2,
760 "IterInverse=1, TolInverse=%s, NIterInverse=%s"
761 % (self.config.astInversionTolerance / 2.0, self.config.astInversionMaxIter),
762 )
763 forwardDetCoeffs = self.make_ast_polymap_coeffs(
764 fpDegree, pixToFocalPlane[d]["x"], pixToFocalPlane[d]["y"]
765 )
766 ccdToFP = ast.PolyMap(
767 forwardDetCoeffs,
768 2,
769 "IterInverse=1, TolInverse=%s, NIterInverse=%s"
770 % (self.config.astInversionTolerance / 2.0, self.config.astInversionMaxIter),
771 )
772 zoomMap = ast.ZoomMap(2, scaleConvert)
773 fullDetMap = normMap.then(ccdToFP).then(rotateAndFlip).then(zoomMap)
774 ccdToFPMapping = TransformPoint2ToPoint2(fullDetMap)
775 detectorBuilder.setTransformFromPixelsTo(FOCAL_PLANE, ccdToFPMapping)
777 outputCamera = cameraBuilder.finish()
779 return outputCamera