Merge pull request #80 from lsst/tickets/DM-10237

DM-10237: fix errors in blendedness calculation

include/lsst/meas/base/Blendedness.h

@@ -87,6 +87,41 @@ class BlendednessAlgorithm : public SimpleAlgorithm {
 
     BlendednessAlgorithm(Control const & ctrl, std::string const & name, afw::table::Schema & schema);
 
+    /**
+     *  Compute the posterior expectation value of the true flux in a pixel
+     *  from its (Gaussian) likelihood and a flat nonnegative prior.
+     *
+     *  This computes
+     *  @f[
+     *       \frac{\int_0^\infty f \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(f-z)^2}{2\sigma^2}} df}
+                  {\int_0^\infty \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(f-z)^2}{2\sigma^2}} df}
+     *  @f]
+     *  where @f$z@f$ is the (noisy) pixel value and @f$\sigma^2@f$ is the pixel variance.  This
+     *  approaches @f$z@f$ when @f$z \gg \sigma@f$ and @f$0@f$ when @f$z \ll -\sigma@f$.
+     *
+     *  We use single precision here for performance reasons; this function is called in a loop
+     *  over single-precision pixels, and relies on a number of calls to exp and erfc, which are
+     *  much faster in single precision.
+     */
+    static float computeAbsExpectation(float data, float variance);
+
+    /**
+     *  Compute the bias induced by using the absolute value of a pixel instead of its value.
+     *
+     *  The computation assumes the true distribution for the pixel is a Gaussian
+     *  with mean mu and the given variance.  To compute mu from the data, use
+     *  computeAbsExpectation.
+     *
+     *  This computes
+     *  @f[
+     *      \sqrt{\frac{2}{\pi}}\sigma e^{-\frac{\mu^2}{2\sigma^2}}
+     *          - \mu\,\mathrm{erfc}\left(\frac{\mu}{\sqrt{2}\sigma}\right)
+     *  @f]
+     *  where @f$\mu@f$ is the mean of the underlying distribution and @f$\sigma^2@f$
+     *  is its variance.
+     */
+    static float computeAbsBias(float mu, float variance);
+
     void measureChildPixels(
         afw::image::MaskedImage<float> const & image,
         afw::table::SourceRecord & child

python/lsst/meas/base/blendedness.cc

@@ -65,6 +65,10 @@ PyBlendenessAlgorithm declareBlendednessAlgorithm(py::module &mod) {
     cls.attr("NO_CENTROID") = py::cast(BlendednessAlgorithm::NO_CENTROID);
     cls.attr("NO_SHAPE") = py::cast(BlendednessAlgorithm::NO_SHAPE);
 
+    cls.def_static("computeAbsExpectation", &BlendednessAlgorithm::computeAbsExpectation,
+                   "data"_a, "variance"_a);
+    cls.def_static("computeAbsBias", &BlendednessAlgorithm::computeAbsBias,
+                   "mu"_a, "variance"_a);
     cls.def("measureChildPixels", &BlendednessAlgorithm::measureChildPixels, "image"_a, "child"_a);
     cls.def("measureParentPixels", &BlendednessAlgorithm::measureParentPixels, "image"_a, "child"_a);
     cls.def("measure", &BlendednessAlgorithm::measure, "measRecord"_a, "exposure"_a);

src/Blendedness.cc

@@ -23,8 +23,7 @@
 
 #include <cmath>
 
-#include "boost/math/special_functions/erf.hpp"
-#include <boost/math/constants/constants.hpp>
+#include "boost/math/constants/constants.hpp"
 
 #include "lsst/meas/base/Blendedness.h"
 #include "lsst/afw/detection/HeavyFootprint.h"
@@ -154,7 +153,6 @@ class ShapeAccumulator : public FluxAccumulator {
     double _wdxy;
 };
 
-
 template <typename Accumulator>
 void computeMoments(
     afw::image::MaskedImage<float> const & image,
@@ -209,13 +207,9 @@ void computeMoments(
             float data = pixelIter.image();
             accumulatorRaw(d.getX(), d.getY(), weight, data);
             float variance = pixelIter.variance();
-            float mu = (std::sqrt(variance/(2.0f/boost::math::constants::pi<float>()))*
-                        std::exp(-0.5f*(data*data)/variance)) +
-                0.5f*data*boost::math::erfc(-data/std::sqrt(2.0f*variance));
-            float bias = (std::sqrt(2.0f*variance/boost::math::constants::pi<float>())*
-                          std::exp(-0.5f*(mu*mu)/variance)) -
-                mu*boost::math::erfc(mu/std::sqrt(2.0f*variance));
-            accumulatorAbs(d.getX(), d.getY(), weight, std::max(std::abs(data) - bias, 0.0f));
+            float mu = BlendednessAlgorithm::computeAbsExpectation(data, variance);
+            float bias = BlendednessAlgorithm::computeAbsBias(mu, variance);
+            accumulatorAbs(d.getX(), d.getY(), weight, std::abs(data) - bias);
         }
     }
 }
@@ -284,6 +278,22 @@ BlendednessAlgorithm::BlendednessAlgorithm(Control const & ctrl,  std::string co
     }
 }
 
+float BlendednessAlgorithm::computeAbsExpectation(float data, float variance) {
+    float normalization = 0.5f*std::erfc(-data/std::sqrt(2.0f*variance));
+    if (!(normalization > 0)) {
+        // avoid division by zero; we know the limit at data << -sigma -> 0.
+        return 0.0;
+    }
+    return data + (std::sqrt(0.5f*variance/boost::math::constants::pi<float>()) *
+                   std::exp(-0.5f*(data*data)/variance) / normalization);
+}
+
+float BlendednessAlgorithm::computeAbsBias(float mu, float variance) {
+    return (std::sqrt(2.0f*variance/boost::math::constants::pi<float>()) *
+                      std::exp(-0.5f*(mu*mu)/variance)) -
+            mu*std::erfc(mu/std::sqrt(2.0f*variance));
+}
+
 void BlendednessAlgorithm::_measureMoments(
     afw::image::MaskedImage<float> const & image,
     afw::table::SourceRecord & child,
@@ -349,7 +359,7 @@ void BlendednessAlgorithm::_measureMoments(
             );
         if (_ctrl.doFlux) {
             child.set(fluxRawKey, accumulatorRaw.getFlux());
-            child.set(fluxAbsKey, accumulatorAbs.getFlux());
+            child.set(fluxAbsKey, std::max(accumulatorAbs.getFlux(), 0.0));
         }
         _shapeRawKey.set(child, accumulatorRaw.getShape());
         _shapeAbsKey.set(child, accumulatorAbs.getShape());
@@ -365,7 +375,7 @@ void BlendednessAlgorithm::_measureMoments(
             accumulatorAbs
         );
         child.set(fluxRawKey, accumulatorRaw.getFlux());
-        child.set(fluxAbsKey, accumulatorAbs.getFlux());
+        child.set(fluxAbsKey, std::max(accumulatorAbs.getFlux(), 0.0));
     }
 }
 

