FIX barycentric coordinates and output a point

Author: pv

CHANGELOG.md

@@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 * Fixed `compas_libigl` plugins are not detected.
 * Align barycentric coordinates of libigl to COMPAS.
+* Ray mesh intersection now returns a point.
 
 ### Removed
 

docs/examples/example_intersections.py

@@ -17,16 +17,15 @@
 
 trimesh = mesh.copy()
 trimesh.quads_to_triangles()
-
 # ==============================================================================
 # Rays
 # ==============================================================================
 
 base = Point(*mesh.centroid())
 base.z = 0
 
-theta = np.linspace(0, np.pi, 20, endpoint=False)
-phi = np.linspace(0, 2 * np.pi, 20, endpoint=False)
+theta = np.linspace(0, np.pi, 5, endpoint=False)
+phi = np.linspace(0, 2 * np.pi, 5, endpoint=False)
 theta, phi = np.meshgrid(theta, phi)
 theta = theta.ravel()
 phi = phi.ravel()
@@ -39,11 +38,13 @@
 mask = xyz[:, 2] > 0
 hemi = xyz[mask]
 
+lines = []
 rays = []
 for x, y, z in hemi:
     point = Point(x, y, z)
     vector = point - base
     vector.unitize()
+    lines.append(Line(base, base + vector))
     rays.append((base, vector))
 
 # ==============================================================================
@@ -52,18 +53,13 @@
 
 index_face = {index: face for index, face in enumerate(mesh.faces())}
 
-hits_per_ray = intersection_rays_mesh(rays, mesh.to_vertices_and_faces())
+hits_per_rays = intersection_rays_mesh(rays, trimesh.to_vertices_and_faces())
 
-intersections = []
-for ray, hits in zip(rays, hits_per_ray):
-    if hits:
-        idx, u, v, w = hits[0]
-        vertices = mesh.face_vertices(idx)
-        p1 = mesh.vertex_coordinates(vertices[0])
-        p2 = mesh.vertex_coordinates(vertices[1])
-        p3 = mesh.vertex_coordinates(vertices[2])
-        point = barycenter_to_point(u, v, w, p1, p2, p3)
-        intersections.append(point)
+intersection_points = []
+for hit in hits_per_rays:
+    if hit:
+        pt, idx, u, v, w = hit[0]
+        intersection_points.append(pt)
 
 # ==============================================================================
 # Visualisation
@@ -73,7 +69,12 @@
 
 viewer.scene.add(mesh, opacity=0.7, show_points=False)
 
-for intersection in intersections:
-    viewer.scene.add(Line(base, intersection), linecolor=Color.blue(), linewidth=3)
+for point in intersection_points:
+    viewer.scene.add(Line(base, point), linecolor=Color.blue(), linewidth=3)
+    viewer.scene.add(point, pointcolor=Color.red(), pointsize=10)
+
+for line in lines:
+    for i in range(20):
+        viewer.scene.add(line.point_at(i / 20), pointcolor=Color.red(), pointsize=5)
 
 viewer.show()

docs/examples/example_intersections_barycentric.py

@@ -18,20 +18,23 @@
 ray = (p0, compas.geometry.Vector(0, 0, 1))
 hits_per_ray = intersection_ray_mesh(ray, mesh.to_vertices_and_faces())
 
-idx, u, v, w = hits_per_ray[0][0], hits_per_ray[0][1], hits_per_ray[0][2], hits_per_ray[0][3]
+point, idx, u, v, w = hits_per_ray[0][0], hits_per_ray[0][1], hits_per_ray[0][2], hits_per_ray[0][3], hits_per_ray[0][4]
 
 index_face = {index: face for index, face in enumerate(mesh.faces())}
 
 intersections = []
 for hit in hits_per_ray:
-    idx, u, v, w = hit
+    point, idx, w, u, v = hit
     point = barycenter_to_point(u, v, w, p1, p2, p3)
     intersections.append(point)
 
 bary_coords = compas.geometry.barycentric_coordinates(intersections[0], [p1, p2, p3])
