Detector: add crosstalk matrix

Contains coefficients for intra-CCD crosstalk correction.

include/lsst/afw/cameraGeom/Detector.h

@@ -63,6 +63,8 @@ enum DetectorType {
  */
 class Detector {
 public:
+    typedef ndarray::Array<float const, 2> CrosstalkMatrix;
+
     /**
      * Make a Detector
      *
@@ -77,6 +79,7 @@ class Detector {
      * @param transforms map of CameraSys: afw::geom::Transform, where each
      *                   Transform's forward transform transforms from PIXELS
      *                   to the specified camera system
+     * @param crosstalk matrix of crosstalk coefficients
      *
      * @throws lsst::pex::exceptions::InvalidParameterError if:
      * - any amplifier names are not unique
@@ -89,7 +92,8 @@ class Detector {
     explicit Detector(std::string const &name, int id, DetectorType type, std::string const &serial,
                       geom::Box2I const &bbox, lsst::afw::table::AmpInfoCatalog const &ampInfoCatalog,
                       Orientation const &orientation, geom::Extent2D const &pixelSize,
-                      TransformMap::Transforms const &transforms);
+                      TransformMap::Transforms const &transforms,
+                      CrosstalkMatrix const &crosstalk=CrosstalkMatrix());
 
     ~Detector() {}
 
@@ -131,6 +135,14 @@ class Detector {
     /** Get the transform registry */
     TransformMap const getTransformMap() const { return _transformMap; }
 
+    /** Have we got crosstalk coefficients? */
+    bool hasCrosstalk() const {
+        return !(_crosstalk.isEmpty() || _crosstalk.getShape() == ndarray::makeVector(0, 0));
+    }
+
+    /** Get the crosstalk coefficients */
+    CrosstalkMatrix const getCrosstalk() const { return _crosstalk; }
+
     /** Get iterator to beginning of amplifier list */
     lsst::afw::table::AmpInfoCatalog::const_iterator begin() const { return _ampInfoCatalog.begin(); }
 
@@ -304,6 +316,7 @@ class Detector {
     geom::Extent2D _pixelSize;              ///< pixel size (mm)
     CameraSys _nativeSys;                   ///< native coordinate system of this detector
     TransformMap _transformMap;             ///< registry of coordinate transforms
+    CrosstalkMatrix _crosstalk;             ///< crosstalk coefficients
 };
 }  // namespace cameraGeom
 }  // namespace afw

python/lsst/afw/cameraGeom/cameraConfig.py

@@ -1,5 +1,6 @@
 from __future__ import absolute_import, division, print_function
 
+import numpy as np
 import lsst.pex.config as pexConfig
 from .transformConfig import TransformMapConfig
 
@@ -45,6 +46,22 @@ class DetectorConfig(pexConfig.Config):
     transposeDetector = pexConfig.Field(
         "Transpose the pixel grid before orienting in focal plane?", bool)
 
