Clearer documentation for intersection_ray_mesh

[yiqiaowang-arch]

Feature Request

As a developer for CityEnergyAnalyst, I want clearer documentation for intersection_ray_mesh so that I could understand the relationship between u, v, w and the vertices of the mesh. For example, is p_hit = up1 + vp2 + wp3 or wp1 + vp2 + up3? What's the relationship between this barycentric coordinate and another function provided by compas, barycentric_coordinates?

Details

In the following example code, I want to reconstruct the coordinates of the intersection point between a ray vertically extruded from a point at XY plane and a triangular mesh.

the code is listed as below:

import compas.geometry
import compas.datastructures


p0 = compas.geometry.Point(2, 0, 0)
p1 = compas.geometry.Point(3, 0, 13)
p2 = compas.geometry.Point(0, 0, 10)
p3 = compas.geometry.Point(0, 2, 10)

mesh = compas.datastructures.Mesh.from_points([[p1.x, p1.y, p1.z],
                                               [p2.x, p2.y, p2.z],
                                               [p3.x, p3.y, p3.z]])

ray = (p0, compas.geometry.Vector(0, 0, 1))
hits = compas.geometry.intersection_ray_mesh(ray, mesh.to_vertices_and_faces())
idx, u, v, t = hits[0]
print("intersection length:", t)

w = 1 - u - v # barycentric coordinates
print("barycentric coordinates:", u, v, w)

phit = compas.geometry.Point(
    w * p1.x + v * p2.x + u * p3.x,
    w * p1.y + v * p2.y + u * p3.y,
    w * p1.z + v * p2.z + u * p3.z
)

print("intersection point", phit.x, phit.y, phit.z)

bary_coords = compas.geometry.barycentric_coordinates([p0.x, p0.y, p0.z], [p1, p2, p3])
print("reconstructed barycentric coordinates:", *bary_coords)

And the output is as such:

intersection length: 12.0
barycentric coordinates: -0.0 0.3333333432674408 0.6666666567325592
intersection point 1.9999999701976776 0.0 11.999999970197678
reconstructed barycentric coordinates: -1.3333333333333335 2.3333333333333335 0.0

It seems that w is corresponding to p1, v to p2 and u to p3. However, in the documentation, this is not clearly listed.

Further, the documentation listed the second and the third output to both the u coordinate. I believe that the third output should be v, instead of u.

Finally, the u, v, w is not identical as the reconstructed barycentric coordinates provided by compas.geometry.barycentric_coordinates. Instead of [0, 0.33, 0.67], the compas function outputs [-1.33, 2.33, 0]. Could you help me understand the difference here?

Describe the solution you'd like
It would be nice if clear documentation for the role of u, v, w is listed; also, the typo in documentation should be fixed. Further, it would be beneficial to unify the two barycentric calculations in compas and in compas_libigl.

[tomvanmele]

@petrasvestartas could you follow up?

[petrasvestartas]

Difference:

• libigl

  • Uses the Möller-Trumbore ray-triangle intersection algorithm
  • Formula: P = (1-u-v)v0 + uv1 + v*v2
  • Returns (face_id, u, v, distance)

• COMPAS compas.geometry.barycentric_coordinates(point, triangle):

  • Uses standard barycentric coordinate calculation with dot products
  • Formula: P = uv0 + vv1 + w*v2 where u + v + w = 1

@tomvanmele
Should we update compas_libigl to align to compas notation, or compas to libigl?

This is a correct doc-string for libigl:

@plugin(category="intersections")
def intersection_ray_mesh(ray, M):
    """Compute the intersection(s) between a ray and a mesh.

    Parameters
    ----------
    ray : tuple[list[float], list[float]]
        A ray represented by a point and a direction vector.
    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]]
        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
        
        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
    """
    point, vector = ray
    vertices, faces = M
    P = np.asarray(point, dtype=np.float64)
    D = np.asarray(vector, dtype=np.float64)
    V = np.asarray(vertices, dtype=np.float64)
    F = np.asarray(faces, dtype=np.int32)
    return _intersections.intersection_ray_mesh(P, D, V, F)

To match compas barycentric coordinate calculation you need to do this:

def convert_libigl_barycentric_to_compas(u_libigl, v_libigl):
    """Convert libigl barycentric coordinates to COMPAS convention.
    
    Parameters
    ----------
    u_libigl : float
        The u coordinate from libigl ray intersection
    v_libigl : float  
        The v coordinate from libigl ray intersection
        
    Returns
    -------
    tuple[float, float, float]
        Barycentric coordinates (w, u, v) where w = 1-u-v
        Compatible with 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
    """
    w = 1.0 - u_libigl - v_libigl
    return (w, v_libigl, u_libigl)

[petrasvestartas]

@yiqiaowang-arch

This is the pull request to align libigl notation to compas:
#39

It would be good that you could check this locally on your pc.

[yiqiaowang-arch]

@yiqiaowang-arch

This is the pull request to align libigl notation to compas: #39

It would be good that you could check this locally on your pc.

Thanks! I will check this weekend. :)

[yiqiaowang-arch]

