geometry_utils

Utilities to process geometry (validation, triangulation, fusion, etc).

Functions:

Name	Description
fix_polygon_2d	Transform an input polygon into one or multiple valid polygons.
flatten_trimesh	Flatten the mesh on the Z-axis, and triangulate it.
merge_trimeshes	Merge the given Trimesh objects into a single Trimesh.
orient_polygons_z_up	Orient the given mesh so that the normal of all its triangles point up.
triangulate_linear_ring_2d	Triangulate a 2D LinearRing.
triangulate_polygon_2d	Triangulate a 2D Polygon.
triangulate_surface_3d	Triangulate a 3D Surface.

fix_polygon_2d(polygon)

Transform an input polygon into one or multiple valid polygons.

Parameters:

Name	Type	Description	Default
polygon	shapely.Polygon	The potentially invalid 2D polygon.	required

Returns:

Type	Description
list[shapely.Polygon]	The valid polygons equivalent to the input polygon.

Raises:

Type	Description
RuntimeError	If the operations don't succeed to make valid polygons.

Source code in python/src/data_pipeline/utils/geometry_utils.py

def fix_polygon_2d(polygon: Polygon) -> list[Polygon]:
    """
    Transform an input polygon into one or multiple valid polygons.

    Parameters
    ----------
    polygon : Polygon
        The potentially invalid 2D polygon.

    Returns
    -------
    list[Polygon]
        The valid polygons equivalent to the input polygon.

    Raises
    ------
    RuntimeError
        If the operations don't succeed to make valid polygons.
    """
    polygon = shapely.remove_repeated_points(polygon)
    if not shapely.is_valid(polygon):
        valid_geometry = shapely.make_valid(polygon, method="structure")
    else:
        valid_geometry = polygon
    if isinstance(valid_geometry, MultiPolygon):
        return list(valid_geometry.geoms)
    elif isinstance(valid_geometry, Polygon):
        return [valid_geometry]
    elif isinstance(valid_geometry, LineString):
        # if valid_geometry.is_closed:
        #     return [Polygon(shell=valid_geometry)]
        # else:
        #     return []
        return []
    else:
        raise RuntimeError(
            f"Unexpected output after making the polygon valid: {valid_geometry}"
        )

flatten_trimesh(mesh, z_value=None)

Flatten the mesh on the Z-axis, and triangulate it.

Warning

This function is a bit buggy, and the triangulating algorithm sometimes crashes without giving any error message. This seems to be often related to problems in the input that are not well handled by the algorithm, such as duplicate vertices.

Parameters:

Name	Type	Description	Default
mesh	trimesh.Trimesh	The Trimesh to flatten.	required
z_value	numpy.float64 | None	The Z value to assign to all vertices. If None, it is the minimum of the Z values of all vertices in the mesh. By default None.	None

Returns:

Type	Description
trimesh.Trimesh	The flattened Trimesh.

Source code in python/src/data_pipeline/utils/geometry_utils.py

def flatten_trimesh(
    mesh: trimesh.Trimesh, z_value: np.float64 | None = None
) -> trimesh.Trimesh:
    """
    Flatten the mesh on the Z-axis, and triangulate it.

    Warning
    -------
    This function is a bit buggy, and the triangulating algorithm sometimes crashes without giving any error message.
    This seems to be often related to problems in the input that are not well handled by the algorithm, such as duplicate vertices.

    Parameters
    ----------
    mesh : trimesh.Trimesh
        The Trimesh to flatten.
    z_value : np.float64 | None, optional
        The Z value to assign to all vertices.
        If None, it is the minimum of the Z values of all vertices in the mesh.
        By default None.

    Returns
    -------
    trimesh.Trimesh
        The flattened Trimesh.
    """
    try:
        logging.debug("Start flattening")
        z_value = np.min(mesh.vertices[:, 2]) if z_value is None else z_value
        # apad is set to 1 mm to prevent errors from araising with triangulation
        # TODO: understand the error better and fix it in a better way?
        flatten_path_2D: Path2D = mesh.projected(normal=(0, 0, 1), apad=0.001)
        logging.debug("Start triangulating")
        vertices, faces = flatten_path_2D.triangulate(
            **{"engine": "triangle", "triangle_args": "p"}
        )
        flat_mesh = trimesh.Trimesh(vertices=vertices, faces=faces)
        flat_mesh.vertices = np.hstack(
            (flat_mesh.vertices, np.full((flat_mesh.vertices.shape[0], 1), z_value))
        )
        orient_polygons_z_up(flat_mesh)
        return flat_mesh
    except Exception as e:
        logging.error(f"Error during flattening the Trimesh: {e}")
        raise e

