DM-54397: Add task for computing a rough shapelet expansion of isolated stars.

include/lsst/meas/algorithms/SpanSetMoments.h

@@ -0,0 +1,132 @@
+// -*- lsst-c++ -*-
+/*
+ * LSST Data Management System
+ * Copyright 2008-2013 LSST Corporation.
+ *
+ * This product includes software developed by the
+ * LSST Project (http://www.lsst.org/).
+ *
+ * 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 LSST License Statement and
+ * the GNU General Public License along with this program.  If not,
+ * see <http://www.lsstcorp.org/LegalNotices/>.
+ */
+
+#ifndef LSST_MEAS_ALGORITHMS_SpanSetMoments_h_INCLUDED
+#define LSST_MEAS_ALGORITHMS_SpanSetMoments_h_INCLUDED
+
+#include "ndarray.h"
+#include "lsst/geom/Point.h"
+#include "lsst/afw/image/MaskedImage.h"
+#include "lsst/afw/geom/ellipses/Quadrupole.h"
+#include "lsst/afw/geom/SpanSet.h"
+#include "lsst/shapelet/ShapeletFunction.h"
+
+namespace lsst {
+namespace meas {
+namespace algorithms {
+
+/**
+ *  A struct that computes the unweighted moments of the pixels in an
+ * `lsst.afw.geom.SpanSet`.
+ *
+ * This class provides low-level pixel-processing code for
+ * `ComputeRoughPsfShapeletsTask`.
+ */
+struct SpanSetMoments {
+
+    /// Total flux within the SpanSet.
+    double flux;
+
+    /// Total variance within the SpanSet.
+    double variance;
+
+    /// Center derived from the unweighted first moment of the image.
+    geom::Point2D center;
+
+    /// Shape derived from the unweighted second moment of the image.
+    afw::geom::ellipses::Quadrupole shape;
+
+    /// The pixels actually used to compute the moments.
+    std::shared_ptr<afw::geom::SpanSet> spans;
+
+    /// Flag set if there were too many bad pixels to compute the moments.
+    bool too_many_bad_pixels = false;
+
+    /// Flag set if the center did not lie within the SpanSet.
+    bool center_out_of_bounds = false;
+
+    /// Flag set if there was a bad pixel too close to the center.
+    bool bad_pixel_in_center = false;
+
+    /// Flag set if the second moments did not resolve to an ellipse.
+    bool singular_second_moments = false;
+
+    /// Test whether any failure flag is set.
+    bool any_flags_set() const {
+        return too_many_bad_pixels || center_out_of_bounds || bad_pixel_in_center || singular_second_moments;
+    }
+
+    /// Return a flattened array of the x coordinates in `spans`.
+    ndarray::Array<float, 1, 1> get_x_array() const;
+
+    /// Return a flattened array of the y coordinates in `spans`.
+    ndarray::Array<float, 1, 1> get_y_array() const;
+
+    /**
+     * Compute the unweighted moments of an image within a SpanSet.
+     *
+     * @param[in] spans          Pixel region to use.
+     * @param[in] masked_image   Image to measure.
+     * @param[in] bad_bitmask    Mask of bad pixels to remove from `spans`
+     *                           before computing the moments.
+     * @param[in] bad_pixel_mask_fraction     Maximum fraction of the pixels
+     *                                        in `spans` that can be bad
+     *                                        before giving up.
+     * @param[in] bad_pixel_exlusion_radius   Radius around the estimated
+     *                                        center where the present of a
+     *                                        bad pixels will cause the
+     *                                        algorithm to give up.
+     */
+    static std::shared_ptr<SpanSetMoments> compute(afw::geom::SpanSet const& spans,
+                                                   afw::image::MaskedImage<float> const& masked_image,
+                                                   afw::image::MaskPixel bad_bitmask,
+                                                   double bad_pixel_max_fraction,
+                                                   double bad_pixel_exclusion_radius);
+
+    /**
+     * Fit a common shapelet expansion to multiple sources whose moments have
+     * already been computed.
+     *
+     * The expected use case is a sample of stars with similar but not
+     * identical moments, which might plausibly have an approximately common
+     * shapelet representation while benefitting from allowing the ellipse
+     * used for the expansion to vary (i.e. to partially capture spatial
+     * variation in the PSF).
+     *
+     * @param[in] masked_image    Image ot measure.
+     * @param[in] moments         Moments already measured on stars.
+     * @param[in] order           Order of the shapelet expansion.
+     * @param[in] scale           Factor to scale the moments ellipses by to
+     *                            yield the ellipse for the shapelet fit.
+     */
+    static shapelet::ShapeletFunction fit_shapelets(
+            afw::image::MaskedImage<float> const& masked_image,
+            std::vector<std::shared_ptr<SpanSetMoments>> const& moments,
+            int order, double scale);
+};
+
+}  // namespace algorithms
+}  // namespace meas
+}  // namespace lsst
+
+#endif  // !LSST_MEAS_ALGORITHMS_SpanSetMoments_h_INCLUDED