Hi, I have checked everything again, first there was an error in my example code, I should get the barycentric coordinate using phit, not p0 in my original code. I updated the code, in compas_libigl = 0.7.4 I was able to recreate the mismatch between compas.geometry.barycentric_coordinates and compas_libigl.intersections.interaction_ray_mesh:

import compas.geometry
import compas.datastructures
from compas_libigl.intersections import intersection_ray_mesh


p0 = compas.geometry.Point(2, 0, 0)
p1 = compas.geometry.Point(3, 0, 13)
p2 = compas.geometry.Point(0, 0, 10)
p3 = compas.geometry.Point(0, 2, 10)

mesh = compas.datastructures.Mesh.from_points([[p1.x, p1.y, p1.z],
                                               [p2.x, p2.y, p2.z],
                                               [p3.x, p3.y, p3.z]])

ray = (p0, compas.geometry.Vector(0, 0, 1))
hits = intersection_ray_mesh(ray, mesh.to_vertices_and_faces())
idx, u, v, t = hits[0]
print("intersection length:", t)

w = 1 - u - v # barycentric coordinates
print("barycentric coordinates:", u, v, w)

phit = compas.geometry.Point(
    w * p1.x + v * p2.x + u * p3.x,
    w * p1.y + v * p2.y + u * p3.y,
    w * p1.z + v * p2.z + u * p3.z
)

print("intersection point", phit.x, phit.y, phit.z)

bary_coords = compas.geometry.barycentric_coordinates(phit, [p1, p2, p3])  # change from p0 to phit
print("reconstructed barycentric coordinates:", *bary_coords)
print("compas_libigl version: ", compas_libigl.__version__)

output (notice how 0, 0.33, 0.66 turns into 0.66, 0.33, 0):

intersection length: 12.0
barycentric coordinates: -0.0 0.3333333432674408 0.6666666567325592
intersection point 1.9999999701976776 0.0 11.999999970197678
reconstructed barycentric coordinates: 0.6666666567325592 0.3333333432674408 0.0
compas_libigl version:  0.7.4

Then in compas_libigl = 0.7.6, same behaviour is observed:

intersection length: 12.0
barycentric coordinates: -0.0 0.3333333432674408 0.6666666567325592
intersection point 1.9999999701976776 0.0 11.999999970197678
reconstructed barycentric coordinates: 0.6666666567325592 0.3333333432674408 0.0
compas_libigl version:  0.7.6

Finally in PR #39 (the output is modified so I have to update my code):

import compas.geometry
import compas.datastructures
import compas_libigl
from compas_libigl.intersections import intersection_ray_mesh


p0 = compas.geometry.Point(2, 0, 0)
p1 = compas.geometry.Point(3, 0, 13)
p2 = compas.geometry.Point(0, 0, 10)
p3 = compas.geometry.Point(0, 2, 10)

mesh = compas.datastructures.Mesh.from_points([[p1.x, p1.y, p1.z],
                                               [p2.x, p2.y, p2.z],
                                               [p3.x, p3.y, p3.z]])

ray = (p0, compas.geometry.Vector(0, 0, 1))
hits = intersection_ray_mesh(ray, mesh.to_vertices_and_faces())
point, idx, u, v, w = hits[0]
print("intersection length:", point.z)

# w = 1 - u - v # barycentric coordinates
print("barycentric coordinates from compas_libigl:", u, v, w)

phit_wvu = compas.geometry.Point(
    w * p1.x + v * p2.x + u * p3.x,
    w * p1.y + v * p2.y + u * p3.y,
    w * p1.z + v * p2.z + u * p3.z
)
phit_uvw = compas.geometry.Point(
    u * p1.x + v * p2.x + w * p3.x,
    u * p1.y + v * p2.y + w * p3.y,
    u * p1.z + v * p2.z + w * p3.z
)

print("reconstructed intersection point with wp1 + vp2 + up3", phit_wvu.x, phit_wvu.y, phit_wvu.z)
print("reconstructed intersection point with up1 + vp2 + wp3", phit_uvw.x, phit_uvw.y, phit_uvw.z)

bary_coords = compas.geometry.barycentric_coordinates(phit_wvu, [p1, p2, p3])  # change from p0 to phit
print("reconstructed barycentric coordinates:", *bary_coords)
print("compas_libigl version: ", compas_libigl.__version__)

However, the output seems still not aligned with the compas barycentric calculation:

intersection length: 11.999999970197678
barycentric coordinates from compas_libigl: -0.0 0.3333333432674408 0.6666666567325592
reconstructed intersection point with wp1 + vp2 + up3 1.9999999701976776 0.0 11.999999970197678
reconstructed intersection point with up1 + vp2 + wp3 0.0 1.3333333134651184 10.0
reconstructed barycentric coordinates: 0.6666666567325592 0.3333333432674408 0.0
compas_libigl version:  0.7.6

Maybe I am wrong somewhere, but from the image above, p1 should have 0.67 contribution to the phit, p2 0.33 and p3 should have no contribution. the barycentric coordinates from compas_libigl results is however still reversed.

Also, in the newest doc, the output typehint is: list[tuple[list[float], int, float, float, float]], where the first one of the tuple is a list[float].
However, the first item of the tuple is actually a compas.geometry.point.Point.