Merge branch 'tickets/DM-23624'

include/lsst/ip/isr/isr.h

@@ -97,12 +97,11 @@ namespace isr {
         );
 
 
-    template<typename ImagePixelT, typename FunctionT>
-    void fitOverscanImage(
-        std::shared_ptr<lsst::afw::math::Function1<FunctionT> > &overscanFunction,
+    template<typename ImagePixelT>
+    std::vector<double> fitOverscanImage(
         lsst::afw::image::MaskedImage<ImagePixelT> const& overscan,
-        double ssize=1.,
-        int sigma=1
+        std::vector<std::string> badPixelMask,
+        bool isTransposed
         );
 
 }}} // namespace lsst::ip::isr

python/lsst/ip/isr/isr.cc

@@ -20,6 +20,7 @@
  * see <https://www.lsstcorp.org/LegalNotices/>.
  */
 #include "pybind11/pybind11.h"
+#include "pybind11/stl.h"
 
 #include <memory>
 
@@ -59,15 +60,16 @@ static void declareAll(py::module& mod, std::string const& suffix) {
     declareCountMaskedPixels<PixelT>(mod, suffix);
 
     mod.def("maskNans", &maskNans<PixelT>, "maskedImage"_a, "maskVal"_a, "allow"_a = 0);
-    mod.def("fitOverscanImage", &fitOverscanImage<PixelT, double>, "overscanFunction"_a, "overscan"_a,
-            "stepSize"_a = 1.1, "sigma"_a = 1);
+    mod.def("fitOverscanImage", &fitOverscanImage<PixelT>,
+            "maskedImage"_a, "badPixelMask"_a, "isTransposed"_a);
 }
 
 }  // namespace lsst::ip::isr::<anonymous>
 
 PYBIND11_MODULE(isr, mod) {
     declareAll<float>(mod, "F");
     declareAll<double>(mod, "D");
+    declareAll<int>(mod, "I");
 }
 
 }  // isr

python/lsst/ip/isr/overscan.py

@@ -20,11 +20,14 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import numpy as np
+import time
 import lsst.afw.math as afwMath
 import lsst.afw.image as afwImage
 import lsst.pipe.base as pipeBase
 import lsst.pex.config as pexConfig
 
+from .isr import fitOverscanImage
+
 __all__ = ["OverscanCorrectionTaskConfig", "OverscanCorrectionTask"]
 
 
