Tickets/DM-13806

include/lsst/meas/modelfit/RegularizedMoments.h

@@ -0,0 +1,132 @@
+// -*- lsst-c++ -*-
+/*
+ * This file is part of package_name.
+ *
+ * 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/>.
+ */
+
+#ifndef LSST_MEAS_MODELFIT_REGULARIZEDMOMENTS_H
+#define LSST_MEAS_MODELFIT_REGULARIZEDMOMENTS_H
+
+#include <Eigen/Dense>
+
+namespace lsst {
+namespace meas {
+namespace modelfit {
+
+/**
+ *  @brief A model for calculating observed shape moments and their derivatives
+ *
+ *  This class models the biased weighted moments of an astronomical source given a vector of
+ *  moments intrinsic to an object, and a vector of moments representing a weighting function.
+ *
+ *  The constructor takes the moments of weight function, and the moments of the object
+ *  are set by calling the setParameters method. This allows a single instance to be reused to evaluate
+ *  multiple positions in object moment space. This structure is computationally more efficient
+ *  as well since the weight function would remain constant for many evaluations of image moments
+ *
+ *  An instance of this class is designed to be used as an objective function in an optimizer
+ *  which will find the maximally likely image moments given moments measured from an image.
+ *  As such the class provides two methods in addition to the setParameters method already mentioned.
+ *  One method computes the biased weighted moments given the specified weight and "true" moments,
+ *  and the other computes the derivative of the objective function along each of the image moment
+ *  variables.
+ */
+class MomentsModel {
+public:
+    using Element = double;
+    using Moments = Eigen::Matrix<Element, 6, 1>;
+    using Jacobian = Eigen::Matrix<Element, 6, 6>;
+    using FirstMoment = Eigen::Matrix<Element, 2, 1>;
+    using SecondMoment = Eigen::Matrix<Element, 2, 2>;
+
+    /**
+     * Construct a MomentsModel object with a given vector of weight function moments
+     *
+     * @param[in]   W       An eigen vector of the moments of the elliptical Gaussian weighting function.
+     *                      Most likely the first element (the zeroth moment) should be one, representing
+     *                      a normalized weight function.
+     */
+    explicit MomentsModel(Moments const &W) : W(W) {}
+
+    /*
+     * default copy and move constructors
+     */
+    MomentsModel(const MomentsModel &) = default;
+    MomentsModel(MomentsModel &&) = default;
+
+    /**
+     *  Sets the moments corresponding to the astronomical object's intrinsic moments
+     *
+     *  When this function is called, it sets the astronomical object's intrinsic moments for all subsequent
+     *  calculations. Additionally after setting the moment, this function calculates and caches the value
+     *  of the biased weighted moments given the weights and intrinsic moments.
+     *
+     *  @param[in]  Q       An eigen vector of the intrinsic moments of the astronomical object
+     */
+    void setParameters(Moments const &Q);
+
+    /**
+     *  Returns the value of the biased weighted moments given the input weight and object's intrinsic moments
+     *
+     *  @param[out]         Eigen vector of the model biased weighted moments given the input weight and
+                            intrinsic moments
+     */
+    Moments computeValues() const;
+
+    /**
+     *  Computes and returns the gradient of the model biased weighted moments along each of the intrinsic
+     *  moment dimensions.
+     *
+     *  @param[out]         A six by six Eigen matrix, with the biased weighted moments down the column, and
+     *                      the derivative of the column moment with respect to the intrinsic moments along
+     *                      the rows.
+     */
+    Jacobian computeJacobian() const;
+
+private:
+    // Calculates the Model evaluated at the input intrinsic moments given the specified weighted moments
+    void evaluateModel();
+
+    // Storage variables
+    Moments W, Q, values;
+
+    Element scalar;
+
+    FirstMoment alpha;
+
+    SecondMoment beta;
+};
+
+// Tests for classes in anonymous name spaces
+double constexpr DEFAULT_TEST_TOLERANCE = 1.0e-6;
+bool testScalar(double tol = DEFAULT_TEST_TOLERANCE);
+
+bool testAlphaX(double tol = DEFAULT_TEST_TOLERANCE);
+bool testAlphaY(double tol = DEFAULT_TEST_TOLERANCE);
+
+bool testBetaX(double tol = DEFAULT_TEST_TOLERANCE);
+bool testBetaY(double tol = DEFAULT_TEST_TOLERANCE);
+bool testBetaXY(double tol = DEFAULT_TEST_TOLERANCE);
+}  // namespace modelfit
+}  // namespace meas
+}  // namespace lsst
+
+#endif  // LSST_MEAS_MODELFIT_REGULARIZEDMOMENTS_H

