yamlCamera.py

Go to the documentation of this file.

# This file is part of obs_base.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.

import yaml

import numpy as np
import lsst.afw.cameraGeom as cameraGeom
import lsst.geom as geom
import lsst.afw.geom as afwGeom
from lsst.afw.cameraGeom import Amplifier, Camera, ReadoutCorner


__all__ = ["makeCamera"]


def makeCamera(cameraFile):
    """An imaging camera (e.g. the LSST 3Gpix camera)

    Parameters
    ----------
    cameraFile : `str`
        Camera description YAML file.

    Returns
    -------
    camera : `lsst.afw.cameraGeom.Camera`
        The desired Camera
    """

    with open(cameraFile) as fd:
        cameraParams = yaml.load(fd, Loader=yaml.CLoader)

    cameraName = cameraParams["name"]

    #
    # Handle distortion models.
    #
    plateScale = geom.Angle(cameraParams["plateScale"], geom.arcseconds)
    nativeSys = cameraGeom.CameraSys(cameraParams["transforms"].pop("nativeSys"))
    transforms = makeTransformDict(nativeSys, cameraParams["transforms"], plateScale)

    ccdParams = cameraParams["CCDs"]
    detectorConfigList = makeDetectorConfigList(ccdParams)

    amplifierDict = {}
    for ccdName, ccdValues in ccdParams.items():
        amplifierDict[ccdName] = makeAmplifierList(ccdValues)

    return makeCameraFromCatalogs(cameraName, detectorConfigList, nativeSys, transforms, amplifierDict)


def makeDetectorConfigList(ccdParams):
    """Make a list of detector configs

    Returns
    -------
    detectorConfig : `list` of `lsst.afw.cameraGeom.DetectorConfig`
        A list of detector configs.
    """
    detectorConfigs = []
    for name, ccd in ccdParams.items():
        detectorConfig = cameraGeom.DetectorConfig()
        detectorConfigs.append(detectorConfig)

        detectorConfig.name = name
        detectorConfig.id = ccd['id']
        detectorConfig.serial = ccd['serial']
        detectorConfig.detectorType = ccd['detectorType']
        if 'physicalType' in ccd:
            detectorConfig.physicalType = ccd['physicalType']
        # This is the orientation we need to put the serial direction along the x-axis
        detectorConfig.bbox_x0, detectorConfig.bbox_y0 = ccd['bbox'][0]
        detectorConfig.bbox_x1, detectorConfig.bbox_y1 = ccd['bbox'][1]
        detectorConfig.pixelSize_x, detectorConfig.pixelSize_y = ccd['pixelSize']
        detectorConfig.transformDict.nativeSys = ccd['transformDict']['nativeSys']
        transforms = ccd['transformDict']['transforms']
        detectorConfig.transformDict.transforms = None if transforms == 'None' else transforms
        detectorConfig.refpos_x, detectorConfig.refpos_y = ccd['refpos']
        detectorConfig.offset_x, detectorConfig.offset_y = ccd['offset']
        detectorConfig.transposeDetector = ccd['transposeDetector']
        detectorConfig.pitchDeg = ccd['pitch']
        detectorConfig.yawDeg = ccd['yaw']
        detectorConfig.rollDeg = ccd['roll']
        if 'crosstalk' in ccd:
            detectorConfig.crosstalk = ccd['crosstalk']

    return detectorConfigs