+    crosstalk = pexConfig.ListField(
+        dtype=float,
+        doc=("Flattened crosstalk coefficient matrix; should have nAmps x nAmps entries. "
+             "Once 'reshape'-ed, ``coeffs[i][j]`` is the fraction of the j-th amp present on the i-th amp."),
+        optional=True
+    )
+
+    def getCrosstalk(self, numAmps):
+        """Return a 2-D numpy array of crosstalk coefficients of the proper shape"""
+        if not self.crosstalk:
+            return None
+        try:
+            return np.array(self.crosstalk, dtype=np.float32).reshape((numAmps, numAmps))
+        except:
+            raise RuntimeError("Cannot reshape 'crosstalk' coefficients to square matrix")
+
 
 class CameraConfig(pexConfig.Config):
     """!A configuration that represents (and can be used to construct) a Camera

python/lsst/afw/cameraGeom/cameraFactory.py

@@ -42,7 +42,7 @@ def makeDetector(detectorConfig, ampInfoCatalog, focalPlaneToField):
         pixelSizeMm = pixelSizeMm,
     )
 
-    return Detector(
+    args = [
         detectorConfig.name,
         detectorConfig.id,
         DetectorType(detectorConfig.detectorType),
@@ -52,7 +52,11 @@ def makeDetector(detectorConfig, ampInfoCatalog, focalPlaneToField):
         orientation,
         pixelSizeMm,
         transforms,
-    )
+    ]
+    crosstalk = detectorConfig.getCrosstalk(len(ampInfoCatalog))
+    if crosstalk:
+        args.append(crosstalk)
+    return Detector(*args)
 
 
 def copyDetector(detector, ampInfoCatalog=None):
@@ -76,7 +80,7 @@ def copyDetector(detector, ampInfoCatalog=None):
     return Detector(detector.getName(), detector.getId(), detector.getType(),
                     detector.getSerial(), detector.getBBox(),
                     ampInfoCatalog, detector.getOrientation(), detector.getPixelSize(),
-                    transformDict)
+                    transformDict, detector.getCrosstalk())
 
 
 def makeOrientation(detectorConfig):

python/lsst/afw/cameraGeom/detector.cc

@@ -26,6 +26,9 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 
+#include "numpy/arrayobject.h"
+#include "ndarray/pybind11.h"
+
 #include "lsst/afw/cameraGeom/CameraPoint.h"
 #include "lsst/afw/cameraGeom/CameraSys.h"
 #include "lsst/afw/cameraGeom/Orientation.h"
@@ -46,6 +49,12 @@ namespace cameraGeom {
 PYBIND11_PLUGIN(_detector) {
     py::module mod("_detector", "Python wrapper for afw _detector library");
 
+    // Need to import numpy for ndarray and eigen conversions
+    if (_import_array() < 0) {
+        PyErr_SetString(PyExc_ImportError, "numpy.core.multiarray failed to import");
+        return nullptr;
+    }
+
     /* Module level */
     py::class_<Detector, std::shared_ptr<Detector>> cls(mod, "Detector");
 
@@ -60,9 +69,9 @@ PYBIND11_PLUGIN(_detector) {
     /* Constructors */
     cls.def(py::init<std::string const &, int, DetectorType, std::string const &, geom::Box2I const &,
                      table::AmpInfoCatalog const &, Orientation const &, geom::Extent2D const &,
-                     TransformMap::Transforms const &>(),
+                     TransformMap::Transforms const &, Detector::CrosstalkMatrix const &>(),
             "name"_a, "id"_a, "type"_a, "serial"_a, "bbox"_a, "ampInfoCatalog"_a, "orientation"_a,
-            "pixelSize"_a, "transforms"_a);
+            "pixelSize"_a, "transforms"_a, "crosstalk"_a=Detector::CrosstalkMatrix());
 
     /* Operators */
     cls.def("__getitem__",
@@ -92,6 +101,8 @@ PYBIND11_PLUGIN(_detector) {
     cls.def("getAmpInfoCatalog", &Detector::getAmpInfoCatalog);
     cls.def("getOrientation", &Detector::getOrientation);
     cls.def("getPixelSize", &Detector::getPixelSize);
+    cls.def("hasCrosstalk", &Detector::hasCrosstalk);
+    cls.def("getCrosstalk", &Detector::getCrosstalk);
     cls.def("getTransformMap", &Detector::getTransformMap);
     cls.def("hasTransform", (bool (Detector::*)(CameraSys const &) const) & Detector::hasTransform,
             "cameraSys"_a);

python/lsst/afw/cameraGeom/testUtils.py

@@ -39,6 +39,7 @@ def __init__(self,
                  orientation=Orientation(),
                  plateScale=20.0,
                  radialDistortion=0.925,
+                 crosstalk=None,
                  modFunc=None,
                  ):
         """!Construct a DetectorWrapper