-print("libigl barycentric coordinates: ", u, v, w)
+print("libigl barycentric coordinates: ", w, u, v)
 print("compas barycentric coordinates: ", *bary_coords)
 
+print(barycenter_to_point(bary_coords[0], bary_coords[1], bary_coords[2], p1, p2, p3))
+print(barycenter_to_point(w, u, v, p1, p2, p3))
+
 # ==============================================================================
 # Visualisation
 # ==============================================================================
@@ -43,4 +46,6 @@
 for intersection in intersections:
     viewer.scene.add(Line(p0, intersection), linecolor=Color.blue(), linewidth=3)
 
+viewer.scene.add(point, pointcolor=Color.red(), pointsize=10)
+
 viewer.show()

src/compas_libigl/intersections.py

@@ -1,70 +1,86 @@
 import numpy as np
+from compas.geometry import Point
 from compas.plugins import plugin
 
 from compas_libigl import _intersections
 
 
-def _conversion_libigl_to_compas(hits_per_ray):
+def _conversion_libigl_to_compas(hits_per_ray, M):
     """Convert libigl barycentric coordinates to COMPAS barycentric coordinates.
 
     Parameters
     ----------
     hits_per_ray : list[tuple[int, float, float, float]]
         Tuples of (face_index, u, v, distance) from libigl ray intersection
+    M : tuple[list[list[float]], list[list[int]]]
+        A mesh represented by a tuple of (vertices, faces)
+        where vertices are 3D points and faces are triangles
 
     Returns
     -------
-    list[tuple[int, float, float, float]]
-        Tuples of (face_index, w, u, v) in COMPAS barycentric coordinate ordering
+    list[tuple[list[float], int, float, float, float]]
+        Tuples of (point, face_index, u, v, w) in COMPAS barycentric coordinate ordering
 
     Note
     ----
     libigl uses: P = (1-u-v)*v0 + u*v1 + v*v2
-    This function returns [w, u, v] = [1-u-v, u, v] to match COMPAS ordering
+    COMPAS uses: P = u*v0 + v*v1 + w*v2 where u + v + w = 1
+    This function converts libigl coordinates to match COMPAS barycentric coordinate ordering
     """
+    vertices = M[0]
+    faces = M[1]
 
     hits_compas = []
     for h in hits_per_ray:
-        idx, u, v, _ = h
-        w = 1.0 - u - v
-        hits_compas.append([idx, w, v, u])
+        idx, u_libigl, v_libigl, _ = h
+        w = 1.0 - u_libigl - v_libigl  # libigl's (1-u-v) coefficient
+        u = u_libigl  # libigl's u coefficient  
+        v = v_libigl  # libigl's v coefficient
+
+        face = faces[idx]
+        p1, p2, p3 = vertices[face[0]], vertices[face[1]], vertices[face[2]]
+        point = barycenter_to_point(u, v, w, p1, p2, p3)
+        
+        # To match COMPAS barycentric coordinates exactly:
+        # COMPAS expects coordinates in order [p1_weight, p2_weight, p3_weight]
+        # Our formula is P = w*p1 + u*p2 + v*p3, so COMPAS order should be [w, u, v]
+        hits_compas.append([point, idx, u, v, w])
     return hits_compas
 
 
 def barycenter_to_point(u, v, w, p1, p2, p3):
-    """Convert COMPAS barycentric coordinates to a point.
+    """Convert barycentric coordinates to a point using the working interpolation formula.
 
     Parameters
     ----------
     u : float
-        The u coordinate
+        The u coordinate (weight for p2)
     v : float
-        The v coordinate
+        The v coordinate (weight for p3)
     w : float
-        The w coordinate
+        The w coordinate (weight for p1)
     p1 : tuple[float, float, float]
-        The first point
+        The first vertex
     p2 : tuple[float, float, float]
-        The second point
+        The second vertex
     p3 : tuple[float, float, float]
-        The third point
-
+        The third vertex
 
     Returns
     -------
-    list[float]
-        The point at the intersection of the ray and the mesh
+    Point
+        The interpolated point
 
     Note
     ----
