not yet ready, working at #5

.gitignore

@@ -16,3 +16,4 @@ i18n/latest_run
 .coverage
 test/gisquestion224812.py
 qgis-auth.db
+test/gisquestion224812-2.py

ghini_tree_position.py

@@ -581,6 +581,13 @@ def place_initial_three_points(points, distances, gps):
         v13 = R.dot([Cx, -Cy])
     # now rotate as of v12
     points[P3]['coordinates'] = tuple(v13 + gps[P1]['coordinates'])
+    #
+    # and compute connectivity to these three points for all others
+    for p in points.values():
+        p.setdefault('prio', 0)
+    for p in [P1, P2, P3]:
+        for q in distances[p]:
+            points[q]['prio'] += 1
 
 
 def rigid_transform_points(points, x, y, theta):
@@ -590,15 +597,13 @@ def rigid_transform_points(points, x, y, theta):
     import numpy as np
     from math import cos, sin, pi
     theta = theta / 180.0 * pi
-    R = np.array([[cos(theta), -sin(theta)],
-                  [sin(theta), cos(theta)]])
+    R = np.array([[cos(theta), -sin(theta), x],
+                  [sin(theta), cos(theta), y]])
     result = {}
     for k, p in points.items():
         result[k] = dict(p)
-        pt = np.array(p['coordinates'])
-        if theta != 0:
-            pt = R.dot(pt)
-        result[k]['coordinates'] = tuple(pt + (x, y))
+        pt = np.array(tuple(p['coordinates'][:2]) + (1, ))
+        result[k]['coordinates'] = tuple(R.dot(pt))
     return result
 
 
@@ -615,14 +620,18 @@ def distance_between_homonyms(p, q):
 
 def compute_minimal_distance_transformation(p, q):
     """computes the x, y, theta rigid transformation that minimizes the SSD
+
+    the rigid transformation applied to the frame 'p' that approximates it
+    to the points 'q'
+
     """
 
     def target(x):
         return distance_between_homonyms(rigid_transform_points(p, *x), q)
 
     import numpy as np
     import scipy.optimize
-    optres = scipy.optimize.minimize(target, (0, 0, 0))
+    optres = scipy.optimize.minimize(target, (0, 0, 0), method='Powell')
     return optres.x
 
 

test/test_initial_three_points.py

@@ -290,6 +290,15 @@ def test_rigid_transform_points_rotate_90(self):
         assert_almost_equal(p['p2']['coordinates'], (-5.0, 4.0))
         assert_almost_equal(p['p3']['coordinates'], (-5.0, 1.0))
 
