Geometric distance can now return the points used for computing the distance

src/geomtree/affine_transform.py

@@ -233,7 +233,7 @@ def belongs(tf,x,alpha=0.95):
     b = T @ tf.translation
     return norm(T @ x -b)#<np.sqrt(chi2.ppf(alpha, 3))
 
-def geometric_distance(tf1,tf2,alpha=0.95,epsilon=1e-6,iterations_max=float('Inf'),view=True): #wip
+def geometric_distance(tf1,tf2,alpha=0.95,epsilon=1e-6,iterations_max=float('Inf'),output=False,view=False): #wip
     # Algo from "On the distance between two ellipsoids" by A. Lin and S.-P. Han
 
     import sys
@@ -276,9 +276,15 @@ def geometric_distance(tf1,tf2,alpha=0.95,epsilon=1e-6,iterations_max=float('Inf
         t1 = max(s1) if len(s1)>0 else 1
         t2 = min(s2) if len(s2)>0 else 0
 
+        xbar = c1 + t1 * (c2 - c1)
+        ybar = c1 + t2 * (c2 - c1)
+
         if t2<t1:
             #print('\n')
-            return 0
+            if output:
+                return 0, xbar,ybar
+            else:
+                return 0
 
         if view:
             vec = lambda t: (1 - t) * c1 + t * c2
@@ -293,9 +299,6 @@ def geometric_distance(tf1,tf2,alpha=0.95,epsilon=1e-6,iterations_max=float('Inf
             plt.ylim(0, chi2.ppf(alpha, 3))
             plt.show()
 
-        xbar = c1 + t1 * (c2 - c1)
-        ybar = c1 + t2 * (c2 - c1)
-
         if view:
             pc = PointCloud([c1, c2, xbar, ybar])
             actors = tf1.get_actor() + tf2.get_actor() + pc.get_actor(color='green', glyph='sphere', scale=0.05) + [segment_primitive(c1, c2, color='#808080', lw=1)]
@@ -305,7 +308,10 @@ def geometric_distance(tf1,tf2,alpha=0.95,epsilon=1e-6,iterations_max=float('Inf
 
         if theta1<epsilon and theta2<epsilon:
             #print('\n')
-            return norm(xbar-ybar)
+            if output:
+                return norm(xbar-ybar),xbar,ybar
+            else:
+                return norm(xbar-ybar)
 
         gamma1 = 1 / norm(inv(
             tf1.linear @ tf1.linear.transpose()))  # < 1/max(abs(np.linalg.eigvals(inv(tf1.linear @ tf1.linear.transpose()))))
@@ -330,7 +336,10 @@ def geometric_distance(tf1,tf2,alpha=0.95,epsilon=1e-6,iterations_max=float('Inf
         j+=1
 
     # out of the loop
-    return norm(xbar - ybar)
+    if output:
+        return norm(xbar - ybar), xbar, ybar
+    else:
+        return norm(xbar - ybar)