diff --git a/include/schaapcommon/fitters/gaussianfitter.h b/include/schaapcommon/fitters/gaussianfitter.h
index 8dcf774262f7f6af689f57cb81a1c771a434130d..160f2511d5c87cddb6e0ba88b724daf77d4bf511 100644
--- a/include/schaapcommon/fitters/gaussianfitter.h
+++ b/include/schaapcommon/fitters/gaussianfitter.h
@@ -25,105 +25,41 @@ struct Ellipse {
  */
 Ellipse DeconvolveGaussian(const Ellipse& largest, const Ellipse& smallest);
 
-class GaussianFitter {
- public:
-  void Fit2DGaussianCentred(const float* image, size_t width, size_t height,
-                            double beam_est, double& beam_major,
-                            double& beam_minor, double& beam_pa,
-                            double box_scale_factor = 10.0,
-                            bool verbose = false);
-
-  void Fit2DCircularGaussianCentred(const float* image, size_t width,
-                                    size_t height, double& beam_size,
-                                    double box_scale_factor = 10.0);
-
-  void Fit2DGaussianFull(const float* image, size_t width, size_t height,
-                         double& val, double& pos_x, double& pos_y,
-                         double& beam_major, double& beam_minor,
-                         double& beam_pa, double* floor_level = nullptr);
-
-  const float* Image() const { return image_; }
-  size_t Width() const { return width_; }
-  size_t Height() const { return height_; }
-  size_t ScaleFactor() const { return scale_factor_; }
-  double XInit() const { return x_init_; };
-  double YInit() const { return y_init_; };
-  double PosConstrained() const { return pos_constrained_; };
-
- private:
-  void Fit2DGaussianCentredInBox(const float* image, size_t width,
-                                 size_t height, double beam_est,
-                                 double& beam_major, double& beam_minor,
-                                 double& beam_pa, size_t box_width,
-                                 size_t box_height, bool verbose);
-
-  void Fit2DCircularGaussianCentredInBox(const float* image, size_t width,
-                                         size_t height, double& beam_size,
-                                         size_t box_width, size_t box_height);
-
-  /**
-   * This function performs a single fit of a Gaussian. The position of the
-   * Gaussian is constrained to be in the centre of the image. The Gaussian is
-   * fitted such that the squared residuals (data - model) are minimal.
-   *
-   * This function is typically used to find the beam-shape of the point-spread
-   * function. The beam estimate is used as initial value for the minor and
-   * major shape.
-   */
-  void SingleFit2DGaussianCentred(const float* image, size_t width,
-                                  size_t height, double beam_est,
-                                  double& beam_major, double& beam_minor,
-                                  double& beam_pa, bool verbose);
-
-  void SingleFit2DCircularGaussianCentred(const float* image, size_t width,
-                                          size_t height, double& beam_size);
-
-  void Fit2DGaussianWithAmplitudeInBox(const float* image, size_t width,
-                                       size_t height, double& val,
-                                       double& pos_x, double& pos_y,
-                                       double& beam_major, double& beam_minor,
-                                       double& beam_pa, double* floor_level,
-                                       size_t x_start, size_t x_end,
-                                       size_t y_start, size_t y_end);
-
-  /**
-   * Fits the position, size and amplitude of a Gaussian. If floor_level is not
-   * a nullptr, the floor (background level, or zero level) is fitted too.
-   */
-  void Fit2DGaussianWithAmplitude(const float* image, size_t width,
-                                  size_t height, double& val, double& pos_x,
-                                  double& pos_y, double& beam_major,
-                                  double& beam_minor, double& beam_pa,
-                                  double* floor_level);
+/**
+ * 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);
 
-  /**
-   * Like SingleFit2DGaussianCentred(), but includes Gaussian centre X and Y
-   * position and amplitude in the fitted parameters.
-   *
-   * This function can typically be used for source fitting.
-   */
-  void Fit2DGaussianWithAmplitude(double& val, double& pos_x, double& pos_y,
-                                  double& beam_major, double& beam_minor,
-                                  double& beam_pa);
+/**
+ * Fits a circular (one parameter) Gaussian to the shape in the centre of the
+ * image.
+ * @param [in,out] beam_size on input, the initial value for beam size. On
+ * output, the fit result.
+ */
+void Fit2DCircularGaussianCentred(const float* image, size_t width,
+                                  size_t height, double& beam_size,
+                                  double box_scale_factor = 10.0);
 
-  /**
-   * Like Fit2DGaussianWithAmplitude(), but includes floor_level as fitted
-   * parameter. Floor is the background/zero level on which the Gaussian
-   * resides.
-   */
-  void Fit2DGaussianWithAmplitudeWithFloor(double& val, double& pos_x,
-                                           double& pos_y, double& beam_major,
-                                           double& beam_minor, double& beam_pa,
-                                           double& floor_level);
+/**
+ * Fit all Gaussian parameters to a shape in an image. Parameters
+ * @p pos_x, @p pos_y and @p beam_major are input and output function
+ * parameters: on input they provide initial values for the fitting,
+ * and on output they provide the result. Parameters @p val,
+ * @p beam_minor, @p beam_pa and @p floor_level are only used for outputting
+ * the result.
+ * @param [out] val amplitude of Gaussian.
+ * @param [in,out] floor_level If not nullptr, will be used to store the floor
+ * level parameter. If nullptr, this value is not fitted and assumed to be zero.
+ */
+void Fit2DGaussianFull(const float* image, size_t width, size_t height,
+                       double& val, double& pos_x, double& pos_y,
+                       double& beam_major, double& beam_minor, double& beam_pa,
+                       double* floor_level = nullptr);
 
-  const float* image_ = nullptr;
-  size_t width_ = 0;
-  size_t height_ = 0;
-  size_t scale_factor_ = 0;
-  double x_init_ = 0.0;
-  double y_init_ = 0.0;
-  double pos_constrained_ = 0.0;
-};
 }  // namespace fitters
 }  // namespace schaapcommon
 #endif
diff --git a/src/fitters/gaussianfitter.cc b/src/fitters/gaussianfitter.cc
index f94a219b3cdbba58ab5f3ac02059ace8e2a66b19..1a6750c6c9275f34e3d668ac0bd88ac31a0e8506 100644
--- a/src/fitters/gaussianfitter.cc
+++ b/src/fitters/gaussianfitter.cc
@@ -3,6 +3,7 @@
 
 #include "gaussianfitter.h"
 
+#include <algorithm>
 #include <cmath>
 #include <iostream>
 #include <limits>
@@ -15,13 +16,18 @@
 
 using aocommon::Matrix2x2;
 
