Add uncalibrateImage to PhotoCalib

include/lsst/afw/image/PhotoCalib.h

@@ -350,6 +350,20 @@ class PhotoCalib final : public table::io::PersistableFacade<PhotoCalib>, public
     MaskedImage<float> calibrateImage(MaskedImage<float> const &maskedImage,
                                       bool includeScaleUncertainty = true) const;
 
+    /**
+     * Return a un-calibrated image, with pixel values in ADU (or whatever the original input to
+     * this photoCalib was).
+     *
+     * Mask pixels are propagated directly from the input image.
+     *
+     * @param maskedImage The masked image with pixel units of nJy to uncalibrate.
+     * @param includeScaleUncertainty Include the uncertainty on the calibration in the resulting variance?
+     *
+     * @return The uncalibrated masked image.
+     */
+    MaskedImage<float> uncalibrateImage(MaskedImage<float> const &maskedImage,
+                                        bool includeScaleUncertainty = true) const;
+
     /**
      * Return a flux calibrated catalog, with new `_flux`, `_fluxErr`, `_mag`, and `_magErr` fields.
      *

python/lsst/afw/image/_photoCalib.cc

@@ -171,6 +171,8 @@ void declarePhotoCalib(lsst::cpputils::python::WrapperCollection &wrappers) {
 
                 cls.def("calibrateImage", &PhotoCalib::calibrateImage, "maskedImage"_a,
                         "includeScaleUncertainty"_a = true);
+                cls.def("uncalibrateImage", &PhotoCalib::uncalibrateImage, "maskedImage"_a,
+                        "includeScaleUncertainty"_a = true);
 
                 cls.def("calibrateCatalog",
                         py::overload_cast<afw::table::SourceCatalog const &,

src/image/PhotoCalib.cc

@@ -94,6 +94,31 @@ void toNanojanskyVariance(ndarray::Array<float const, 2, 1> const &instFlux,
                (eigenInstFluxVar / eigenInstFlux.square() + (scaleErr / eigenFlux * eigenInstFlux).square());
 }
 
+/**
+ * Compute the variance of an array of fluxes, for calculations on nJy calibrated MaskedImages.
+ *
+ * MaskedImage stores the variance instead of the standard deviation, so we can skip a sqrt().
+ * Usage in calibrateImage() is to compute instFlux (ADU) directly, so that the calibration scale is
+ * implictly `flux/instFlux` (and thus not passed as an argument).
+ *
+ * @param flux[in] The flux in nJy.
+ * @param fluxVar[in] The variance of the fluxes.
+ * @param scaleErr[in] The error on the calibration scale.
+ * @param instFlux[in] The instrumental fluxes calculated from the fluxes.
+ * @param out[out] The output array to fill with the instrumental variance values.
+ */
+void fromNanojanskyVariance(ndarray::Array<float const, 2, 1> const &flux,
+                            ndarray::Array<float const, 2, 1> const &fluxVar, float scaleErr,
+                            ndarray::Array<float const, 2, 1> const &instFlux,
+                            ndarray::Array<float, 2, 1> out) {
+    auto eigenFlux = ndarray::asEigen<Eigen::ArrayXpr>(flux);
+    auto eigenFluxVar = ndarray::asEigen<Eigen::ArrayXpr>(fluxVar);
+    auto eigenInstFlux = ndarray::asEigen<Eigen::ArrayXpr>(instFlux);
+    auto eigenOut = ndarray::asEigen<Eigen::ArrayXpr>(out);
+    eigenOut = (eigenFluxVar / eigenFlux.square() - (scaleErr / eigenFlux * eigenInstFlux).square()) *
+               eigenInstFlux.square();
+}
+
 double toMagnitudeErr(double instFlux, double instFluxErr, double scale, double scaleErr) {
     return 2.5 / std::log(10.0) * hypot(instFluxErr / instFlux, scaleErr / scale);
 }
@@ -280,6 +305,28 @@ MaskedImage<float> PhotoCalib::calibrateImage(MaskedImage<float> const &maskedIm
     return result;
 }
 