def makeAmplifierList(ccd):
    """Construct a list of AmplifierBuilder objects
    """
    # Much of this will need to be filled in when we know it.
    assert len(ccd) > 0
    amp = list(ccd['amplifiers'].values())[0]

    rawBBox = makeBBoxFromList(amp['rawBBox'])  # total in file
    xRawExtent, yRawExtent = rawBBox.getDimensions()

    readCorners = {"LL": ReadoutCorner.LL,
                   "LR": ReadoutCorner.LR,
                   "UL": ReadoutCorner.UL,
                   "UR": ReadoutCorner.UR}

    amplifierList = []
    for name, amp in sorted(ccd['amplifiers'].items(), key=lambda x: x[1]['hdu']):
        amplifier = Amplifier.Builder()
        amplifier.setName(name)

        ix, iy = amp['ixy']
        perAmpData = amp['perAmpData']
        if perAmpData:
            x0, y0 = 0, 0           # origin of data within each amp image
        else:
            x0, y0 = ix*xRawExtent, iy*yRawExtent

        rawDataBBox = makeBBoxFromList(amp['rawDataBBox'])  # Photosensitive area
        xDataExtent, yDataExtent = rawDataBBox.getDimensions()
        amplifier.setBBox(geom.BoxI(
            geom.PointI(ix*xDataExtent, iy*yDataExtent), rawDataBBox.getDimensions()))

        rawBBox = makeBBoxFromList(amp['rawBBox'])
        rawBBox.shift(geom.ExtentI(x0, y0))
        amplifier.setRawBBox(rawBBox)

        rawDataBBox = makeBBoxFromList(amp['rawDataBBox'])
        rawDataBBox.shift(geom.ExtentI(x0, y0))
        amplifier.setRawDataBBox(rawDataBBox)

        rawSerialOverscanBBox = makeBBoxFromList(amp['rawSerialOverscanBBox'])
        rawSerialOverscanBBox.shift(geom.ExtentI(x0, y0))
        amplifier.setRawHorizontalOverscanBBox(rawSerialOverscanBBox)

        rawParallelOverscanBBox = makeBBoxFromList(amp['rawParallelOverscanBBox'])
        rawParallelOverscanBBox.shift(geom.ExtentI(x0, y0))
        amplifier.setRawVerticalOverscanBBox(rawParallelOverscanBBox)

        rawSerialPrescanBBox = makeBBoxFromList(amp['rawSerialPrescanBBox'])
        rawSerialPrescanBBox.shift(geom.ExtentI(x0, y0))
        amplifier.setRawPrescanBBox(rawSerialPrescanBBox)

        if perAmpData:
            amplifier.setRawXYOffset(geom.Extent2I(ix*xRawExtent, iy*yRawExtent))
        else:
            amplifier.setRawXYOffset(geom.Extent2I(0, 0))

        amplifier.setReadoutCorner(readCorners[amp['readCorner']])
        amplifier.setGain(amp['gain'])
        amplifier.setReadNoise(amp['readNoise'])
        amplifier.setSaturation(amp['saturation'])
        amplifier.setSuspectLevel(amp.get('suspect', np.nan))

        # flip data when assembling if needs be (e.g. data from the serial at the top of a CCD)
        flipX, flipY = amp.get("flipXY")

        amplifier.setRawFlipX(flipX)
        amplifier.setRawFlipY(flipY)
        # linearity placeholder stuff
        amplifier.setLinearityCoeffs([float(val) for val in amp['linearityCoeffs']])
        amplifier.setLinearityType(amp['linearityType'])
        amplifier.setLinearityThreshold(float(amp['linearityThreshold']))
        amplifier.setLinearityMaximum(float(amp['linearityMax']))
        amplifier.setLinearityUnits("DN")
        amplifierList.append(amplifier)
    return amplifierList


def makeAmpInfoCatalog(ccd):
    """Backward compatible name.
    """
    return makeAmplifierList(ccd)


def makeBBoxFromList(ylist):
    """Given a list [(x0, y0), (xsize, ysize)], probably from a yaml file,
    return a BoxI
    """
    (x0, y0), (xsize, ysize) = ylist
    return geom.BoxI(geom.PointI(x0, y0), geom.ExtentI(xsize, ysize))


def makeTransformDict(nativeSys, transformDict, plateScale):
    """Make a dictionary of TransformPoint2ToPoint2s from yaml, mapping from nativeSys

    Parameters
    ----------
    nativeSys : `lsst.afw.cameraGeom.CameraSys`
    transformDict : `dict`
        A dict specifying parameters of transforms; keys are camera system names.
    plateScale : `lsst.geom.Angle`
        The size of a pixel in angular units/mm (e.g. 20 arcsec/mm for LSST)

    Returns
    -------
    transforms : `dict`
        A dict of `lsst.afw.cameraGeom.CameraSys` : `lsst.afw.geom.TransformPoint2ToPoint2`

    The resulting dict's keys are `~lsst.afw.cameraGeom.CameraSys`,
    and the values are Transforms *from* NativeSys *to* CameraSys
    """
    # As other comments note this is required, and this is one function where it's assumed
    assert nativeSys == cameraGeom.FOCAL_PLANE, "Cameras with nativeSys != FOCAL_PLANE are not supported."

    resMap = dict()

    for key, transform in transformDict.items():
        transformType = transform["transformType"]
        knownTransformTypes = ["affine", "radial"]
        if transformType not in knownTransformTypes:
            raise RuntimeError("Saw unknown transform type for %s: %s (known types are: [%s])" % (
                key, transform["transformType"], ", ".join(knownTransformTypes)))

        if transformType == "affine":
            affine = geom.AffineTransform(np.array(transform["linear"]),
                                          np.array(transform["translation"]))

            transform = afwGeom.makeTransform(affine)
        elif transformType == "radial":
            # radial coefficients of the form [0, 1 (no units), C2 (rad), usually 0, C3 (rad^2), ...]
            # Radial distortion is modeled as a radial polynomial that converts from focal plane radius
            # (in mm) to field angle (in radians). The provided coefficients are divided by the plate
            # scale (in radians/mm) meaning that C1 is always 1.
            radialCoeffs = np.array(transform["coeffs"])

            radialCoeffs *= plateScale.asRadians()
            transform = afwGeom.makeRadialTransform(radialCoeffs)
        else:
            raise RuntimeError("Impossible condition \"%s\" is not in: [%s])" % (
                transform["transformType"], ", ".join(knownTransformTypes)))

        resMap[cameraGeom.CameraSys(key)] = transform

    return resMap


