Merge branch 'move-draw-gaussian' into 'master'

Move Gaussian drawing function from WSClean

See merge request !155

CMakeLists.txt

@@ -4,7 +4,7 @@
 # Top level CMakeLists.txt file for SchaapCommon
 cmake_minimum_required(VERSION 3.13) # CMP0079 needs CMake 3.13
 
-set(SCHAAPCOMMON_ALL_MODULES facets fft fitters h5parm ducc0)
+set(SCHAAPCOMMON_ALL_MODULES facets fitters h5parm math ducc0)
 
 # Define a couple of dedicated variables, to prevent that standard CMAKE
 # variables get screwed up when schaapcommon included in another package

include/schaapcommon/fitters/gaussianfitter.h

@@ -6,15 +6,11 @@
 
 #include <cstring>
 
+#include "../math/ellipse.h"
+
 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,
@@ -23,16 +19,18 @@ struct Ellipse {
  * 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);
+math::Ellipse DeconvolveGaussian(const math::Ellipse& largest,
+                                 const math::Ellipse& smallest);
 
 /**
  * Fit Gaussian parameters to the shape in the center of the image.
  * @param beam_est provides the initial value for the beam's size. It must be
  * somewhat close and must be non-zero for the algorithm to converge (quickly).
  */
-Ellipse Fit2DGaussianCentred(const float* image, size_t width, size_t height,
-                             double beam_est, double box_scale_factor = 10.0,
-                             bool verbose = false);
+math::Ellipse Fit2DGaussianCentred(const float* image, size_t width,
+                                   size_t height, double beam_est,
+                                   double box_scale_factor = 10.0,
+                                   bool verbose = false);
 
 /**
  * Fits a circular (one parameter) Gaussian to the shape in the centre of the

include/schaapcommon/fft/compositefft.h → include/schaapcommon/math/compositefft.h

@@ -7,8 +7,7 @@
 #include <fftw3.h>
 #include <aocommon/staticfor.h>
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 void FftR2CComposite(fftwf_plan plan_r2c, fftwf_plan plan_c2c,
                      size_t image_height, size_t image_width, const float* in,
                      fftwf_complex* out, aocommon::StaticFor<size_t>& loop);
@@ -17,7 +16,6 @@ void FftC2RComposite(fftwf_plan plan_c2c, fftwf_plan plan_c2r,
                      size_t image_height, size_t image_width,
                      const fftwf_complex* in, float* out,
                      aocommon::StaticFor<size_t>& loop);
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math
 
-#endif
\ No newline at end of file
+#endif

include/schaapcommon/fft/convolution.h → include/schaapcommon/math/convolution.h

@@ -6,8 +6,7 @@
 
 #include <cstring>
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 
 /**
  * @brief Make the FFTW float planner thread safe.
@@ -56,7 +55,6 @@ void PrepareConvolutionKernel(float* dest, const float* source,
  */
 void Convolve(float* image, const float* kernel, size_t image_width,
               size_t image_height);
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math
 
 #endif

include/schaapcommon/math/drawgaussian.h

+#ifndef SCHAAPCOMMON_MATH_DRAW_GAUSSIAN_H_
+#define SCHAAPCOMMON_MATH_DRAW_GAUSSIAN_H_
+
+#include <cstring>
+
+#include "ellipse.h"
+
+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);
+
+}  // namespace schaapcommon::math
+
+#endif

include/schaapcommon/math/ellipse.h

+#ifndef SCHAAPCOMMON_MATH_ELLIPSE_H_
+#define SCHAAPCOMMON_MATH_ELLIPSE_H_
+
+namespace schaapcommon::math {
+
+struct Ellipse {
+  double major;
+  double minor;
+  double position_angle;
+};
+
+}  // namespace schaapcommon::math
+
+#endif

include/schaapcommon/fft/resampler.h → include/schaapcommon/math/resampler.h

@@ -14,8 +14,7 @@
 
 #include <fftw3.h>
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 
 class Resampler {
  private:
@@ -107,7 +106,7 @@ class Resampler {
   aocommon::Lane<Task> tasks_;
   std::vector<std::thread> threads_;
 };
-}  // namespace fft
-}  // namespace schaapcommon
+
+}  // namespace schaapcommon::math
 
 #endif

include/schaapcommon/fft/restoreimage.h → include/schaapcommon/math/restoreimage.h

