diff --git a/CHANGELOG.md b/CHANGELOG.md
index a805765..169062b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -30,6 +30,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 
 * Fixed `compas_libigl` plugins are not detected.
+* Align barycentric coordinates of libigl to COMPAS.
+* Ray mesh intersection now returns a point.
 * Add project dependency groups in pyproject.toml.
 * Change requirements.txt to pyproject.toml.
 
diff --git a/docs/examples/example_intersections.py b/docs/examples/example_intersections.py
index 301384d..df735b8 100644
--- a/docs/examples/example_intersections.py
+++ b/docs/examples/example_intersections.py
@@ -7,6 +7,7 @@
 from compas_viewer import Viewer
 
 from compas_libigl.intersections import intersection_rays_mesh
+from compas_libigl.intersections import barycenter_to_point
 
 # ==============================================================================
 # Input geometry
@@ -16,7 +17,6 @@
 
 trimesh = mesh.copy()
 trimesh.quads_to_triangles()
-
 # ==============================================================================
 # Rays
 # ==============================================================================
@@ -24,8 +24,8 @@
 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()
@@ -38,30 +38,27 @@
 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))
 
 # ==============================================================================
 # Intersections
 # ==============================================================================
 
-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:
-        base, vector = ray
-        index = hits[0][0]
-        distance = hits[0][3]
-        face = index_face[index]
-        point = base + vector * distance
-        intersections.append(point)
+intersection_points = []
+for hits_per_ray in hits_per_rays:
+    if hits_per_ray:
+        for hit in hits_per_ray:
+            pt, idx, u, v, w = hit
+            intersection_points.append(pt)
 
 # ==============================================================================
 # Visualisation
@@ -71,7 +68,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()
diff --git a/docs/examples/example_intersections_barycentric.py b/docs/examples/example_intersections_barycentric.py
new file mode 100644
index 0000000..903cc8e
--- /dev/null
+++ b/docs/examples/example_intersections_barycentric.py
@@ -0,0 +1,49 @@
+import compas.geometry
+import compas.datastructures
+from compas_libigl.intersections import intersection_ray_mesh
+from compas_libigl.intersections import barycenter_to_point
+from compas_viewer import Viewer
+from compas.colors import Color
+from compas.geometry import Line
+import compas
+
+
+p0 = compas.geometry.Point(2, 0, 0)
+p1 = compas.geometry.Point(3 + 2, 0 - 2, 13)
+p2 = compas.geometry.Point(0 - 2, 0 - 2, 10)
+p3 = compas.geometry.Point(0 - 2, 2 + 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_per_ray = intersection_ray_mesh(ray, mesh.to_vertices_and_faces())
+
+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]
+
+intersections = []
+for hit in hits_per_ray:
+    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: ", 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
+# ==============================================================================
+
+viewer = Viewer(width=1600, height=900)
+
+viewer.scene.add(mesh, opacity=0.7, show_points=False)
+
+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()
diff --git a/docs/examples/example_mapping.py b/docs/examples/example_mapping.py
index 3a080a8..4bc585a 100644
--- a/docs/examples/example_mapping.py
+++ b/docs/examples/example_mapping.py
@@ -61,8 +61,8 @@
 # Offset mesh by normals, normals are interpolated from the original mesh.
 # ==============================================================================
 mesh_mapped_offset = mesh_mapped.copy()
-for i in range(mesh_mapped.number_of_vertices()):   
-    mesh_mapped_offset.vertex_attributes(i, "xyz", mesh_mapped.vertex_attributes(i, "xyz") - mn[i]*0.001)
+for i in range(mesh_mapped.number_of_vertices()):
+    mesh_mapped_offset.vertex_attributes(i, "xyz", mesh_mapped.vertex_attributes(i, "xyz") - mn[i] * 0.001)
 
 # ==============================================================================
 # Get Boundary Polylines
@@ -74,7 +74,7 @@
     points = []
     for j in range(len(mf[i])):
         id = mf[i][j]