python/lsst/meas/modelfit/SConscript

@@ -16,6 +16,7 @@ scripts.BasicSConscript.pybind11(
      'integrals',
      'unitTransformedLikelihood',
      'cmodel/cmodel',
+     'regularizedMoments/regularizedMoments',
      ],
     addUnderscore=False
 )

python/lsst/meas/modelfit/regularizedMoments/__init__.py

@@ -0,0 +1,2 @@
+from .regularizedMoments import *
+from .regularizedMomentsContinued import *

python/lsst/meas/modelfit/regularizedMoments/regularizedMoments.cc

@@ -0,0 +1,71 @@
+/*
+ * This file is part of meas_modelfit.
+ *
+ * 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/>.
+ */
+
+#include "pybind11/pybind11.h"
+#include "pybind11/stl.h"
+
+#include "numpy/arrayobject.h"
+#include "ndarray/pybind11.h"
+#include "ndarray/eigen.h"
+
+#include "lsst/meas/modelfit/RegularizedMoments.h"
+
+namespace py = pybind11;
+using namespace py::literals;
+
+namespace lsst {
+namespace meas {
+namespace modelfit {
+
+namespace {
+using PyMomentsModel = py::class_<MomentsModel, std::shared_ptr<MomentsModel>>;
+
+PYBIND11_PLUGIN(regularizedMoments) {
+    if (_import_array() < 0) {
+        PyErr_SetString(PyExc_ImportError, "numpy.core.multiarray failed to import");
+        return nullptr;
+    }
+
+    py::module mod("regularizedMoments");
+    PyMomentsModel cls(mod, "MomentsModel");
+
+    cls.def(py::init<MomentsModel::Moments const &>(), "weights"_a);
+
+    cls.def("setParameters", &MomentsModel::setParameters);
+    cls.def("computeValues", &MomentsModel::computeValues);
+    cls.def("computeJacobian", &MomentsModel::computeJacobian);
+
+    // Wrap the test functions for the anonymous classes
+    mod.def("testScalar", &testScalar, "tol"_a = DEFAULT_TEST_TOLERANCE);
+    mod.def("testAlphaX", &testAlphaX, "tol"_a = DEFAULT_TEST_TOLERANCE);
+    mod.def("testAlphaY", &testAlphaY, "tol"_a = DEFAULT_TEST_TOLERANCE);
+    mod.def("testBetaX", &testBetaX, "tol"_a = DEFAULT_TEST_TOLERANCE);
+    mod.def("testBetaXY", &testBetaXY, "tol"_a = DEFAULT_TEST_TOLERANCE);
+    mod.def("testBetaY", &testBetaY, "tol"_a = DEFAULT_TEST_TOLERANCE);
+
+    return mod.ptr();
+}
+}  // namespace
+}  // namespace modelfit
+}  // namespace meas
+}  // namespace lsst

python/lsst/meas/modelfit/regularizedMoments/regularizedMomentsContinued.py

