Skip to content

Commit

Permalink
Add option to build camera in gbdes task
Browse files Browse the repository at this point in the history
  • Loading branch information
@cmsaunders
cmsaunders committed Sep 5, 2024
1 parent 61e55fd commit 7b757ef
Show file tree
Hide file tree
Showing 3 changed files with 160 additions and 76 deletions.
123 changes: 60 additions & 63 deletions python/lsst/drp/tasks/build_camera.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,37 +20,23 @@
# see <https://www.lsstcorp.org/LegalNotices/>.
#
import re
import time
from collections import defaultdict

import astshim as ast
import lsst.pex.config
import numpy as np
import scipy.optimize
from lsst.afw.cameraGeom import (
ACTUAL_PIXELS,
FIELD_ANGLE,
FOCAL_PLANE,
PIXELS,
TAN_PIXELS,
Camera,
CameraSys,
DetectorConfig,
)
from lsst.afw.geom import TransformPoint2ToPoint2, transformConfig
from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, Camera
from lsst.afw.geom import TransformPoint2ToPoint2
from lsst.geom import degrees
from lsst.pipe.base import Task

from .gbdesAstrometricFit import _degreeFromNCoeffs, _nCoeffsFromDegree

__all__ = [
"BuildCameraConfig",
"BuildCameraTask",
"BuildCameraFromAstrometryConfig",
"BuildCameraFromAstrometryTask",
]

cameraSysList = [FIELD_ANGLE, FOCAL_PLANE, PIXELS, TAN_PIXELS, ACTUAL_PIXELS]
cameraSysMap = dict((sys.getSysName(), sys) for sys in cameraSysList)

# Set up global variable to use in minimization callback.
Nfeval = 0


Expand Down Expand Up @@ -149,7 +135,7 @@ def _z_function_dy(params, x, y, order=4):
return z


class BuildCameraConfig(lsst.pex.config.Config):
class BuildCameraFromAstrometryConfig(lsst.pex.config.Config):
"""Configuration for BuildCameraTask."""