def makeCameraFromCatalogs(cameraName, detectorConfigList, nativeSys, transformDict, amplifierDict,
                           pupilFactoryClass=cameraGeom.pupil.PupilFactory):
    """Construct a Camera instance from a dictionary of
       detector name : `lsst.afw.cameraGeom.amplifier`

    Parameters
    ----------
    cameraName : `str`
        The name of the camera
    detectorConfigList : `list`
        A list of `lsst.afw.cameraGeom.cameraConfig.DetectorConfig`
    nativeSys : `lsst.afw.cameraGeom.CameraSys`
        The native transformation type; must be `lsst.afw.cameraGeom.FOCAL_PLANE`
    transformDict : `dict`
        A dict of lsst.afw.cameraGeom.CameraSys : `lsst.afw.geom.TransformPoint2ToPoint2`
    amplifierDict : `dict`
        A dictionary of detector name :
                           `lsst.afw.cameraGeom.Amplifier.Builder`
    pupilFactoryClass : `type`, optional
        Class to attach to camera;
             `lsst.default afw.cameraGeom.PupilFactory`

    Returns
    -------
    camera : `lsst.afw.cameraGeom.Camera`
        New Camera instance.

    Notes
    ------
    Copied from `lsst.afw.cameraGeom.cameraFactory` with permission and encouragement
    from Jim Bosch
    """

    # nativeSys=FOCAL_PLANE seems to be assumed in various places in this file
    # (e.g. the definition of TAN_PIXELS), despite CameraConfig providing the
    # illusion that it's configurable.
    # Note that we can't actually get rid of the nativeSys config option
    # without breaking lots of on-disk camera configs.
    assert nativeSys == cameraGeom.FOCAL_PLANE, "Cameras with nativeSys != FOCAL_PLANE are not supported."

    focalPlaneToField = transformDict[cameraGeom.FIELD_ANGLE]

    cameraBuilder = Camera.Builder(cameraName)
    cameraBuilder.setPupilFactoryClass(pupilFactoryClass)

    # Ensure all transforms in the camera transform dict are included.
    for toSys, transform in transformDict.items():
        cameraBuilder.setTransformFromFocalPlaneTo(toSys, transform)

    for detectorConfig in detectorConfigList:
        # This should build all detector pixel -> focalPlane transforms.
        cameraGeom.addDetectorBuilderFromConfig(cameraBuilder, detectorConfig,
                                                amplifierDict[detectorConfig.name],
                                                focalPlaneToField)

        # For reasons I don't understand, some obs_ packages (e.g. HSC) set
        # nativeSys to None for their detectors (which doesn't seem to be
        # permitted by the config class!), but they really mean PIXELS. For
        # backwards compatibility we use that as the default...
        detectorNativeSys = detectorConfig.transformDict.nativeSys
        detectorNativeSys = (cameraGeom.PIXELS if detectorNativeSys is None else
                             cameraGeom.CameraSysPrefix(detectorNativeSys))

        # ...well, actually, it seems that we've always assumed down in C++
        # that the answer is always PIXELS without ever checking that it is.
        # So let's assert that it is, since there are hints all over this file
        # (e.g. the definition of TAN_PIXELS) that other parts of the codebase
        # have regularly made that assumption as well.  Note that we can't
        # actually get rid of the nativeSys config option without breaking
        # lots of on-disk camera configs.
        assert detectorNativeSys == cameraGeom.PIXELS, \
            "Detectors with nativeSys != PIXELS are not supported."
        detectorNativeSys = cameraGeom.CameraSys(detectorNativeSys, detectorConfig.name)

    return cameraBuilder.finish()