merge_trimeshes(meshes, fix_geometry)

Merge the given Trimesh objects into a single Trimesh. Optionally fix the geometry of the final mesh.

Parameters:

Name	Type	Description	Default
meshes	list[trimesh.Trimesh]	List of meshes to merge.	required
fix_geometry	bool	Whether to apply different operations to fix the mesh.	required

Returns:

Type	Description
trimesh.Trimesh	The merged mesh.

Raises:

Type	Description
RuntimeError	If the object obtained from merging is not a Trimesh.

Source code in python/src/data_pipeline/utils/geometry_utils.py

def merge_trimeshes(
    meshes: list[trimesh.Trimesh], fix_geometry: bool
) -> trimesh.Trimesh:
    """
    Merge the given Trimesh objects into a single Trimesh.
    Optionally fix the geometry of the final mesh.

    Parameters
    ----------
    meshes : list[trimesh.Trimesh]
        List of meshes to merge.
    fix_geometry : bool
        Whether to apply different operations to fix the mesh.

    Returns
    -------
    trimesh.Trimesh
        The merged mesh.

    Raises
    ------
    RuntimeError
        If the object obtained from merging is not a Trimesh.
    """
    full_mesh = trimesh.util.concatenate(meshes)
    if not isinstance(full_mesh, trimesh.Trimesh):
        raise RuntimeError("The combination of the meshes isn't a Trimesh.")
    if full_mesh.is_empty:
        logging.debug("Empty!")
        return full_mesh
    if fix_geometry:
        full_mesh.process(validate=True)
        full_mesh.update_faces(full_mesh.unique_faces())
        full_mesh.fix_normals()
        full_mesh.fill_holes()
    else:
        full_mesh.process()
    return full_mesh

orient_polygons_z_up(mesh)

Orient the given mesh so that the normal of all its triangles point up. The input mesh is expected to be flat, parallel to the x,y plane.

Warning

This function modifies the mesh directly.

Parameters:

Name	Type	Description	Default
mesh	trimesh.Trimesh	The input flat mesh.	required

Raises:

Type	Description
RuntimeError	If the mesh has a triangle that is not parallel to the x,y plane, with a difference larger than the threshold used.

Source code in python/src/data_pipeline/utils/geometry_utils.py

def orient_polygons_z_up(mesh: trimesh.Trimesh) -> None:
    """
    Orient the given mesh so that the normal of all its triangles point up.
    The input mesh is expected to be flat, parallel to the x,y plane.

    Warning
    -------
    This function modifies the mesh directly.

    Parameters
    ----------
    mesh : trimesh.Trimesh
        The input flat mesh.

    Raises
    ------
    RuntimeError
        If the mesh has a triangle that is not parallel to the x,y plane, with a difference larger than the threshold used.
    """
    faces = mesh.faces.copy()
    for i, normal in enumerate(mesh.face_normals):
        # Check if the face is actually in the x,y plane
        if not np.all(np.abs(normal[:2]) < 1e-3):
            raise RuntimeError(
                f"A value is larger than expected in absolute value: {np.max(np.abs(normal[:2]))}"
            )

        if normal[2] < 0:
            faces[i, 1], faces[i, 2] = faces[i, 2], faces[i, 1]

    mesh.faces = faces

triangulate_linear_ring_2d(lring_2d, holes_lrings_2d=None)

Triangulate a 2D LinearRing.

Parameters:

Name	Type	Description	Default
lring_2d	shapely.geometry.LinearRing	The outer boundary of the 2D LinearRing to triangulate.	required
holes_lrings_2d	list[shapely.geometry.LinearRing] | None	The list of 2D LinearRings representing the holes. By default None.	None

Returns:

Name	Type	Description
vertices	numpy.typing.NDArray[numpy.float64]	The vertices of the output triangulation.
triangles	numpy.typing.NDArray[numpy.int64]	The triangles of the output triangulation, refering to vertices.

Source code in python/src/data_pipeline/utils/geometry_utils.py