python/lsst/meas/algorithms/SConscript

@@ -15,5 +15,6 @@ scripts.BasicSConscript.python(['_algorithmsLib'], [
     'pcaPsf.cc',
     'psfCandidate/psfCandidate.cc',
     'singleGaussianPsf.cc',
+    'spanSetMoments.cc',
     'spatialModelPsf.cc',
     'warpedPsf.cc'])

python/lsst/meas/algorithms/_algorithmsLib.cc

@@ -45,6 +45,7 @@ void wrapCoaddPsf(WrapperCollection &wrappers);
 void wrapCoaddBoundedField(WrapperCollection &wrappers);
 void wrapCr(WrapperCollection &wrappers);
 void wrapCloughTocher2DInterpolatorUtils(WrapperCollection &wrappers);
+void wrapSpanSetMoments(WrapperCollection & wrappers);
 
 PYBIND11_MODULE(_algorithmsLib, mod) {
 WrapperCollection wrappers(mod, "lsst.meas.algorithms");
@@ -67,9 +68,10 @@ WrapperCollection wrappers(mod, "lsst.meas.algorithms");
     wrapCoaddBoundedField(wrappers);
     wrapCr(wrappers);
     wrapCloughTocher2DInterpolatorUtils(wrappers);
+    wrapSpanSetMoments(wrappers);
     wrappers.finish();
 }
 
 }  // algorithms
 }  // meas
-}  // lsst
+}  // lsst

python/lsst/meas/algorithms/computePsfChiSq.py

@@ -0,0 +1,90 @@
+# This file is part of meas_algorithms.
+#
+# 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/>.
+
+__all__ = ["ComputePsfChiSqTask", "ComputePsfChiSqConfig"]
+
+import numpy as np
+
+from lsst.afw.detection import Psf
+from lsst.afw.image import ImageF, MaskedImageF
+from lsst.geom import Point2D
+from lsst.pex.config import Config, Field, ListField
+from lsst.pipe.base import Struct, Task
+
+
+class ComputePsfChiSqConfig(Config):
+    padding = Field(
+        "Number of pixels to grow the PSF model bbox by on all sides.",
+        dtype=int,
+        default=10,
+    )
+    bad_mask_planes = ListField(
+        "Mask planes to identify pixels to drop from the calculation.",
+        dtype=str,
+        default=["SAT", "SUSPECT"],
+    )
+
+
+class ComputePsfChiSqTask(Task):
+    ConfigClass = ComputePsfChiSqConfig
+    _DefaultName = "computePsfChiSq"
+    config: ComputePsfChiSqConfig
+
+    def run(
+        self, *, masked_image: MaskedImageF, psf: Psf, x: np.ndarray, y: np.ndarray
+    ) -> Struct:
+        bitmask = masked_image.mask.getPlaneBitMask(self.config.bad_mask_planes)
+        result = Struct(
+            chisq=np.zeros(x.shape, dtype=np.float32),
+            snr=np.zeros(x.shape, dtype=np.float32),
+            n_missing_pixels=np.zeros(x.shape, dtype=np.uint16),
+        )
+        for i, (sx, sy) in enumerate(zip(x, y)):
+            inner_model = psf.computeImage(Point2D(sx, sy))
+            inner_bbox = inner_model.getBBox()
+            outer_bbox = inner_bbox.dilatedBy(self.config.padding)
+            n_clipped_pixels = 0
+            if not masked_image.getBBox().contains(outer_bbox):
+                original_bbox_area = outer_bbox.getArea()
+                outer_bbox.clip(masked_image.getBBox())
+                n_clipped_pixels = original_bbox_area - outer_bbox.getArea()
+            outer_model = ImageF(outer_bbox)
+            outer_model.array[:, :] = 0
+            outer_model[inner_bbox] = inner_model.convertF()
+            outer_data = masked_image.image[outer_bbox]
+            outer_variance = masked_image.variance[outer_bbox]
+            mask = np.logical_not(masked_image.mask[outer_bbox].array & bitmask)
+            flat_model = outer_model.array[mask]
+            flat_data = outer_data.array[mask]
+            flat_variance = outer_variance.array[mask]
+            sigma = flat_variance**0.5
+            weighted_model = flat_model / sigma
+            weighted_data = flat_data / sigma
+            inv_var_alpha = np.dot(weighted_model, weighted_model)
+            alpha = np.dot(weighted_model, weighted_data) / inv_var_alpha
+            weighted_residuals = weighted_model * alpha
+            weighted_residuals -= weighted_data
+            result.chisq[i] = sum(weighted_residuals**2) / np.sum(mask)
+            result.snr[i] = alpha * inv_var_alpha**0.5
+            result.n_missing_pixels[i] = (
+                outer_bbox.getArea() - np.sum(mask) + n_clipped_pixels
+            )
+        return result