Merge branch 'deconvolve-gaussians' into 'master'

Add function to deconvolve two 2D-Gaussian shapes

See merge request !153

include/schaapcommon/fitters/gaussianfitter.h

@@ -8,6 +8,23 @@
 
 namespace schaapcommon {
 namespace fitters {
+
+struct Ellipse {
+  double major;
+  double minor;
+  double position_angle;
+};
+
+/**
+ * Deconvolve two Gaussian kernels. The result convolved with @p smallest
+ * results again in @p largest. If @p smallest is larger than @p largest,
+ * an ellipse is returned with all parameters set to NaN.
+ *
+ * As the output axes are returned with the same units as the input axes,
+ * the input may be both in FWHM or both in sigma.
+ */
+Ellipse DeconvolveGaussian(const Ellipse& largest, const Ellipse& smallest);
+
 class GaussianFitter {
  public:
   void Fit2DGaussianCentred(const float* image, size_t width, size_t height,

src/fitters/gaussianfitter.cc

@@ -5,6 +5,7 @@
 
 #include <cmath>
 #include <iostream>
+#include <limits>
 
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_multifit_nlin.h>
@@ -22,7 +23,7 @@ namespace {
 const long double kSigmaToBeam = 2.0L * std::sqrt(2.0L * std::log(2.0L));
 
 void ToAnglesAndFwhm(double sx, double sy, double beta, double& ellipse_major,
-                     double& ellipse_minor, double& ellipse_phase_angle) {
+                     double& ellipse_minor, double& ellipse_position_angle) {
   const double beta_factor = 1.0 - beta * beta;
   double cov[4];
   cov[0] = sx * sx / beta_factor;
@@ -40,11 +41,11 @@ void ToAnglesAndFwhm(double sx, double sy, double beta, double& ellipse_major,
       vec1[0] = vec2[0];
       vec1[1] = vec2[1];
     }
-    ellipse_phase_angle = -std::atan2(vec1[0], vec1[1]);
+    ellipse_position_angle = -std::atan2(vec1[0], vec1[1]);
   } else {
     ellipse_major = std::sqrt(std::fabs(sx)) * kSigmaToBeam;
     ellipse_minor = std::sqrt(std::fabs(sx)) * kSigmaToBeam;
-    ellipse_phase_angle = 0.0;
+    ellipse_position_angle = 0.0;
   }
 }
 
@@ -54,7 +55,7 @@ void ToAnglesAndFwhm(double sx, double sy, double beta, double& ellipse_major,
  * @param y Y coordinate.
  * @param sx Sigma value for the x direction.
  * @param sy Sigma value for the y direction.
- * @param beta Beta value.
+ * @param beta Beta value; beta = 2 cov(x, y) / (sx * sy)
  */
 double GaussCentered(double x, double y, double sx, double sy, double beta) {
   return std::exp(-x * x / (2.0 * sx * sx) + beta * x * y / (sx * sy) -
@@ -787,5 +788,58 @@ void GaussianFitter::Fit2DGaussianWithAmplitudeWithFloor(
   beam_major *= scale_factor_;
   beam_minor *= scale_factor_;
 }
+
+Ellipse DeconvolveGaussian(const Ellipse& a, const Ellipse& b) {
+  /**
+   * These calculations can be derived from the following three steps necessary
+   * to deconvolve two Gaussian kernels. First, calculating the covariance
+   * matrix of both ellipses, then applying the formula for the variance of the
+   * sum of two Gaussians a = b + c (solving for the covariance of c), followed
+   * by converting the covariance matrix back to ellipse parameters.
+   * These steps can then be compacted to avoid having to first calculate the
+   * full parameters of a and b. The covariance can be calculated as:
+   * var_x = semiMajorAxis² * cos(phi)² + semiMinorAxis² * sin(phi)²
+   * var_y = semiMajorAxis² * sin(phi)² + semiMinorAxis² * cos(phi)²
+   * cov_xy = (semiMajorAxis² - semiMinorAxis²) * sin(phi) * cos(phi)
+   * beta = r = cov_xy / (sxsy)
+   */
+  const double cos_phi_a = std::cos(a.position_angle);
+  const double cos_square_a = cos_phi_a * cos_phi_a;
+  const double sin_phi_a = std::sin(a.position_angle);
+  const double sin_square_a = sin_phi_a * sin_phi_a;
+  const double cos_phi_b = std::cos(b.position_angle);
+  const double cos_square_b = cos_phi_b * cos_phi_b;
+  const double sin_phi_b = std::sin(b.position_angle);
+  const double sin_square_b = sin_phi_b * sin_phi_b;
+  const double a_maj_squared = a.major * a.major;
+  const double a_min_squared = a.minor * a.minor;
+  const double b_maj_squared = b.major * b.major;
+  const double b_min_squared = b.minor * b.minor;
+  const double cov_x =
+      (a_maj_squared * cos_square_a) + (a_min_squared * sin_square_a) -
+      (b_maj_squared * cos_square_b) - (b_min_squared * sin_square_b);
+  const double cov_y =
+      (a_maj_squared * sin_square_a) + (a_min_squared * cos_square_a) -
+      (b_maj_squared * sin_square_b) - (b_min_squared * cos_square_b);
+  const double cov_xy =
+      2.0 * ((a_min_squared - a_maj_squared) * sin_phi_a * cos_phi_a -
+             (b_min_squared - b_maj_squared) * sin_phi_b * cos_phi_b);
+  const double d = cov_x + cov_y;
+  const double r =
+      std::sqrt((cov_x - cov_y) * (cov_x - cov_y) + cov_xy * cov_xy);
+  Ellipse result;
+  const double epsilon = (std::fabs(cov_x) + std::fabs(cov_y)) * 1e-5;
+  if (cov_x < -epsilon || cov_y < -epsilon || d < r) {
+    result.major = std::numeric_limits<double>::quiet_NaN();
+    result.minor = std::numeric_limits<double>::quiet_NaN();
+    result.position_angle = std::numeric_limits<double>::quiet_NaN();
+  } else {
+    result.major = std::sqrt(0.5 * (d + r));
+    result.minor = std::sqrt(0.5 * (d - r));
+    result.position_angle = 0.5 * std::atan2(-cov_xy, cov_x - cov_y);
+  }
+  return result;
+}
+
 }  // namespace fitters
 }  // namespace schaapcommon

src/fitters/test/tgaussianfitter.cc

@@ -143,4 +143,80 @@ BOOST_AUTO_TEST_CASE(fit_small_beam) {
       1.0e-4);
 }
 
+BOOST_AUTO_TEST_CASE(deconvolve) {
+  using schaapcommon::fitters::DeconvolveGaussian;
+  using schaapcommon::fitters::Ellipse;
+  Ellipse result;
+
+  // Evaluate all four kwadrants.
+  for (double phi = 0.0; phi < 2.0 * M_PI; phi += M_PI / 4.0) {
+    // Simple case: equal PAs and zero minor axis for second ellipse, basically
+    // a 1D deconvolution.
+    const double input_phi = 0.5 + phi;
+    result = DeconvolveGaussian(Ellipse{2.0, 1.0, input_phi},
+                                Ellipse{0.2, 0.0, input_phi});
+    BOOST_CHECK_CLOSE_FRACTION(result.major, std::sqrt(2.0 * 2.0 - 0.2 * 0.2),
+                               1e-5);
+    BOOST_CHECK_CLOSE_FRACTION(result.minor, 1.0, 1e-5);
+    double expected = input_phi;
+    while (expected > M_PI * 0.5) expected -= M_PI;
+    BOOST_CHECK_CLOSE_FRACTION(result.position_angle, expected, 1e-5);
+
+    // Same as above, but rotate second ellipse by 180 deg (which should not
+    // change anything).
+    result = DeconvolveGaussian(Ellipse{2.0, 1.0, input_phi},
+                                Ellipse{0.2, 0.0, input_phi + M_PI});
+    BOOST_CHECK_CLOSE_FRACTION(result.major, std::sqrt(2.0 * 2.0 - 0.2 * 0.2),
+                               1e-5);
+    BOOST_CHECK_CLOSE_FRACTION(result.minor, 1.0, 1e-5);
+    BOOST_CHECK_CLOSE_FRACTION(result.position_angle, expected, 1e-5);
+  }
+
+  // Deconvolve with zero size ellipse: should return first ellipse.
+  result = DeconvolveGaussian(Ellipse{2.0, 1.0, 0.5}, Ellipse{0.0, 0.0, 0.1});
+  BOOST_CHECK_CLOSE_FRACTION(result.major, std::sqrt(2.0 * 2.0 - 0.0 * 0.0),
+                             1e-5);
+  BOOST_CHECK_CLOSE_FRACTION(result.minor, 1.0, 1e-5);
+  BOOST_CHECK_CLOSE_FRACTION(result.position_angle, 0.5, 1e-5);
+
+  // Equal sizes: should return zero-sized ellipse.
+  result = DeconvolveGaussian(Ellipse{2.0, 1.0, 0.5}, Ellipse{2.0, 1.0, 0.5});
+  BOOST_CHECK_LT(std::abs(result.major), 1e-5);
+  BOOST_CHECK_LT(std::abs(result.minor), 1e-5);
+  // position angle is somewhat undefined, but should still be finite...
+  BOOST_CHECK(std::isfinite(result.position_angle));
+  // ...and within bounds.
+  BOOST_CHECK_LE(result.position_angle, 2.0 * M_PI);
+  BOOST_CHECK_GE(result.position_angle, -2.0 * M_PI);
+
+  // 90 degree different position angle.
+  result =
+      DeconvolveGaussian(Ellipse{5.0, 2.0, 0.0}, Ellipse{1.0, 0.0, M_PI * 0.5});
+  BOOST_CHECK_CLOSE_FRACTION(result.major, 5.0, 1e-5);
+  BOOST_CHECK_CLOSE_FRACTION(result.minor, std::sqrt(2.0 * 2.0 - 1.0 * 1.0),
+                             1e-5);
+  BOOST_CHECK_LT(std::abs(result.position_angle), 1e-5);
+
+  // Circular deconvolution.
+  result = DeconvolveGaussian(Ellipse{5.0, 5.0, -0.3}, Ellipse{3.0, 3.0, 0.7});
+  BOOST_CHECK_CLOSE_FRACTION(result.major, 4.0, 1e-5);
+  BOOST_CHECK_CLOSE_FRACTION(result.minor, 4.0, 1e-5);
+  BOOST_CHECK(std::isfinite(result.position_angle));
+  BOOST_CHECK_LE(result.position_angle, 2.0 * M_PI);
+  BOOST_CHECK_GE(result.position_angle, -2.0 * M_PI);
+
+  // A complex case (numbers were calculated using the code, assuming it is
+  // correct).
+  result = DeconvolveGaussian(Ellipse{10.0, 8.0, 0.3}, Ellipse{7.0, 5.0, 1.4});
+  BOOST_CHECK_CLOSE_FRACTION(result.major, 8.477876113222349, 1e-5);
+  BOOST_CHECK_CLOSE_FRACTION(result.minor, 4.2574190079030156, 1e-5);
+  BOOST_CHECK_CLOSE_FRACTION(result.position_angle, 0.11532393547063115, 1e-5);
+
+  // Overflow situation.
+  result = DeconvolveGaussian(Ellipse{3.0, 3.0, 0.0}, Ellipse{4.0, 0.0, 0.0});
+  BOOST_CHECK(!std::isfinite(result.major));
+  BOOST_CHECK(!std::isfinite(result.major));
+  BOOST_CHECK(!std::isfinite(result.major));
+}
+
 BOOST_AUTO_TEST_SUITE_END()