Merge 5783c0e into 0664209

clifford/tools/g3c/object_fitting.py

@@ -0,0 +1,132 @@
+
+from . import *
+
+
+@numba.njit
+def val_fit_circle(point_list):
+    """
+    Performs Leo Dorsts circle fitting technique
+    """
+    # Loop over our points and construct the matrix
+    accumulator_matrix = np.zeros((32, 32))
+    for i in range(point_list.shape[0]):
+        # Get the point as a left gmt matrix
+        P_i_l = get_left_gmt_matrix(point_list[i,:])
+        # Multiply and add
+        accumulator_matrix += P_i_l @ mask0 @ P_i_l
+    accumulator_matrix = accumulator_matrix @ mask1
+    # Find the eigenvalues of the matrix
+    e_vals, e_vecs = np.linalg.eig(accumulator_matrix)
+    # Find the smallest and second smallest non negative eigenvalues
+    min_eval = np.inf
+    min_eval_index = -1
+    min_eval_index2 = -1
+    for i in range(len(e_vals)):
+        this_e_val = e_vals[i]
+        if this_e_val < min_eval and this_e_val > 0:
+            min_eval = this_e_val
+            min_eval_index2 = min_eval_index
+            min_eval_index = i
+    best_sphere = val_normalised(mask1@np.real(e_vecs[:, min_eval_index]))
+    second_best_sphere = val_normalised(mask1@np.real(e_vecs[:, min_eval_index2]))
+    best_circle = val_normalised(mask3@dual_func(omt_func(best_sphere,second_best_sphere)))
+    return best_circle
+
+
+def fit_circle(point_list):
+    """
+    Performs Leo Dorsts circle fitting technique
+    """
+    return layout.MultiVector(value=val_fit_circle(np.array(point_list)))
+
+
+@numba.njit
+def val_fit_line(point_list):
+    """
+    Does line fitting with combo J.Lasenbys method and L. Dorsts
+    """
+    accumulator_matrix = np.zeros((32,32))
+    for i in range(point_list.shape[0]):
+        P_i_l = get_left_gmt_matrix(point_list[i,:])
+        P_i_r = get_right_gmt_matrix(point_list[i,:])
+        accumulator_matrix += mask3@P_i_l@P_i_r
+    # Find the eigenvalues of the matrix
+    e_vals, e_vecs = np.linalg.eig(accumulator_matrix)
+    # Find the smallest non negative eigenvalue
+    min_eval = np.inf
+    min_eval_index = -1
+    for i in range(len(e_vals)):
+        if e_vals[i] < min_eval and e_vals[i] > 0:
+            min_eval = e_vals[i]
+            min_eval_index = i
+    best_line = mask3@omt_func(dual_func(e_vecs[:, min_eval_index]),ninf_val)
+    return val_normalised(best_line)
+
+
+def fit_line(point_list):
+    """
+    Does line fitting with combo J.Lasenbys method and L. Dorsts
+    """
+    return layout.MultiVector(value=val_fit_line(np.array(point_list)))
+
+
+@numba.njit
+def val_fit_sphere(point_list):
+    """
+    Performs Leo Dorsts sphere fitting technique
+    """
+    # Loop over our points and construct the matrix
+    accumulator_matrix = np.zeros((32, 32))
+    for i in range(point_list.shape[0]):
+        # Get the point as a left gmt matrix
+        P_i_l = get_left_gmt_matrix(point_list[i,:])
+        # Multiply and add
+        accumulator_matrix += P_i_l @ mask0 @ P_i_l
+    accumulator_matrix = accumulator_matrix @ mask1
+    # Find the eigenvalues of the matrix
+    e_vals, e_vecs = np.linalg.eig(accumulator_matrix)
+    # Find the smallest non negative eigenvalues
+    min_eval = np.inf
+    min_eval_index = -1
+    for i in range(len(e_vals)):
+        if e_vals[i] < min_eval and e_vals[i] > 0:
+            min_eval = e_vals[i]
+            min_eval_index = i
+    best_sphere = val_normalised(mask4@dual_func(np.real(e_vecs[:, min_eval_index])))
+    return best_sphere
+
+
+def fit_sphere(point_list):
+    """
+    Performs Leo Dorsts sphere fitting technique
+    """
+    return layout.MultiVector(value=val_fit_sphere(np.array(point_list)))
+
+
+@numba.njit
+def val_fit_plane(point_list):
+    """
+    Does plane fitting with combo J.Lasenbys method and L. Dorsts
+    """
+    accumulator_matrix = np.zeros((32,32))
+    for i in range(point_list.shape[0]):
+        P_i_l = get_left_gmt_matrix(point_list[i,:])
+        P_i_r = get_right_gmt_matrix(point_list[i,:])
+        accumulator_matrix += mask4@P_i_l@P_i_r
+    e_vals, e_vecs = np.linalg.eig(accumulator_matrix)
+    # Find the smallest non negative eigenvalues
+    min_eval = np.inf
+    min_eval_index = -1
+    for i in range(len(e_vals)):
+        if e_vals[i] < min_eval and e_vals[i] > 0:
+            min_eval = e_vals[i]
+            min_eval_index = i
+    best_plane = val_normalised(e_vecs[:, min_eval_index])
+    return best_plane
+
+
+def fit_plane(point_list):
+    """
+    Does plane fitting with combo J.Lasenbys method and L. Dorsts
+    """
+    return layout.MultiVector(value=val_fit_plane(np.array(point_list)))

