Include feedback from review

* Docstrings and comments
* Fixed tests that behave a bit differently now due to the figure of merit calculation

src/libertem/analysis/fullmatch.py

@@ -3,8 +3,6 @@
 
 import libertem.analysis.gridmatching as grm
 from libertem.utils import make_polar
-# import (Match, PointSelection, CorrelationResult,
-# make_polar, NotFoundException, make_polar_vectors)
 
 
 def full_match(centers, zero=None, parameters={}):
@@ -42,7 +40,7 @@ def full_match(cls, correlation_result: grm.CorrelationResult, zero, cand=[], pa
         parameters
             Parameters for the matching.
             min_angle: Minimum angle between two vectors to be considered candidates
-            tolerance: Relative position tolerance for peaks to be considered matches
+            tolerance: Absolute position tolerance in px for peaks to be considered matches
             min_points: Minimum points to try clustering matching. Otherwise match directly
             min_match: Minimum matched clusters from clustering matching to be considered successful
             min_cluster_size_fraction: Tuning parameter for clustering matching. Larger values allow
@@ -143,7 +141,7 @@ def hdbscan_candidates(points, parameters):
     In the end we return the shortest matches.
 
     '''
-    cutoff = parameters["tolerance"] * parameters["max_delta"]
+    cutoff = parameters["tolerance"]
     clusterer = hdbscan.HDBSCAN(**make_hdbscan_config(points, parameters))
     vectors = grm.make_polar_vectors(points, parameters)
     clusterer.fit(vectors)

src/libertem/analysis/gridmatching.py

@@ -75,6 +75,12 @@ def __init__(self, correlation_result: CorrelationResult,
         super().__init__(correlation_result, selector)
         assert len(indices) == len(self)
 
+    def __str__(self):
+        result = "zero: %s\n"\
+            "a: %s\n"\
+            "b: %s\n"
+        return result % (str(self.zero), str(self.a), str(self.b))
+
     @classmethod
     def from_point_selection(cls, point_selection: PointSelection,
             zero, a, b, indices, selector=None, parameters={}):
@@ -332,33 +338,47 @@ def _do_match(cls, point_selection: PointSelection, zero, polar_vectors, paramet
     @classmethod
     def _find_best_vector_match(cls, point_selection: PointSelection, zero, candidates, parameters):
         '''
-        Return the match that matches most points
+        Return the match that matches with the best figure of merit
+
+        Good properties for vectors are
+        * Matching many points in the result
+        * Orthogonal
+        * Equal length
+        * Short
+
+        The function implements a heuristic to calculate a figure of merit that boosts
+        candidates for each of the criteria that they fulfill or nearly fulfill.
+
+        FIXME the figure of merit is currently determined heuristically based on discrete
+        thresholds. A smooth function would be more elegant.
+
+        FIXME improve this on more real-world examples; define test cases.
         '''
         def fom(d):
             m = d[1]
             na = np.linalg.norm(m.a)
             nb = np.linalg.norm(m.b)
+            # Matching many points is good
             res = len(m)
-            # Orthogonal
+            # Boost for orthogonality
             if np.abs(np.dot(m.a, m.b)) < 0.1 * na * nb:
                 res *= 5
             elif np.abs(np.dot(m.a, m.b)) < 0.3 * na * nb:
                 res *= 2
-            # nearly equal length
-            if (na - nb) < 0.1 * max(na, nb):
+            # Boost fo nearly equal length
+            if np.abs(na - nb) < 0.1 * max(na, nb):
                 res *= 5
-            if (na - nb) < 0.3 * max(na, nb):
+            if np.abs(na - nb) < 0.3 * max(na, nb):
                 res *= 2
+            # The division favors short vectors
             res /= na
             res /= nb
             return res
 
         match_matrix = cls._do_match(point_selection, zero, candidates, parameters)
         if match_matrix:
-            # we select the entry with highest number of matches
-            candidate_index, match = sorted(
-                match_matrix.items(), key=fom, reverse=True
-            )[0]
+            # we select the entry with highest figure of merit (fom)
+            candidate_index, match = max(match_matrix.items(), key=fom)
             return match
         else:
             raise NotFoundException

tests/test_gridmatching.py

@@ -232,32 +232,36 @@ def test_fullmatch_three_residual(zero, a, b):
     grid_2 = _fullgrid(zero, aa, bb, 4, skip_zero=True)
 
     random = np.array([
-        (0.3, 0.5),
-        (-0.3, 12.5),
-        (-0.3, -17.5),
+        (0.3, 0.4),
+        (-0.33, 12.5),
+        (-0.37, -17.4),
     ])
 
     grid = np.vstack((grid_1, grid_2, random))
 
     parameters = {
         'min_delta': 0.3,
         'max_delta': 3,
+        'tolerance': 0.05
     }
 
     (matches, unmatched, weak) = fm.full_match(grid, zero=zero, parameters=parameters)
 
+    print(matches[0])
+    print(matches[1])
+
     assert(len(matches) == 2)
 
     assert(len(unmatched) == len(random))
     assert(len(weak) == 0)
 
     match1 = matches[0]
 
-    assert(np.allclose(zero, match1.zero))
+    assert np.allclose(zero, match1.zero)
     assert(np.allclose(a, match1.a) or np.allclose(b, match1.a)
-           or np.allclose(-a, match1.a) or np.allclose(-b, match1.a))
+          or np.allclose(-a, match1.a) or np.allclose(-b, match1.a))
     assert(np.allclose(a, match1.b) or np.allclose(b, match1.b)
-           or np.allclose(-a, match1.b) or np.allclose(-b, match1.b))
+          or np.allclose(-a, match1.b) or np.allclose(-b, match1.b))
     assert(len(match1) == len(grid_1))
     assert(np.allclose(match1.calculated_refineds, grid_1))
 
@@ -297,6 +301,7 @@ def test_fullmatch_one_residual(zero, a, b):
     parameters = {
         'min_delta': 0.3,
         'max_delta': 3,
+        'tolerance': 0.1
     }
     (matches, unmatched, weak) = fm.full_match(grid, zero=zero, parameters=parameters)
 
@@ -320,9 +325,13 @@ def test_fullmatch_no_residual(zero, a, b):
     parameters = {
         'min_delta': 0.3,
         'max_delta': 3,
+        'tolerance': 0.1
     }
     (matches, unmatched, weak) = fm.full_match(grid, zero=zero, parameters=parameters)
 
+    print(matches[0])
+    print(matches[1])
+
     assert(len(matches) == 2)
 
     assert(len(unmatched) == 0)
@@ -347,6 +356,7 @@ def test_fullmatch_weak(zero, a, b):
     parameters = {
         'min_delta': 0.3,
         'max_delta': 3,
+        'tolerance': 0.05
     }
 
     values = np.ones(len(grid))