@@ -6,8 +6,7 @@
 
 #include <cstddef>
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 
 /**
  * @brief Restore a diffuse image (e.g. produced with a multi-scale clean) using
@@ -34,6 +33,5 @@ void RestoreImage(float* image_data, const float* model_data,
                   long double beam_position_angle, long double pixel_scale_l,
                   long double pixel_scale_m);
 
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math
 #endif

src/fitters/gaussianfitter.cc

@@ -16,6 +16,8 @@
 
 using aocommon::Matrix2x2;
 
+using schaapcommon::math::Ellipse;
+
 namespace schaapcommon::fitters {
 namespace {
 

src/fitters/test/tgaussianfitter.cc

@@ -8,8 +8,8 @@
 #include <aocommon/image.h>
 #include <aocommon/threadpool.h>
 
-#include "../../../include/schaapcommon/fft/convolution.h"
-#include "../../../include/schaapcommon/fft/restoreimage.h"
+#include "../../../include/schaapcommon/math/convolution.h"
+#include "../../../include/schaapcommon/math/restoreimage.h"
 
 namespace {
 constexpr size_t kThreadCount = 2;
@@ -18,7 +18,7 @@ constexpr size_t kHeight = 64;
 constexpr long double kPixelSize = 1 /*amin*/ * (M_PI / 180.0 / 60.0);
 }  // namespace
 
-using schaapcommon::fitters::Ellipse;
+using schaapcommon::math::Ellipse;
 
 BOOST_AUTO_TEST_SUITE(gaussian_fitter)
 