test/test_g3c_tools.py

@@ -11,6 +11,7 @@
 from clifford.tools.g3c.rotor_estimation import *
 from clifford.tools.g3c.object_clustering import *
 from clifford.tools.g3c.scene_simplification import *
+from clifford.tools.g3c.object_fitting import *
 from clifford.tools.g3c.model_matching import *
 from clifford.tools.g3 import random_euc_mv
 from clifford.tools.g3c.GAOnline import draw_objects, GAScene, GanjaScene
@@ -29,6 +30,47 @@
 import random
 
 
+class TestFitObjects(unittest.TestCase):
+    def test_fit_circle(self):
+        noise = 0.1
+        trueP = random_circle()
+        point_list = project_points_to_circle([random_conformal_point() for i in range(100)], trueP)
+        point_list = [up(down(P) + noise * random_euc_mv()) for P in point_list]
+        print(trueP)
+        circle = fit_circle(point_list)
+        print(circle)
+        #draw(point_list + [circle], static=False, scale=0.1)
+
+    def test_fit_line(self):
+        noise = 0.1
+        trueP = random_line()
+        point_list = project_points_to_line([random_conformal_point() for i in range(100)], trueP)
+        point_list = [up(down(P) + noise * random_euc_mv()) for P in point_list]
+        print(trueP)
+        line = fit_line(point_list)
+        print(line)
+        #draw(point_list + [line], static=False, scale=0.1)
+
+    def test_fit_sphere(self):
+        noise = 0.1
+        trueP = random_sphere()
+        point_list = project_points_to_sphere([random_conformal_point() for i in range(100)], trueP)
+        point_list = [up(down(P) + noise * random_euc_mv()) for P in point_list]
+        print(trueP)
+        sphere = fit_sphere(point_list)
+        print(sphere)
+        #draw([sphere] + point_list, static=False, scale=0.1)
+
+    def test_fit_plane(self):
+        noise = 0.1
+        trueP = random_plane()
+        point_list = project_points_to_plane([random_conformal_point() for i in range(100)], trueP)
+        point_list = [up(down(P) + noise * random_euc_mv()) for P in point_list]
+        print(trueP)
+        plane = fit_plane(point_list)
+        print(plane)
+        #draw(point_list + [plane], static=False, scale=0.1)
+
 
 class TestGeneralLogarithm(unittest.TestCase):