tests/testBlendedness.py

@@ -64,6 +64,26 @@ def tearDown(self):
         del self.bbox
         del self.dataset
 
+    def testAbsExpectation(self):
+        f = lsst.meas.base.BlendednessAlgorithm.computeAbsExpectation
+        # comparison values computed with Mathematica
+        self.assertFloatsAlmostEqual(f(-1.0, 1.5**2), 0.897767011, rtol=1E-5)
+        self.assertFloatsAlmostEqual(f(0.0, 1.5**2), 1.19682684, rtol=1E-5)
+        self.assertFloatsAlmostEqual(f(1.0, 1.5**2), 1.64102639, rtol=1E-5)
+        self.assertFloatsAlmostEqual(f(-1.0, 0.3**2), 0.0783651228, rtol=1E-5)
+        self.assertFloatsAlmostEqual(f(0.0, 0.3**2), 0.239365368, rtol=1E-5)
+        self.assertFloatsAlmostEqual(f(1.0, 0.3**2), 1.00046288, rtol=1E-5)
+
+    def testAbsBias(self):
+        f = lsst.meas.base.BlendednessAlgorithm.computeAbsBias
+        # comparison values computed with Mathematica
+        self.assertFloatsAlmostEqual(f(0.0, 1.5**2), 1.19682684, rtol=1E-5)
+        self.assertFloatsAlmostEqual(f(0.5, 1.5**2), 0.762708343, rtol=1E-5)
+        self.assertFloatsAlmostEqual(f(1.0, 1.5**2), 0.453358941, rtol=1E-5)
+        self.assertFloatsAlmostEqual(f(0.0, 0.3**2), 0.239365368, rtol=1E-5)
+        self.assertFloatsAlmostEqual(f(0.5, 0.3**2), 0.011895931, rtol=1E-5)
+        self.assertFloatsAlmostEqual(f(1.0, 0.3**2), 0.0000672467314, rtol=1E-5)
+
     def testBlendedness(self):
         """
         Check that we measure a positive blendedness for two overlapping sources