@@ -36,8 +36,8 @@ BOOST_AUTO_TEST_CASE(fit) {
     const long double beam_position_angle =
         beam_position_angle_index * M_PI / 10.0;
 
-    schaapcommon::fft::MakeFftwfPlannerThreadSafe();
-    schaapcommon::fft::RestoreImage(
+    schaapcommon::math::MakeFftwfPlannerThreadSafe();
+    schaapcommon::math::RestoreImage(
         restored.Data(), model.Data(), width, height, beam_major, beam_minor,
         beam_position_angle, kPixelSize, kPixelSize);
 
@@ -68,8 +68,8 @@ BOOST_AUTO_TEST_CASE(fit_full) {
     long double beam_minor = 5.0L * kPixelSize;
     long double beam_position_angle = beam_position_angle_index * M_PI / 10.0;
 
-    schaapcommon::fft::MakeFftwfPlannerThreadSafe();
-    schaapcommon::fft::RestoreImage(
+    schaapcommon::math::MakeFftwfPlannerThreadSafe();
+    schaapcommon::math::RestoreImage(
         restored.Data(), model.Data(), width, height, beam_major, beam_minor,
         beam_position_angle, kPixelSize, kPixelSize);
 
@@ -114,10 +114,10 @@ BOOST_AUTO_TEST_CASE(fit_circular) {
   const long double beam_minor = 4.0L * kPixelSize;
   const long double beam_position_angle = 0.0;
   const long double estimated_beam_pixel = 1.0;  // this is on purpose way off
-  schaapcommon::fft::MakeFftwfPlannerThreadSafe();
-  schaapcommon::fft::RestoreImage(restored.Data(), model.Data(), kWidth,
-                                  kHeight, beam_major, beam_minor,
-                                  beam_position_angle, kPixelSize, kPixelSize);
+  schaapcommon::math::MakeFftwfPlannerThreadSafe();
+  schaapcommon::math::RestoreImage(restored.Data(), model.Data(), kWidth,
+                                   kHeight, beam_major, beam_minor,
+                                   beam_position_angle, kPixelSize, kPixelSize);
 
   const Ellipse result = schaapcommon::fitters::Fit2DGaussianCentred(
       restored.Data(), restored.Width(), restored.Height(),
@@ -152,10 +152,10 @@ BOOST_AUTO_TEST_CASE(fit_small_beam) {
   const long double beam_position_angle = 0.0;
   const long double estimated_beam_pixel = 1.0;  // this is on purpose way off
 
-  schaapcommon::fft::MakeFftwfPlannerThreadSafe();
-  schaapcommon::fft::RestoreImage(restored.Data(), model.Data(), kWidth,
-                                  kHeight, beam_major, beam_minor,
-                                  beam_position_angle, kPixelSize, kPixelSize);
+  schaapcommon::math::MakeFftwfPlannerThreadSafe();
+  schaapcommon::math::RestoreImage(restored.Data(), model.Data(), kWidth,
+                                   kHeight, beam_major, beam_minor,
+                                   beam_position_angle, kPixelSize, kPixelSize);
 
   Ellipse ellipse = schaapcommon::fitters::Fit2DGaussianCentred(
       restored.Data(), restored.Width(), restored.Height(),

src/fft/CMakeLists.txt → src/math/CMakeLists.txt

@@ -19,13 +19,15 @@ get_filename_component(MODULE ${CMAKE_CURRENT_SOURCE_DIR} NAME)
 set(PUBLIC_HEADER_DIR ${SCHAAPCOMMON_SOURCE_DIR}/include/schaapcommon/${MODULE})
 
 set(PUBLIC_HEADERS
-    ${PUBLIC_HEADER_DIR}/convolution.h ${PUBLIC_HEADER_DIR}/resampler.h
+    ${PUBLIC_HEADER_DIR}/convolution.h ${PUBLIC_HEADER_DIR}/drawgaussian.h
+    ${PUBLIC_HEADER_DIR}/ellipse.h ${PUBLIC_HEADER_DIR}/resampler.h
     ${PUBLIC_HEADER_DIR}/restoreimage.h)
 
 target_sources(
   ${SCHAAPCOMMON_PROJECT_NAME}
   PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/convolution.cc
           ${CMAKE_CURRENT_SOURCE_DIR}/compositefft.cc
+          ${CMAKE_CURRENT_SOURCE_DIR}/drawgaussian.cc
           ${CMAKE_CURRENT_SOURCE_DIR}/resampler.cc
           ${CMAKE_CURRENT_SOURCE_DIR}/restoreimage.cc)
 

src/fft/compositefft.cc → src/math/compositefft.cc

@@ -19,8 +19,7 @@ size_t RoundUp(size_t a, size_t b) { return ((a + b) / b) * b; }
 
 }  // namespace
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 void FftR2CComposite(fftwf_plan plan_r2c, fftwf_plan plan_c2c,
                      size_t image_height, size_t image_width, const float* in,
                      fftwf_complex* out, aocommon::StaticFor<size_t>& loop) {
@@ -162,5 +161,4 @@ void FftC2RComposite(fftwf_plan plan_c2c, fftwf_plan plan_c2r,
 
   fftwf_free(temp1);
 }
-}  // namespace fft
-}  // namespace schaapcommon
\ No newline at end of file
+}  // namespace schaapcommon::math

src/fft/convolution.cc → src/math/convolution.cc

@@ -15,8 +15,7 @@
 
 #include <iostream>
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 
 void MakeFftwfPlannerThreadSafe() { fftwf_make_planner_thread_safe(); }
 
@@ -167,5 +166,4 @@ void Convolve(float* image, const float* kernel, size_t image_width,
   fftwf_destroy_plan(plan_c2r);
 }
 
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math

src/math/drawgaussian.cc

+#include "drawgaussian.h"
+
+#include <cmath>
+
+#include <aocommon/imagecoordinates.h>
+
+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);
+}
+}  // namespace
+
+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) {
+  // 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;
+  const long double cos_angle = std::cos(angle);
+  const long double sin_angle = std::sin(angle);
+
+  // Make rotation matrix
+  long double transf[4];
+  transf[0] = cos_angle;
+  transf[1] = -sin_angle;
+  transf[2] = sin_angle;
+  transf[3] = cos_angle;
+
+  const double sigma_max = std::max(std::fabs(sigma_major_axis * transf[0]),
+                                    std::fabs(sigma_major_axis * transf[1]));
+  // Multiply with scaling matrix to make variance 1.
+  transf[0] /= sigma_major_axis;
+  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(position_ra, position_dec, ra, 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 x_left =
+      std::clamp(source_x - bounding_box_size, 0, int(image_width));
+  const int x_right =
+      std::clamp(source_x + bounding_box_size, x_left, int(image_width));
+  const int y_top =
+      std::clamp(source_y - bounding_box_size, 0, int(image_height));
+  const int y_bottom =
+      std::clamp(source_y + bounding_box_size, y_top, int(image_height));
+
+  std::vector<double> values;
+  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 dist =
+          std::sqrt(l_transf * l_transf + m_transf * m_transf);
+      long double v = Gaussian(dist, 1.0L);
+      flux_sum += static_cast<double>(v);
+      values.emplace_back(v);
+    }
+  }
+  const double* iter = values.data();
+  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;
+    }
+  }
+}
+
+}  // namespace schaapcommon::math

src/fft/resampler.cc → src/math/resampler.cc

@@ -3,8 +3,7 @@
 #include <complex>
 
 using aocommon::WindowFunction;
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 
 Resampler::Resampler(size_t input_width, size_t input_height,
                      size_t output_width, size_t output_height,
@@ -248,5 +247,4 @@ void Resampler::UnapplyWindow(float* data) const {
     data[i] /= window_out_[i];
   }
 }
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math

src/fft/restoreimage.cc → src/math/restoreimage.cc

@@ -12,14 +12,13 @@
 
 #include "convolution.h"
 