+MaskedImage<float> PhotoCalib::uncalibrateImage(MaskedImage<float> const &maskedImage,
+                                                bool includeScaleUncertainty) const {
+    // Deep copy construct, as we're mutiplying in-place.
+    auto result = MaskedImage<float>(maskedImage, true);
+
+    if (_isConstant) {
+        *(result.getImage()) /= _calibrationMean;
+    } else {
+        _calibration->divideImage(*(result.getImage()), true);  // only in the overlap region
+    }
+    if (includeScaleUncertainty) {
+        fromNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(),
+                               _calibrationErr, result.getImage()->getArray(),
+                               result.getVariance()->getArray());
+    } else {
+        fromNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(), 0,
+                               result.getImage()->getArray(), result.getVariance()->getArray());
+    }
+
+    return result;
+}
+
 afw::table::SourceCatalog PhotoCalib::calibrateCatalog(afw::table::SourceCatalog const &catalog,
                                                        std::vector<std::string> const &instFluxFields) const {
     auto const &inSchema = catalog.getSchema();

tests/test_photoCalib.py

@@ -594,11 +594,17 @@ def testCalibrateImageConstant(self):
         photoCalib = lsst.afw.image.PhotoCalib(self.calibration, self.calibrationErr)
         result = photoCalib.calibrateImage(maskedImage)
         self.assertMaskedImagesAlmostEqual(expect, result)
+        uncalibrated = photoCalib.uncalibrateImage(result)
+        # High rtol because our specific choice of calib values results in a
+        # very-close-to-but-not-1 value, which results in poor round-tripping.
+        self.assertMaskedImagesAlmostEqual(maskedImage, uncalibrated, rtol=4e-2)
 
         # same test, but without using the calibration error
         expect = makeCalibratedMaskedImageNoCalibrationError(image, mask, variance, outImage)
         result = photoCalib.calibrateImage(maskedImage, includeScaleUncertainty=False)
         self.assertMaskedImagesAlmostEqual(expect, result)
+        uncalibrated = photoCalib.uncalibrateImage(result, includeScaleUncertainty=False)
+        self.assertMaskedImagesAlmostEqual(maskedImage, uncalibrated)
 
     def testCalibrateImageNonConstant(self):
         """Test a spatially-varying calibration."""
@@ -612,11 +618,17 @@ def testCalibrateImageNonConstant(self):
         photoCalib = lsst.afw.image.PhotoCalib(self.linearXCalibration, self.calibrationErr)
         result = photoCalib.calibrateImage(maskedImage)
         self.assertMaskedImagesAlmostEqual(expect, result)
+        uncalibrated = photoCalib.uncalibrateImage(result)
+        # High rtol because our specific choice of calib values results in a
+        # very-close-to-but-not-1 value, which results in poor round-tripping.
+        self.assertMaskedImagesAlmostEqual(maskedImage, uncalibrated, rtol=4e-2)
 
         # same test, but without using the calibration error
         expect = makeCalibratedMaskedImageNoCalibrationError(image, mask, variance, outImage)
         result = photoCalib.calibrateImage(maskedImage, includeScaleUncertainty=False)
         self.assertMaskedImagesAlmostEqual(expect, result)
+        uncalibrated = photoCalib.uncalibrateImage(result, includeScaleUncertainty=False)
+        self.assertMaskedImagesAlmostEqual(maskedImage, uncalibrated)
 
     def testCalibrateImageNonConstantSubimage(self):
         """Test a non-constant calibration on a sub-image, to ensure we're
@@ -634,11 +646,17 @@ def testCalibrateImageNonConstantSubimage(self):
         photoCalib = lsst.afw.image.PhotoCalib(self.linearXCalibration, self.calibrationErr)
         result = photoCalib.calibrateImage(subImage)
         self.assertMaskedImagesAlmostEqual(expect.subset(subBox), result)
+        uncalibrated = photoCalib.uncalibrateImage(result)
+        # High rtol because our specific choice of calib values results in a
+        # very-close-to-but-not-1 value, which results in poor round-tripping.
+        self.assertMaskedImagesAlmostEqual(subImage, uncalibrated, rtol=2e-3)
 
         # same test, but without using the calibration error
         expect = makeCalibratedMaskedImageNoCalibrationError(image, mask, variance, outImage)
         result = photoCalib.calibrateImage(maskedImage, includeScaleUncertainty=False)
         self.assertMaskedImagesAlmostEqual(expect, result)
+        uncalibrated = photoCalib.uncalibrateImage(result, includeScaleUncertainty=False)
+        self.assertMaskedImagesAlmostEqual(subImage, uncalibrated)
 
     def testNonPositiveMeans(self):
         # no negative calibrations