@@ -0,0 +1,139 @@
+# This file is part of meas modelfit.
+#
+# 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 numpy as np
+
+
+def makeGaussian(x, y, scale, muX, muY, varX, varXY, varY):
+    """Create an elliptical Gaussian.
+
+    Parameters
+    ----------
+    x : 2D numpy array
+        An array containing the x coordinates the Gaussian will be evaluated on.
+        most likely the result of a numpy.indices call
+    y : 2D numpy array
+        An array containing the y coordinates the Gaussian will be evaluated on.
+        most likely the result of a numpy.indices call
+    scale : `float`
+        The value the resulting Gaussian will have when summed over all pixels.
+    muX : `float`
+        The central position of the Gaussian in the x direction
+    muY : `float`
+        The central position of the Gaussian in the y direction
+    varX : `float`
+        The variance of the Gaussian about the muX position
+    varXY : `float`
+        The covariance of the Gaussian in x and y
+    varY : `float`
+        The variance of the Gaussian about the muY position
+
+    Returns
+    -------
+    Gaussian : 2D numpy array
+        The Gaussian array generated from the input values
+    """
+
+    rho = varXY/(varX**0.5*varY**0.5)
+
+    psf = np.exp(-1/(2*(1-rho**2)) *
+                 ((x-muX)**2/varX+(y - muY)**2/varY -
+                 2*rho*(x-muX)*(y-muY)/(varX**0.5*varY**0.5)))
+
+    psf /= psf.sum()
+    return scale*psf
+
+
+def buildUncertanty(imShape, W, uncertanty):
+    """ Propagate pixel uncertainties to uncertainties in weighted moments
+
+    Parameters
+    ----------
+    imShape : tuple(float, float)
+        The shape of image for which weighted moments have been calculated
+    W : iterable
+        An iterable object with six elements corresponding to the moments used
+        in the weighted moment calculation, scale, mean in x, mean in y, variance
+        in x, covariance between x and y, and variance in y.
+    uncertanty : `float`
+        Uncertainty in the pixel value. This is a single value, as this routine
+        assumes errors are background dominated, and uncorrelated
+
+    Returns
+    -------
+    covarianceMatrix : 2D 6x6 numpy array of floats
+        This is the covariance matrix on the measured moments with uncertainties
+        propagated from pixel uncertainties
+    """
+    yInd, xInd = np.indices(imShape)
+    weightImage = makeGaussian(xInd, yInd, *W)
+    sigmaImage = np.eye(weightImage.size)*uncertanty
+    MomentWeightMatrix = np.zeros((6, weightImage.size))
+
+    weightImageFlat = weightImage.flatten()
+    xIndFlat = xInd.flatten()
+    yIndFlat = yInd.flatten()
+
+    MomentWeightMatrix[0] = weightImageFlat
+    MomentWeightMatrix[1] = weightImageFlat*xIndFlat
+    MomentWeightMatrix[2] = weightImageFlat*yIndFlat
+    MomentWeightMatrix[3] = weightImageFlat*xIndFlat**2
+    MomentWeightMatrix[4] = weightImageFlat*xIndFlat*yIndFlat
+    MomentWeightMatrix[5] = weightImageFlat*yIndFlat**2
+
+    return np.dot(MomentWeightMatrix, np.dot(sigmaImage, np.transpose(MomentWeightMatrix)))
+
+
+def measureMoments(image, W):
+    """ Calculate weighted moments of the input image with the given weight array
+
+    Parameters
+    ----------
+    image : 2D numpy array of floats
+        This is the input postage stamp of a source for which the weighted moments are
+        to be measured
+    W : 2D numpy array of floats
+        Array of floats that are used as weights when calculating moments on the input
+        image. Array must be the same shape image
+
+    Returns
+    -------
+    moments : 6 element numpy array
+        These are the weighted moments as measured from the input image in the order of
+        0th, 1st X, 1st Y, 2nd X, 2nd XY, 2nd Y
+
+    Raises
+    ------
+    AssertionError: Raises if the input arrays are not the same shape
+    """
+    assert image.shape == W.shape, "Input image and weight array must be the same shape"
+
+    yInd, xInd = np.indices(image.shape)
+    weightImage = makeGaussian(xInd, yInd, *W)
+
+    zero = np.sum(image*weightImage)
+    oneX = np.sum(image*weightImage*xInd)
+    oneY = np.sum(image*weightImage*yInd)
+    twoX = np.sum(image*weightImage*xInd**2)
+    twoXY = np.sum(image*weightImage*xInd*yInd)
+    twoY = np.sum(image*weightImage*yInd**2)
+
+    return np.array((zero, oneX, oneY, twoX, twoXY, twoY)), weightImage