diff --git a/include/schaapcommon/fitters/gaussianfitter.h b/include/schaapcommon/fitters/gaussianfitter.h
index 8e088d035efc4baa9dd04da443ca01f4f1581a1b..8dcf774262f7f6af689f57cb81a1c771a434130d 100644
--- a/include/schaapcommon/fitters/gaussianfitter.h
+++ b/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,
diff --git a/src/fitters/gaussianfitter.cc b/src/fitters/gaussianfitter.cc
index 573e05375f441867083aff4500db9375347f4e66..f94a219b3cdbba58ab5f3ac02059ace8e2a66b19 100644
--- a/src/fitters/gaussianfitter.cc
+++ b/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
diff --git a/src/fitters/test/tgaussianfitter.cc b/src/fitters/test/tgaussianfitter.cc
index 73d99272adb4bfa2af7cbfe2d4a4890c97785607..0f8e3a6916cd84d27e61ea0074be406b87e2e00c 100644
--- a/src/fitters/test/tgaussianfitter.cc
+++ b/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()