-namespace schaapcommon {
-namespace fft {
+namespace schaapcommon::math {
 void RestoreImage(float* image_data, const float* model_data,
                   size_t image_width, size_t image_height,
                   long double beam_major_axis, long double beam_minor_axis,
                   long double beam_position_angle, long double pixel_scale_l,
                   long double pixel_scale_m) {
-  if (beam_major_axis == 0.0 && beam_minor_axis == 0.0) {
+  if (beam_major_axis == 0.0L && beam_minor_axis == 0.0L) {
     for (size_t j = 0; j != image_width * image_height; ++j) {
       image_data[j] += model_data[j];
     }
@@ -83,13 +82,12 @@ void RestoreImage(float* image_data, const float* model_data,
     aocommon::UVector<float> convolved_model(
         model_data, model_data + image_width * image_height);
 
-    schaapcommon::fft::ResizeAndConvolve(convolved_model.data(), image_width,
-                                         image_height, kernel.data(),
-                                         bounding_box_size);
+    schaapcommon::math::ResizeAndConvolve(convolved_model.data(), image_width,
+                                          image_height, kernel.data(),
+                                          bounding_box_size);
     for (size_t j = 0; j != image_width * image_height; ++j) {
       image_data[j] += convolved_model[j];
     }
   }
 }
-}  // namespace fft
-}  // namespace schaapcommon
+}  // namespace schaapcommon::math

src/fft/test/CMakeLists.txt → src/math/test/CMakeLists.txt

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

src/fft/test/tconvolution.cc → src/math/test/tconvolution.cc

@@ -25,7 +25,7 @@ BOOST_AUTO_TEST_CASE(prepare_kernel) {
                             {10.0, 11.0, 8.0, 9.0, 14.0, 15.0, 12.0, 13.0, 2.0,
                              3.0, 0.0, 1.0, 6.0, 7.0, 4.0, 5.0});
   aocommon::Image kernel_out(kWidth, kHeight);
-  schaapcommon::fft::PrepareConvolutionKernel(
+  schaapcommon::math::PrepareConvolutionKernel(
       kernel_out.Data(), kernel_in.Data(), kWidth, kHeight);
   for (size_t i = 0; i != kernel_out.Size(); ++i) {
     BOOST_CHECK_CLOSE(kernel_out[i], static_cast<float>(i), 1e-4);
@@ -38,11 +38,11 @@ BOOST_AUTO_TEST_CASE(prepare_small_kernel) {
   aocommon::Image kernel_in(kWidth / 2, kHeight / 2, {16.0, 13.0, 4.0, 1.0});
   aocommon::Image kernel_out(kWidth, kHeight, 0.0);
   BOOST_CHECK_THROW(
-      schaapcommon::fft::PrepareSmallConvolutionKernel(
+      schaapcommon::math::PrepareSmallConvolutionKernel(
           kernel_out.Data(), kWidth, kHeight, kernel_in.Data(), kWidth * 2),
       std::runtime_error);
 
-  schaapcommon::fft::PrepareSmallConvolutionKernel(
+  schaapcommon::math::PrepareSmallConvolutionKernel(
       kernel_out.Data(), kWidth, kHeight, kernel_in.Data(), kWidth / 2);
 
   // kernel_out should have non-zero corner values only
@@ -68,14 +68,14 @@ BOOST_AUTO_TEST_CASE(convolve) {
 
   kernel[kWidth * kHeight / 2 + kWidth / 2] = dirac_scale * 1.0f;
 
-  schaapcommon::fft::MakeFftwfPlannerThreadSafe();
+  schaapcommon::math::MakeFftwfPlannerThreadSafe();
   BOOST_CHECK_THROW(
-      schaapcommon::fft::ResizeAndConvolve(image.Data(), kWidth, kHeight,
-                                           kernel.Data(), kWidth * 2),
+      schaapcommon::math::ResizeAndConvolve(image.Data(), kWidth, kHeight,
+                                            kernel.Data(), kWidth * 2),
       std::runtime_error);
 
-  schaapcommon::fft::ResizeAndConvolve(image.Data(), kWidth, kHeight,
-                                       kernel.Data(), kWidth);
+  schaapcommon::math::ResizeAndConvolve(image.Data(), kWidth, kHeight,
+                                        kernel.Data(), kWidth);
 
   // Convolution with (scaled) dirac kernel should return same signal,
   // scaled by a constant factor.

src/fft/test/tresampler.cc → src/math/test/tresampler.cc

@@ -9,7 +9,7 @@
 
 #include <aocommon/image.h>
 
-using schaapcommon::fft::Resampler;
+using schaapcommon::math::Resampler;
 
 namespace {
 constexpr size_t kInputWidth = 16;
@@ -137,4 +137,4 @@ BOOST_AUTO_TEST_CASE(upsample_regular_window) {
   CheckResampledImage(input_image, output_image, scale_output);
 }
 
-BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
+BOOST_AUTO_TEST_SUITE_END()