tangentPlaneDegree = lsst.pex.config.Field(
Expand Down Expand Up @@ -184,7 +170,7 @@ class BuildCameraConfig(lsst.pex.config.Config):
plateScale = lsst.pex.config.Field(
dtype=float,
doc=("Scaling between camera coordinates in mm and angle on the sky in" " arcsec."),
default=1.0,
default=20.005867576692737,
)
astInversionTolerance = lsst.pex.config.Field(
dtype=float,
Expand All @@ -201,29 +187,20 @@ class BuildCameraConfig(lsst.pex.config.Config):
)


class BuildCameraTask(Task):
class BuildCameraFromAstrometryTask(Task):
"""Build an `lsst.afw.cameraGeom.Camera` object out of the `gbdes`
polynomials mapping from pixels to the tangent plane.
Parameters
----------
camera : `lsst.afw.cameraGeom.Camera`
Camera object from which to take pupil function, name, and other
properties.
detectorList : `list` [`int`]
List of detector ids.
visitList : `list` [`int`]
List of ids for visits that were used to train the input model.
"""

ConfigClass = BuildCameraConfig
ConfigClass = BuildCameraFromAstrometryConfig
_DefaultName = "buildCamera"

def __init__(self, camera, detectorList, visitList, **kwargs):
def __init__(self, **kwargs):
super().__init__(**kwargs)
self.camera = camera
self.detectorList = detectorList
self.visitList = visitList

# The gbdes model normalizes the pixel positions to the range -1 - +1.
X = np.arange(-1, 1, 0.1)
Expand All @@ -232,7 +209,7 @@ def __init__(self, camera, detectorList, visitList, **kwargs):
self.x = x.ravel()
self.y = y.ravel()

def run(self, mapParams, mapTemplate):
def run(self, mapParams, mapTemplate, detectorList, visitList, inputCamera, rotationAngle):
"""Convert the model parameters into a Camera object.
Parameters
Expand All @@ -243,6 +220,15 @@ def run(self, mapParams, mapTemplate):
mapTemplate : `dict`
Dictionary describing the format of the astrometric model,
following the convention in `gbdes`.
detectorList : `list` [`int`]
List of detector ids.
visitList : `list` [`int`]
List of ids for visits that were used to train the input model.
camera : `lsst.afw.cameraGeom.Camera`
Camera object from which to take pupil function, name, and other
properties.
rotationAngle : `float`
Value in radians of the average rotation angle of the input visits.
Returns
-------
Expand All @@ -251,27 +237,29 @@ def run(self, mapParams, mapTemplate):
"""

# Normalize the model.
newParams, newIntX, newIntY = self._normalizeModel(mapParams, mapTemplate)
newParams, newIntX, newIntY = self._normalizeModel(mapParams, mapTemplate, detectorList, visitList)

if self.config.modelSplitting == "basic":
# Put all of the camera distortion into the pixels->focal plane
# part of the distortion model, with the focal plane->tangent plane
# part only used for scaling between the focal plane and sky.
pixToFocalPlane, focalPlaneToTangentPlane = self._basicModel(newParams)
pixToFocalPlane, focalPlaneToTangentPlane = self._basicModel(newParams, detectorList)

else:
# Fit two polynomials, such that the first describes the pixel->
# focal plane part of the model, and the second describes the focal
# plane->tangent plane part of the model, with the goal of
# generating a more physically-motivated distortion model.
pixToFocalPlane, focalPlaneToTangentPlane = self._splitModel(newIntX, newIntY)
pixToFocalPlane, focalPlaneToTangentPlane = self._splitModel(newIntX, newIntY, detectorList)

# Turn the mappings into a Camera object.
camera = self._translateToAfw(pixToFocalPlane, focalPlaneToTangentPlane)
camera = self._translateToAfw(
pixToFocalPlane, focalPlaneToTangentPlane, detectorList, inputCamera, rotationAngle
)

return camera

def _normalizeModel(self, mapParams, mapTemplate):
def _normalizeModel(self, mapParams, mapTemplate, detectorList, visitList):
"""Normalize the camera mappings, such that they correspond to the
average visit.
Expand Down Expand Up @@ -306,7 +294,7 @@ def _normalizeModel(self, mapParams, mapTemplate):
deviceParams = []
visitParams = []
for element in mapTemplate["BAND/DEVICE"]["Elements"]:
for detector in self.detectorList:
for detector in detectorList:
detectorTemplate = element.replace("DEVICE", str(detector))
detectorTemplate = detectorTemplate.replace("BAND", ".+")
for k, params in mapParams.items():
Expand All @@ -315,7 +303,7 @@ def _normalizeModel(self, mapParams, mapTemplate):
deviceArray = np.vstack(deviceParams)

for element in mapTemplate["EXPOSURE"]["Elements"]:
for visit in self.visitList:
for visit in visitList:
visitTemplate = element.replace("EXPOSURE", str(visit))
for k, params in mapParams.items():
if re.fullmatch(visitTemplate, k):
Expand All @@ -342,7 +330,7 @@ def _normalizeModel(self, mapParams, mapTemplate):
newIntY = []
for deviceMap in newDeviceArray:
nCoeffsDev = len(deviceMap) // 2
deviceDegree = _degreeFromNCoeffs(nCoeffsDev)
deviceDegree = int(-1.5 + 0.5 * (1 + 8 * nCoeffsDev) ** 0.5)

intX = _z_function(deviceMap[:nCoeffsDev], self.x, self.y, order=deviceDegree)
intY = _z_function(deviceMap[nCoeffsDev:], self.x, self.y, order=deviceDegree)
Expand All @@ -354,7 +342,7 @@ def _normalizeModel(self, mapParams, mapTemplate):

return newDeviceArray, newIntX, newIntY

def _basicModel(self, modelParameters):
def _basicModel(self, modelParameters, detectorList):
"""This will just convert the pix->fp parameters into the right format,
and return an identity mapping for the fp->tp part.
Expand All @@ -376,15 +364,15 @@ def _basicModel(self, modelParameters):

nCoeffsFP = modelParameters.shape[1] // 2
pixelsToFocalPlane = defaultdict(dict)
for d, det in enumerate(self.detectorList):
for d, det in enumerate(detectorList):
pixelsToFocalPlane[det]["x"] = modelParameters[d][:nCoeffsFP]
pixelsToFocalPlane[det]["y"] = modelParameters[d][nCoeffsFP:]

focalPlaneToTangentPlane = {"x": [0, 1, 0], "y": [0, 0, 1]}

return pixelsToFocalPlane, focalPlaneToTangentPlane

def _splitModel(self, targetX, targetY, pixToFPGuess=None, fpToTpGuess=None):
def _splitModel(self, targetX, targetY, detectorList, pixToFPGuess=None, fpToTpGuess=None):
"""Fit a two-step model, with one polynomial per detector fitting the
pixels to focal plane part, followed by a polynomial fitting the focal
plane to tangent plane part.
Expand Down Expand Up @@ -419,10 +407,10 @@ def _splitModel(self, targetX, targetY, pixToFPGuess=None, fpToTpGuess=None):
tpDegree = self.config.tangentPlaneDegree
fpDegree = self.config.focalPlaneDegree

nCoeffsFP = _nCoeffsFromDegree(fpDegree)
nCoeffsFP_tot = len(self.detectorList) * nCoeffsFP * 2
nCoeffsFP = int((fpDegree + 2) * (fpDegree + 1) / 2)
nCoeffsFP_tot = len(detectorList) * nCoeffsFP * 2

nCoeffsTP = _nCoeffsFromDegree(tpDegree)
nCoeffsTP = int((tpDegree + 2) * (tpDegree + 1) / 2)
# The constant and linear terms will be fixed to remove degeneracy with
# the pix->fp part of the model
nCoeffsFixed = 3
Expand All @@ -447,7 +435,7 @@ def two_part_function(params):
# The function giving the split model.
intX_tot = []
intY_tot = []
for i in range(len(self.detectorList)):
for i in range(len(detectorList)):
intX = _z_function(
params[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree
)
Expand Down Expand Up @@ -487,7 +475,7 @@ def jac(params):
intX_tot = []
intY_tot = []
# Fill in the derivatives for the pix->fp terms.
for i in range(len(self.detectorList)):
for i in range(len(detectorList)):
dX_i = dX[nX * i : nX * (i + 1)]
dY_i = dY[nX * i : nX * (i + 1)]
intX = _z_function(
Expand Down Expand Up @@ -537,7 +525,7 @@ def hessian(params):

# Loop over the detectors to fill in the d(pix->fp)**2 and
# d(pix->fp)d(fp->tp) cross terms.
for i in range(len(self.detectorList)):
for i in range(len(detectorList)):
intX = _z_function(
fp_params[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree
)
Expand Down Expand Up @@ -579,7 +567,7 @@ def hessian(params):

intX_tot = np.array(intX_tot).ravel()
intY_tot = np.array(intY_tot).ravel()
tmpZ_array = np.zeros((nCoeffsFree, nX * len(self.detectorList)))
tmpZ_array = np.zeros((nCoeffsFree, nX * len(detectorList)))
# Finally, get the d(fp->tp)**2 terms
for j in range(nCoeffsFree):
tParams = np.zeros(nCoeffsTP)
Expand All @@ -604,7 +592,7 @@ def callbackMF(params):
initGuess = np.zeros(nCoeffsFP_tot + 2 * nCoeffsFree)
if pixToFPGuess:
useVar = min(nCoeffsFP, len(pixToFPGuess[0]["x"]))
for i, det in enumerate(self.detectorList):
for i, det in enumerate(detectorList):
initGuess[2 * nCoeffsFP * i : 2 * nCoeffsFP * i + useVar] = pixToFPGuess[det]["x"][:useVar]
initGuess[(2 * i + 1) * nCoeffsFP : (2 * i + 1) * nCoeffsFP + useVar] = pixToFPGuess[det][
"y"
Expand All @@ -631,7 +619,7 @@ def callbackMF(params):

# Convert parameters to a dictionary.
pixToFP = {}
for i, det in enumerate(self.detectorList):
for i, det in enumerate(detectorList):
pixToFP[det] = {
"x": res.x[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP],
"y": res.x[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)],
Expand Down Expand Up @@ -676,7 +664,9 @@ def makeAstPolyMapCoeffs(self, order, xCoeffs, yCoeffs):

return polyArray

def _translateToAfw(self, pixToFocalPlane, focalPlaneToTangentPlane):
def _translateToAfw(
self, pixToFocalPlane, focalPlaneToTangentPlane, detectorList, inputCamera, rotationAngle
):
"""Convert the model parameters to a Camera object.
Parameters
Expand All @@ -687,6 +677,8 @@ def _translateToAfw(self, pixToFocalPlane, focalPlaneToTangentPlane):
polynomials.
focalPlaneToTangentPlane : `dict` [`string`, `float`]
Dictionary giving the focal plane to tangent plane mapping.
rotationAngle : `float`
Value in radians of the average rotation angle of the input visits.
Returns
-------
Expand All @@ -695,21 +687,25 @@ def _translateToAfw(self, pixToFocalPlane, focalPlaneToTangentPlane):
"""
if self.config.modelSplitting == "basic":
tpDegree = 1
fpDegree = _degreeFromNCoeffs(len(pixToFocalPlane[self.detectorList[0]]["x"]))
nCoeffsFP = len(pixToFocalPlane[detectorList[0]]["x"])
fpDegree = int(-1.5 + 0.5 * (1 + 8 * nCoeffsFP) ** 0.5)
else:
tpDegree = self.config.tangentPlaneDegree
fpDegree = self.config.focalPlaneDegree

scaleConvert = (1 * degrees).asArcseconds() / self.config.plateScale

cameraBuilder = Camera.Builder(self.camera.getName())
cameraBuilder.setPupilFactoryName(self.camera.getPupilFactoryName())
cameraBuilder = Camera.Builder(inputCamera.getName())
cameraBuilder.setPupilFactoryName(inputCamera.getPupilFactoryName())

# Convert fp->tp to AST format:
forwardCoeffs = self.makeAstPolyMapCoeffs(
tpDegree, focalPlaneToTangentPlane["x"], focalPlaneToTangentPlane["y"]
)
rotateAndFlipCoeffs = self.makeAstPolyMapCoeffs(1, [0, 0, -1], [0, -1, 0])
# Reverse rotation from input visits and flip x-axis.
cosRot = np.cos(rotationAngle)
sinRot = np.sin(rotationAngle)
rotateAndFlipCoeffs = self.makeAstPolyMapCoeffs(1, [0, cosRot, -sinRot], [0, -sinRot, cosRot])

ccdZoom = ast.ZoomMap(2, 1 / scaleConvert)
ccdToSky = ast.PolyMap(
Expand All @@ -731,14 +727,15 @@ def _translateToAfw(self, pixToFocalPlane, focalPlaneToTangentPlane):
cameraBuilder.setTransformFromFocalPlaneTo(FIELD_ANGLE, focalPlaneToTPMapping)

# Convert pix->fp to AST format:
for detector in self.camera:
for detector in inputCamera:
d = detector.getId()
if d not in pixToFocalPlane:
# TODO: just grab detector map from the obs_ package.
# This camera will not include detectors that were not used in astrometric fit.
continue

detectorBuilder = cameraBuilder.add(detector.getName(), detector.getId())
# TODO: Add all the detector details like pixel size, etc.
detectorBuilder.setBBox(detector.getBBox())
detectorBuilder.setPixelSize(detector.getPixelSize())
for amp in detector.getAmplifiers():
detectorBuilder.append(amp.rebuild())

Expand Down
Loading

0 comments on commit 7b757ef

Please sign in to comment.