-    libigl uses: P = (1-u-v)*v0 + u*v1 + v*v2
-    This function returns [w, u, v] = [1-u-v, u, v] to match COMPAS ordering
+    Uses barycentric interpolation: P = w*p1 + u*p2 + v*p3
+    where w + u + v = 1
     """
-    w = 1 - u - v  # barycentric coordinates
-
-    phit = [u * p1[0] + v * p2[0] + w * p3[0], u * p1[1] + v * p2[1] + w * p3[1], u * p1[2] + v * p2[2] + w * p3[2]]
+    phit = [w * p1[0] + u * p2[0] + v * p3[0], 
+            w * p1[1] + u * p2[1] + v * p3[1], 
+            w * p1[2] + u * p2[2] + v * p3[2]]
 
-    return phit
+    return Point(*phit)
 
 
 @plugin(category="intersections")
@@ -81,21 +97,23 @@ def intersection_ray_mesh(ray, M):
 
     Returns
     -------
-    list[tuple[int, float, float, float]]
+    list[tuple[list[float], int, float, float, float]]
         The array contains a tuple per intersection of the ray with the mesh.
         Each tuple contains:
 
-        0. the index of the intersected face
-        1. the u coordinate of the intersection in the barycentric coordinates of the face
-        2. the v coordinate of the intersection in the barycentric coordinates of the face
-        3. the distance between the ray origin and the hit
+        0. the point of intersection
+        1. the index of the intersected face
+        2. the u coordinate of the intersection in COMPAS barycentric coordinates
+        3. the v coordinate of the intersection in COMPAS barycentric coordinates
+        4. the w coordinate of the intersection in COMPAS barycentric coordinates
+        
 
         Note
         ----
-        The barycentric coordinates (u, v) follow the libigl convention where:
-        - For a triangle with vertices (v0, v1, v2) at face indices F[face_id]
-        - The intersection point P = (1-u-v)*v0 + u*v1 + v*v2
-        - This differs from COMPAS barycentric_coordinates which uses a different vertex ordering
+        The returned barycentric coordinates follow COMPAS convention where:
+        - For a triangle with vertices (p1, p2, p3) at face indices F[face_id]
+        - The intersection point P = u*p1 + v*p2 + w*p3 where u + v + w = 1
+        - These coordinates match those returned by compas.geometry.barycentric_coordinates
     """
     point, vector = ray
     vertices, faces = M
@@ -107,7 +125,7 @@ def intersection_ray_mesh(ray, M):
     hits_per_ray = _intersections.intersection_ray_mesh(P, D, V, F)
 
     # Convert libigl barycentric coordinates to COMPAS convention
-    hits_compas = _conversion_libigl_to_compas(hits_per_ray)
+    hits_compas = _conversion_libigl_to_compas(hits_per_ray, M)
 
     return hits_compas
 
@@ -125,14 +143,16 @@ def intersection_rays_mesh(rays, M):
 
     Returns
     -------
-    list[list[tuple[int, float, float, float]]]
+    list[list[tuple[list[float], int, float, float, float]]]
         List of intersection results, one per ray.
         Each intersection result contains tuples with:
 
-        0. the index of the intersected face
-        1. the u coordinate of the intersection in the barycentric coordinates of the face
-        2. the v coordinate of the intersection in the barycentric coordinates of the face
-        3. the distance between the ray origin and the hit
+        0. the point of intersection
+        1. the index of the intersected face
+        2. the u coordinate of the intersection in COMPAS barycentric coordinates
+        3. the v coordinate of the intersection in COMPAS barycentric coordinates
+        4. the w coordinate of the intersection in COMPAS barycentric coordinates
+        
     """
     points, vectors = zip(*rays)
     vertices, faces = M
@@ -146,6 +166,6 @@ def intersection_rays_mesh(rays, M):
     # Convert libigl barycentric coordinates to COMPAS convention
     hits_per_ray_compas = []
     for hit in hits_per_ray:
-        hits_per_ray_compas.append(_conversion_libigl_to_compas(hit))
+        hits_per_ray_compas.append(_conversion_libigl_to_compas(hit, M))
 
     return hits_per_ray_compas