-namespace schaapcommon {
-namespace fitters {
-
+namespace schaapcommon::fitters {
 namespace {
 
 const long double kSigmaToBeam = 2.0L * std::sqrt(2.0L * std::log(2.0L));
 
+struct FitData {
+  const float* image = nullptr;
+  size_t width = 0;
+  size_t height = 0;
+  size_t scale_factor = 0;
+};
+
 void ToAnglesAndFwhm(double sx, double sy, double beta, double& ellipse_major,
                      double& ellipse_minor, double& ellipse_position_angle) {
   const double beta_factor = 1.0 - beta * beta;
@@ -49,6 +55,13 @@ void ToAnglesAndFwhm(double sx, double sy, double beta, double& ellipse_major,
   }
 }
 
+Ellipse ToAnglesAndFwhm(double sx, double sy, double beta) {
+  Ellipse ellipse;
+  ToAnglesAndFwhm(sx, sy, beta, ellipse.major, ellipse.minor,
+                  ellipse.position_angle);
+  return ellipse;
+}
+
 /**
  * Calculates a two dimensional gaussian with the specified parameters.
  * @param x X coordinate.
@@ -77,22 +90,23 @@ double GaussCircularCentered(double x, double y, double s) {
  * squared errors(/residuals).
  */
 int FittingCentered(const gsl_vector* xvec, void* data, gsl_vector* f) {
-  const GaussianFitter& fitter = *static_cast<const GaussianFitter*>(data);
+  const FitData& fitter = *static_cast<const FitData*>(data);
   const double sx = gsl_vector_get(xvec, 0);
   const double sy = gsl_vector_get(xvec, 1);
   const double beta = gsl_vector_get(xvec, 2);
-  const size_t width = fitter.Width();
-  const size_t height = fitter.Height();
+  const size_t width = fitter.width;
+  const size_t height = fitter.height;
   const int x_mid = width / 2;
   const int y_mid = height / 2;
-  const double scale = 1.0 / fitter.ScaleFactor();
+  const double scale = 1.0 / fitter.scale_factor;
 
   size_t data_index = 0;
   for (size_t yi = 0; yi != height; ++yi) {
-    double y = (yi - y_mid) * scale;
+    const double y = (yi - y_mid) * scale;
     for (size_t xi = 0; xi != width; ++xi) {
-      double x = (xi - x_mid) * scale;
-      double e = GaussCentered(x, y, sx, sy, beta) - fitter.Image()[data_index];
+      const double x = (xi - x_mid) * scale;
+      const double e =
+          GaussCentered(x, y, sx, sy, beta) - fitter.image[data_index];
       gsl_vector_set(f, data_index, e);
       ++data_index;
     }
@@ -101,20 +115,21 @@ int FittingCentered(const gsl_vector* xvec, void* data, gsl_vector* f) {
 }
 
 int FittingCircularCentered(const gsl_vector* xvec, void* data, gsl_vector* f) {
-  const GaussianFitter& fitter = *static_cast<const GaussianFitter*>(data);
+  const FitData& fitter = *static_cast<const FitData*>(data);
   const double s = gsl_vector_get(xvec, 0);
-  const size_t width = fitter.Width();
-  const size_t height = fitter.Height();
+  const size_t width = fitter.width;
+  const size_t height = fitter.height;
   const int x_mid = width / 2;
   const int y_mid = height / 2;
-  const double scale = 1.0 / fitter.ScaleFactor();
+  const double scale = 1.0 / fitter.scale_factor;
 
   size_t data_index = 0;
   for (size_t yi = 0; yi != height; ++yi) {
-    double y = (yi - y_mid) * scale;
+    const double y = (yi - y_mid) * scale;
     for (size_t xi = 0; xi != width; ++xi) {
-      double x = (xi - x_mid) * scale;
-      double e = GaussCircularCentered(x, y, s) - fitter.Image()[data_index];
+      const double x = (xi - x_mid) * scale;
+      const double e =
+          GaussCircularCentered(x, y, s) - fitter.image[data_index];
       gsl_vector_set(f, data_index, e);
       ++data_index;
     }
@@ -127,27 +142,27 @@ int FittingCircularCentered(const gsl_vector* xvec, void* data, gsl_vector* f) {
  */
 int FittingDerivativeCentered(const gsl_vector* xvec, void* data,
                               gsl_matrix* J) {
-  const GaussianFitter& fitter = *static_cast<const GaussianFitter*>(data);
+  const FitData& fitter = *static_cast<const FitData*>(data);
   const double sx = gsl_vector_get(xvec, 0);
   const double sy = gsl_vector_get(xvec, 1);
   const double beta = gsl_vector_get(xvec, 2);
-  const size_t width = fitter.Width();
-  const size_t height = fitter.Height();
+  const size_t width = fitter.width;
+  const size_t height = fitter.height;
   const int x_mid = width / 2;
   const int y_mid = height / 2;
-  const double scale = 1.0 / fitter.ScaleFactor();
+  const double scale = 1.0 / fitter.scale_factor;
 
   size_t data_index = 0;
   for (size_t yi = 0; yi != height; ++yi) {
-    double y = (yi - y_mid) * scale;
+    const double y = (yi - y_mid) * scale;
     for (size_t xi = 0; xi != width; ++xi) {
-      double x = (xi - x_mid) * scale;
-      double exp_term = GaussCentered(x, y, sx, sy, beta);
-      double dsx =
+      const double x = (xi - x_mid) * scale;
+      const double exp_term = GaussCentered(x, y, sx, sy, beta);
+      const double dsx =
           (beta * x * y / (sx * sx * sy) + x * x / (sx * sx * sx)) * exp_term;
-      double dsy =
+      const double dsy =
           (beta * x * y / (sy * sy * sx) + y * y / (sy * sy * sy)) * exp_term;
-      double dbeta = x * y / (sx * sy) * exp_term;
+      const double dbeta = x * y / (sx * sy) * exp_term;
       gsl_matrix_set(J, data_index, 0, dsx);
       gsl_matrix_set(J, data_index, 1, dsy);
       gsl_matrix_set(J, data_index, 2, dbeta);
@@ -159,24 +174,24 @@ int FittingDerivativeCentered(const gsl_vector* xvec, void* data,
 
 int FittingDerivativeCircularCentered(const gsl_vector* xvec, void* data,
                                       gsl_matrix* J) {
-  const GaussianFitter& fitter = *static_cast<const GaussianFitter*>(data);
+  const FitData& fitter = *static_cast<const FitData*>(data);
   const double s = gsl_vector_get(xvec, 0);
-  const size_t width = fitter.Width();
-  const size_t height = fitter.Height();
+  const size_t width = fitter.width;
+  const size_t height = fitter.height;
   const int x_mid = width / 2, y_mid = height / 2;
-  const double scale = 1.0 / fitter.ScaleFactor();
+  const double scale = 1.0 / fitter.scale_factor;
 
   size_t data_index = 0;
   for (size_t yi = 0; yi != height; ++yi) {
-    double y = (yi - y_mid) * scale;
+    const double y = (yi - y_mid) * scale;
     for (size_t xi = 0; xi != width; ++xi) {
-      double x = (xi - x_mid) * scale;
-      double exp_term = GaussCircularCentered(x, y, s);
+      const double x = (xi - x_mid) * scale;
+      const double exp_term = GaussCircularCentered(x, y, s);
       // derivative of exp((-x*x - y*y)/(2.0*s*s)) to s
       // = (-x*x - y*y)/2.0*-2/(s*s*s)
       // = (-x*x - y*y)/(-s*s*s)
       // = (x*x + y*y)/(s*s*s)
-      double ds = ((x * x + y * y) / (s * s * s)) * exp_term;
+      const double ds = ((x * x + y * y) / (s * s * s)) * exp_term;
       gsl_matrix_set(J, data_index, 0, ds);
       ++data_index;
     }
@@ -202,26 +217,26 @@ int FittingBothCircularCentered(const gsl_vector* x, void* data, gsl_vector* f,
 }
 
 int FittingWithAmplitude(const gsl_vector* xvec, void* data, gsl_vector* f) {
-  const GaussianFitter& fitter = *static_cast<const GaussianFitter*>(data);
+  const FitData& fitter = *static_cast<const FitData*>(data);
   const double v = gsl_vector_get(xvec, 0);
   const double xc = gsl_vector_get(xvec, 1);
   const double yc = gsl_vector_get(xvec, 2);
   const double sx = gsl_vector_get(xvec, 3);
   const double sy = gsl_vector_get(xvec, 4);
   const double beta = gsl_vector_get(xvec, 5);
-  const size_t width = fitter.Width();
-  const size_t height = fitter.Height();
+  const size_t width = fitter.width;
+  const size_t height = fitter.height;
   const int x_mid = width / 2;
   const int y_mid = height / 2;
-  const double scale = 1.0 / fitter.ScaleFactor();
+  const double scale = 1.0 / fitter.scale_factor;
 
   size_t data_index = 0;
   for (int yi = 0; yi != int(height); ++yi) {
-    double yS = yc + (yi - y_mid) * scale;
+    const double yS = yc + (yi - y_mid) * scale;
     for (int xi = 0; xi != int(width); ++xi) {
-      double xS = xc + (xi - x_mid) * scale;
-      double e =
-          GaussCentered(xS, yS, sx, sy, beta) * v - fitter.Image()[data_index];
+      const double xS = xc + (xi - x_mid) * scale;
+      const double e =
+          GaussCentered(xS, yS, sx, sy, beta) * v - fitter.image[data_index];
       gsl_vector_set(f, data_index, e);
       ++data_index;
     }
@@ -231,41 +246,34 @@ int FittingWithAmplitude(const gsl_vector* xvec, void* data, gsl_vector* f) {
 
 int FittingDerivativeWithAmplitude(const gsl_vector* xvec, void* data,
                                    gsl_matrix* J) {
-  const GaussianFitter& fitter = *static_cast<const GaussianFitter*>(data);
-  const double scale = 1.0 / fitter.ScaleFactor();
+  const FitData& fitter = *static_cast<const FitData*>(data);
+  const double scale = 1.0 / fitter.scale_factor;
   const double v = gsl_vector_get(xvec, 0);
   const double xc = gsl_vector_get(xvec, 1);
   const double yc = gsl_vector_get(xvec, 2);
   const double sx = gsl_vector_get(xvec, 3);
   const double sy = gsl_vector_get(xvec, 4);
   const double beta = gsl_vector_get(xvec, 5);
-  if (fitter.PosConstrained() != 0.0 &&
-      (std::fabs(xc - fitter.XInit()) > fitter.PosConstrained() * scale ||
-       std::fabs(yc - fitter.YInit()) > fitter.PosConstrained() * scale)) {
-    std::cout << "GSL_EDOM\n";
-    return GSL_EDOM;
-  }
-  const size_t width = fitter.Width();
-  const size_t height = fitter.Height();
+  const size_t width = fitter.width;
+  const size_t height = fitter.height;
   const int x_mid = width / 2;
   const int y_mid = height / 2;
 
   size_t data_index = 0;
   for (int yi = 0; yi != int(height); ++yi) {
-    double y = yc + (yi - y_mid) * scale;
+    const double y = yc + (yi - y_mid) * scale;
     for (int xi = 0; xi != int(width); ++xi) {
       // TODO I need to go over the signs -- ds, dy, dsx, dsy in particular
-      double x = xc + (xi - x_mid) * scale;
-      double exp_term = GaussCentered(x, y, sx, sy, beta);
-      double dv = exp_term;
-      exp_term *= v;
-      double dx = (-beta * y / (sx * sy) - x / (sx * sx)) * exp_term;
-      double dy = (-beta * x / (sy * sx) - y / (sy * sy)) * exp_term;
-      double dsx =
+      const double x = xc + (xi - x_mid) * scale;
+      const double dv = GaussCentered(x, y, sx, sy, beta);
+      const double exp_term = dv * v;
+      const double dx = (-beta * y / (sx * sy) - x / (sx * sx)) * exp_term;
+      const double dy = (-beta * x / (sy * sx) - y / (sy * sy)) * exp_term;
+      const double dsx =
           (beta * x * y / (sx * sx * sy) + x * x / (sx * sx * sx)) * exp_term;
-      double dsy =
+      const double dsy =
           (beta * x * y / (sy * sy * sx) + y * y / (sy * sy * sy)) * exp_term;
-      double dbeta = x * y / (sx * sy) * exp_term;
+      const double dbeta = x * y / (sx * sy) * exp_term;
       gsl_matrix_set(J, data_index, 0, dv);
       gsl_matrix_set(J, data_index, 1, dx);
       gsl_matrix_set(J, data_index, 2, dy);
@@ -287,8 +295,8 @@ int FittingBothWithAmplitude(const gsl_vector* x, void* data, gsl_vector* f,
 
 int FittingWithAmplitudeAndFloor(const gsl_vector* xvec, void* data,
                                  gsl_vector* f) {
-  const GaussianFitter& fitter = *static_cast<const GaussianFitter*>(data);
-  const double scale = 1.0 / fitter.ScaleFactor();
+  const FitData& fitter = *static_cast<const FitData*>(data);
+  const double scale = 1.0 / fitter.scale_factor;
   const double v = gsl_vector_get(xvec, 0);
   const double xc = gsl_vector_get(xvec, 1);
   const double yc = gsl_vector_get(xvec, 2);
@@ -296,23 +304,18 @@ int FittingWithAmplitudeAndFloor(const gsl_vector* xvec, void* data,
   const double sy = gsl_vector_get(xvec, 4);
   const double beta = gsl_vector_get(xvec, 5);
   const double fl = gsl_vector_get(xvec, 6);
-  if (fitter.PosConstrained() != 0.0 &&
-      (std::fabs(xc - fitter.XInit()) > fitter.PosConstrained() * scale ||
-       std::fabs(yc - fitter.YInit()) > fitter.PosConstrained() * scale)) {
-    return GSL_EDOM;
-  }
-  const size_t width = fitter.Width();
-  const size_t height = fitter.Height();
+  const size_t width = fitter.width;
+  const size_t height = fitter.height;
   const int x_mid = width / 2;
   const int y_mid = height / 2;
 
   size_t data_index = 0;
   for (int yi = 0; yi != int(height); ++yi) {
-    double yS = yc + (yi - y_mid) * scale;
+    const double yS = yc + (yi - y_mid) * scale;
     for (int xi = 0; xi != int(width); ++xi) {
-      double xS = xc + (xi - x_mid) * scale;
-      double e = GaussCentered(xS, yS, sx, sy, beta) * v -
-                 fitter.Image()[data_index] + fl;
+      const double xS = xc + (xi - x_mid) * scale;
+      const double e = GaussCentered(xS, yS, sx, sy, beta) * v -
+                       fitter.image[data_index] + fl;
       gsl_vector_set(f, data_index, e);
       ++data_index;
     }
@@ -322,35 +325,34 @@ int FittingWithAmplitudeAndFloor(const gsl_vector* xvec, void* data,
 
 int FittingDerivativeWithAmplitudeAndFloor(const gsl_vector* xvec, void* data,
                                            gsl_matrix* J) {
-  const GaussianFitter& fitter = *static_cast<const GaussianFitter*>(data);
+  const FitData& fitter = *static_cast<const FitData*>(data);
   const double v = gsl_vector_get(xvec, 0);
   const double xc = gsl_vector_get(xvec, 1);
   const double yc = gsl_vector_get(xvec, 2);
   const double sx = gsl_vector_get(xvec, 3);
   const double sy = gsl_vector_get(xvec, 4);
   const double beta = gsl_vector_get(xvec, 5);
-  const size_t width = fitter.Width();
-  const size_t height = fitter.Height();
+  const size_t width = fitter.width;
+  const size_t height = fitter.height;
   const int x_mid = width / 2;
   const int y_mid = height / 2;
-  const double scale = 1.0 / fitter.ScaleFactor();
+  const double scale = 1.0 / fitter.scale_factor;
 
   size_t data_index = 0;
   for (int yi = 0; yi != int(height); ++yi) {
-    double y = yc + (yi - y_mid) * scale;
+    const double y = yc + (yi - y_mid) * scale;
     for (int xi = 0; xi != int(width); ++xi) {
-      double x = xc + (xi - x_mid) * scale;
-      double exp_term = GaussCentered(x, y, sx, sy, beta);
-      double dv = exp_term;
-      exp_term *= v;
-      double dx = (-beta * y / (sx * sy) - x / (sx * sx)) * exp_term;
-      double dy = (-beta * x / (sy * sx) - y / (sy * sy)) * exp_term;
-      double dsx =
+      const double x = xc + (xi - x_mid) * scale;
+      const double dv = GaussCentered(x, y, sx, sy, beta);
+      const double exp_term = dv * v;
+      const double dx = (-beta * y / (sx * sy) - x / (sx * sx)) * exp_term;
+      const double dy = (-beta * x / (sy * sx) - y / (sy * sy)) * exp_term;
+      const double dsx =
           (beta * x * y / (sx * sx * sy) + x * x / (sx * sx * sx)) * exp_term;
-      double dsy =
+      const double dsy =
           (beta * x * y / (sy * sy * sx) + y * y / (sy * sy * sy)) * exp_term;
-      double dbeta = x * y / (sx * sy) * exp_term;
-      double dfl = 1.0;
+      const double dbeta = x * y / (sx * sy) * exp_term;
+      const double dfl = 1.0;
       gsl_matrix_set(J, data_index, 0, dv);
       gsl_matrix_set(J, data_index, 1, dx);
       gsl_matrix_set(J, data_index, 2, dy);
@@ -371,193 +373,46 @@ int FittingBothWithAmplitudeAndFloor(const gsl_vector* x, void* data,
   return GSL_SUCCESS;
 }
 
-}  // namespace
-
-void GaussianFitter::Fit2DGaussianCentred(
-    const float* image, size_t width, size_t height, double beam_estimate,
-    double& beam_major, double& beam_minor, double& beam_phase_angle,
-    double box_scale_factor, bool verbose) {
-  size_t preferred_size = std::max<size_t>(
-      std::ceil(box_scale_factor), std::ceil(beam_estimate * box_scale_factor));
-  if (preferred_size % 2 != 0) ++preferred_size;
-  if (preferred_size < width || preferred_size < height) {
-    size_t n_iterations = 0;
-    bool box_was_large_enough;
-    do {
-      size_t box_width = std::min(preferred_size, width);
-      size_t box_height = std::min(preferred_size, height);
-      if (verbose) std::cout << "Fit initial value:" << beam_estimate << "\n";
-      Fit2DGaussianCentredInBox(image, width, height, beam_estimate, beam_major,
-                                beam_minor, beam_phase_angle, box_width,
-                                box_height, verbose);
-      if (verbose) {
-        std::cout << "Fit result:" << beam_major << " x " << beam_minor
-                  << " px, " << beam_phase_angle << " (box was " << box_width
-                  << " x " << box_height << ")\n";
-      }
-
-      box_was_large_enough =
-          (beam_major * box_scale_factor * 0.8 < box_width ||
-           box_width >= width) &&
-          (beam_major * box_scale_factor * 0.8 < box_height ||
-           box_height >= height);
-      if (!box_was_large_enough) {
-        preferred_size =
-            std::max<size_t>(std::ceil(box_scale_factor),
-                             std::ceil(beam_major * box_scale_factor));
-        if (preferred_size % 2 != 0) ++preferred_size;
-        beam_estimate = std::max(beam_major, beam_estimate);
-      }
-      ++n_iterations;
-    } while (!box_was_large_enough && n_iterations < 5);
-  } else {
-    if (verbose) std::cout << "Image is as large as the fitting box.\n";
-    SingleFit2DGaussianCentred(image, width, height, beam_estimate, beam_major,
-                               beam_minor, beam_phase_angle, verbose);
-  }
-}
-
-void GaussianFitter::Fit2DCircularGaussianCentred(const float* image,
-                                                  size_t width, size_t height,
-                                                  double& beam_size,
-                                                  double box_scale_factor) {
-  double initial_value = beam_size;
-  size_t preferred_size = std::max<size_t>(
-      std::ceil(box_scale_factor), std::ceil(beam_size * box_scale_factor));
-  if (preferred_size % 2 != 0) ++preferred_size;
-  if (preferred_size < width || preferred_size < height) {
-    size_t box_width = std::min(preferred_size, width);
-    size_t box_height = std::min(preferred_size, height);
-    size_t n_iterations = 0;
-    bool box_was_large_enough;
-    do {
-      Fit2DCircularGaussianCentredInBox(image, width, height, beam_size,
-                                        box_width, box_height);
-
-      box_was_large_enough = (beam_size * box_scale_factor * 0.8 < box_width ||
-                              width >= box_width) &&
-                             (beam_size * box_scale_factor * 0.8 < box_height ||
-                              height >= box_height);
-      if (!box_was_large_enough) {
-        preferred_size =
-            std::max<size_t>(std::ceil(box_scale_factor),
-                             std::ceil(beam_size * box_scale_factor));
-        if (preferred_size % 2 != 0) ++preferred_size;
-        beam_size = std::max(initial_value, beam_size);
-      }
-      ++n_iterations;
-    } while (!box_was_large_enough && n_iterations < 5);
-  } else {
-    SingleFit2DCircularGaussianCentred(image, width, height, beam_size);
-  }
-}
-
-void GaussianFitter::Fit2DGaussianFull(const float* image, size_t width,
-                                       size_t height, double& val,
-                                       double& pos_x, double& pos_y,
-                                       double& beam_major, double& beam_minor,
-                                       double& beam_phase_angle,
-                                       double* floor_level) {
-  size_t preferred_size = std::max<size_t>(10, std::ceil(beam_major * 10.0));
-  if (preferred_size % 2 != 0) ++preferred_size;
-  if (preferred_size < width || preferred_size < height) {
-    size_t x_start =
-        std::max<int>(0, int(std::round(pos_x)) - int(preferred_size) / 2);
-    size_t x_end =
-        std::min(width, size_t(std::round(pos_x)) + preferred_size / 2);
-    size_t y_start =
-        std::max<int>(0, int(std::round(pos_y)) - int(preferred_size) / 2);
-    size_t y_end =
-        std::min(height, size_t(std::round(pos_y)) + preferred_size / 2);
-    size_t n_iterations = 0;
-    bool box_was_large_enough;
-    do {
-      Fit2DGaussianWithAmplitudeInBox(
-          image, width, height, val, pos_x, pos_y, beam_major, beam_minor,
-          beam_phase_angle, floor_level, x_start, x_end, y_start, y_end);
-
-      size_t box_width = x_end - x_start;
-      size_t box_height = y_end - y_start;
-      box_was_large_enough =
-          (beam_major * 4.0 < box_width || width >= box_width) &&
-          (beam_major * 4.0 < box_height || height >= box_height);
-      if (!box_was_large_enough) {
-        preferred_size = std::max<size_t>(10, std::ceil(beam_major * 10.0));
-        if (preferred_size % 2 != 0) ++preferred_size;
-      }
-      ++n_iterations;
-    } while (!box_was_large_enough && n_iterations < 5);
-  } else {
-    Fit2DGaussianWithAmplitude(image, width, height, val, pos_x, pos_y,
-                               beam_major, beam_minor, beam_phase_angle,
-                               floor_level);
-  }
-}
-
-void GaussianFitter::Fit2DGaussianCentredInBox(
-    const float* image, size_t width, size_t height, double beam_estimate,
-    double& beam_major, double& beam_minor, double& beam_phase_angle,
-    size_t box_width, size_t box_height, bool verbose) {
-  const size_t x_start = (width - box_width) / 2;
-  const size_t y_start = (height - box_height) / 2;
-  aocommon::UVector<float> small_image(box_width * box_height);
-  for (size_t y = y_start; y != (height + box_height) / 2; ++y) {
-    std::copy_n(&image[y * width + x_start], box_width,
-                &small_image[(y - y_start) * box_width]);
-  }
-
-  SingleFit2DGaussianCentred(small_image.data(), box_width, box_height,
-                             beam_estimate, beam_major, beam_minor,
-                             beam_phase_angle, verbose);
-}
-
-void GaussianFitter::Fit2DCircularGaussianCentredInBox(
-    const float* image, size_t width, size_t height, double& beam_size,
-    size_t box_width, size_t box_height) {
-  const size_t x_start = (width - box_width) / 2;
-  const size_t y_start = (height - box_height) / 2;
-  aocommon::UVector<float> small_image(box_width * box_height);
-  for (size_t y = y_start; y != (height + box_height) / 2; ++y) {
-    std::copy_n(&image[y * width + x_start], box_width,
-                &small_image[(y - y_start) * box_width]);
-  }
-
-  SingleFit2DCircularGaussianCentred(small_image.data(), box_width, box_height,
-                                     beam_size);
-}
-
-void GaussianFitter::SingleFit2DGaussianCentred(
-    const float* image, size_t width, size_t height, double beam_estimate,
-    double& beam_major, double& beam_minor, double& beam_phase_angle,
-    bool verbose) {
+/**
+ * This function performs a single fit of a Gaussian. The position of the
+ * Gaussian is constrained to be in the centre of the image. The Gaussian is
+ * fitted such that the squared residuals (data - model) are minimal.
+ *
+ * This function is typically used to find the beam-shape of the point-spread
+ * function. The beam estimate is used as initial value for the minor and
+ * major shape.
+ */
+Ellipse SingleFit2DGaussianCentred(const float* image, size_t width,
+                                   size_t height, double beam_estimate,
+                                   bool verbose) {
   if (width * height < 3) {
-    beam_major = std::numeric_limits<double>::quiet_NaN();
-    beam_minor = std::numeric_limits<double>::quiet_NaN();
-    beam_phase_angle = std::numeric_limits<double>::quiet_NaN();
-    return;
+    return {std::numeric_limits<double>::quiet_NaN(),
+            std::numeric_limits<double>::quiet_NaN(),
+            std::numeric_limits<double>::quiet_NaN()};
   }
 
-  width_ = width;
-  height_ = height;
-  image_ = image;
-  scale_factor_ = (width + height) / 2;
+  FitData data;
+  data.width = width;
+  data.height = height;
+  data.image = image;
+  data.scale_factor = (width + height) / 2;
 
   const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmsder;
   gsl_multifit_fdfsolver* solver =
-      gsl_multifit_fdfsolver_alloc(T, width_ * height_, 3);
+      gsl_multifit_fdfsolver_alloc(T, width * height, 3);
 
   gsl_multifit_function_fdf fdf;
   fdf.f = &FittingCentered;
   fdf.df = &FittingDerivativeCentered;
   fdf.fdf = &FittingBothCentered;
-  fdf.n = width_ * height_;
+  fdf.n = width * height;
   fdf.p = 3;
-  fdf.params = this;
+  fdf.params = &data;
 
   // Using the FWHM formula for a Gaussian:
   double initial_values_array[3] = {
-      beam_estimate / (scale_factor_ * double(kSigmaToBeam)),
-      beam_estimate / (scale_factor_ * double(kSigmaToBeam)), 0.0};
+      beam_estimate / (data.scale_factor * double(kSigmaToBeam)),
+      beam_estimate / (data.scale_factor * double(kSigmaToBeam)), 0.0};
   gsl_vector_view initial_values =
       gsl_vector_view_array(initial_values_array, 3);
   gsl_multifit_fdfsolver_set(solver, &fdf, &initial_values.vector);
@@ -580,35 +435,51 @@ void GaussianFitter::SingleFit2DGaussianCentred(
 
   gsl_multifit_fdfsolver_free(solver);
 
-  ToAnglesAndFwhm(sx, sy, beta, beam_major, beam_minor, beam_phase_angle);
-  beam_major *= scale_factor_;
-  beam_minor *= scale_factor_;
+  Ellipse ellipse = ToAnglesAndFwhm(sx, sy, beta);
+  ellipse.major *= data.scale_factor;
+  ellipse.minor *= data.scale_factor;
+  return ellipse;
 }
 
-void GaussianFitter::SingleFit2DCircularGaussianCentred(const float* image,
-                                                        size_t width,
-                                                        size_t height,
-                                                        double& beam_size) {
-  width_ = width;
-  height_ = height;
-  image_ = image;
-  scale_factor_ = (width + height) / 2;
+Ellipse Fit2DGaussianCentredInBox(const float* image, size_t width,
+                                  size_t height, double beam_estimate,
+                                  size_t box_width, size_t box_height,
+                                  bool verbose) {
+  const size_t x_start = (width - box_width) / 2;
+  const size_t y_start = (height - box_height) / 2;
+  aocommon::UVector<float> small_image(box_width * box_height);
+  for (size_t y = y_start; y != (height + box_height) / 2; ++y) {
+    std::copy_n(&image[y * width + x_start], box_width,
+                &small_image[(y - y_start) * box_width]);
+  }
+
+  return SingleFit2DGaussianCentred(small_image.data(), box_width, box_height,
+                                    beam_estimate, verbose);
+}
+
+void SingleFit2DCircularGaussianCentred(const float* image, size_t width,
+                                        size_t height, double& beam_size) {
+  FitData data;
+  data.width = width;
+  data.height = height;
+  data.image = image;
+  data.scale_factor = (width + height) / 2;
 
   const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmsder;
   gsl_multifit_fdfsolver* solver =
-      gsl_multifit_fdfsolver_alloc(T, width_ * height_, 1);
+      gsl_multifit_fdfsolver_alloc(T, width * height, 1);
 
   gsl_multifit_function_fdf fdf;
   fdf.f = &FittingCircularCentered;
   fdf.df = &FittingDerivativeCircularCentered;
   fdf.fdf = &FittingBothCircularCentered;
-  fdf.n = width_ * height_;
+  fdf.n = width * height;
   fdf.p = 1;
-  fdf.params = this;
+  fdf.params = &data;
 
   // Using the FWHM formula for a Gaussian:
   double initial_values_array[1] = {beam_size /
-                                    (scale_factor_ * double(kSigmaToBeam))};
+                                    (data.scale_factor * double(kSigmaToBeam))};
   gsl_vector_view initial_values =
       gsl_vector_view_array(initial_values_array, 1);
   gsl_multifit_fdfsolver_set(solver, &fdf, &initial_values.vector);
@@ -628,80 +499,56 @@ void GaussianFitter::SingleFit2DCircularGaussianCentred(const float* image,
   const double s = gsl_vector_get(solver->x, 0);
   gsl_multifit_fdfsolver_free(solver);
 
-  beam_size = s * kSigmaToBeam * scale_factor_;
+  beam_size = s * kSigmaToBeam * data.scale_factor;
 }
 
-void GaussianFitter::Fit2DGaussianWithAmplitudeInBox(
-    const float* image, size_t width, size_t /*height*/, double& val,
-    double& pos_x, double& pos_y, double& beam_major, double& beam_minor,
-    double& beam_phase_angle, double* floor_level, size_t x_start, size_t x_end,
-    size_t y_start, size_t y_end) {
-  const size_t box_width = x_end - x_start;
-  const size_t box_height = y_end - y_start;
+void Fit2DCircularGaussianCentredInBox(const float* image, size_t width,
+                                       size_t height, double& beam_size,
+                                       size_t box_width, size_t box_height) {
+  const size_t x_start = (width - box_width) / 2;
+  const size_t y_start = (height - box_height) / 2;
   aocommon::UVector<float> small_image(box_width * box_height);
-  for (size_t y = y_start; y != y_end; ++y) {
+  for (size_t y = y_start; y != (height + box_height) / 2; ++y) {
     std::copy_n(&image[y * width + x_start], box_width,
                 &small_image[(y - y_start) * box_width]);
   }
 
-  pos_x -= x_start;
-  pos_y -= y_start;
-  Fit2DGaussianWithAmplitude(small_image.data(), box_width, box_height, val,
-                             pos_x, pos_y, beam_major, beam_minor,
-                             beam_phase_angle, floor_level);
-  pos_x += x_start;
-  pos_y += y_start;
+  SingleFit2DCircularGaussianCentred(small_image.data(), box_width, box_height,
+                                     beam_size);
 }
 
 /**
- * Fits the position, size and amplitude of a Gaussian. If floor_level is not
- * a nullptr, the floor (background level, or zero level) is fitted too.
+ * Like SingleFit2DGaussianCentred(), but includes Gaussian centre X and Y
+ * position and amplitude in the fitted parameters.
+ *
+ * This function can typically be used for source fitting.
  */
-void GaussianFitter::Fit2DGaussianWithAmplitude(
-    const float* image, size_t width, size_t height, double& val, double& pos_x,
-    double& pos_y, double& beam_major, double& beam_minor,
-    double& beam_phase_angle, double* floor_level) {
-  width_ = width;
-  height_ = height;
-  image_ = image;
-  scale_factor_ = (width + height) / 2;
-
-  if (floor_level == nullptr) {
-    Fit2DGaussianWithAmplitude(val, pos_x, pos_y, beam_major, beam_minor,
-                               beam_phase_angle);
-  } else {
-    Fit2DGaussianWithAmplitudeWithFloor(val, pos_x, pos_y, beam_major,
-                                        beam_minor, beam_phase_angle,
-                                        *floor_level);
-  }
-}
-
-void GaussianFitter::Fit2DGaussianWithAmplitude(double& val, double& pos_x,
-                                                double& pos_y,
-                                                double& beam_major,
-                                                double& beam_minor,
-                                                double& beam_phase_angle) {
+void Fit2DGaussianWithAmplitudeNoFloor(const FitData& data, double& val,
+                                       double& pos_x, double& pos_y,
+                                       double& beam_major, double& beam_minor,
+                                       double& beam_position_angle) {
   const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmsder;
   gsl_multifit_fdfsolver* solver =
-      gsl_multifit_fdfsolver_alloc(T, width_ * height_, 6);
+      gsl_multifit_fdfsolver_alloc(T, data.width * data.height, 6);
 
   gsl_multifit_function_fdf fdf;
   fdf.f = &FittingWithAmplitude;
   fdf.df = &FittingDerivativeWithAmplitude;
   fdf.fdf = &FittingBothWithAmplitude;
-  fdf.n = width_ * height_;
+  fdf.n = data.width * data.height;
   fdf.p = 6;
-  fdf.params = this;
+  fdf.params = const_cast<FitData*>(&data);
 
-  // Using the FWHM formula for a Gaussian:
-  x_init_ = -(pos_x - width_ / 2) / scale_factor_;
-  y_init_ = -(pos_y - height_ / 2) / scale_factor_;
+  double x_init = -(pos_x - data.width / 2) / data.scale_factor;
+  double y_init = -(pos_y - data.height / 2) / data.scale_factor;
   double initial_values_array[6] = {
       val,
-      x_init_,
-      y_init_,
-      beam_major / (scale_factor_ * double(kSigmaToBeam)),
-      beam_major / (scale_factor_ * double(kSigmaToBeam)),
+      x_init,
+      y_init,
+      beam_major /
+          (data.scale_factor *
+           double(kSigmaToBeam)),  // Uses the FWHM formula for a Gaussian
+      beam_major / (data.scale_factor * double(kSigmaToBeam)),
       0.0};
   gsl_vector_view initial_values =
       gsl_vector_view_array(initial_values_array, 6);
@@ -720,42 +567,51 @@ void GaussianFitter::Fit2DGaussianWithAmplitude(double& val, double& pos_x,
   } while (status == GSL_CONTINUE && iter < 500);
 
   val = gsl_vector_get(solver->x, 0);
-  pos_x = -1.0 * gsl_vector_get(solver->x, 1) * scale_factor_ + width_ / 2;
-  pos_y = -1.0 * gsl_vector_get(solver->x, 2) * scale_factor_ + height_ / 2;
+  pos_x =
+      -1.0 * gsl_vector_get(solver->x, 1) * data.scale_factor + data.width / 2;
+  pos_y =
+      -1.0 * gsl_vector_get(solver->x, 2) * data.scale_factor + data.height / 2;
   double sx = gsl_vector_get(solver->x, 3), sy = gsl_vector_get(solver->x, 4),
          beta = gsl_vector_get(solver->x, 5);
 
   gsl_multifit_fdfsolver_free(solver);
 
-  ToAnglesAndFwhm(sx, sy, beta, beam_major, beam_minor, beam_phase_angle);
-  beam_major *= scale_factor_;
-  beam_minor *= scale_factor_;
+  ToAnglesAndFwhm(sx, sy, beta, beam_major, beam_minor, beam_position_angle);
+  beam_major *= data.scale_factor;
+  beam_minor *= data.scale_factor;
 }
 
-void GaussianFitter::Fit2DGaussianWithAmplitudeWithFloor(
-    double& val, double& pos_x, double& pos_y, double& beam_major,
-    double& beam_minor, double& beam_phase_angle, double& floor_level) {
+/**
+ * Like Fit2DGaussianWithAmplitude(), but includes floor_level as fitted
+ * parameter. Floor is the background/zero level on which the Gaussian
+ * resides.
+ */
+void Fit2DGaussianWithAmplitudeWithFloor(const FitData& data, double& val,
+                                         double& pos_x, double& pos_y,
+                                         double& beam_major, double& beam_minor,
+                                         double& beam_position_angle,
+                                         double& floor_level) {
   const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmsder;
   gsl_multifit_fdfsolver* solver =
-      gsl_multifit_fdfsolver_alloc(T, width_ * height_, 7);
+      gsl_multifit_fdfsolver_alloc(T, data.width * data.height, 7);
 
   gsl_multifit_function_fdf fdf;
   fdf.f = &FittingWithAmplitudeAndFloor;
   fdf.df = &FittingDerivativeWithAmplitudeAndFloor;
   fdf.fdf = &FittingBothWithAmplitudeAndFloor;
-  fdf.n = width_ * height_;
+  fdf.n = data.width * data.height;
   fdf.p = 7;
-  fdf.params = this;
+  fdf.params = const_cast<FitData*>(&data);
 
   // Using the FWHM formula for a Gaussian:
-  x_init_ = -(pos_x - width_ / 2) / scale_factor_;
-  y_init_ = -(pos_y - height_ / 2) / scale_factor_;
+  const double x_init = -(pos_x - data.width / 2) / data.scale_factor;
+  const double y_init = -(pos_y - data.height / 2) / data.scale_factor;
   double initial_values_array[7] = {
       val,
-      x_init_,
-      y_init_,
-      beam_major / (scale_factor_ * double(kSigmaToBeam)),
-      beam_major / (scale_factor_ * double(kSigmaToBeam)),
+      x_init,
+      y_init,
+      beam_major / (data.scale_factor * double(kSigmaToBeam)),
+      beam_major / (data.scale_factor * double(kSigmaToBeam)),
       0.0,
       0.0};
   gsl_vector_view initial_values =
@@ -775,8 +631,10 @@ void GaussianFitter::Fit2DGaussianWithAmplitudeWithFloor(
   } while (status == GSL_CONTINUE && iter < 500);
 
   val = gsl_vector_get(solver->x, 0);
-  pos_x = -1.0 * gsl_vector_get(solver->x, 1) * scale_factor_ + width_ / 2;
-  pos_y = -1.0 * gsl_vector_get(solver->x, 2) * scale_factor_ + height_ / 2;
+  pos_x =
+      -1.0 * gsl_vector_get(solver->x, 1) * data.scale_factor + data.width / 2;
+  pos_y =
+      -1.0 * gsl_vector_get(solver->x, 2) * data.scale_factor + data.height / 2;
   double sx = gsl_vector_get(solver->x, 3);
   double sy = gsl_vector_get(solver->x, 4);
   double beta = gsl_vector_get(solver->x, 5);
@@ -784,9 +642,176 @@ void GaussianFitter::Fit2DGaussianWithAmplitudeWithFloor(
 
   gsl_multifit_fdfsolver_free(solver);
 
-  ToAnglesAndFwhm(sx, sy, beta, beam_major, beam_minor, beam_phase_angle);
-  beam_major *= scale_factor_;
-  beam_minor *= scale_factor_;
+  ToAnglesAndFwhm(sx, sy, beta, beam_major, beam_minor, beam_position_angle);
+  beam_major *= data.scale_factor;
+  beam_minor *= data.scale_factor;
+}
+
+/**
+ * Fits the position, size and amplitude of a Gaussian. If floor_level is not
+ * a nullptr, the floor (background level, or zero level) is fitted too.
+ */
+void Fit2DGaussianWithAmplitude(const float* image, size_t width, size_t height,
+                                double& val, double& pos_x, double& pos_y,
+                                double& beam_major, double& beam_minor,
+                                double& beam_position_angle,
+                                double* floor_level) {
+  FitData data;
+  data.width = width;
+  data.height = height;
+  data.image = image;
+  data.scale_factor = (width + height) / 2;
+
+  if (floor_level == nullptr) {
+    Fit2DGaussianWithAmplitudeNoFloor(data, val, pos_x, pos_y, beam_major,
+                                      beam_minor, beam_position_angle);
+  } else {
+    Fit2DGaussianWithAmplitudeWithFloor(data, val, pos_x, pos_y, beam_major,
+                                        beam_minor, beam_position_angle,
+                                        *floor_level);
+  }
+}
+
+void Fit2DGaussianWithAmplitudeInBox(
+    const float* image, size_t width, size_t /*height*/, double& val,
+    double& pos_x, double& pos_y, double& beam_major, double& beam_minor,
+    double& beam_position_angle, double* floor_level, size_t x_start,
+    size_t x_end, size_t y_start, size_t y_end) {
+  const size_t box_width = x_end - x_start;
+  const size_t box_height = y_end - y_start;
+  aocommon::UVector<float> small_image(box_width * box_height);
+  for (size_t y = y_start; y != y_end; ++y) {
+    std::copy_n(&image[y * width + x_start], box_width,
+                &small_image[(y - y_start) * box_width]);
+  }
+
+  pos_x -= x_start;
+  pos_y -= y_start;
+  Fit2DGaussianWithAmplitude(small_image.data(), box_width, box_height, val,
+                             pos_x, pos_y, beam_major, beam_minor,
+                             beam_position_angle, floor_level);
+  pos_x += x_start;
+  pos_y += y_start;
+}
+
+}  // namespace
+
+Ellipse Fit2DGaussianCentred(const float* image, size_t width, size_t height,
+                             double beam_estimate, double box_scale_factor,
+                             bool verbose) {
+  size_t preferred_size = std::max<size_t>(
+      std::ceil(box_scale_factor), std::ceil(beam_estimate * box_scale_factor));
+  if (preferred_size % 2 != 0) ++preferred_size;
+  Ellipse result;
+  if (preferred_size < width || preferred_size < height) {
+    size_t n_iterations = 0;
+    bool box_was_large_enough;
+    do {
+      size_t box_width = std::min(preferred_size, width);
+      size_t box_height = std::min(preferred_size, height);
+      if (verbose) std::cout << "Fit initial value:" << beam_estimate << "\n";
+      result = Fit2DGaussianCentredInBox(image, width, height, beam_estimate,
+                                         box_width, box_height, verbose);
+      if (verbose) {
+        std::cout << "Fit result:" << result.major << " x " << result.minor
+                  << " px, " << result.position_angle << " (box was "
+                  << box_width << " x " << box_height << ")\n";
+      }
+
+      box_was_large_enough =
+          (result.major * box_scale_factor * 0.8 < box_width ||
+           box_width >= width) &&
+          (result.major * box_scale_factor * 0.8 < box_height ||
+           box_height >= height);
+      if (!box_was_large_enough) {
+        preferred_size =
+            std::max<size_t>(std::ceil(box_scale_factor),
+                             std::ceil(result.major * box_scale_factor));
+        if (preferred_size % 2 != 0) ++preferred_size;
+        beam_estimate = std::max(result.major, beam_estimate);
+      }
+      ++n_iterations;
+    } while (!box_was_large_enough && n_iterations < 5);
+  } else {
+    if (verbose) std::cout << "Image is as large as the fitting box.\n";
+    result = SingleFit2DGaussianCentred(image, width, height, beam_estimate,
+                                        verbose);
+  }
+  return result;
+}
+
+void Fit2DCircularGaussianCentred(const float* image, size_t width,
+                                  size_t height, double& beam_size,
+                                  double box_scale_factor) {
+  const double initial_value = beam_size;
+  size_t preferred_size = std::max<size_t>(
+      std::ceil(box_scale_factor), std::ceil(beam_size * box_scale_factor));
+  if (preferred_size % 2 != 0) ++preferred_size;
+  if (preferred_size < width || preferred_size < height) {
+    size_t box_width = std::min(preferred_size, width);
+    size_t box_height = std::min(preferred_size, height);
+    size_t n_iterations = 0;
+    bool box_was_large_enough;
+    do {
+      Fit2DCircularGaussianCentredInBox(image, width, height, beam_size,
+                                        box_width, box_height);
+
+      box_was_large_enough = (beam_size * box_scale_factor * 0.8 < box_width ||
+                              width >= box_width) &&
+                             (beam_size * box_scale_factor * 0.8 < box_height ||
+                              height >= box_height);
+      if (!box_was_large_enough) {
+        preferred_size =
+            std::max<size_t>(std::ceil(box_scale_factor),
+                             std::ceil(beam_size * box_scale_factor));
+        if (preferred_size % 2 != 0) ++preferred_size;
+        beam_size = std::max(initial_value, beam_size);
+      }
+      ++n_iterations;
+    } while (!box_was_large_enough && n_iterations < 5);
+  } else {
+    SingleFit2DCircularGaussianCentred(image, width, height, beam_size);
+  }
+}
+
+void Fit2DGaussianFull(const float* image, size_t width, size_t height,
+                       double& val, double& pos_x, double& pos_y,
+                       double& beam_major, double& beam_minor,
+                       double& beam_position_angle, double* floor_level) {
+  size_t preferred_size = std::max<size_t>(10, std::ceil(beam_major * 10.0));
+  if (preferred_size % 2 != 0) ++preferred_size;
+  if (preferred_size < width || preferred_size < height) {
+    const size_t x_start =
+        std::max<int>(0, int(std::round(pos_x)) - int(preferred_size) / 2);
+    const size_t x_end =
+        std::min(width, size_t(std::round(pos_x)) + preferred_size / 2);
+    const size_t y_start =
+        std::max<int>(0, int(std::round(pos_y)) - int(preferred_size) / 2);
+    const size_t y_end =
+        std::min(height, size_t(std::round(pos_y)) + preferred_size / 2);
+    size_t n_iterations = 0;
+    bool box_was_large_enough;
+    do {
+      Fit2DGaussianWithAmplitudeInBox(
+          image, width, height, val, pos_x, pos_y, beam_major, beam_minor,
+          beam_position_angle, floor_level, x_start, x_end, y_start, y_end);
+
+      const size_t box_width = x_end - x_start;
+      const size_t box_height = y_end - y_start;
+      box_was_large_enough =
+          (beam_major * 4.0 < box_width || width >= box_width) &&
+          (beam_major * 4.0 < box_height || height >= box_height);
+      if (!box_was_large_enough) {
+        preferred_size = std::max<size_t>(10, std::ceil(beam_major * 10.0));
+        if (preferred_size % 2 != 0) ++preferred_size;
+      }
+      ++n_iterations;
+    } while (!box_was_large_enough && n_iterations < 5);
+  } else {
+    Fit2DGaussianWithAmplitude(image, width, height, val, pos_x, pos_y,
+                               beam_major, beam_minor, beam_position_angle,
+                               floor_level);
+  }
 }
 
 Ellipse DeconvolveGaussian(const Ellipse& a, const Ellipse& b) {
@@ -841,5 +866,4 @@ Ellipse DeconvolveGaussian(const Ellipse& a, const Ellipse& b) {
   return result;
 }
 
-}  // namespace fitters
-}  // namespace schaapcommon
+}  // namespace schaapcommon::fitters
diff --git a/src/fitters/test/tgaussianfitter.cc b/src/fitters/test/tgaussianfitter.cc
index 0f8e3a6916cd84d27e61ea0074be406b87e2e00c..52ebb200d3cd50da866f7dcfe4b9c0be9ca165c1 100644
--- a/src/fitters/test/tgaussianfitter.cc
+++ b/src/fitters/test/tgaussianfitter.cc
@@ -18,12 +18,14 @@ constexpr size_t kHeight = 64;
 constexpr long double kPixelSize = 1 /*amin*/ * (M_PI / 180.0 / 60.0);
 }  // namespace
 
+using schaapcommon::fitters::Ellipse;
+
 BOOST_AUTO_TEST_SUITE(gaussian_fitter)
 
 BOOST_AUTO_TEST_CASE(fit) {
   aocommon::ThreadPool::GetInstance().SetNThreads(kThreadCount);
-  for (size_t beam_phase_angle_index = 0; beam_phase_angle_index != 10;
-       ++beam_phase_angle_index) {
+  for (size_t beam_position_angle_index = 0; beam_position_angle_index != 10;
+       ++beam_position_angle_index) {
     const size_t width = 512, height = 512;
     aocommon::Image model(width, height, 0.0);
     aocommon::Image restored(width, height, 0.0);
@@ -31,26 +33,65 @@ BOOST_AUTO_TEST_CASE(fit) {
     const long double kPixelSize = 1.0L /*amin*/ * (M_PI / 180.0 / 60.0);
     const long double beam_major = 20.0L * kPixelSize;
     const long double beam_minor = 5.0L * kPixelSize;
-    const long double beam_phase_angle = beam_phase_angle_index * M_PI / 10.0;
+    const long double beam_position_angle =
+        beam_position_angle_index * M_PI / 10.0;
+
+    schaapcommon::fft::MakeFftwfPlannerThreadSafe();
+    schaapcommon::fft::RestoreImage(
+        restored.Data(), model.Data(), width, height, beam_major, beam_minor,
+        beam_position_angle, kPixelSize, kPixelSize);
+
+    const Ellipse ellipse = schaapcommon::fitters::Fit2DGaussianCentred(
+        restored.Data(), width, height, 5.0, 10.0, false);
+
+    BOOST_CHECK_CLOSE_FRACTION(ellipse.major, 20.0, 1.0e-3);
+    BOOST_CHECK_CLOSE_FRACTION(ellipse.minor, 5.0, 1.0e-3);
+    const double fit_position_angle =
+        std::fmod((ellipse.position_angle + 2.0 * M_PI), M_PI);
+    BOOST_CHECK_CLOSE_FRACTION(fit_position_angle, beam_position_angle, 1.0e-3);
+  }
+}
+
+BOOST_AUTO_TEST_CASE(fit_full) {
+  aocommon::ThreadPool::GetInstance().SetNThreads(kThreadCount);
+  const size_t width = 512;
+  const size_t height = 512;
+  const size_t source_x = 200;
+  const size_t source_y = 300;
+  for (size_t beam_position_angle_index = 0; beam_position_angle_index != 10;
+       ++beam_position_angle_index) {
+    aocommon::Image model(width, height, 0.0);
+    aocommon::Image restored(width, height, 0.0);
+    model[(source_y * width) + source_x] = 2.0;
+    const long double kPixelSize = 1.0L /*amin*/ * (M_PI / 180.0 / 60.0);
+    long double beam_major = 20.0L * kPixelSize;
+    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(restored.Data(), model.Data(), width,
-                                    height, beam_major, beam_minor,
-                                    beam_phase_angle, kPixelSize, kPixelSize);
-
-    schaapcommon::fitters::GaussianFitter fitter;
-    // Check that Fit2DGaussianCentred updates these variables.
-    double fit_major = 0.0;  // Should become 20.0.
-    double fit_minor = 0.0;  // Should become 5.0.
-    double fit_phase_angle =
-        -1.0;  // Should become beam_phase_angle, which can be 0.0.
-    fitter.Fit2DGaussianCentred(restored.Data(), width, height, 5.0, fit_major,
-                                fit_minor, fit_phase_angle, 10.0, false);
-    fit_phase_angle = std::fmod((fit_phase_angle + 2.0 * M_PI), M_PI);
-
-    BOOST_CHECK_CLOSE_FRACTION(fit_major, 20.0, 1.0e-3);
-    BOOST_CHECK_CLOSE_FRACTION(fit_minor, 5.0, 1.0e-3);
-    BOOST_CHECK_CLOSE_FRACTION(fit_phase_angle, beam_phase_angle, 1.0e-3);
+    schaapcommon::fft::RestoreImage(
+        restored.Data(), model.Data(), width, height, beam_major, beam_minor,
+        beam_position_angle, kPixelSize, kPixelSize);
+
+    Ellipse fit_ellipse;
+    // The following four parameters are used as initial values.
+    fit_ellipse.major = 10.0;
+    double fit_amplitude = 1.0;
+    double fit_x = source_x + 3;
+    double fit_y = source_y - 3;
+    schaapcommon::fitters::Fit2DGaussianFull(
+        restored.Data(), width, height, fit_amplitude, fit_x, fit_y,
+        fit_ellipse.major, fit_ellipse.minor, fit_ellipse.position_angle,
+        nullptr);
+
+    BOOST_CHECK_CLOSE_FRACTION(fit_amplitude, 2.0, 1.0e-3);
+    BOOST_CHECK_CLOSE_FRACTION(fit_x, source_x, 1.0e-3);
+    BOOST_CHECK_CLOSE_FRACTION(fit_y, source_y, 1.0e-3);
+    BOOST_CHECK_CLOSE_FRACTION(fit_ellipse.major, 20.0, 1.0e-3);
+    BOOST_CHECK_CLOSE_FRACTION(fit_ellipse.minor, 5.0, 1.0e-3);
+    const double fit_position_angle =
+        std::fmod((fit_ellipse.position_angle + 2.0 * M_PI), M_PI);
+    BOOST_CHECK_CLOSE_FRACTION(fit_position_angle, beam_position_angle, 1.0e-3);
   }
 }
 
@@ -59,13 +100,8 @@ BOOST_AUTO_TEST_CASE(insufficient_data_fit) {
   aocommon::Image model(width, height, 0.0);
   aocommon::Image restored(width, height, 0.0);
   model[((height / 2) * width) + (width / 2)] = 1.0;
-  schaapcommon::fitters::GaussianFitter fitter;
-  double fit_major = 0.0;
-  double fit_minor = 0.0;
-  double fit_phase_angle = 1.0;
-  BOOST_CHECK_NO_THROW(fitter.Fit2DGaussianCentred(
-      restored.Data(), width, height, 1.0, fit_major, fit_minor,
-      fit_phase_angle, 10.0, false));
+  BOOST_CHECK_NO_THROW(schaapcommon::fitters::Fit2DGaussianCentred(
+      restored.Data(), width, height, 1.0, 10.0, false));
 }
 
 BOOST_AUTO_TEST_CASE(fit_circular) {
@@ -76,27 +112,22 @@ BOOST_AUTO_TEST_CASE(fit_circular) {
 
   const long double beam_major = 4.0L * kPixelSize;
   const long double beam_minor = 4.0L * kPixelSize;
-  const long double beam_phase_angle = 0.0;
+  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_phase_angle, kPixelSize, kPixelSize);
-
-  schaapcommon::fitters::GaussianFitter fitter;
-  // Check that Fit2DGaussianCentred updates these variables.
-  double fit_major = 0.0;         // Should become 4.0.
-  double fit_minor = 0.0;         // Should become 4.0.
-  double fit_phase_angle = -1.0;  // Should become beam_phase_angle (0.0).
-  fitter.Fit2DGaussianCentred(
+                                  beam_position_angle, kPixelSize, kPixelSize);
+
+  const Ellipse result = schaapcommon::fitters::Fit2DGaussianCentred(
       restored.Data(), restored.Width(), restored.Height(),
-      estimated_beam_pixel, fit_major, fit_minor, fit_phase_angle, 10.0, false);
+      estimated_beam_pixel, 10.0, false);
 
-  BOOST_CHECK_CLOSE_FRACTION(fit_major, 4.0, 1.0e-4);
-  BOOST_CHECK_CLOSE_FRACTION(fit_minor, 4.0, 1.0e-4);
-  BOOST_CHECK_SMALL(
-      std::abs(fit_phase_angle - static_cast<double>(beam_phase_angle)),
-      1.0e-4);
+  BOOST_CHECK_CLOSE_FRACTION(result.major, 4.0, 1.0e-4);
+  BOOST_CHECK_CLOSE_FRACTION(result.minor, 4.0, 1.0e-4);
+  BOOST_CHECK_SMALL(std::abs(result.position_angle -
+                             static_cast<double>(beam_position_angle)),
+                    1.0e-4);
 }
 
 BOOST_AUTO_TEST_CASE(little_data_circular_fit) {
@@ -104,10 +135,9 @@ BOOST_AUTO_TEST_CASE(little_data_circular_fit) {
   aocommon::Image model(width, height, 0.0);
   aocommon::Image restored(width, height, 0.0);
   model[((height / 2) * width) + (width / 2)] = 1.0;
-  schaapcommon::fitters::GaussianFitter fitter;
   double fit = 0.0;
-  BOOST_CHECK_NO_THROW(
-      fitter.Fit2DCircularGaussianCentred(restored.Data(), width, height, fit));
+  BOOST_CHECK_NO_THROW(schaapcommon::fitters::Fit2DCircularGaussianCentred(
+      restored.Data(), width, height, fit));
 }
 
 BOOST_AUTO_TEST_CASE(fit_small_beam) {
@@ -119,33 +149,27 @@ BOOST_AUTO_TEST_CASE(fit_small_beam) {
 
   const long double beam_major = 4.0L * kPixelSize;
   const long double beam_minor = 0.5L * kPixelSize;
-  const long double beam_phase_angle = 0.0;
+  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_phase_angle, kPixelSize, kPixelSize);
-
-  schaapcommon::fitters::GaussianFitter fitter;
-  // Check that Fit2DGaussianCentred updates these variables.
-  double fit_major = 0.0;         // Should become 4.0.
-  double fit_minor = 0.0;         // Should become 0.5.
-  double fit_phase_angle = -1.0;  // Should become beam_phase_angle (0.0).
-  fitter.Fit2DGaussianCentred(
+                                  beam_position_angle, kPixelSize, kPixelSize);
+
+  Ellipse ellipse = schaapcommon::fitters::Fit2DGaussianCentred(
       restored.Data(), restored.Width(), restored.Height(),
-      estimated_beam_pixel, fit_major, fit_minor, fit_phase_angle, 10.0, false);
+      estimated_beam_pixel, 10.0, false);
 
-  BOOST_CHECK_CLOSE_FRACTION(fit_major, 4.0, 1.0e-4);
-  BOOST_CHECK_CLOSE_FRACTION(fit_minor, 0.5, 1.0e-4);
-  BOOST_CHECK_SMALL(
-      std::abs(fit_phase_angle - static_cast<double>(beam_phase_angle)),
-      1.0e-4);
+  BOOST_CHECK_CLOSE_FRACTION(ellipse.major, 4.0, 1.0e-4);
+  BOOST_CHECK_CLOSE_FRACTION(ellipse.minor, 0.5, 1.0e-4);
+  BOOST_CHECK_SMALL(std::abs(ellipse.position_angle -
+                             static_cast<double>(beam_position_angle)),
+                    1.0e-4);
 }
 
 BOOST_AUTO_TEST_CASE(deconvolve) {
   using schaapcommon::fitters::DeconvolveGaussian;
-  using schaapcommon::fitters::Ellipse;
   Ellipse result;
 
   // Evaluate all four kwadrants.