def triangulate_linear_ring_2d(
    lring_2d: sg.LinearRing,
    holes_lrings_2d: list[sg.LinearRing] | None = None,
) -> Tuple[NDArray[np.float64], NDArray[np.int64]]:
    """
    Triangulate a 2D LinearRing.

    Parameters
    ----------
    lring_2d : sg.LinearRing
        The outer boundary of the 2D LinearRing to triangulate.
    holes_lrings_2d : list[sg.LinearRing] | None, optional
        The list of 2D LinearRings representing the holes.
        By default None.

    Returns
    -------
    vertices: NDArray[np.float64]
        The vertices of the output triangulation.
    triangles: NDArray[np.int64]
        The triangles of the output triangulation, refering to `vertices`.
    """
    # Check that the holes inputs are correct
    if holes_lrings_2d is None:
        holes_lrings_2d = []

    # Prepare the triangulation
    points_2d = np.array(lring_2d.coords[:-1], dtype=np.float64)[::-1]
    ring_segments = [(i, (i + 1) % len(points_2d)) for i in range(len(points_2d))]

    for hole_lring in holes_lrings_2d:
        offset = points_2d.shape[0]
        hole_points_2d = np.array(hole_lring.coords[:-1], dtype=np.float64)[::-1]
        ring_segments.extend(
            [
                (offset + i, offset + (i + 1) % len(hole_points_2d))
                for i in range(len(hole_points_2d))
            ]
        )
        points_2d = np.concat((points_2d, hole_points_2d))

    # Remove duplicate vertices
    tol = 1e-6
    scale = 1.0 / tol
    rounded = np.round(points_2d * scale).astype(np.int64)

    # Use a structured array so np.unique can work on rows efficiently
    dtype = np.dtype([("x", np.int64), ("y", np.int64)])
    structured = rounded.view(dtype).reshape(-1)

    _, uniq_idx, inv_map = np.unique(structured, return_index=True, return_inverse=True)

    previous_len = points_2d.shape[0]
    points_2d = points_2d[uniq_idx]
    if points_2d.shape[0] != previous_len:
        logging.debug(f"{previous_len} -> {points_2d.shape[0]}")
    ring_segments = [inv_map[np.array(segment)] for segment in ring_segments]

    triangulation_input = {"vertices": points_2d, "segments": ring_segments}
    # logging.debug(triangulation_input)
    polygon_2d = sg.Polygon(lring_2d, holes=holes_lrings_2d)
    # logging.debug(polygon_2d)

    triangulation = tri.triangulate(triangulation_input, "p")
    if "triangles" not in triangulation:
        logging.debug("No triangle!")
        # logging.debug(f"{np.array(lring_2d.coords) = }")
        # logging.debug(f"{holes_lrings_2d = }")
        return np.empty((0, 2), dtype=np.float64), np.empty((0, 3), dtype=np.int64)
    vertices = np.array(triangulation["vertices"], dtype=np.float64)
    triangles = np.array(triangulation["triangles"], dtype=np.int64)

    # Remove triangles that are outside the polygon (especially holes)
    triangles_to_keep = []
    polygon_2d = sg.Polygon(lring_2d, holes=holes_lrings_2d)
    for i, triangle in enumerate(triangles):
        # Check if the centroid of the triangle is in the polygon
        v0, v1, v2 = vertices[triangle]
        centroid = sg.Point((v0 + v1 + v2) / 3)
        if polygon_2d.contains(centroid):
            triangles_to_keep.append(i)

    triangles = triangles[triangles_to_keep]

    return vertices, triangles

triangulate_polygon_2d(polygon)

Triangulate a 2D Polygon.

Parameters:

Name	Type	Description	Default
polygon	shapely.Polygon	The 2D Polygon to triangulate.	required

Returns:

Name	Type	Description
vertices	numpy.typing.NDArray[numpy.float64]	The vertices of the output triangulation.
triangles	numpy.typing.NDArray[numpy.int64]	The triangles of the output triangulation, refering to vertices.
valid	bool	Whether the output is a valid triangulation.

Source code in python/src/data_pipeline/utils/geometry_utils.py

