Add a test for Gaussian draw function

include/schaapcommon/math/drawgaussian.h

@@ -11,11 +11,12 @@ namespace schaapcommon::math {
  * Draw a two-dimensional Gaussian onto an lm-space image.
  */
 void DrawGaussianToLm(float* image_data, size_t image_width,
-                      size_t image_height, long double ra, long double dec,
-                      long double pixel_scale_l, long double pixel_scale_m,
-                      long double l_shift, long double m_shift,
-                      long double position_ra, long double position_dec,
-                      const Ellipse& ellipse, long double integrated_flux);
+                      size_t image_height, long double phase_centre_ra,
+                      long double phase_centre_dec, long double pixel_scale_l,
+                      long double pixel_scale_m, long double l_shift,
+                      long double m_shift, long double source_ra,
+                      long double source_dec, const Ellipse& ellipse,
+                      long double integrated_flux);
 
 }  // namespace schaapcommon::math
 

src/math/drawgaussian.cc

@@ -17,11 +17,12 @@ inline long double Gaussian(long double x, long double sigma) {
 namespace schaapcommon::math {
 
 void DrawGaussianToLm(float* image_data, size_t image_width,
-                      size_t image_height, long double ra, long double dec,
-                      long double pixel_scale_l, long double pixel_scale_m,
-                      long double l_shift, long double m_shift,
-                      long double position_ra, long double position_dec,
-                      const Ellipse& ellipse, long double integrated_flux) {
+                      size_t image_height, long double phase_centre_ra,
+                      long double phase_centre_dec, long double pixel_scale_l,
+                      long double pixel_scale_m, long double l_shift,
+                      long double m_shift, long double source_ra,
+                      long double source_dec, const Ellipse& ellipse,
+                      long double integrated_flux) {
   // Using the FWHM formula for a Gaussian:
   const long double fwhm_constant = 2.0L * std::sqrt(2.0L * M_LN2);
   const long double sigma_major_axis = ellipse.major / fwhm_constant;
@@ -51,8 +52,8 @@ void DrawGaussianToLm(float* image_data, size_t image_width,
   const int bounding_box_size = std::ceil(sigma_max * 20.0 / min_pixel_scale);
   long double source_l;
   long double source_m;
-  ImageCoordinates::RaDecToLM(position_ra, position_dec, ra, dec, source_l,
-                              source_m);
+  ImageCoordinates::RaDecToLM(source_ra, source_dec, phase_centre_ra,
+                              phase_centre_dec, source_l, source_m);
 
   // Calculate the bounding box
   int source_x;

src/math/test/CMakeLists.txt

@@ -3,4 +3,5 @@
 
 include(${SCHAAPCOMMON_SOURCE_DIR}/cmake/unittest.cmake)
 
-add_unittest(math runtests.cc tconvolution.cc tresampler.cc trestoreimage.cc)
+add_unittest(math runtests.cc tconvolution.cc tdrawgaussian.cc tresampler.cc
+             trestoreimage.cc)

src/math/test/tdrawgaussian.cc

+#include <boost/test/unit_test.hpp>
+
+#include <aocommon/image.h>
+
+#include "drawgaussian.h"
+
+using schaapcommon::math::Ellipse;
+
+BOOST_AUTO_TEST_SUITE(draw_gaussian)
+
+BOOST_AUTO_TEST_CASE(to_lm) {
+  constexpr size_t kHeight = 200;
+  constexpr size_t kWidth = 500;
+  aocommon::Image image(kWidth, kHeight, 0.0f);
+  constexpr long double kRa = 20.0 * (M_PI / 180.0);   // 20 degrees
+  constexpr long double kDec = 60.0 * (M_PI / 180.0);  // 60 degrees
+  constexpr long double kPixelScaleL = 1.0 * (M_PI / 180.0 / 60.0);  // 1 amin
+  constexpr long double kPixelScaleM = 1.0 * (M_PI / 180.0 / 60.0);  // 1 amin
+  constexpr long double kLShift = 0.0;
+  constexpr long double kMShift = 0.0;
+  constexpr long double kSourceRa = kRa;
+  constexpr long double kSourceDec = kDec;
+  Ellipse ellipse;
+  ellipse.major = kPixelScaleL * 20.0;
+  ellipse.minor = kPixelScaleM * 5.0;
+  // Rotate Gaussian ten degrees to the east
+  ellipse.position_angle = 10.0 * (M_PI / 180.0);
+  constexpr long double kFlux = 150.0;
+
+  DrawGaussianToLm(image.Data(), kWidth, kHeight, kRa, kDec, kPixelScaleL,
+                   kPixelScaleM, kLShift, kMShift, kSourceRa, kSourceDec,
+                   ellipse, kFlux);
+
+  BOOST_CHECK_CLOSE_FRACTION(image.Sum(), kFlux, 1e-5);
+
+  constexpr size_t kCentralPixel = (kHeight / 2) * kWidth + kWidth / 2;
+  BOOST_CHECK_CLOSE_FRACTION(image[kCentralPixel], 1.32381356, 1e-5);
+
+  constexpr size_t kLeftOfCenterPixel = 90 * kWidth + 247;
+  BOOST_CHECK_CLOSE_FRACTION(image[kLeftOfCenterPixel], 0.0631099567, 1e-5);
+
+  constexpr size_t kRightOfCenterPixel = 90 * kWidth + 253;
+  BOOST_CHECK_CLOSE_FRACTION(image[kRightOfCenterPixel], 0.532994568, 1e-5);
+}
+
+BOOST_AUTO_TEST_SUITE_END()