@@ -62,7 +65,7 @@ class OverscanCorrectionTaskConfig(pexConfig.Config):
     maskPlanes = pexConfig.ListField(
         dtype=str,
         doc="Mask planes to reject when measuring overscan",
-        default=['SAT'],
+        default=['BAD', 'SAT'],
     )
     overscanIsInt = pexConfig.Field(
         dtype=bool,
@@ -475,8 +478,19 @@ def measureVectorOverscan(self, image):
         calcImage, isTransposed = self.transpose(calcImage)
         masked = self.maskOutliers(calcImage)
 
+        startTime = time.perf_counter()
+
         if self.config.fitType == 'MEDIAN_PER_ROW':
-            overscanVector = self.collapseArrayMedian(masked)
+            mi = afwImage.MaskedImageI(image.getBBox())
+            masked = masked.astype(int)
+            if isTransposed:
+                masked = masked.transpose()
+
+            mi.image.array[:, :] = masked.data[:, :]
+            if bool(masked.mask.shape):
+                mi.mask.array[:, :] = masked.mask[:, :]
+
+            overscanVector = fitOverscanImage(mi, self.config.maskPlanes, isTransposed)
             maskArray = self.maskExtrapolated(overscanVector)
         else:
             collapsed = self.collapseArray(masked)
@@ -507,6 +521,8 @@ def measureVectorOverscan(self, image):
                 # Otherwise we can just use things as normal.
                 overscanVector = evaler(indices, coeffs)
                 maskArray = self.maskExtrapolated(collapsed)
+        endTime = time.perf_counter()
+        self.log.info(f"Overscan measurement took {endTime - startTime}s for {self.config.fitType}")
         return pipeBase.Struct(overscanValue=np.array(overscanVector),
                                maskArray=maskArray,
                                isTransposed=isTransposed)

src/Isr.cc

@@ -49,52 +49,50 @@ size_t maskNans(afw::image::MaskedImage<PixelT> const& mi, afw::image::MaskPixel
     return nPix;
 }
 
-template<typename ImagePixelT, typename FunctionT>
-void fitOverscanImage(
-    std::shared_ptr< afw::math::Function1<FunctionT> > &overscanFunction,
+template<typename ImagePixelT>
+std::vector<double> fitOverscanImage(
     afw::image::MaskedImage<ImagePixelT> const& overscan,
-    double ssize,
-    int sigma
+    std::vector<std::string> badPixelMask,
+    bool isTransposed
 ) {
     typedef afw::image::MaskedImage<ImagePixelT> MaskedImage;
 
-
+    /**
+    This is transposed here to match the existing numpy-array ordering.
+    This effectively transposes the image for us.
+    **/
     const int height = overscan.getHeight();
     const int width  = overscan.getWidth();
-    std::vector<double> values(height);
-    std::vector<double> errors(height);
-    std::vector<double> positions(height);
-
-    std::vector<double> parameters(overscanFunction->getNParameters(), 0.);
-    std::vector<double> stepsize(overscanFunction->getNParameters(), ssize);
-
-    for (int y = 0; y < height; ++y) {
-        /**
-        geom::Box2I bbox       = geom::Box2I( geom::Point2I(0, y),
-                                              geom::Point2I(0, width) );
-        The above was how this was defined before ticket #1556.  As I understand it
-        the following is the new way to do this
-        **/
-        geom::Box2I bbox = geom::Box2I(geom::Point2I(0,y), geom::Point2I(width,y));
-        MaskedImage mi         = MaskedImage(overscan, bbox);
-        afw::math::Statistics stats = afw::math::makeStatistics(*(mi.getImage()), afw::math::MEAN | afw::math::STDEV);
-
-        values[y]    = stats.getValue(afw::math::MEAN);
-        errors[y]    = stats.getValue(afw::math::STDEV);
-        positions[y] = y;
 
+    int length = height;
+    if (isTransposed) {
+        length = width;
+    }
+
+    std::vector<double> values(length);
+
+    afw::math::StatisticsControl statControl;
+    statControl.setAndMask(overscan.getMask()->getPlaneBitMask(badPixelMask));
+
+    const int x0 = overscan.getX0();
+    const int y0 = overscan.getY0();
+    auto origin = geom::Point2I(x0, y0);
+    geom::Extent2I shifter;
+    geom::Extent2I extents;
+    if (isTransposed) {
+        shifter = geom::Extent2I(1, 0);
+        extents = geom::Extent2I(1, height);
+    } else {
+        shifter = geom::Extent2I(0, 1);
+        extents = geom::Extent2I(width, 1);
     }
-    afw::math::FitResults fitResults = afw::math::minimize(
-        *overscanFunction,
-        parameters,
-        stepsize,
-        values,
-        errors,
-        positions,
-        sigma
-        );
-
-    overscanFunction->setParameters(fitResults.parameterList);
+
+    for (int x = 0; x < length; ++x) {
+        MaskedImage mi = MaskedImage(overscan, geom::Box2I(origin, extents));
+        values[x] = afw::math::makeStatistics(mi, afw::math::MEDIAN, statControl).getValue();
+        origin.shift(shifter);
+    }
+    return values;
 }
 
 std::string between(std::string &s, char ldelim, char rdelim) {
@@ -115,26 +113,25 @@ std::string between(std::string &s, char ldelim, char rdelim) {
 // Explicit instantiations
 
 template
-void fitOverscanImage(
-     std::shared_ptr<afw::math::Function1<double> > &overscanFunction,
-    afw::image::MaskedImage<float> const& overscan,
-    double ssize,
-    int sigma);
-
+std::vector<double> fitOverscanImage<int>(
+    afw::image::MaskedImage<int> const&, std::vector<std::string> badPixelMask, bool isTransposed);
+template
+std::vector<double> fitOverscanImage<float>(
+    afw::image::MaskedImage<float> const&, std::vector<std::string> badPixelMask, bool isTransposed);
 template
-void fitOverscanImage(
-     std::shared_ptr<afw::math::Function1<double> > &overscanFunction,
-    afw::image::MaskedImage<double> const& overscan,
-    double ssize,
-    int sigma);
+std::vector<double> fitOverscanImage<double>(
+    afw::image::MaskedImage<double> const&, std::vector<std::string> badPixelMask, bool isTransposed);
 
 template class CountMaskedPixels<float>;
 template class CountMaskedPixels<double>;
+template class CountMaskedPixels<int>;
 
 // Function to mask nans in a masked image
 template size_t maskNans<float>(afw::image::MaskedImage<float> const&, afw::image::MaskPixel,
                                 afw::image::MaskPixel);
 template size_t maskNans<double>(afw::image::MaskedImage<double> const&, afw::image::MaskPixel,
                                  afw::image::MaskPixel);
+template size_t maskNans<int>(afw::image::MaskedImage<int> const&, afw::image::MaskPixel,
+                              afw::image::MaskPixel);
 
 }}} // namespace lsst::ip::isr