Repository: transientskp/pyse

Issue:
max_pix_variance_factor does not make a measurable difference
When using moments, the peak spectral brightness is the maximum pixel values across an island times a correction factor called [fudge_max_pix](https://github.com/transientskp/pyse/blob/a9c8d875a1efdaaccbd6baadac22473f2254dc04/sourcefinder/image.py#L100) to take account of the fact that the maximum pixel is an underestimate of the true peak spectral brightness since the peak will typically not be centered on a pixel.

This correction introduced its own variance, called [max_pix_variance_factor](https://github.com/transientskp/pyse/blob/a9c8d875a1efdaaccbd6baadac22473f2254dc04/sourcefinder/image.py#L101C14-L101C37). 
However, this additional variance is so small compared to the rms noise that we might as well disregard it.
To clarify: `max_pix_variance_factor` is very close to 1.

Discussion:

Patch:
diff --git a/sourcefinder/extract.py b/sourcefinder/extract.py
index a586ee6..65ac0bd 100644
--- a/sourcefinder/extract.py
+++ b/sourcefinder/extract.py
@@ -204,16 +204,16 @@ class Island(object):
         """Deviation"""
         return (self.data / self.rms_orig).max()
 
-    def fit(self, fudge_max_pix_factor, max_pix_variance_factor, beamsize,
-            correlation_lengths, fixed=None):
+    def fit(self, fudge_max_pix_factor, beamsize, correlation_lengths,
+            fixed=None):
         """Fit the position"""
         try:
             measurement, gauss_island, gauss_residual = \
                 source_profile_and_errors(self.data, self.threshold(),
                                           self.rms, self.noise(), self.beam,
                                           fudge_max_pix_factor,
-                                          max_pix_variance_factor, beamsize,
-                                          correlation_lengths, fixed=fixed)
+                                          beamsize, correlation_lengths,
+                                          fixed=fixed)
         except ValueError:
             # Fitting failed
             logger.error("Moments & Gaussian fitting failed at %s" % (
@@ -326,55 +326,54 @@ class ParamSet(MutableMapping):
         boolean entry is used to loosen a bound when a fit becomes unfeasible,
         see the documentation on scipy.optimize.Bounds.
         """
-        if self.moments:
-            if hasattr(self["peak"], "value"):
-                self.bounds["peak"] = (0.5 * self["peak"].value,
-                                       1.5 * self["peak"].value, False)
-            else:
-                self.bounds["peak"] = (0.5 * self["peak"], 1.5 * self["peak"],
-                                       False)
-            # Accommodate for negative heights. The "bounds" argument in
-            # scipy.optimize.least-squares demands that the lower bound is
-            # smaller than the upper bound and with the (0.5, 1.5) fractions
-            # this does not work out well for negative heights. Background:
-            # only to be expected in images from visibilities from polarization
-            # products other than Stokes I.
-            if self.bounds["peak"][0] > self.bounds["peak"][1]:
-                true_upper = self.bounds["peak"][0]
-                true_lower = self.bounds["peak"][1]
-                self.bounds["peak"] = (true_lower, true_upper, False)
-
-            if hasattr(self["semimajor"], "value"):
-                self.bounds["semimajor"] = (0.5 * self["semimajor"].value,
-                                            1.5 * self["semimajor"].value,
-                                            False)
-            else:
-                self.bounds["semimajor"] = (0.5 * self["semimajor"],
-                                            1.5 * self["semimajor"], False)
+        if hasattr(self["peak"], "value"):
+            self.bounds["peak"] = (0.5 * self["peak"].value,
+                                   1.5 * self["peak"].value, False)
+        else:
+            self.bounds["peak"] = (0.5 * self["peak"], 1.5 * self["peak"],
+                                   False)
+        # Accommodate for negative heights. The "bounds" argument in
+        # scipy.optimize.least-squares demands that the lower bound is
+        # smaller than the upper bound and with the (0.5, 1.5) fractions
+        # this does not work out well for negative heights. Background:
+        # only to be expected in images from visibilities from polarization
+        # products other than Stokes I.
+        if self.bounds["peak"][0] > self.bounds["peak"][1]:
+            true_upper = self.bounds["peak"][0]
+            true_lower = self.bounds["peak"][1]
+            self.bounds["peak"] = (true_lower, true_upper, False)
+
+        if hasattr(self["semimajor"], "value"):
+            self.bounds["semimajor"] = (0.5 * self["semimajor"].value,
+                                        1.5 * self["semimajor"].value,
+                                        False)
+        else:
+            self.bounds["semimajor"] = (0.5 * self["semimajor"],
+                                        1.5 * self["semimajor"], False)
 
-            if hasattr(self["semiminor"], "value"):
-                self.bounds["semiminor"] = (0.5 * self["semiminor"].value,
-                                            1.5 * self["semiminor"].value,
-                                            False)
-            else:
-                self.bounds["semiminor"] = (0.5 * self["semiminor"],
-                                            1.5 * self["semiminor"], False)
-
-            self.bounds["xbar"] = (0., data_shape[0], False)
-            self.bounds["ybar"] = (0., data_shape[1], False)
-            # The upper bound for theta is a bit odd, one would expect
-            # np.pi/2 here, but that will yield imperfect fits in
-            # cases where the axes are aligned with the coordinate axes,
-            # which, in turn, breaks the AxesSwapGaussTest.testFitHeight
-            # and AxesSwapGaussTest.testFitSize unit tests. Thus, some
-            # margin needs to be applied.
-            # keep_feasibly is True here to accommodate for a fitting process
-            # where the margin is exceeded.
-            self.bounds["theta"] = (-np.pi/2, np.pi, True)
+        if hasattr(self["semiminor"], "value"):
+            self.bounds["semiminor"] = (0.5 * self["semiminor"].value,
+                                        1.5 * self["semiminor"].value,
+                                        False)
+        else:
+            self.bounds["semiminor"] = (0.5 * self["semiminor"],
+                                        1.5 * self["semiminor"], False)
+
+        self.bounds["xbar"] = (0., data_shape[0], False)
+        self.bounds["ybar"] = (0., data_shape[1], False)
+        # The upper bound for theta is a bit odd, one would expect
+        # np.pi/2 here, but that will yield imperfect fits in
+        # cases where the axes are aligned with the coordinate axes,
+        # which, in turn, breaks the AxesSwapGaussTest.testFitHeight
+        # and AxesSwapGaussTest.testFitSize unit tests. Thus, some
+        # margin needs to be applied.
+        # keep_feasibly is True here to accommodate for a fitting process
+        # where the margin is exceeded.
+        self.bounds["theta"] = (-np.pi/2, np.pi, True)
 
         return self
 
-    def calculate_errors(self, noise, max_pix_variance_factor, correlation_lengths, threshold):
+    def calculate_errors(self, noise, correlation_lengths, threshold):
         """Calculate positional errors
 
         Uses _condon_formulae() if this object is based on a Gaussian fit,
@@ -383,17 +382,16 @@ class ParamSet(MutableMapping):
 
         if self.gaussian:
             return self._condon_formulae(noise, correlation_lengths)
-        elif self.moments or self.thin_detection:
-            # If thin_detection == True, error_bars_from_moments probably
+        else:
+            # If thin_detection == True, _error_bars_from_moments probably
             # gives a better estimate of the true errors than errors from
-            # Gaussian fits, since they are generally larger. May still be
-            # underestimating the true errors, though.
+            # Gaussian fits, i.e. from _condon_formulae, since they are
+            # generally larger. _error_bars_from_moments could still be
+            # underestimating the true errors in case of thin detections.
             if not threshold:
                 threshold = 0
-            return self._error_bars_from_moments(noise, max_pix_variance_factor,
-                                                 correlation_lengths, threshold)
-        else:
-            return False
+            return self._error_bars_from_moments(noise, correlation_lengths,
+                                                 threshold)
 
     def _condon_formulae(self, noise, correlation_lengths):
         """Returns the errors on parameters from Gaussian fits according to
@@ -428,7 +426,6 @@ class ParamSet(MutableMapping):
 
         rho1 = np.sqrt(rho_sq1)
         rho2 = np.sqrt(rho_sq2)
-        rho3 = np.sqrt(rho_sq3)
 
         denom1 = np.sqrt(2. * np.log(2.)) * rho1
         denom2 = np.sqrt(2. * np.log(2.)) * rho2
@@ -462,8 +459,14 @@ class ParamSet(MutableMapping):
         if errortheta > np.pi:
             errortheta = np.pi
 
-        peak += -noise ** 2 / peak + self.clean_bias
+        # Formula (34) of the NVSS paper, without the general fitting bias,
+        # since we don't observe that. On the contrary, there may be peak
+        # fitting bias, but that has the opposite sign, i.e. fitted peaks are
+        # systematically biased too low, but more investigation is needed to
+        # quantify that effect.
+        peak += self.clean_bias
 
+        # Formula (37) of the NVSS paper.
         errorpeaksq = ((self.frac_flux_cal_error * peak) ** 2 +
                        self.clean_bias_error ** 2 +
                        2. * peak ** 2 / rho_sq3)
@@ -473,11 +476,17 @@ class ParamSet(MutableMapping):
         help1 = (errorsmaj / smaj) ** 2
         help2 = (errorsmin / smin) ** 2
         help3 = theta_B * theta_b / (4. * smaj * smin)
+
+        # Since we recomputed the peak spectral brightness following formula
+        # (34) from the NVSS paper, we also need to recompute the flux density
+        # following formula (35) of that paper.
+        flux *= peak / self['peak'].value
+
         errorflux = np.abs(flux) * np.sqrt(
             errorpeaksq / peak ** 2 + help3 * (help1 + help2))
 
         self['peak'] = Uncertain(peak, errorpeak)
-        self['flux'].error = errorflux
+        self['flux'] = Uncertain(flux, errorflux)
         self['xbar'].error = errorx
         self['ybar'].error = errory
         self['semimajor'].error = errorsmaj
@@ -486,7 +495,7 @@ class ParamSet(MutableMapping):
 
         return self
 
-    def _error_bars_from_moments(self, noise, max_pix_variance_factor, correlation_lengths,
+    def _error_bars_from_moments(self, noise, correlation_lengths,
                                  threshold):
         """Provide reasonable error estimates from the moments"""
 
@@ -559,8 +568,7 @@ class ParamSet(MutableMapping):
         # the error from the (corrected) maximum pixel method,
         # while it is part of the expression for rho_sq above.
         errorpeaksq = ((frac_flux_cal_error * peak) ** 2 +
-                       clean_bias_error ** 2 + noise ** 2 +
-                       max_pix_variance_factor * peak ** 2)
+                       clean_bias_error ** 2 + noise ** 2)
         errorpeak = np.sqrt(errorpeaksq)
 
         help1 = (errorsmaj / smaj) ** 2
@@ -692,9 +700,9 @@ class ParamSet(MutableMapping):
         return self
 
 
-def source_profile_and_errors(data, threshold, rms, noise,
-                              beam, fudge_max_pix_factor, max_pix_variance_factor,
-                              beamsize, correlation_lengths, fixed=None):
+def source_profile_and_errors(data, threshold, rms, noise, beam,
+                              fudge_max_pix_factor, beamsize,
+                              correlation_lengths, fixed=None):
     """Return a number of measurable properties with errorbars
 
     Given an island of pixels it will return a number of measurable
@@ -708,7 +716,8 @@ def source_profile_and_errors(data, threshold, rms, noise,
 
     Args:
 
-        data (np.ndarray): array of pixel values, can be a masked
+        data: np.ma.Maskedarray or np.ndarray
+            Array of pixel values, can be a masked
             array, which is necessary for proper Gauss fitting,
             because the pixels below the threshold in the corners and
             along the edges should not be included in the fitting
@@ -728,10 +737,6 @@ def source_profile_and_errors(data, threshold, rms, noise,
         fudge_max_pix_factor(float): Correct for the underestimation of the peak
                                      by taking the maximum pixel value.
 
-        max_pix_variance_factor(float): Take account of additional variance
-                                        induced by the maximum pixel method,
-                                        on top of the background noise.
-
         beamsize(float): The FWHM size of the clean beam
 
         correlation_lengths(tuple): Tuple of two floats describing the distance
@@ -770,12 +775,15 @@ def source_profile_and_errors(data, threshold, rms, noise,
     # We can always find the maximum pixel value and derive the barycenter
     # position. The other three Gaussian parameters we can copy from the clean
     # beam.
+
+    peak = data.max()
     # Are we fitting a -ve or +ve Gaussian?
-    if data.mean() >= 0:
-        # The peak is always underestimated when you take the highest pixel.
-        peak = data.max() * fudge_max_pix_factor
-    else:
+    if peak <= 0:
         peak = data.min()
+    # The peak is always underestimated when you take the highest (or lowest)
+    # pixel.
+    peak *= fudge_max_pix_factor
+
     total = data.sum()
     x, y = np.indices(data.shape)
     xbar = float((x * data).sum() / total)
@@ -783,7 +791,7 @@ def source_profile_and_errors(data, threshold, rms, noise,
 
     param.update({
         "peak": peak,
-        "flux": peak,
+        "flux": np.pi * peak * beam[0] * beam[1] / beamsize,
         "xbar": xbar,
         "ybar": ybar,
         "semimajor": beam[0],
@@ -791,6 +799,9 @@ def source_profile_and_errors(data, threshold, rms, noise,
         "theta": beam[2]
     })
 
+    if fixed:
+        param.update(fixed)
+
     # data_as_ones is constructed to help determine if the island has enough
     # width for Gauss fitting.
     try:
@@ -805,20 +816,22 @@ def source_profile_and_errors(data, threshold, rms, noise,
     minimum_width = min(max_along_x, max_along_y)
 
     if minimum_width > 2:
+        # We can realistically compute moments and perform Gaussian fits if
+        # the island has a thickness of more than 2 in both dimensions.
         param.thin_detection = False
-        # If the island or subisland has thickness of more than 2 in both
-        # dimensions we can properly compute moments and fit a Gaussian with
-        # six free parameters.
-        try:
-            # We can compute moments and possibly try fitting if the island has a
-            # thickness of more than 2 in both dimensions.
-            param.update(fitting.moments(data, fudge_max_pix_factor, beamsize,
-                                         moments_threshold))
-            param.moments = True
-        except ValueError:
-            logger.warning('Moments computations failed, use defaults.')
+        if not fixed:
+            # Moments can only be computed if no parameters are fixed.
+            try:
+                param.update(fitting.moments(data, fudge_max_pix_factor, beam,
+                                             beamsize, moments_threshold))
+                if moments_threshold:
+                    # A complete moments computation is only possible if we have
+                    # imposed a non-zero threshold (and no parameters are fixed).
+                    param.moments = True
+                    param.compute_bounds(data.shape)
+            except ValueError:
+                logger.warning('Moments computations failed, use defaults.')
         try:
-            param.compute_bounds(data.shape)
             gaussian_soln = fitting.fitgaussian(data, param, fixed=fixed,
                                                 bounds=param.bounds)
             param.update(gaussian_soln)
@@ -831,13 +844,9 @@ def source_profile_and_errors(data, threshold, rms, noise,
                      " Proceeding with defaults %s""",
                      str(param))
 
-    if fixed and not param.gaussian:
-        # moments can't handle fixed params
-        raise ValueError("fit failed with given fixed parameters")
-
     param["flux"] = (np.pi * param["peak"] * param["semimajor"] *
                      param["semiminor"] / beamsize)
-    param.calculate_errors(noise, max_pix_variance_factor, correlation_lengths, threshold)
+    param.calculate_errors(noise, correlation_lengths, threshold)
     param.deconvolve_from_clean_beam(beam)
 
     # Calculate residuals
@@ -1425,7 +1434,6 @@ def source_measurements_pixels_and_celestial_vectorised(num_islands, npixs,
                                                         indices,
                                                         labelled_data, labels, wcs,
                                                         fudge_max_pix_factor,
-                                                        max_pix_variance_factor,
                                                         beam, beamsize,
                                                         correlation_lengths,
                                                         eps_ra, eps_dec):
@@ -1524,10 +1532,6 @@ def source_measurements_pixels_and_celestial_vectorised(num_islands, npixs,
                  ensemble of unresolved sources, when a circular restoring beam
                  is appropriate.
 
-    :param max_pix_variance_factor: float to take account of the additional
-             uncertainty introduced by fudge_max_pix. Can probably be removed
-             since its effect is negligible in all sensible cases.
-
     :param beam: tuple of 3 floats describing the restoring beam in terms of
                  its semi-major and semi-minor axes (in pixels) and the position
                  angle of the semi-major axis, east from local north, in
@@ -1733,8 +1737,7 @@ def source_measurements_pixels_and_celestial_vectorised(num_islands, npixs,
     fitting.moments_enhanced(sources, noises, chunk_positions, xpositions,
                              ypositions, minimum_widths, npixs, thresholds,
                              local_noise_levels, maxposs, maxis,
-                             fudge_max_pix_factor, max_pix_variance_factor,
-                             np.array(beam), beamsize,
+                             fudge_max_pix_factor, np.array(beam), beamsize,
                              np.array(correlation_lengths), 0, 0,
                              Gaussian_islands, Gaussian_residuals, dummy,
                              moments_of_sources, sig, chisq, reduced_chisq)
diff --git a/sourcefinder/fitting.py b/sourcefinder/fitting.py
index 75ea2b3..35a80bb 100644
--- a/sourcefinder/fitting.py
+++ b/sourcefinder/fitting.py
@@ -26,7 +26,7 @@ def find_true_peak(peak, T, epsilon, msq, maxpix):
                                                       (peak / T - 1))
 
 
-def moments(data, fudge_max_pix_factor, beamsize, threshold=0):
+def moments(data, fudge_max_pix_factor, beam, beamsize, threshold=0):
     """
     Calculate source positional values using moments.
 
@@ -38,11 +38,19 @@ def moments(data, fudge_max_pix_factor, beamsize, threshold=0):
     ----------
     data : np.ma.MaskedArray or np.ndarray
         Actual 2D image data.
+
     fudge_max_pix_factor : float
         Correct for the underestimation of the peak by taking the maximum
         pixel value.
+
+    beam : 3-tuple
+        tuple of 3 floats: (semimajor axis, semiminor axis, theta).
+        The axes are in units of pixels and theta, the position angle of the
+        major axis wrt the positive y-axis, is in radians.
+
     beamsize : float
         The FWHM size of the clean beam.
+
     threshold : float, default: 0
         Source parameters like the semimajor and semiminor axes derived
         from moments can be underestimated if one does not take account of
@@ -74,119 +82,116 @@ def moments(data, fudge_max_pix_factor, beamsize, threshold=0):
             or np.isnan(working2)):
         raise ValueError("Unable to estimate Gauss shape")
 
-    maxpos = np.unravel_index(np.abs(data).argmax(), data.shape)
+    maxpos = np.unravel_index(data.argmax(), data.shape)
 
     # Are we fitting a -ve or +ve Gaussian?
-    if data.mean() >= 0:
+    if data[maxpos] > 0:
         peak = data[maxpos]
-        # The peak is always underestimated when you take the highest or lowest
-        # - for images other than Stokes I, where that may apply - pixel.
-        peak *= fudge_max_pix_factor
     else:
         peak = data.min()
+    # The peak is always underestimated when you take the highest or lowest
+    # - for images other than Stokes I, where that may apply - pixel.
+    peak *= fudge_max_pix_factor
 
     # Some problems arise with the sqrt of (working1-working2) when they are
     # equal, this happens with islands that have a thickness of only one pixel
     # in at least one dimension.  Due to rounding errors this difference
     # becomes negative--->math domain error in sqrt.
-    if len(data.nonzero()[0]) > 1:
-        semimajor_tmp = (working1 + working2) * 2.0 * math.log(2.0)
-        semiminor_tmp = (working1 - working2) * 2.0 * math.log(2.0)
-
+    semimajor_tmp = (working1 + working2) * 2.0 * math.log(2.0)
+    semiminor_tmp = (working1 - working2) * 2.0 * math.log(2.0)
+
+    if semiminor_tmp <= 0:
+        return {
+               "peak": peak,
+               "flux": total,
+               "xbar": xbar,
+               "ybar": ybar,
+               "semimajor": beam[0],
+               "semiminor": beam[1],
+               "theta": beam[2]
+        }
+    if abs(semimajor_tmp - semiminor_tmp) < 0.01:
+        # In this case the island is very thin in at least one dimension.
+        # We set the position angle to zero.
+        theta = 0.0
+    else:
         # Theta is not affected by the cut-off at the threshold, see Spreeuw's
         # thesis (2010), page 45, so we can easily compute the position angle
         # first.
-        if abs(semimajor_tmp - semiminor_tmp) < 0.01:
-            # short circuit!
-            theta = 0.
+        if xxbar != yybar:
+            theta = math.atan(2. * xybar / (xxbar - yybar)) / 2.
         else:
-            if xxbar != yybar:
-                theta = math.atan(2. * xybar / (xxbar - yybar)) / 2.
-            else:
-                theta = np.sign(xybar) * math.pi / 4.0
+            theta = np.sign(xybar) * math.pi / 4.0
 
-            if theta * xybar > 0.:
-                if theta < 0.:
-                    theta += math.pi / 2.0
-                else:
-                    theta -= math.pi / 2.0
-
-        if semiminor_tmp > 0:
-            rounded_barycenter = int(round(xbar)), int(round(ybar))
-            # basevalue and basepos will be needed for "tweaked moments".
-            try:
-                if not data.mask[rounded_barycenter]:
-                    basepos = rounded_barycenter
-                else:
-                    basepos = maxpos
-            except IndexError:
-                basepos = maxpos
-            except AttributeError:
-                # If the island is not masked at all, we can safely set basepos to
-                # the rounded barycenter position.
-                basepos = rounded_barycenter
-            basevalue = data[basepos]
-
-            if np.sign(threshold) == np.sign(basevalue):
-                # Implementation of "tweaked moments", equation 2.67 from
-                # Spreeuw's thesis. In that formula the "base position" was the
-                # maximum pixel position, though. Here that is the rounded
-                # barycenter position, unless it's masked. If it's masked, it
-                # will be the maximum pixel position.
-                deltax, deltay = xbar - basepos[0], ybar - basepos[1]
-
-                epsilon = np.log(2.) * ((np.cos(theta) * deltax +
-                                         np.sin(theta) * deltay) ** 2
-                                        + (np.cos(theta) * deltay -
-                                           np.sin(theta) * deltax) ** 2
-                                        * semiminor_tmp / semimajor_tmp)
-
-                # Set limits for the root finder similar to the bounds for
-                # Gaussian fits.
-                if basevalue > 0:
-                    low_bound = 0.5 * basevalue
-                    upp_bound = 1.5 * basevalue
-                else:
-                    low_bound = 1.5 * basevalue
-                    upp_bound = 0.5 * basevalue
-
-                # The number of iterations used for the root finder is also
-                # returned, but not used here.
-                peak, _ = newton_raphson_root_finder(find_true_peak, basevalue,
-                                                     low_bound, upp_bound,
-                                                     1e-8, 100,
-                                                     threshold, epsilon,
-                                                     semiminor_tmp, basevalue)
-                # The corrections below for the semi-major and semi-minor axes are
-                # to compensate for the underestimate of these quantities
-                # due to the cutoff at the threshold.
-                ratio = threshold / peak
-                semimajor_tmp /= (1.0 + math.log(ratio) * ratio / (1.0 - ratio))
-                semiminor_tmp /= (1.0 + math.log(ratio) * ratio / (1.0 - ratio))
+        if theta * xybar > 0.:
+            if theta < 0.:
+                theta += math.pi / 2.0
+            else:
+                theta -= math.pi / 2.0
 
-        elif np.sign(threshold) == np.sign(peak):
-            # A semi-minor axis exactly zero gives all kinds of problems.
-            # For instance wrt conversion to celestial coordinates.
-            # This is a quick fix.
+    rounded_barycenter = int(round(xbar)), int(round(ybar))
+    # basevalue and basepos will be needed for "tweaked moments".
+    try:
+        if not data.mask[rounded_barycenter]:
+            basepos = rounded_barycenter
+        else:
+            basepos = maxpos
+    except IndexError:
+        basepos = maxpos
+    except AttributeError:
+        # If the island is not masked at all, we can safely set basepos to
+        # the rounded barycenter position.
+        basepos = rounded_barycenter
+    basevalue = data[basepos]
+
+    if np.sign(threshold) == np.sign(basevalue):
+        # Implementation of "tweaked moments", equation 2.67 from
+        # Spreeuw's thesis. In that formula the "base position" was the
+        # maximum pixel position, though. Here that is the rounded
+        # barycenter position, unless it's masked. If it's masked, it
+        # will be the maximum pixel position.
+        deltax, deltay = xbar - basepos[0], ybar - basepos[1]
+
+        epsilon = np.log(2.) * ((np.cos(theta) * deltax +
+                                 np.sin(theta) * deltay) ** 2
+                                + (np.cos(theta) * deltay -
+                                   np.sin(theta) * deltax) ** 2
+                                * semiminor_tmp / semimajor_tmp)
+
+        # Set limits for the root finder similar to the bounds for
+        # Gaussian fits.
+        if basevalue > 0:
+            low_bound = 0.5 * basevalue
+            upp_bound = 1.5 * basevalue
+        else:
+            low_bound = 1.5 * basevalue
+            upp_bound = 0.5 * basevalue
+
+        # The number of iterations used for the root finder is also
+        # returned, but not used here.
+        peak, _ = newton_raphson_root_finder(find_true_peak, basevalue,
+                                             low_bound, upp_bound,
+                                             1e-8, 100,
+                                             threshold, epsilon,
+                                             semiminor_tmp, basevalue)
+
+    if np.sign(threshold) == np.sign(peak):
+        # The corrections below for the semi-major and semi-minor axes are
+        # to compensate for the underestimate of these quantities
+        # due to the cutoff at the threshold.
+        if peak > threshold:
             ratio = threshold / peak
-            semimajor_tmp /= (1.0 + np.log(ratio) * ratio / (1.0 - ratio))
-            semiminor_tmp = beamsize ** 2 / (np.pi ** 2 * semimajor_tmp)
-
-        semimajor = math.sqrt(semimajor_tmp)
-        semiminor = math.sqrt(semiminor_tmp)
+            semimajor_tmp /= (1.0 + math.log(ratio) * ratio / (1.0 - ratio))
+            semiminor_tmp /= (1.0 + math.log(ratio) * ratio / (1.0 - ratio))
 
-    else:
-        # This is the case when the island (or more likely subisland) has
-        # a size of only one pixel.
-        theta = 0.
-        semiminor = np.sqrt(beamsize / np.pi)
-        semimajor = np.sqrt(beamsize / np.pi)
+    semimajor = math.sqrt(semimajor_tmp)
+    semiminor = math.sqrt(semiminor_tmp)
 
     # NB: a dict should give us a bit more flexibility about arguments;
     # however, all those here are ***REQUIRED***.
     return {
         "peak": peak,
-        "flux": total,
+        "flux": np.pi * peak * semimajor * semiminor / beamsize,
         "xbar": xbar,
         "ybar": ybar,
         "semimajor": semimajor,
@@ -196,20 +201,19 @@ def moments(data, fudge_max_pix_factor, beamsize, threshold=0):
 
 
 @guvectorize([(float32[:], float32[:], int32[:], int32[:], int32[:], int32,
-               int32, float32, float32, int32[:], float32, float64, float64,
+               int32, float32, float32, int32[:], float32, float64,
                float64[:], float64, float64[:], float64, float64, float32[:, :],
                float32[:, :], float32[:, :], float32[:, :], float32[:],
                float32[:], float32[:])],
-             ('(n), (n), (m), (n), (n), (), (), (), (), (m), (), (), (), (k),' +
+             ('(n), (n), (m), (n), (n), (), (), (), (), (m), (), (), (k),' +
               '(), (m), (), (), (q, r), (q, r), (l, p) -> (l, p), (), (), ()'),
-             nopython=True)
+             nopython=True, target='parallel')
 def moments_enhanced(source_island, noise_island, chunkpos, posx, posy,
                      min_width, no_pixels, threshold, noise, maxpos, maxi,
-                     fudge_max_pix_factor, max_pix_variance_factor, beam,
-                     beamsize, correlation_lengths, clean_bias_error,
-                     frac_flux_cal_error, Gaussian_islands_map,
-                     Gaussian_residuals_map, dummy, computed_moments,
-                     significance, chisq, reduced_chisq):
+                     fudge_max_pix_factor, beam, beamsize, correlation_lengths,
+                     clean_bias_error, frac_flux_cal_error,
+                     Gaussian_islands_map, Gaussian_residuals_map, dummy,
+                     computed_moments, significance, chisq, reduced_chisq):
     """
     Calculate source properties using moments.
 
@@ -280,13 +284,10 @@ def moments_enhanced(source_island, noise_island, chunkpos, posx, posy,
         Correction factor for underestimation of the peak by considering the 
         maximum pixel value.
 
-    max_pix_variance_factor : float
-        Additional variance induced by the maximum pixel method, beyond the 
-        background noise.
-
     beam : np.ndarray
-        Array of three floats: [semimajor axis, semiminor axis, theta]. 
-        Units: pixels.
+        Array of three floats: (semimajor axis, semiminor axis, theta).
+        The axes are in units of pixels and theta, the position angle of the
+        major axis wrt the positive y-axis, is in radians.
 
     beamsize : float
         FWHM size of the clean beam. Units: pixels.
@@ -396,29 +397,34 @@ def moments_enhanced(source_island, noise_island, chunkpos, posx, posy,
         raise ValueError("Unable to estimate Gauss shape")
 
     # Are we fitting a -ve or +ve Gaussian?
-    if source_island.mean() >= 0:
+    if maxi > 0:
         # The peak is always underestimated when you take the highest pixel.
-        peak = maxi * fudge_max_pix_factor
+        peak = maxi
     else:
         peak = source_island.min()
+    # The peak is always underestimated when you take the highest or lowest
+    # - for images other than Stokes I, where that may apply - pixel.
+    peak *=  fudge_max_pix_factor
 
-    # We need this width to determine Gaussian shape parameters in a meaningful
-    # way.
-    if min_width > 2:
-        # The idea of the try except here is that, even though we require a
-        # minimum width of the island, there may still be occasions where
-        # working2 can be slightly higher than working1, perhaps due to rounding
-        # errors.
-        semimajor_tmp = (working1 + working2) * 2.0 * math.log(2.0)
-        semiminor_tmp = (working1 - working2) * 2.0 * math.log(2.0)
+    semimajor_tmp = (working1 + working2) * 2.0 * math.log(2.0)
+    semiminor_tmp = (working1 - working2) * 2.0 * math.log(2.0)
 
-        # Theta is not affected by the cut-off at the threshold, see Spreeuw's
-        # thesis (2010), page 45, so we can easily compute the position angle
-        # first.
+    if semiminor_tmp <= 0 or min_width <= 2:
+        # Short circuit in one way or another!
+        # We need a minimal width to determine Gaussian shape parameters in
+        # a meaningful way.
+        smaj = beam[0]
+        smin = beam[1]
+        theta = beam[2]
+    else:
         if abs(semimajor_tmp - semiminor_tmp) < 0.01:
-            # short circuit!
-            theta = 0.
+            # In this case the island is very thin in at least one dimension.
+            # We set the position angle to zero.
+            theta = 0.0
         else:
+            # Theta is not affected by the cut-off at the threshold, see Spreeuw's
+            # thesis (2010), page 45, so we can easily compute the position angle
+            # first.
             if xxbar != yybar:
                 theta = math.atan(2. * xybar / (xxbar - yybar)) / 2.
             else:
@@ -430,99 +436,72 @@ def moments_enhanced(source_island, noise_island, chunkpos, posx, posy,
                 else:
                     theta -= math.pi / 2.0
 
-        if semiminor_tmp > 0:
-            # We will be extrapolating from a pixel centred position to the real
-            # position of the Gaussian peak.
-            rounded_barycenter = int(round(xbar)), int(round(ybar))
-            # First we try the pixel centred position closest to the barycenter
-            # position to extrapolate from.
-            if rounded_barycenter[0] in posx and rounded_barycenter[1] in posy:
-                basepos = rounded_barycenter
-                # We need to loop over all possible positions to find the
-                # source_island index corresponding to the rounded barycenter
-                # position.
-                for i in (posx==rounded_barycenter[0]).nonzero()[0]:
-                    if posy[i] == rounded_barycenter[1]:
-                        basevalue = source_island[i]
-                        break
+        # We will be extrapolating from a pixel centred position to the real
+        # position of the Gaussian peak.
+        rounded_barycenter = int(round(xbar)), int(round(ybar))
+        # First we try the pixel centred position closest to the barycenter
+        # position to extrapolate from.
+        if rounded_barycenter[0] in posx and rounded_barycenter[1] in posy:
+            basepos = rounded_barycenter
+            # We need to loop over all possible positions to find the
+            # source_island index corresponding to the rounded barycenter
+            # position.
+            for i in (posx==rounded_barycenter[0]).nonzero()[0]:
+                if posy[i] == rounded_barycenter[1]:
+                    basevalue = source_island[i]
+                    break
+        else:
+            # The rounded barycenter position is not in source_island, so we
+            # revert to the maximum pixel position, which is always included.
+            basepos = maxpos[0] - chunkpos[0], maxpos[1] - chunkpos[1]
+            # In this case we do not need to figure out the source_island
+            # index corresponding to the maximum pixel position, because
+            # maxi, the maximum pixel value has already been provided, as an
+            # argument to this function.
+            basevalue = maxi
+
+        if np.sign(threshold) == np.sign(basevalue):
+            # Implementation of "tweaked moments", equation 2.67 from
+            # Spreeuw's thesis. In that formula the "base position" was the
+            # maximum pixel position, though. Here that is the rounded
+            # barycenter position, unless it's masked. If it's masked, it
+            # will be the maximum pixel position.
+            deltax, deltay = xbar - basepos[0], ybar - basepos[1]
+
+            epsilon = np.log(2.) * ((np.cos(theta) * deltax +
+                                     np.sin(theta) * deltay) ** 2
+                                    + (np.cos(theta) * deltay -
+                                       np.sin(theta) * deltax) ** 2
+                                    * semiminor_tmp / semimajor_tmp)
+
+            # Set limits for the root finder similar to the bounds for
+            # Gaussian fits.
+            if basevalue > 0:
+                low_bound = 0.5 * basevalue
+                upp_bound = 1.5 * basevalue
             else:
-                # The rounded barycenter position is not in source_island, so we
-                # revert to the maximum pixel position, which is always included.
-                basepos = maxpos[0] - chunkpos[0], maxpos[1] - chunkpos[1]
-                # In this case we do not need to figure out the source_island
-                # index corresponding to the maximum pixel position, because
-                # maxi, the maximum pixel value has already been provided, as an
-                # argument to this function.
-                basevalue = maxi
-
-            if np.sign(threshold) == np.sign(basevalue):
-                # Implementation of "tweaked moments", equation 2.67 from
-                # Spreeuw's thesis. In that formula the "base position" was the
-                # maximum pixel position, though. Here that is the rounded
-                # barycenter position, unless it's masked. If it's masked, it
-                # will be the maximum pixel position.
-                deltax, deltay = xbar - basepos[0], ybar - basepos[1]
-
-                epsilon = np.log(2.) * ((np.cos(theta) * deltax +
-                                         np.sin(theta) * deltay) ** 2
-                                        + (np.cos(theta) * deltay -
-                                           np.sin(theta) * deltax) ** 2
-                                        * semiminor_tmp / semimajor_tmp)
-
-                # Set limits for the root finder similar to the bounds for
-                # Gaussian fits.
-                if basevalue > 0:
-                    low_bound = 0.5 * basevalue
-                    upp_bound = 1.5 * basevalue
-                else:
-                    low_bound = 1.5 * basevalue
-                    upp_bound = 0.5 * basevalue
-
-                # The number of iterations used for the root finder is also
-                # returned, but not used here.
-                peak, _ = newton_raphson_root_finder(find_true_peak, basevalue,
-                                                     low_bound, upp_bound,
-                                                     1e-8, 100,
-                                                     threshold, epsilon,
-                                                     semiminor_tmp, basevalue)
-                # The corrections below for the semi-major and semi-minor axes are
-                # to compensate for the underestimate of these quantities
-                # due to the cutoff at the threshold.
+                low_bound = 1.5 * basevalue
+                upp_bound = 0.5 * basevalue
+
+            # The number of iterations used for the root finder is also
+            # returned, but not used here.
+            peak, _ = newton_raphson_root_finder(find_true_peak, basevalue,
+                                                 low_bound, upp_bound,
+                                                 1e-8, 100,
+                                                 threshold, epsilon,
+                                                 semiminor_tmp, basevalue)
+
+        if np.sign(threshold) == np.sign(peak):
+            # The corrections below for the semi-major and semi-minor axes are
+            # to compensate for the underestimate of these quantities
+            # due to the cutoff at the threshold.
+            if peak > threshold:
                 ratio = threshold / peak
                 semimajor_tmp /= (1.0 + math.log(ratio) * ratio / (1.0 - ratio))
                 semiminor_tmp /= (1.0 + math.log(ratio) * ratio / (1.0 - ratio))
 
-        elif np.sign(threshold) == np.sign(peak):
-            # A semi-minor axis exactly zero gives all kinds of problems.
-            # For instance wrt conversion to celestial coordinates.
-            # This is a quick fix.
-            ratio = threshold / peak
-            semimajor_tmp /= (1.0 + np.log(ratio) * ratio / (1.0 - ratio))
-            semiminor_tmp = beamsize ** 2 / (np.pi ** 2 * semimajor_tmp)
-
-        try:
-            smaj = math.sqrt(semimajor_tmp)
-            smin = math.sqrt(semiminor_tmp)
-            if smin == 0:
-                # A semi-minor axis exactly zero gives all kinds of problems.
-                # For instance wrt conversion to celestial coordinates.
-                # This is a quick fix.
-                smin = beamsize / (np.pi * smaj)
-
-        # In all cases where we hit a math domain error, we can use the clean
-        # beam parameters to derive reasonable estimates for the Gaussian shape
-        # parameters.
-        except Exception:
-            smaj = beam[0]
-            smin = beam[1]
-            theta = beam[2]
-    else:
-        # In the case that the island has insufficient width, we can also derive
-        # reasonable estimates for Gaussian shape parameters using the clean
-        # beam.
-        smaj = beam[0]
-        smin = beam[1]
-        theta = beam[2]
+        smaj = math.sqrt(semimajor_tmp)
+        smin = math.sqrt(semiminor_tmp)
 
     # Reconstruct the Gaussian profile we just derived at the pixel positions
     # of the island. Initialise a flat array to hold these values.
@@ -626,8 +605,7 @@ def moments_enhanced(source_island, noise_island, chunkpos, posx, posy,
     #                utils.maximum_pixel_method_variance(
     #                    beam[0], beam[1], beam[2]) * peak ** 2)
     errorpeaksq = ((frac_flux_cal_error * peak) ** 2 +
-                   clean_bias_error ** 2 + noise ** 2 +
-                   max_pix_variance_factor * peak ** 2)
+                   clean_bias_error ** 2 + noise ** 2)
     errorpeak = np.sqrt(errorpeaksq)
 
     help1 = (errorsmaj / smaj) ** 2
@@ -764,7 +742,7 @@ def fitgaussian(pixels, params, fixed=None, max_nfev=None, bounds={}):
 
     Parameters
     ----------
-    pixels : np.ma.MaskedArray
+    pixels : np.ma.MaskedArray or np.ndarray
         Pixel values (with bad pixels masked)
     params : dict
         Initial fit parameters (possibly estimated using the moments() function)
diff --git a/sourcefinder/image.py b/sourcefinder/image.py
index e9faf18..62db4a4 100644
--- a/sourcefinder/image.py
+++ b/sourcefinder/image.py
@@ -19,10 +19,7 @@ from multiprocessing import Pool
 from functools import cached_property
 from functools import partial
 
-try:
-    import ndimage
-except ImportError:
-    from scipy import ndimage
+from scipy import ndimage
 from numba import guvectorize, float32, int32
 import os
 
@@ -102,8 +99,6 @@ class ImageData(object):
         # These three quantities are only dependent on the beam, so should be calculated
         # once the beam is known and not for each source separately.
         self.fudge_max_pix_factor = utils.fudge_max_pix(beam[0], beam[1], beam[2])
-        self.max_pix_variance_factor = utils.maximum_pixel_method_variance(
-                           beam[0], beam[1], beam[2])
         self.beamsize = utils.calculate_beamsize(beam[0], beam[1])
         self.correlation_lengths = utils.calculate_correlation_lengths(beam[0], beam[1])
         self.clip = {}
@@ -748,7 +743,6 @@ class ImageData(object):
                 fitme, threshold_at_pixel, self.rmsmap[chunk],
                 self.rmsmap[int(x), int(y)],
                 self.beam, self.fudge_max_pix_factor,
-                self.max_pix_variance_factor,
                 self.beamsize, self.correlation_lengths, fixed=fixed)
         except ValueError:
             # Fit failed to converge
@@ -778,7 +772,7 @@ class ImageData(object):
 
         Parameters
         ----------
-        positions : tuple
+        positions : list of tuples
             List of (RA, Dec) tuples. Positions to be fit, in decimal degrees.
         boxsize : int
             Length of the square section of the image to use for the fit.
@@ -987,12 +981,12 @@ class ImageData(object):
                 slices)
 
     @staticmethod
-    def fit_islands(fudge_max_pix_factor, max_pix_variance_factor, beamsize,
-                    correlation_lengths, fixed, island):
+    def fit_islands(fudge_max_pix_factor, beamsize, correlation_lengths,
+                    fixed, island):
         """This function was created to enable the use of 'partial' such that
         we can parallellize source measurements"""
-        return island.fit(fudge_max_pix_factor, max_pix_variance_factor, beamsize,
-                          correlation_lengths, fixed=fixed)
+        return island.fit(fudge_max_pix_factor, beamsize, correlation_lengths,
+                          fixed=fixed)
 
     # “The sliced_to_indices function has been accelerated by a factor 2.7
     # using ChatGPT 4.0.  Its AI-output has been verified for correctness,
@@ -1208,7 +1202,6 @@ class ImageData(object):
             with Pool(psutil.cpu_count()) as p:
                 fit_islands_partial = partial(ImageData.fit_islands,
                                               self.fudge_max_pix_factor,
-                                              self.max_pix_variance_factor,
                                               self.beamsize,
                                               self.correlation_lengths, fixed)
                 fit_results = p.map(fit_islands_partial, island_list)
@@ -1253,8 +1246,8 @@ class ImageData(object):
                     num_islands, npixs, maxposs, maxis, self.data_bgsubbed.data,
                     self.rmsmap.data, analysisthresholdmap.data, indices,
                     labelled_data, labels, self.wcs, self.fudge_max_pix_factor,
-                    self.max_pix_variance_factor,  self.beam, self.beamsize,
-                    self.correlation_lengths, self.eps_ra, self.eps_dec)
+                    self.beam, self.beamsize, self.correlation_lengths,
+                    self.eps_ra, self.eps_dec)
 
             if self.islands:
                 self.Gaussian_islands = Gaussian_islands
diff --git a/sourcefinder/utils.py b/sourcefinder/utils.py
index 54e9271..cb4f471 100644
--- a/sourcefinder/utils.py
+++ b/sourcefinder/utils.py
@@ -194,47 +194,6 @@ def fudge_max_pix(semimajor, semiminor, theta):
     return correction
 
 
-def maximum_pixel_method_variance(semimajor, semiminor, theta):
-    """Estimate variance for peak flux at pixel position of maximum
-
-    When we use the maximum pixel method, with a correction
-    fudge_max_pix, there should be no bias, unless the peaks of the
-    Gaussians are not randomly distributed, but relatively close to
-    the centres of the pixels due to selection effects from detection
-    thresholds.
-
-    Disregarding the latter effect and noise, we can compute the
-    variance of the maximum pixel method by integrating (the true
-    flux-the average true flux)^2 = (the true flux-fudge_max_pix)^2
-    over the pixel area and dividing by the pixel area ( = 1).  This
-    is just equal to integral of the true flux^2 over the pixel area
-    - fudge_max_pix^2.
-    """
-
-    # scipy.integrate.dblquad: Computes a double integral
-    # from the scipy docs:
-    #   Return the double (definite) integral of f1(y,x) from x=a..b
-    #   and y=f2(x)..f3(x).
-
-    log20 = np.log(2.0)
-    cos_theta = np.cos(theta)
-    sin_theta = np.sin(theta)
-
-    def landscape(y, x):
-        return np.exp(2.0 * log20 * (
-                      math.pow(((cos_theta * x + sin_theta * y) / semiminor),
-                                  2) +
-                      math.pow(((cos_theta * y - sin_theta * x) / semimajor),
-                                  2)))
-
-    (result, abserr) = scipy.integrate.dblquad(landscape, -0.5, 0.5,
-                                               lambda ymin: -0.5,
-                                               lambda ymax: 0.5)
-    variance = result - math.pow(fudge_max_pix(semimajor, semiminor, theta), 2)
-
-    return variance
-
-
 def flatten(nested_list):
     """Flatten a nested list
 


Test patch:
diff --git a/test/test_errors.py b/test/test_errors.py
index 79d3249..b3ed53b 100644
--- a/test/test_errors.py
+++ b/test/test_errors.py
@@ -5,7 +5,7 @@ from sourcefinder.extract import Detection
 from sourcefinder.extract import ParamSet
 from sourcefinder.utility.coordinates import WCS
 from sourcefinder.utility.uncertain import Uncertain
-from sourcefinder.utils import maximum_pixel_method_variance, calculate_correlation_lengths
+from sourcefinder.utils import calculate_correlation_lengths
 
 
 class DummyImage(object):
@@ -37,7 +37,6 @@ class TestFluxErrors(unittest.TestCase):
         self.beam = (1.0, 1.0, 0.0)
         self.noise = 1.0  # RMS at position of source
         self.threshold = 3.0  # significance * rms at position of source
-        self.max_pix_variance_factor = maximum_pixel_method_variance(*self.beam)
         self.correlation_lengths = calculate_correlation_lengths(self.beam[0], self.beam[1])
 
     def test_positive_flux_condon(self):
@@ -46,8 +45,8 @@ class TestFluxErrors(unittest.TestCase):
         self.assertGreaterEqual(self.p['flux'].error, 0)
 
     def test_positive_flux_moments(self):
-        self.p._error_bars_from_moments(self.noise, self.max_pix_variance_factor,
-                                        self.correlation_lengths, self.threshold)
+        self.p._error_bars_from_moments(self.noise, self.correlation_lengths,
+                                        self.threshold)
         self.assertGreaterEqual(self.p['peak'].error, 0)
         self.assertGreaterEqual(self.p['flux'].error, 0)
 
@@ -65,8 +64,7 @@ class TestFluxErrors(unittest.TestCase):
         # impossible. Cancel any further processing of the image.
         self.p['peak'] *= -1
         self.assertRaises(ValueError, self.p._error_bars_from_moments,
-                          self.noise, self.max_pix_variance_factor,
-                          self.correlation_lengths, self.threshold)
+                          self.noise, self.correlation_lengths, self.threshold)
 
 
 class TestPositionErrors(unittest.TestCase):
diff --git a/test/test_gaussian.py b/test/test_gaussian.py
index a8c8d85..4fc0e26 100644
--- a/test/test_gaussian.py
+++ b/test/test_gaussian.py
@@ -8,7 +8,7 @@ import numpy as np
 from sourcefinder.extract import source_profile_and_errors
 from sourcefinder.fitting import moments, fitgaussian, FIT_PARAMS
 from sourcefinder.gaussian import gaussian
-from sourcefinder.utils import fudge_max_pix, calculate_beamsize, maximum_pixel_method_variance, \
+from sourcefinder.utils import fudge_max_pix, calculate_beamsize, \
     calculate_correlation_lengths
 
 # The units that are tested often require information about the resolution element:
@@ -25,15 +25,18 @@ class SimpleGaussTest(unittest.TestCase):
         self.height = 10
         self.x = 250
         self.y = 250
-        self.maj = 40
-        self.min = 20
+        self.smaj = 40
+        self.smin = 20
         self.theta = 0
         self.mygauss = np.ma.array(
-            gaussian(self.height, self.x, self.y, self.maj, self.min,
+            gaussian(self.height, self.x, self.y, self.smaj, self.smin,
                      self.theta)(Xin, Yin))
-        self.beamsize = calculate_beamsize(self.maj, self.min)
-        self.fudge_max_pix_factor = fudge_max_pix(self.maj, self.min, self.theta)
-        self.moments = moments(self.mygauss, self.fudge_max_pix_factor, self.beamsize, 0)
+        self.beam = (self.smaj, self.smin, self.theta)
+        self.beamsize = calculate_beamsize(self.smaj, self.smin)
+        self.fudge_max_pix_factor = fudge_max_pix(self.smaj, self.smin,
+                                                  self.theta)
+        self.moments = moments(self.mygauss, self.fudge_max_pix_factor,
+                               self.beam, self.beamsize, 0)
         self.fit = fitgaussian(self.mygauss, self.moments)
 
     def testHeight(self):
@@ -51,12 +54,12 @@ class SimpleGaussTest(unittest.TestCase):
         self.assertAlmostEqual(self.fit["ybar"], self.y)
 
     def testMomentSize(self):
-        self.assertAlmostEqual(self.moments["semimajor"], self.maj, 3)
-        self.assertAlmostEqual(self.moments["semiminor"], self.min, 3)
+        self.assertAlmostEqual(self.moments["semimajor"], self.smaj, 3)
+        self.assertAlmostEqual(self.moments["semiminor"], self.smin, 3)
 
     def testFitSize(self):
-        self.assertAlmostEqual(self.fit["semimajor"], self.maj)
-        self.assertAlmostEqual(self.fit["semiminor"], self.min)
+        self.assertAlmostEqual(self.fit["semimajor"], self.smaj)
+        self.assertAlmostEqual(self.fit["semiminor"], self.smin)
 
     def testMomentAngle(self):
         self.assertAlmostEqual(self.moments["theta"], self.theta)
@@ -73,15 +76,18 @@ class NegativeGaussTest(SimpleGaussTest):
         self.height = -10
         self.x = 250
         self.y = 250
-        self.maj = 40
-        self.min = 20
+        self.smaj = 40
+        self.smin = 20
         self.theta = 0
         self.mygauss = np.ma.array(
-            gaussian(self.height, self.x, self.y, self.maj, self.min,
+            gaussian(self.height, self.x, self.y, self.smaj, self.smin,
                      self.theta)(Xin, Yin))
-        self.beamsize = calculate_beamsize(self.maj, self.min)
-        self.fudge_max_pix_factor = fudge_max_pix(self.maj, self.min, self.theta)
-        self.moments = moments(self.mygauss, self.fudge_max_pix_factor, self.beamsize, 0)
+        self.beam = (self.smaj, self.smin, self.theta)
+        self.beamsize = calculate_beamsize(self.smaj, self.smin)
+        self.fudge_max_pix_factor = fudge_max_pix(self.smaj, self.smin,
+                                                  self.theta)
+        self.moments = moments(self.mygauss, self.fudge_max_pix_factor, self.beam,
+                               self.beamsize, 0)
         self.fit = fitgaussian(self.mygauss, self.moments)
 
     def testHeight(self):
@@ -96,15 +102,18 @@ class CircularGaussTest(SimpleGaussTest):
         self.height = 10
         self.x = 250
         self.y = 250
-        self.maj = 40
-        self.min = 40
+        self.smaj = 40
+        self.smin = 40
         self.theta = 0
         self.mygauss = np.ma.array(
-            gaussian(self.height, self.x, self.y, self.maj, self.min,
+            gaussian(self.height, self.x, self.y, self.smaj, self.smin,
                      self.theta)(Xin, Yin))
-        self.beamsize = calculate_beamsize(self.maj, self.min)
-        self.fudge_max_pix_factor = fudge_max_pix(self.maj, self.min, self.theta)
-        self.moments = moments(self.mygauss, self.fudge_max_pix_factor, self.beamsize, 0)
+        self.beam = (self.smaj, self.smin, self.theta)
+        self.beamsize = calculate_beamsize(self.smaj, self.smin)
+        self.fudge_max_pix_factor = fudge_max_pix(self.smaj, self.smin,
+                                                  self.theta)
+        self.moments = moments(self.mygauss, self.fudge_max_pix_factor,
+                               self.beam, self.beamsize, 0)
         self.fit = fitgaussian(self.mygauss, self.moments)
 
     def testMomentAngle(self):
@@ -122,15 +131,18 @@ class NarrowGaussTest(SimpleGaussTest):
         self.height = 10
         self.x = 250
         self.y = 250
-        self.maj = 40
-        self.min = 1
+        self.smaj = 40
+        self.smin = 1
         self.theta = 0
         self.mygauss = np.ma.array(gaussian(
-            self.height, self.x, self.y, self.maj, self.min, self.theta)(Xin,
+            self.height, self.x, self.y, self.smaj, self.smin, self.theta)(Xin,
                                                                          Yin))
-        self.beamsize = calculate_beamsize(self.maj, self.min)
-        self.fudge_max_pix_factor = fudge_max_pix(self.maj, self.min, self.theta)
-        self.moments = moments(self.mygauss, self.fudge_max_pix_factor, self.beamsize, 0)
+        self.beam = (self.smaj, self.smin, self.theta)
+        self.beamsize = calculate_beamsize(self.smaj, self.smin)
+        self.fudge_max_pix_factor = fudge_max_pix(self.smaj, self.smin,
+                                                  self.theta)
+        self.moments = moments(self.mygauss, self.fudge_max_pix_factor,
+                               self.beam, self.beamsize, 0)
         self.fit = fitgaussian(self.mygauss, self.moments)
 
 
@@ -142,15 +154,17 @@ class RotatedGaussTest(SimpleGaussTest):
         self.height = 10
         self.x = 250
         self.y = 250
-        self.maj = 40
-        self.min = 20
+        self.smaj = 40
+        self.smin = 20
         self.theta = np.pi / 4
         self.mygauss = np.ma.array(gaussian(
-            self.height, self.x, self.y, self.maj, self.min, self.theta)(Xin,
+            self.height, self.x, self.y, self.smaj, self.smin, self.theta)(Xin,
                                                                          Yin))
-        self.beamsize = calculate_beamsize(self.maj, self.min)
-        self.fudge_max_pix_factor = fudge_max_pix(self.maj, self.min, self.theta)
-        self.moments = moments(self.mygauss, self.fudge_max_pix_factor, self.beamsize, 0)
+        self.beam = (self.smaj, self.smin, self.theta)
+        self.beamsize = calculate_beamsize(self.smaj, self.smin)
+        self.fudge_max_pix_factor = fudge_max_pix(self.smaj, self.smin, self.theta)
+        self.moments = moments(self.mygauss, self.fudge_max_pix_factor,
+                               self.beam, self.beamsize, 0)
         self.fit = fitgaussian(self.mygauss, self.moments)
 
     def testMomentAngle(self):
@@ -165,15 +179,17 @@ class RotatedGaussTest2(SimpleGaussTest):
         self.height = 10
         self.x = 250
         self.y = 250
-        self.maj = 40
-        self.min = 20
+        self.smaj = 40
+        self.smin = 20
         self.theta = 3 * np.pi / 4
         self.mygauss = np.ma.array(gaussian(
-            self.height, self.x, self.y, self.maj, self.min, self.theta)(Xin,
+            self.height, self.x, self.y, self.smaj, self.smin, self.theta)(Xin,
                                                                          Yin))
-        self.beamsize = calculate_beamsize(self.maj, self.min)
-        self.fudge_max_pix_factor = fudge_max_pix(self.maj, self.min, self.theta)
-        self.moments = moments(self.mygauss, self.fudge_max_pix_factor, self.beamsize, 0)
+        self.beam = (self.smaj, self.smin, self.theta)
+        self.beamsize = calculate_beamsize(self.smaj, self.smin)
+        self.fudge_max_pix_factor = fudge_max_pix(self.smaj, self.smin, self.theta)
+        self.moments = moments(self.mygauss, self.fudge_max_pix_factor,
+                               self.beam, self.beamsize, 0)
         self.fit = fitgaussian(self.mygauss, self.moments)
 
     def testMomentAngle(self):
@@ -192,15 +208,18 @@ class AxesSwapGaussTest(SimpleGaussTest):
         self.height = 10
         self.x = 250
         self.y = 250
-        self.maj = 20
-        self.min = 40
+        self.smaj = 20
+        self.smin = 40
         self.theta = 0
         self.mygauss = np.ma.array(gaussian(
-            self.height, self.x, self.y, self.maj, self.min, self.theta)(Xin,
+            self.height, self.x, self.y, self.smaj, self.smin, self.theta)(Xin,
                                                                          Yin))
-        self.beamsize = calculate_beamsize(self.maj, self.min)
-        self.fudge_max_pix_factor = fudge_max_pix(self.maj, self.min, self.theta)
-        self.moments = moments(self.mygauss, self.fudge_max_pix_factor, self.beamsize, 0)
+        self.beam = (self.smaj, self.smin, self.theta)
+        self.beamsize = calculate_beamsize(self.smaj, self.smin)
+        self.fudge_max_pix_factor = fudge_max_pix(self.smaj, self.smin,
+                                                  self.theta)
+        self.moments = moments(self.mygauss, self.fudge_max_pix_factor,
+                               self.beam, self.beamsize, 0)
         self.fit = fitgaussian(self.mygauss, self.moments)
 
     def testMomentAngle(self):
@@ -222,12 +241,12 @@ class AxesSwapGaussTest(SimpleGaussTest):
         self.assertAlmostEqual(theta, np.pi / 2)
 
     def testMomentSize(self):
-        self.assertAlmostEqual(self.moments["semiminor"], self.maj, 5)
-        self.assertAlmostEqual(self.moments["semimajor"], self.min, 5)
+        self.assertAlmostEqual(self.moments["semiminor"], self.smaj, 5)
+        self.assertAlmostEqual(self.moments["semimajor"], self.smin, 5)
 
     def testFitSize(self):
-        self.assertAlmostEqual(self.fit["semiminor"], self.maj)
-        self.assertAlmostEqual(self.fit["semimajor"], self.min)
+        self.assertAlmostEqual(self.fit["semiminor"], self.smaj)
+        self.assertAlmostEqual(self.fit["semimajor"], self.smin)
 
 
 class RandomGaussTest(unittest.TestCase):
@@ -237,12 +256,14 @@ class RandomGaussTest(unittest.TestCase):
     def setUp(self):
         Xin, Yin = np.indices((500, 500))
         self.mygauss = np.random.random(Xin.shape)
+        self.beam = (beam[0], beam[1], beam[2])
         self.beamsize = calculate_beamsize(beam[0], beam[1])
         self.fudge_max_pix_factor = fudge_max_pix(beam[0], beam[1], beam[2])
 
     def testMoments(self):
         try:
-            moments(self.mygauss, self.fudge_max_pix_factor, self.beamsize, 0)
+            moments(self.mygauss, self.fudge_max_pix_factor, self.beam,
+                    self.beamsize, 0)
         except:
             self.fail('Moments method failed to run.')
 
@@ -261,9 +282,12 @@ class NoisyGaussTest(unittest.TestCase):
         self.mygauss = np.ma.array(
             gaussian(self.peak, self.xbar, self.ybar,
                      self.semimajor, self.semiminor, self.theta)(Xin, Yin))
+        self.beam = (self.semimajor, self.semiminor, self.theta)
         self.beamsize = calculate_beamsize(self.semimajor, self.semiminor)
-        self.fudge_max_pix_factor = fudge_max_pix(self.semimajor, self.semiminor, self.theta)
-        self.moments = moments(self.mygauss, self.fudge_max_pix_factor, self.beamsize, 0)
+        self.fudge_max_pix_factor = fudge_max_pix(self.semimajor, self.semiminor,
+                                                  self.theta)
+        self.moments = moments(self.mygauss, self.fudge_max_pix_factor,
+                               self.beam, self.beamsize, 0)
 
     def test_tiny_pixel_offset(self):
         pixel_noise = 0.0001
@@ -315,8 +339,6 @@ class NoisyGaussTest(unittest.TestCase):
             beam=(self.semimajor, self.semiminor, self.theta),
             fudge_max_pix_factor=fudge_max_pix(self.semimajor, self.semiminor,
                                                self.theta),
-            max_pix_variance_factor=maximum_pixel_method_variance(
-                self.semimajor, self.semiminor, self.theta),
             correlation_lengths=calculate_correlation_lengths(
                 self.semimajor, self.semiminor),
             beamsize=calculate_beamsize(self.semimajor, self.semiminor)
diff --git a/test/test_image.py b/test/test_image.py
index 4853a8b..b7aabb5 100644
--- a/test/test_image.py
+++ b/test/test_image.py
@@ -120,7 +120,7 @@ class TestFitFixedPositions(unittest.TestCase):
         fit. A cropped image helps to meet GH's disk quota."""
         sample_coord = [[61.42263448, 63.33334492]]
         results = self.cropped_image.fit_fixed_positions(sample_coord, 32.4687)
-        self.assertAlmostEqual(results[0].peak.value, -0.00364, delta=1e-5)
+        self.assertAlmostEqual(results[0].peak.value, -0.00053898)
 
     def testLowFitThreshold(self):
         """
@@ -302,7 +302,7 @@ class TestSimpleImageSourceFind(unittest.TestCase):
         known_result_fit = \
             [1.36896042e+02, 1.40221872e+01,   # RA (deg), DEC (deg)
              5.06084005e-04, 1.29061600e-03,  # Err, err
-             7.24671176e-04, 1.04806706e-04,  # Peak spectral brightness, err
+             7.42220029e-04, 1.07344735e-04,  # Peak spectral brightness, err
              6.03179622e-04, 1.62549622e-04,  # Flux density, err
              6.44646215e+00, 2.55194168e+01,
              # Significance level, beam semimajor-axis width (arcsec)
@@ -310,13 +310,13 @@ class TestSimpleImageSourceFind(unittest.TestCase):
              # Beam semiminor-axis width (arcsec), beam position angle (deg)
              ew_sys_err, ns_sys_err,
              4.97109604e+00, 1.00000000e+00,  # error_radius (arcsec), fit_type
-             6.03417635e-01, 6.67105734e-01]  # chisq, reduced chisq
+             1.47489354e-01, 1.63056225e-01]  # chisq, reduced chisq
 
         known_result_moments = \
             [1.3689603e+02, 1.4022377e+01,  # RA (deg), DEC (deg)
              5.5378844e-04, 1.1825778e-03,  # Err, err
-             7.3612988e-04, 1.1431403e-04,  # Peak spectral brightness, err
-             6.0276804e-04, 1.6508212e-04,  # Flux density, err
+             7.3612988e-04, 1.1412757e-04,  # Peak spectral brightness, err
+             6.0276804e-04, 1.6499560e-04,  # Flux density, err
              6.4464622e+00, 2.4559519e+01,
              # Significance level, beam semimajor-axis width (arcsec)
              1.1146187e+01, 1.7876042e+02,  # Beam semiminor-axis width (arcsec),
diff --git a/test/test_sf_utils.py b/test/test_sf_utils.py
index c19b8a0..202a604 100644
--- a/test/test_sf_utils.py
+++ b/test/test_sf_utils.py
@@ -8,7 +8,7 @@ from sourcefinder.utility.uncertain import Uncertain
 from sourcefinder.utils import circular_mask
 from sourcefinder.utils import generate_result_maps
 from sourcefinder.utils import generate_subthresholds
-from sourcefinder.utils import maximum_pixel_method_variance, fudge_max_pix
+from sourcefinder.utils import fudge_max_pix
 
 
 class TestCircularMask(unittest.TestCase):
@@ -53,27 +53,17 @@ class TestResultMaps(unittest.TestCase):
 
 
 class UtilsTest(unittest.TestCase):
-    # format: semimajor, semiminor, theta, correction, variance
+    # format: semimajor, semiminor, theta, correction
     def setUp(self):
         self.correct_data = (
-            (1.46255645752, 1.30175183614, 1.51619779602, 1.06380431709,
-             0.0017438970306615786),
-            (1.56219151815, 1.25084762573, 1.49080535398, 1.06328376441,
-             0.0017765054643135159),
-            (1.55231285095, 1.24513448079, 1.44660127410, 1.06398549354,
-             0.0018208980607872238)
+            (1.46255645752, 1.30175183614, 1.51619779602, 1.06380431709),
+            (1.56219151815, 1.25084762573, 1.49080535398, 1.06328376441),
+            (1.55231285095, 1.24513448079, 1.44660127410, 1.06398549354)
         )
 
-    def testMaximumPixelMethodVariance(self):
-        for (
-        semimajor, semiminor, theta, correction, variance) in self.correct_data:
-            self.assertEqual(
-                maximum_pixel_method_variance(semimajor, semiminor, theta),
-                variance)
-
     def testFudgeMaxPix(self):
         for (
-        semimajor, semiminor, theta, correction, variance) in self.correct_data:
+        semimajor, semiminor, theta, correction) in self.correct_data:
             self.assertAlmostEqual(fudge_max_pix(semimajor, semiminor, theta),
                                    correction)