def triangulate_polygon_2d(
    polygon: Polygon,
) -> Tuple[NDArray[np.float64], NDArray[np.int64], bool]:
    """
    Triangulate a 2D Polygon.

    Parameters
    ----------
    polygon : Polygon
        The 2D Polygon to triangulate.

    Returns
    -------
    vertices: NDArray[np.float64]
        The vertices of the output triangulation.
    triangles: NDArray[np.int64]
        The triangles of the output triangulation, refering to `vertices`.
    valid: bool
        Whether the output is a valid triangulation.
    """
    # Return nothing if the ring doesn't have enough points to compute a normal
    all_points = list(polygon.exterior.coords[:-1])
    for interior in polygon.interiors:
        all_points.extend(interior.coords[:-1])
    unique_points = np.unique(np.concat(all_points), axis=0)
    if unique_points.shape[0] < 3:
        return (
            np.zeros((0, 2), dtype=np.float64),
            np.zeros((0, 2), dtype=np.int64),
            False,
        )

    # Triangulate the polygon
    holes_lrings_2d = list(polygon.interiors)
    tri_vertices, tri_faces = triangulate_linear_ring_2d(
        polygon.exterior,
        holes_lrings_2d=holes_lrings_2d,
    )

    return tri_vertices, tri_faces, True

triangulate_surface_3d(outer_boundary, holes, vertices)

Triangulate a 3D Surface.

Parameters:

Name	Type	Description	Default
outer_boundary	numpy.typing.NDArray[numpy.int64]	Outer boundary of the Surface, referring to vertices.	required
holes	list[numpy.typing.NDArray[numpy.int64]]	Holes of the Surface, referring to vertices	required
vertices	numpy.typing.NDArray[numpy.float64]	The vertices to extract the boundaries from.	required

Returns:

Name	Type	Description
tri_vertices	numpy.typing.NDArray[numpy.float64]	The vertices of the output triangulation.
tri_triangles	numpy.typing.NDArray[numpy.int64]	The triangles of the output triangulation, refering to tri_vertices.
valid	bool	Whether the output is a valid triangulation.

Source code in python/src/data_pipeline/utils/geometry_utils.py

def triangulate_surface_3d(
    outer_boundary: NDArray[np.int64],
    holes: list[NDArray[np.int64]],
    vertices: NDArray[np.float64],
) -> Tuple[NDArray[np.float64], NDArray[np.int64], bool]:
    """
    Triangulate a 3D Surface.

    Parameters
    ----------
    outer_boundary : NDArray[np.int64]
        Outer boundary of the Surface, referring to `vertices`.
    holes : list[NDArray[np.int64]]
        Holes of the Surface, referring to `vertices`
    vertices : NDArray[np.float64]
        The vertices to extract the boundaries from.

    Returns
    -------
    tri_vertices: NDArray[np.float64]
        The vertices of the output triangulation.
    tri_triangles: NDArray[np.int64]
        The triangles of the output triangulation, refering to `tri_vertices`.
    valid: bool
        Whether the output is a valid triangulation.
    """
    # Do nothing if the object is already a triangle
    if outer_boundary.shape == 3 and len(holes) == 0:
        tri_vertices = vertices[outer_boundary]
        tri_faces = np.array([0, 1, 2]).reshape(1, 3)
        return tri_vertices, tri_faces, True

    # Compute the plane of the surface
    used_ids = np.unique(np.concatenate([outer_boundary] + holes))
    unique_points = np.unique(vertices[used_ids], axis=0)
    plane = Plane3D.from_points(unique_points)
    if not plane.is_valid:
        logging.debug("Invalid plane!")
        return (
            np.zeros((0, 3), dtype=np.float64),
            np.zeros((0, 3), dtype=np.int64),
            False,
        )

    # Compute the 2D points
    outer_boundary_2d = plane.project_points(vertices[outer_boundary])
    holes_2d = [plane.project_points(vertices[hole]) for hole in holes]

    # Create the polygon (repeat the first point for the rings)
    exterior = np.vstack((outer_boundary_2d, outer_boundary_2d[:1]))
    interiors = [np.vstack((hole_2d, hole_2d[:1])) for hole_2d in holes_2d]
    poly = Polygon(shell=exterior, holes=interiors)

    # Fix the polygon in case of problems
    polys = fix_polygon_2d(polygon=poly)

    # Triangulate all the polygons
    tri_vertices_2d_all = np.empty((0, 2), dtype=np.float64)
    tri_faces_all = np.empty((0, 3), dtype=np.int64)
    worked_all = []
    for poly in polys:
        tri_vertices, tri_faces, worked = triangulate_polygon_2d(poly)
        tri_faces_all = np.vstack(
            (tri_faces_all, tri_faces + tri_vertices_2d_all.shape[0])
        )
        tri_vertices_2d_all = np.vstack((tri_vertices_2d_all, tri_vertices))
        worked_all.append(worked)
    worked = any(worked_all)

    # Unproject the points
    tri_vertices_all = plane.unproject_points(tri_vertices_2d_all)

    return tri_vertices_all, tri_faces_all, worked