Add function to draw Gaussian in pixel space

include/schaapcommon/math/drawgaussian.h

@@ -7,8 +7,17 @@
 
 namespace schaapcommon::math {
 
+/**
+ * Draw a two-dimensional Gaussian onto a gridded image. Units are given in
+ * pixels.
+ */
+void DrawGaussianToXy(float* image_data, size_t image_width,
+                      size_t image_height, double source_x, double source_y,
+                      const Ellipse& ellipse, long double integrated_flux);
+
 /**
  * Draw a two-dimensional Gaussian onto an lm-space image.
+ * As of yet, @p pixel_scale_l and @p pixel_scale_m must be equal.
  */
 void DrawGaussianToLm(float* image_data, size_t image_width,
                       size_t image_height, long double phase_centre_ra,

src/math/drawgaussian.cc

 #include "drawgaussian.h"
 
+#include <cassert>
 #include <cmath>
 
 #include <aocommon/imagecoordinates.h>
@@ -7,28 +8,21 @@
 using aocommon::ImageCoordinates;
 
 namespace {
-inline long double Gaussian(long double x, long double sigma) {
-  // Evaluate unnormalized Gaussian (unit valued peak, but not a unit integral.
-  const long double xi = x / sigma;
-  return std::exp(-0.5L * xi * xi);
+inline long double Gaussian(long double x) {
+  // Evaluate unnormalized Gaussian (unit valued peak, but not a unit integral).
+  return std::exp(-0.5L * x * x);
 }
 }  // namespace
 
 namespace schaapcommon::math {
 
-void DrawGaussianToLm(float* image_data, size_t image_width,
-                      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) {
+void DrawGaussianToXy(float* image_data, size_t image_width,
+                      size_t image_height, double source_x, double source_y,
+                      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;
   const long double sigma_minor_axis = ellipse.minor / fwhm_constant;
-  // TODO this won't work for non-equally spaced dimensions
-  const long double min_pixel_scale = std::min(pixel_scale_l, pixel_scale_m);
 
   // Position angle is angle from North:
   const long double angle = ellipse.position_angle + M_PI_2;
@@ -49,57 +43,79 @@ void DrawGaussianToLm(float* image_data, size_t image_width,
   transf[1] /= sigma_major_axis;
   transf[2] /= sigma_minor_axis;
   transf[3] /= sigma_minor_axis;
-  const int bounding_box_size = std::ceil(sigma_max * 20.0 / min_pixel_scale);
-  long double source_l;
-  long double 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;
-  int source_y;
-  ImageCoordinates::LMToXY<long double>(
-      source_l - l_shift, source_m - m_shift, pixel_scale_l, pixel_scale_m,
-      image_width, image_height, source_x, source_y);
+  const int bounding_box_size = std::ceil(sigma_max * 20.0);
   const int x_left =
-      std::clamp(source_x - bounding_box_size, 0, int(image_width));
+      std::clamp<int>(source_x - bounding_box_size, 0, image_width);
   const int x_right =
-      std::clamp(source_x + bounding_box_size, x_left, int(image_width));
+      std::clamp<int>(source_x + bounding_box_size, x_left, image_width);
   const int y_top =
-      std::clamp(source_y - bounding_box_size, 0, int(image_height));
+      std::clamp<int>(source_y - bounding_box_size, 0, image_height);
   const int y_bottom =
-      std::clamp(source_y + bounding_box_size, y_top, int(image_height));
+      std::clamp<int>(source_y + bounding_box_size, y_top, image_height);
 
   std::vector<double> values;
+  values.reserve((x_right - x_left) * (y_bottom - y_top));
   double flux_sum = 0.0;
   for (int y = y_top; y != y_bottom; ++y) {
     for (int x = x_left; x != x_right; ++x) {
-      long double l, m;
-      ImageCoordinates::XYToLM<long double>(x, y, pixel_scale_l, pixel_scale_m,
-                                            image_width, image_height, l, m);
-      l += l_shift;
-      m += m_shift;
-      const long double l_transf =
-          (l - source_l) * transf[0] + (m - source_m) * transf[1];
-      const long double m_transf =
-          (l - source_l) * transf[2] + (m - source_m) * transf[3];
+      const long double x_transf =
+          (x - source_x) * transf[0] + (y - source_y) * transf[1];
+      const long double y_transf =
+          (x - source_x) * transf[2] + (y - source_y) * transf[3];
       const long double dist =
-          std::sqrt(l_transf * l_transf + m_transf * m_transf);
-      long double v = Gaussian(dist, 1.0L);
+          std::sqrt(x_transf * x_transf + y_transf * y_transf);
+      long double v = Gaussian(dist);
       flux_sum += static_cast<double>(v);
       values.emplace_back(v);
     }
   }
   const double* iter = values.data();
+  // While the integral of a continuous Gaussian is known, the Gaussian that is
+  // drawn here is a sampled (gridded) Gaussian. Therefore, it is not accurate
+  // to use the standard Gaussian formula to normalize the Gaussian.
+  // TODO if the Gaussian cuts the side of the image, this leads to unexpected
+  // results.
   const double factor = integrated_flux / flux_sum;
   for (int y = y_top; y != y_bottom; ++y) {
     float* image_data_ptr = image_data + y * image_width + x_left;
     for (int x = x_left; x != x_right; ++x) {
-      (*image_data_ptr) += *iter * factor;
+      *image_data_ptr += *iter * factor;
       ++image_data_ptr;
       ++iter;
     }
   }
 }
 
+void DrawGaussianToLm(float* image_data, size_t image_width,
+                      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) {
+  assert(pixel_scale_l == pixel_scale_m);
+  long double source_l;
+  long double source_m;
+  ImageCoordinates::RaDecToLM(source_ra, source_dec, phase_centre_ra,
+                              phase_centre_dec, source_l, source_m);
+  long double source_x;
+  long double source_y;
+  ImageCoordinates::LMToXYfloat<long double>(
+      source_l - l_shift, source_m - m_shift, pixel_scale_l, pixel_scale_m,
+      image_width, image_height, source_x, source_y);
+
+  Ellipse xy_ellipse;
+  xy_ellipse.major = ellipse.major / pixel_scale_l;
+  xy_ellipse.minor = ellipse.minor / pixel_scale_l;
+
+  // Position angle is North through East. Because l increases to the left
+  // whereas x increases to the right (see e.g. ImageCoordinates::LMToXY()),
+  // the sign is flipped here.
+  xy_ellipse.position_angle = -ellipse.position_angle;
+  DrawGaussianToXy(image_data, image_width, image_height, source_x, source_y,
+                   xy_ellipse, integrated_flux);
+}
+
 }  // namespace schaapcommon::math

src/math/restoreimage.cc

@@ -23,7 +23,8 @@ void RestoreImage(float* image_data, const float* model_data,
       image_data[j] += model_data[j];
     }
   } else {
-    // Using the FWHM formula for a Gaussian:
+    // TODO this can make use of the Gaussian drawing functions in
+    // schaapcommon::math Using the FWHM formula for a Gaussian:
     const long double sigma_major =
         beam_major_axis / (2.0L * std::sqrt(2.0L * std::log(2.0L)));
     const long double sigma_minor =