-        points.append(mesh_mapped.vertex_attributes(id, "xyz") + mn[id]*0.002)
+        points.append(mesh_mapped.vertex_attributes(id, "xyz") + mn[id] * 0.002)
     points.append(points[0])
     polyline = Polyline(points)
     boundaries.append(polyline)
diff --git a/docs/examples/example_mapping_patterns.py b/docs/examples/example_mapping_patterns.py
index 64db92b..01f84ee 100644
--- a/docs/examples/example_mapping_patterns.py
+++ b/docs/examples/example_mapping_patterns.py
@@ -9,7 +9,6 @@
 from compas_libigl.mapping import map_pattern_to_mesh
 
 
-
 # ==============================================================================
 # Input geometry: 3D Mesh
 # ==============================================================================
@@ -17,7 +16,6 @@
 mesh = Mesh.from_obj(Path(__file__).parent.parent.parent / "data" / "minimal_surface.obj")
 
 
-
 for vertex in mesh.vertices():
     x, y, z = mesh.vertex_attributes(vertex, "xyz")  # type: ignore
     mesh.vertex_attributes(vertex, "xyz", [x, -z, y])
@@ -32,7 +30,7 @@
 
 for vertex in mesh.vertices():
     x, y, z = mesh.vertex_attributes(vertex, "xyz")  # type: ignore
-    if abs(z-aabb.zmin) < 1e-3 or abs(z-aabb.zmax) < 1e-3:
+    if abs(z - aabb.zmin) < 1e-3 or abs(z - aabb.zmax) < 1e-3:
         fixed_vertices.append(vertex)
 
 # ==============================================================================
diff --git a/pyproject.toml b/pyproject.toml
index 45601f4..76838a9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -161,4 +161,4 @@ max-doc-length = 179
 
 [tool.ruff.format]
 docstring-code-format = true
-docstring-code-line-length = "dynamic"
+docstring-code-line-length = "dynamic"
\ No newline at end of file
diff --git a/src/compas_libigl/intersections.py b/src/compas_libigl/intersections.py
index 05fb692..b1ab715 100644
--- a/src/compas_libigl/intersections.py
+++ b/src/compas_libigl/intersections.py
@@ -1,9 +1,88 @@
 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, 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[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
+    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_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 barycentric coordinates to a point using the working interpolation formula.
+
+    Parameters
+    ----------
+    u : float
+        The u coordinate (weight for p2)
+    v : float
+        The v coordinate (weight for p3)
+    w : float
+        The w coordinate (weight for p1)
+    p1 : tuple[float, float, float]
+        The first vertex
+    p2 : tuple[float, float, float]
+        The second vertex
+    p3 : tuple[float, float, float]
+        The third vertex
+
+    Returns
+    -------
+    Point
+        The interpolated point
+
+    Note
+    ----
+    Uses barycentric interpolation: P = w*p1 + u*p2 + v*p3
+    where w + u + v = 1
+    """
+    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 Point(*phit)
+
+
 @plugin(category="intersections")
 def intersection_ray_mesh(ray, M):
     """Compute the intersection(s) between a ray and a mesh.
@@ -18,14 +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 u 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 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
@@ -33,7 +121,13 @@ def intersection_ray_mesh(ray, M):
     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)
+
+    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, M)
+
+    return hits_compas
 
 
 def intersection_rays_mesh(rays, M):
@@ -49,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 u 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
@@ -64,4 +160,12 @@ def intersection_rays_mesh(rays, M):
     D = np.asarray(vectors, dtype=np.float64)
     V = np.asarray(vertices, dtype=np.float64)
     F = np.asarray(faces, dtype=np.int32)
-    return _intersections.intersection_rays_mesh(P, D, V, F)
+
+    hits_per_ray = _intersections.intersection_rays_mesh(P, D, V, F)
+
+    # 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, M))
+
+    return hits_per_ray_compas