@@ -93,7 +93,11 @@ def covariance(self):
 
         Returns:
             numpy.ndarray: 2x2 matrix representing covariance matrix in the
-            reference x-y frame.
+            reference x-y frame. Covariance matrix elements are ordered
+            like::
+
+                [[sigma_x**2, rho*sigma_x*sigma_y],
+                 [rho*sigma_x*sigma_y, sigma_y**2]]
         """
         rotated_cov = np.array([[self.semimajor ** 2, 0],
                                 [0, self.semiminor ** 2]],
 
@@ -0,0 +1,74 @@
+import math
+
+import numpy as np
+from scipy.spatial.distance import mahalanobis as scipy_mahalanobis
+
+from fastimgproto.sourcefind.fit import Gaussian2dParams
+
+
+def mahalanobis_distance(gaussian2d_params, x, y):
+    """
+    Compute the Mahalanobis distance from Gaussian-centre to point (x,y).
+
+    Mahalanobis distance is the Euclidean distance normalized by the
+    standard-deviation of the Gaussian:
+    https://en.wikipedia.org/wiki/Mahalanobis_distance
+
+    Args:
+        gaussian2d_params (Gaussian2dParams): Gaussian to compute distance for.
+        x (float): X-coordinate of position
+        y (float): Y-coordinate of position
+
+    Returns:
+        float: Mahalanobis distance
+    """
+    inv_covariance = np.linalg.inv(gaussian2d_params.covariance)
+    g_centre = np.array(
+        [gaussian2d_params.x_centre, gaussian2d_params.y_centre],
+        dtype=np.float_
+    )
+    posn = np.array([x, y], dtype=np.float_)
+    return scipy_mahalanobis(u=g_centre, v=posn, VI=inv_covariance)
+
+
+def pixel_distance(gaussian2d_params, x, y):
+    """
+    Compute the distance from Gaussian-centre to point (x,y).
+
+    (Regular distance measured in units of pixels)
+
+    Args:
+        gaussian2d_params (Gaussian2dParams): Gaussian to compute distance for.
+        x (float): X-coordinate of position
+        y (float): Y-coordinate of position
+
+    Returns:
+        float: Euclidean distance (units of pixels)
+    """
+    xdiff = gaussian2d_params.x_centre - x
+    ydiff = gaussian2d_params.y_centre - y
+    return np.sqrt(xdiff * xdiff + ydiff * ydiff)
+
+
+def shape_difference(g1, g2):
+    """
+    A symmetric shape-difference metric for Gaussians
+
+    Computes the Covariance-difference term of the Bhattacharyya_distance, cf
+    https://en.wikipedia.org/wiki/Bhattacharyya_distance
+
+
+    Args:
+        g1 (Gaussian2dParams):
+        g2 (Gaussian2dParams):
+
+    Returns:
+        float: Dimensionless semipositive value. 0 if g1 and g2 have same shape.
+    """
+    cov1 = g1.covariance
+    cov2 = g2.covariance
+    halfway_cov = (cov1 + cov2) / 2.0
+    ratio = np.linalg.det(halfway_cov) / np.sqrt(
+        np.linalg.det(cov1) * np.linalg.det(cov2)
+    )
+    return 0.5 * math.log(ratio)
@@ -11,8 +11,9 @@
     gaussian_point_source,
     uncorrelated_gaussian_noise_background,
 )
-from fastimgproto.fixtures.sourcefind import \
-    check_single_source_extraction_successful
+from fastimgproto.fixtures.sourcefind import (
+    check_single_source_extraction_successful,
+)
 from fastimgproto.sourcefind.fit import Gaussian2dParams
 from fastimgproto.sourcefind.image import (
     Island,
 
@@ -10,6 +10,38 @@
 from fastimgproto.sourcefind.fit import Gaussian2dParams
 
 
+def test_covariance():
+    """
+    Verify / document the ordering of the the covariance matrix components
+    """
+    # A Gaussian elongated in the x-direction
+    smaj = 2
+    smin = 1.
+    theta = 0
+    x_long = Gaussian2dParams(x_centre=0,
+                              y_centre=0,
+                              amplitude=1,
+                              semimajor=smaj,
+                              semiminor=smin,
+                              theta=theta,
+                              )
+
+    print(x_long.covariance)
+    expected_covariance = np.array(
+        [[smaj ** 2, 0],
+         [0, smin ** 2]],
+        dtype=np.float_)
+    assert_allclose(x_long.covariance, expected_covariance, atol=1e-14)
+
+    # A Gaussian elongated in the y-direction (Note theta)
+    y_long = attr.evolve(x_long, theta=np.pi / 2)
+    expected_covariance = np.array(
+        [[smin ** 2, 0],
+         [0, smaj ** 2]],
+        dtype=np.float_)
+    assert_allclose(y_long.covariance, expected_covariance, atol=1e-14)
+
+
 def test_correlation_coefficient():
     smaj = 1.5
     smin = 1.
@@ -93,18 +125,18 @@ def test_validation_after_evolve_call():
 
 def test_unconstrained_initialization_theta_oob():
     pars = dict(x_centre=48.24, y_centre=52.66,
-                  amplitude=42.,
-                  semimajor=1.5,
-                  semiminor=1.4,
-                  theta=np.pi / 4,
-                  )
+                amplitude=42.,
+                semimajor=1.5,
+                semiminor=1.4,
+                theta=np.pi / 4,
+                )
     g1 = Gaussian2dParams(**pars)
 
     pars.update(theta=g1.theta + np.pi)
     g2 = Gaussian2dParams.from_unconstrained_parameters(**pars)
     assert g1 == g2
 
-    pars.update(theta=g1.theta + 2*np.pi)
+    pars.update(theta=g1.theta + 2 * np.pi)
     g3 = Gaussian2dParams.from_unconstrained_parameters(**pars)
     assert g1 == g3
 
@@ -114,19 +146,18 @@ def test_unconstrained_initialization_theta_oob():
 
     # Check initialization at bounds:
     # Include upper bound
-    pars.update(theta=np.pi/2.)
+    pars.update(theta=np.pi / 2.)
     g5 = Gaussian2dParams(**pars)
 
     # Shouldn't work - exclude lower bound
-    pars.update(theta=-np.pi/2.)
+    pars.update(theta=-np.pi / 2.)
     with pytest.raises(ValueError):
         g6 = Gaussian2dParams(**pars)
     # Should flip lower bound to upper bound:
     g6 = Gaussian2dParams.from_unconstrained_parameters(**pars)
     assert g6.theta == np.pi / 2.
 
 
-
 def test_unconstrained_initialization_flipped_major_minor():
     smaj = 1.5
     smin = 1.3
@@ -142,6 +173,6 @@ def test_unconstrained_initialization_flipped_major_minor():
         g1.x_centre, g1.y_centre, g1.amplitude,
         semimajor=g1.semiminor,
         semiminor=g1.semimajor,
-        theta=g1.theta + np.pi/2.
+        theta=g1.theta + np.pi / 2.
     )
     assert g1 == g2
@@ -0,0 +1,90 @@
+import math
+
+import attr
+import numpy as np
+from pytest import approx
+
+from fastimgproto.sourcefind.fit import Gaussian2dParams
+from fastimgproto.sourcefind.metrics import (
+    mahalanobis_distance,
+    pixel_distance,
+    shape_difference,
+)
+
+
+def test_distance_metrics():
+    x_0, y_0 = 1., 4.
+    g2d_x_long = Gaussian2dParams(x_centre=x_0,
+                                  y_centre=y_0,
+                                  amplitude=1,
+                                  semimajor=2.,
+                                  semiminor=1.,
+                                  theta=0.,
+                                  )
+
+    # A Gaussian elongated in the y-direction (Note theta)
+    g2d_y_long = attr.evolve(g2d_x_long, theta=np.pi / 2)
+
+    x, y = x_0 + 2, y_0
+    assert mahalanobis_distance(g2d_x_long, x=x, y=y) == approx(1.0)
+    assert mahalanobis_distance(g2d_y_long, x=x, y=y) == approx(2.0)
+    assert pixel_distance(g2d_x_long, x, y) == pixel_distance(g2d_y_long, x, y)
+    assert pixel_distance(g2d_x_long, x, y) == 2.
+
+    x, y = x_0, y_0 + 1
+    assert mahalanobis_distance(g2d_x_long, x=x, y=y) == approx(1.0)
+    assert mahalanobis_distance(g2d_y_long, x=x, y=y) == approx(0.5)
+    assert pixel_distance(g2d_x_long, x, y) == pixel_distance(g2d_y_long, x, y)
+    assert pixel_distance(g2d_x_long, x, y) == 1.
+
+    g2d_diag_long = attr.evolve(g2d_x_long, theta=np.pi / 4.)
+    x, y = x_0 + np.cos(np.pi / 4), y_0 + np.sin(np.pi / 4)
+    assert mahalanobis_distance(g2d_diag_long, x=x, y=y) == approx(0.5)
+    assert pixel_distance(g2d_diag_long, x, y) == approx(1.)
+
+
+def test_shape_metric():
+    g_symmetric_1 = Gaussian2dParams(x_centre=0.,
+                                     y_centre=0.,
+                                     amplitude=1,
+                                     semimajor=1.,
+                                     semiminor=1.,
+                                     theta=0.,
+                                     )
+
+    assert shape_difference(g_symmetric_1, g_symmetric_1) == 0
+
+    translated_sym = attr.evolve(g_symmetric_1, x_centre=123,
+                                 y_centre=456, )
+    assert shape_difference(g_symmetric_1, translated_sym) == 0
+
+    # These are all the same 'shape distance' from the symmetric circular
+    # Gaussian:
+    g_x_2 = attr.evolve(g_symmetric_1, semimajor=2., semiminor=1.)
+    g_y_2 = attr.evolve(g_x_2, theta=np.pi / 2.)
+    g_diag_2 = attr.evolve(g_x_2, theta=np.pi / 4.)
+    g_x_half = attr.evolve(g_symmetric_1, semiminor=0.5)
+
+    symmetric_to_elongated_by_2 = shape_difference(g_symmetric_1, g_x_2)
+    for g_pars in (g_x_2, g_y_2, g_diag_2, g_x_half):
+        assert shape_difference(g_symmetric_1,
+                                g_pars) == symmetric_to_elongated_by_2
+
+    # Twice the elongation factor -> Twice the distance metric
+    assert shape_difference(g_x_half, g_x_2) == 2 * symmetric_to_elongated_by_2
+
+    # Elongated in both directions -> Twice the distance metric
+    g_symmetric_2 = attr.evolve(g_symmetric_1, semimajor=2., semiminor=2.)
+    assert (shape_difference(g_symmetric_1, g_symmetric_2) ==
+            approx(2 * symmetric_to_elongated_by_2))
+
+    # Elongated and rotated by 90-degrees:
+    # (Equivalent to shrink in x, stretch in y)
+    # -> Twice the distance metric
+    assert (shape_difference(g_x_2, g_y_2) ==
+            approx(2 * symmetric_to_elongated_by_2))
+
+    # Rotation by 45-degrees:
+    assert (symmetric_to_elongated_by_2 <
+            shape_difference(g_diag_2, g_x_2) <
+            2 * symmetric_to_elongated_by_2)