+    def test_rigid_transform_points_rotate_45(self):
+        p = {'p1': {'coordinates': (1.0, 0.0)},
+             'p2': {'coordinates': (0.0, 1.0)},
+             'p3': {'coordinates': (1.0, 1.0)}, }
+        p = rigid_transform_points(p, x=0, y=0, theta=45)
+        assert_almost_equal(p['p1']['coordinates'], (.70710678, .707106780))
+        assert_almost_equal(p['p2']['coordinates'], (-.70710678, .70710678))
+        assert_almost_equal(p['p3']['coordinates'], (0, 1.41421356))
+
     def test_rigid_transform_points_rotate_m90(self):
         p = {'p1': {'coordinates': (1.0, 1.0)},
              'p2': {'coordinates': (4.0, 5.0)},

test/test_solve_puzzle.py

@@ -5,6 +5,9 @@
 from ghini_tree_position import solve_puzzle, get_distances_from_csv
 from ghini_tree_position import utm_zone_proj4
 from ghini_tree_position import extrapolate_coordinates
+from ghini_tree_position import compute_minimal_distance_transformation
+from ghini_tree_position import place_initial_three_points
+from ghini_tree_position import rigid_transform_points
 
 
 class TestUTMChoice(unittest.TestCase):
@@ -137,27 +140,83 @@ def test_extrapolate_coordinates_two(self):
         assert_almost_equal(tuple(points['t']['coordinates']), (2.0, 4.5))
 
 
+class TestComputeMinimalDistanceTransformation(unittest.TestCase):
+
+    def test_compute_minimal_distance_transformation_2_points(self):
+        p = {'A': {'coordinates': (0, 0)},
+             'B': {'coordinates': (2, 0)}}
+        q = {'A': {'coordinates': (0, 0)},
+             'B': {'coordinates': (4, 0)}}
+        t = compute_minimal_distance_transformation(p, q)
+        assert_almost_equal(t, (1, 0, 0), decimal=6)
+
+    def test_compute_minimal_distance_transformation_345(self):
+        p = {'A': {'coordinates': (0, 0)},
+             'B': {'coordinates': (4, 0)},
+             'C': {'coordinates': (0, 3)},
+             'D': {'coordinates': (4, 3)}}
+        q = {'A': {'coordinates': (0, 0)},
+             'B': {'coordinates': (8, 0)},
+             'C': {'coordinates': (0, 6)},
+             'D': {'coordinates': (8, 6)}}
+        t = compute_minimal_distance_transformation(p, q)
+        assert_almost_equal(t, (2, 1.5, 0), decimal=6)
+
+    def test_compute_minimal_distance_transformation_4_points(self):
+        p = {'A': {'coordinates': (2, 0)},
+             'B': {'coordinates': (1, 0)},
+             'C': {'coordinates': (0, 1)},
+             'D': {'coordinates': (0, 2)}}
+        q = {'A': {'coordinates': (0, 2)},
+             'B': {'coordinates': (0, 1)},
+             'C': {'coordinates': (-1, 0)},
+             'D': {'coordinates': (-2, 0)}}
+        t = compute_minimal_distance_transformation(p, q)
+        assert_almost_equal(t[:2], (0, 0), decimal=6)
+        assert_almost_equal(t[2] / 90, 1.0, decimal=5)
+
+    def test_compute_minimal_distance_transformation_4_points2(self):
+        p = {'A': {'coordinates': (2, 0)},
+             'B': {'coordinates': (1, 0)},
+             'C': {'coordinates': (0, 1)},
+             'D': {'coordinates': (0, 2)}}
+        q = {'A': {'coordinates': (0, 2.1)},
+             'B': {'coordinates': (0, 0.9)},
+             'C': {'coordinates': (-0.9, 0)},
+             'D': {'coordinates': (-2.1, 0)}}
+        t = compute_minimal_distance_transformation(p, q)
+        assert_almost_equal(t[:2], (0, 0), decimal=6)
+        assert_almost_equal(t[2] / 90, 1.0, decimal=5)
+
+
 class TestSolvePuzzle(unittest.TestCase):
 
     def test_solve_puzzle(self):
         points = {
-            'p1': {'id': 'p1'},
-            'p2': {'id': 'p2'},
-            'p3': {'id': 'p3'},
-            'p4': {'id': 'p4'}}
-        proj_gps = {'p1': {'coordinates': (0.0, 0.0)},
-                    'p2': {'coordinates': (4.0, 0.0)},
-                    'p3': {'coordinates': (4.0, 3.0)},
-                    'p4': {'coordinates': (0.0, 3.0)}, }
-        distances_list = [('p1', 'p2', 4.0),
-                          ('p1', 'p3', 5.0),
-                          ('p1', 'p4', 3.0),
-                          ('p2', 'p3', 3.0),
-                          ('p2', 'p4', 5.0),
-                          ('p1', 'p4', 4.0)]
+            'A': {'id': 'A'},
+            'B': {'id': 'B'},
+            'C': {'id': 'C'},
+            'D': {'id': 'D'}}
+        proj_gps = {'A': {'coordinates': (0.0, 0.0)},
+                    'B': {'coordinates': (8.0, 0.0)},
+                    'C': {'coordinates': (8.0, 6.0)},
+                    'D': {'coordinates': (0.0, 6.0)}, }
+        distances_list = [('A', 'B', 4.0),
+                          ('A', 'C', 5.0),
+                          ('A', 'D', 3.0),
+                          ('B', 'C', 3.0),
+                          ('B', 'D', 5.0),
+                          ('C', 'D', 4.0)]
         distances = {}
         for p, q, d in distances_list:
             distances.setdefault(p, {})
             distances.setdefault(q, {})
             distances[p][q] = distances[q][p] = d
-        solution = solve_puzzle(points, distances, proj_gps)
+        place_initial_three_points(points, distances, proj_gps)
+        extrapolate_coordinates(points, distances)
+        t = compute_minimal_distance_transformation(points, proj_gps)
+        p = rigid_transform_points(points, *t)
+        assert_almost_equal(p['A']['coordinates'], (2, 1.5), decimal=6)
+        assert_almost_equal(p['B']['coordinates'], (6, 1.5), decimal=6)
+        assert_almost_equal(p['C']['coordinates'], (6, 4.5), decimal=6)
+        assert_almost_equal(p['D']['coordinates'], (2, 4.5), decimal=6)