@@ -57,6 +58,7 @@ def __init__(self,
             (the r^3 coefficient of the radial distortion polynomial
             that converts FIELD_ANGLE in radians to FOCAL_PLANE in mm);
             0.925 is the value Dave Monet measured for lsstSim data
+        @param[in] crosstalk  crosstalk coefficient matrix
         @param[in] modFunc  a function that can modify attributes just before constructing the detector;
             modFunc receives one argument: a DetectorWrapper with all attributes except detector set.
         """
@@ -105,6 +107,9 @@ def __init__(self,
             CameraSys(TAN_PIXELS, self.name): pixelToTanPixel,
             CameraSys(ACTUAL_PIXELS, self.name): afwGeom.makeRadialTransform([0, 0.95, 0.01]),
         }
+        if crosstalk is None:
+            crosstalk = [[0.0 for _ in range(numAmps)] for _ in range(numAmps)]
+        self.crosstalk = crosstalk
         if modFunc:
             modFunc(self)
         self.detector = Detector(
@@ -117,6 +122,7 @@ def __init__(self,
             self.orientation,
             self.pixelSize,
             self.transMap,
+            np.array(self.crosstalk, dtype=np.float32),
         )
 
 

src/cameraGeom/Detector.cc

@@ -31,7 +31,7 @@ namespace cameraGeom {
 Detector::Detector(std::string const &name, int id, DetectorType type, std::string const &serial,
                    geom::Box2I const &bbox, table::AmpInfoCatalog const &ampInfoCatalog,
                    Orientation const &orientation, geom::Extent2D const &pixelSize,
-                   TransformMap::Transforms const &transforms)
+                   TransformMap::Transforms const &transforms, CrosstalkMatrix const &crosstalk)
         : _name(name),
           _id(id),
           _type(type),
@@ -42,7 +42,8 @@ Detector::Detector(std::string const &name, int id, DetectorType type, std::stri
           _orientation(orientation),
           _pixelSize(pixelSize),
           _nativeSys(CameraSys(PIXELS, name)),
-          _transformMap(_nativeSys, transforms) {
+          _transformMap(_nativeSys, transforms),
+          _crosstalk(crosstalk) {
     _init();
 }
 
@@ -114,6 +115,22 @@ void Detector::_init() {
             throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
         }
     }
+
+    // ensure crosstalk coefficients matrix is square
+    if (hasCrosstalk()) {
+        auto shape = _crosstalk.getShape();
+        assert(shape.size() == 2);  // we've declared this as a 2D array
+        if (shape[0] != shape[1]) {
+            std::ostringstream os;
+            os << "Non-square crosstalk matrix: " << _crosstalk << " for detector \"" << _name << "\"";
+            throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
+        }
+        if (shape[0] != _ampInfoCatalog.size()) {
+            std::ostringstream os;
+            os << "Wrong size crosstalk matrix: " << _crosstalk << " for detector \"" << _name << "\"";
+            throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
+        }
+    }
 }
 
 //

tests/test_detector.py

@@ -28,6 +28,8 @@
 from builtins import range
 from builtins import zip
 
+import numpy as np
+
 import lsst.utils.tests
 import lsst.pex.exceptions
 import lsst.afw.geom as afwGeom
@@ -91,6 +93,18 @@ def addBadCameraSys(dw):
         with self.assertRaises(lsst.pex.exceptions.Exception):
             DetectorWrapper(modFunc=addBadCameraSys)
 
+        # These break in the pybind layer
+        for crosstalk in ([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0],  # 1D and not numpy
+                          np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]),  # 1D, wrong numpy type
+                          np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0], dtype=np.float32),  # 1D
+                          ):
+            self.assertRaises(TypeError, DetectorWrapper, crosstalk=crosstalk)
+        # These break in the Detector ctor: wrong shape
+        self.assertRaises(lsst.pex.exceptions.InvalidParameterError,
+                          DetectorWrapper, crosstalk=np.array([[1.0, 2.0], [3.0, 4.0]]))
+        self.assertRaises(lsst.pex.exceptions.InvalidParameterError,
+                          DetectorWrapper, crosstalk=np.array([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]]))
+
     def testTransform(self):
         """Test the transform method
         """