Search code examples
pythonopen3d

Crop function that slices triangles instead of removing them (open3d)


I have a TriangleMesh in open3d and I would like to crop it using a bounding box.

Open3d has the crop function, which removes triangles if they are fully or partially outside the bounding box.

Is there a function that slices triangles instead if are partially outside the bounding box?

Here is a simple 2D example (see plot below). Given the bounding box and the input triangle, the open3d crop function would simply remove the triangle. I would like a function that takes this triangle that overlaps with the bounding box and slices it. Is there such a function?

2d crop example


Solution

  • I build a solution that can crop AxisAlignedBoundingBoxes (which is what I need):

    Input mesh: enter image description here

    Clean crop with the below code: enter image description here

    Open3d crop: enter image description here

    import open3d as o3d
    import numpy as np
    
    def sliceplane(mesh, axis, value, direction):
        # axis can be 0,1,2 (which corresponds to x,y,z)
        # value where the plane is on that axis
        # direction can be True or False (True means remove everything that is
        # greater, False means less
        # than)
    
        vertices = np.asarray(mesh.vertices)
        triangles = np.asarray(mesh.triangles)
        new_vertices = list(vertices)
        new_triangles = []
    
        # (a, b) -> c
        # c refers to index of new vertex that sits at the intersection between a,b
        # and the boundingbox edge
        # a is always inside and b is always outside
        intersection_edges = dict()
    
        # find axes to compute
        axes_compute = [0,1,2]
        # remove axis that the plane is on
        axes_compute.remove(axis)
    
        def compute_intersection(vertex_in_index, vertex_out_index):
            vertex_in = vertices[vertex_in_index]
            vertex_out = vertices[vertex_out_index]
            if (vertex_in_index, vertex_out_index) in intersection_edges:
                intersection_index = intersection_edges[(vertex_in_index, vertex_out_index)]
                intersection = new_vertices[intersection_index]
            else:
                intersection = [None, None, None]
                intersection[axis] = value
                const_1 = (value - vertex_in[axis])/(vertex_out[axis] - vertex_in[axis])
                c = axes_compute[0]
                intersection[c] = (const_1 * (vertex_out[c] - vertex_in[c])) + vertex_in[c]
                c = axes_compute[1]
                intersection[c] = (const_1 * (vertex_out[c] - vertex_in[c])) + vertex_in[c]
                assert not (None in intersection)
                # save new vertice and remember that this intersection already added an edge
                new_vertices.append(intersection)
                intersection_index = len(new_vertices) - 1
                intersection_edges[(vertex_in_index, vertex_out_index)] = intersection_index
    
            return intersection_index
    
        for t in triangles:
            v1, v2, v3 = t
            if direction:
                v1_out = vertices[v1][axis] > value
                v2_out = vertices[v2][axis] > value
                v3_out = vertices[v3][axis] > value
            else: 
                v1_out = vertices[v1][axis] < value
                v2_out = vertices[v2][axis] < value
                v3_out = vertices[v3][axis] < value
    
            bool_sum = sum([v1_out, v2_out, v3_out])
            # print(f"{v1_out=}, {v2_out=}, {v3_out=}, {bool_sum=}")
    
            if bool_sum == 0:
                # triangle completely inside --> add and continue
                new_triangles.append(t)
            elif bool_sum == 3:
                # triangle completely outside --> skip
                continue
            elif bool_sum == 2:
                # two vertices outside 
                # add triangle using both intersections
                vertex_in_index = v1 if (not v1_out) else (v2 if (not v2_out) else v3)
                vertex_out_1_index = v1 if v1_out else (v2 if v2_out else v3)
                vertex_out_2_index = v3 if v3_out else (v2 if v2_out else v1)
                # print(f"{vertex_in_index=}, {vertex_out_1_index=}, {vertex_out_2_index=}")
                # small sanity check if indices sum matches
                assert sum([vertex_in_index, vertex_out_1_index, vertex_out_2_index]) == sum([v1,v2,v3])
    
                # add new triangle 
                new_triangles.append([vertex_in_index, compute_intersection(vertex_in_index, vertex_out_1_index), 
                    compute_intersection(vertex_in_index, vertex_out_2_index)])
    
            elif bool_sum == 1:
                # one vertice outside
                # add three triangles
                vertex_out_index = v1 if v1_out else (v2 if v2_out else v3)
                vertex_in_1_index = v1 if (not v1_out) else (v2 if (not v2_out) else v3)
                vertex_in_2_index = v3 if (not v3_out) else (v2 if (not v2_out) else v1)
                # print(f"{vertex_out_index=}, {vertex_in_1_index=}, {vertex_in_2_index=}")
                # small sanity check if outdices sum matches
                assert sum([vertex_out_index, vertex_in_1_index, vertex_in_2_index]) == sum([v1,v2,v3])
    
                new_triangles.append([vertex_in_1_index, compute_intersection(vertex_in_1_index, vertex_out_index), vertex_in_2_index])
                new_triangles.append([compute_intersection(vertex_in_1_index, vertex_out_index), 
                    compute_intersection(vertex_in_2_index, vertex_out_index), vertex_in_2_index])
    
            else:
                assert False
    
        # TODO remap indices and remove unused 
    
        mesh = o3d.geometry.TriangleMesh()
        mesh.vertices = o3d.utility.Vector3dVector(np.array(new_vertices))
        mesh.triangles = o3d.utility.Vector3iVector(np.array(new_triangles))
        return mesh
    
    def clean_crop_xy(mesh, min_corner, max_corner):
        min_x = min(min_corner[0], max_corner[0])
        min_y = min(min_corner[1], max_corner[1])
        max_x = max(min_corner[0], max_corner[0])
        max_y = max(min_corner[1], max_corner[1])
    
        # mesh = sliceplane(mesh, 0, min_x, False)
        mesh_sliced = sliceplane(mesh, 0, max_x, True)
        mesh_sliced = sliceplane(mesh_sliced, 0, min_x, False)
        mesh_sliced = sliceplane(mesh_sliced, 1, max_y, True)
        mesh_sliced = sliceplane(mesh_sliced, 1, min_y, False)
        # mesh_sliced = mesh_sliced.paint_uniform_color([0,0,1])
    
        return mesh_sliced
    
    knot_mesh = o3d.data.KnotMesh()
    mesh = o3d.io.read_triangle_mesh(knot_mesh.path)
    v = np.asarray(mesh.vertices)
    
    mesh_cropped = clean_crop_xy(mesh, (-200, 0), (0, -200))
    bounding_box = o3d.geometry.AxisAlignedBoundingBox((-200, 0, -1000), (0,200,1000))
    mesh_crop_bad = mesh.crop(bounding_box)
    
    o3d.visualization.draw_geometries([mesh], mesh_show_back_face=True, mesh_show_wireframe=True)
    o3d.visualization.draw_geometries([mesh_crop_bad], mesh_show_back_face=True, mesh_show_wireframe=True)
    o3d.visualization.draw_geometries([mesh_cropped], mesh_show_back_face=True, mesh_show_wireframe=True)