In this post, I will discuss the process of clipping operations on meshes using the OpenMesh library. For the purpose of this demonstration, we’ll utilize the OpenMesh library to represent our mesh. For those unfamiliar, OpenMesh provides extensive tools for manipulating polygonal meshes. We will be working with a triangulated mesh version for simplicity. Refer to previous post for more details on the OpenMesh library.

Understanding the Plane Equation

The fundamental operation in mesh clipping involves a plane, defined uniquely by a point 𝑃 and a direction 𝑁. The plane equation in 3D space passing through point 𝑃 can be represented as: dot (𝑁, point − 𝑃 ) = 0 where: 𝑁 is a unit normal vector (vector length is always 1), point represents any point on the plane. For instance, consider a point at (5,100,0) with a plane normal pointing in the direction (1,0,0) and P (0,0,0). Substituting these values into the equation gives us 5, indicating that the point is 5 units from the origin along the plane. Positive or negative values will determine which side of the plane the point resides on.

Clipping Operation

The clipping operation involves clipping the input mesh with a defined plane. All triangles on the Positive side of the plane are retained, while those on the Negative are discarded. For triangles that intersect the plane (where some vertices lie on the Positive side and others on the Negative), new “clipped” triangles are created.

Implementing the Clipping Function

I start by writing a function to determine the positional relationship of a triangle relative to the plane. The function takes the mesh, the face, and the plane parameters 𝑁 and 𝑃, and returns whether the triangle falls into the Negative, Mixed, or Positive zone. The algorithm checks all edges of the triangle to ascertain their positions relative to the plane. Here is a basic implementation using OpenMesh:

module;

#include <OpenMesh/Core/IO/MeshIO.hh>
#include <OpenMesh/Core/Mesh/TriMesh_ArrayKernelT.hh>
#include <OpenMesh/Core/Mesh/TriConnectivity.hh>

export module ClipMesh;
import std;

using TriMeshKernel = OpenMesh::TriMesh_ArrayKernelT<OpenMesh::DefaultTraitsDouble>;

enum class Status { Positive, Mixed, Negative };

Status evalClipStatus(const TriMeshKernel& meshData, OpenMesh::FaceHandle fh, OpenMesh::Vec3d P, OpenMesh::Vec3d N) {
	uint32_t positive = 0;
	uint32_t negative = 0;
	uint32_t mixed = 0;

	for (auto heh : meshData.fh_range(fh)) {
		OpenMesh::VertexHandle vh0 = heh.from();
		OpenMesh::VertexHandle vh1 = heh.to();

		auto P0 = meshData.point(vh0);
		auto P1 = meshData.point(vh1);

		double d0 = OpenMesh::dot(N, P0) - OpenMesh::dot(N, P);
		double d1 = OpenMesh::dot(N, P1) - OpenMesh::dot(N, P);
		
		if (d0 <= 0 && d1 <= 0)
			++negative;
		else if (d0 > 0 && d1 > 0)
			++positive;
		else
			++mixed;
	}

	if (mixed > 0)
		return Status::Mixed;

	if (positive == 0 && negative == 3)
		return Status::Negative;

	return Status::Positive;
}

As you can see, the algorithm is quite simple. We loop through all half-edges, take their start and end vertices, plug them into the plane equation, and determine which zone each half-edge vertex falls into. We keep track of these zones using counters for Positive, Negative, and Mixed zones. In the end, if at least one edge is in the Mixed zone, the resulting triangle intersects the plane. If all edges of the triangle are in the Negative zone, the triangle lies behind the plane. The conditions at the end can certainly be simplified, but they are written this way for clarity.

Clip operation

The clip function processes each face of the mesh. Faces in the Negative zone are deleted. Those in the Mixed zone require careful handling to generate new clipped triangles:

export void clip(TriMeshKernel& meshData, OpenMesh::Vec3d P, OpenMesh::Vec3d N) 
{
    for (auto fh : meshData.faces()) {
        Status status = evalClipStatus(meshData, fh, P, N);
        if (status == Status::Negative)
            meshData.delete_face(fh, true);
        else if (status == Status::Mixed) {
         // Detailed logic to handle clipped triangles will be implemented here
        }
    }
}

As we mentioned earlier, in the case of the Negative zone, we completely remove the triangle, i.e., we simply delete it from the mesh. The true parameter means to remove all the vertices of the triangle, not just the face itself. The most challenging part of the clip function remains – generating new clipped triangles that are in the Mixed zone. For this, let’s return to our plane and triangle. Let’s consider the possible scenarios when the plane intersects the triangle.

There are two possible scenarios when a plane intersects a triangle:

Case 1: As a result of the clipping, a triangle remains (P0, Q0, Q1).

Case 2: As a result of the clipping, a convex quadrilateral remains (P0, Q0, Q1, P2), which can be represented by two triangles (P0, Q1, P2) and (P0, Q0, Q1).

This can be implemented in code as follows:

First, we find the vertices of the resulting clipped triangle or convex quadrilateral. Then, we delete the original triangle, but this time we pass the parameter false, which indicates that vertices should not be removed after deleting the face. This is important because some of the vertices of the deleted triangle may be used to create a new clipped triangle (or quadrilateral). For example, some vertices may be in the Positive zone and should not be removed, while others in the Negative zone must be removed. To determine which vertices need to be deleted, we initially store all the vertices of the triangle, and at the very end, we check which of these vertices are isolated. If they are isolated, we then delete them. If there are three vertices, we simply add a new triangle, If there are four vertices, we add two new triangles that form a convex quadrilateral.

Attention, further there will be some inefficient C++ code, but it’s done intentionally for the sake of simplicity in understanding. We’ll discuss how to improve it at the end of the post.

// status == Status::Mixed branch:

std::vector<OpenMesh::VertexHandle> vertices = generateCutVertices(meshData, fh, P, N);
        
std::vector<OpenMesh::VertexHandle> potentialIsolatedVertices;
for (auto vh : meshData.fv_range(fh))
    potentialIsolatedVertices.push_back(vh);

meshData.delete_face(fh, false);

if (std::size(vertices) == 3) 
    meshData.add_face(vertices[0], vertices[1], vertices[2]);

if (std::size(vertices) == 4) {
    meshData.add_face(vertices[0], vertices[1], vertices[2]);
    meshData.add_face(vertices[0], vertices[2], vertices[3]);
}       

for (auto vh : potentialIsolatedVertices)
    if (meshData.is_isolated(vh)) 
        meshData.delete_vertex(vh);

Let’s move on to the most interesting part, the generateCutVertices function, which generates cut vertices. This function uses the same loop as evalClipStatus, but now we will calculate the cut vertices within it. The result will be a vector of 3 or 4 vertices, as mentioned earlier, forming either a triangle or a convex quadrilateral. Some vertices will be newly calculated. The algorithm is quite straightforward: it involves iterating over all the edges of the triangle. If an edge is in the Negative zone, its vertices should not be added to the resulting primitive. If both vertices are in the Positive zone, then both vertices are added to the resulting primitive. If an edge is in the Mixed zone, a new vertex must be calculated based on the positions of the edge’s vertices. If the first vertex is in the Positive zone, then the new vertex should follow the first vertex. Conversely, if the second vertex is in the Positive zone, then the new vertex should precede it. Here’s how you could implement this in code:

std::vector<OpenMesh::VertexHandle> vertexHandles;

for (auto heh : meshData.fh_range(fh)) {
    OpenMesh::VertexHandle vh0 = heh.from();
    OpenMesh::VertexHandle vh1 = heh.to();

    auto P0 = meshData.point(vh0);
    auto P1 = meshData.point(vh1);

    double d0 = OpenMesh::dot(N, P0) - OpenMesh::dot(N, P);
    double d1 = OpenMesh::dot(N, P1) - OpenMesh::dot(N, P);
        
    // The edge is on the negative side, so we don't need vh0, vh1
    if (d0 <= 0 && d1 <= 0)
        continue;

    // We are on the positive side, so we simply add vh0, vh1 to the resulting primitive
    if (d0 > 0 && d1 > 0) {
        if (!std::ranges::contains(vertexHandles, vh0))
            vertexHandles.push_back(vh0);
                
        if (!std::ranges::contains(vertexHandles, vh1))
            vertexHandles.push_back(vh1);
        continue;
    }
    
    // Calculate a new cut vertex between points P0 and P1:
    auto newVh = calculateCutVertex(P0, P1, d0, d1);

    // Vertex vh0 is on the positive side, vh1 on the negative, so we add a cut vertex for vh1
    if (d0 > 0) {
        if (!std::ranges::contains(vertexHandles, vh0))
            vertexHandles.push_back(vh0);
            
        if (!std::ranges::contains(vertexHandles, newVh))
            vertexHandles.push_back(newVh);
    }
    
    // Vertex vh0 is on the negative side, so we add a cut vertex for vh0, vh1 is on the positive side
    if (d1 > 0) {
        if (!std::ranges::contains(vertexHandles, newVh))
            vertexHandles.push_back(newVh);
            
        if (!std::ranges::contains(vertexHandles, vh1))
            vertexHandles.push_back(vh1);
    }
}

Since two edges may contain the same vertex, it is necessary to perform a check using std::ranges::contains.

Finding cut vertex

The question remains how to find newVh.

In Figure , the point P0 is on the positive side of the plane and the point P1 is on the negative side. The point of intersection of the line segment and the plane occurs at Q. To be on the plane we need dot(N, point − P) = 0.

To be on the line segment, we need Q = (1 − t)P0 + tP1 for some t ∈ [0, 1].
Substituting Q into point in the plane equation and solving for t leads to:
t = (dot(N, P) – dot(N, P0)) / (dot(N, P1) – dot(N, P0)).
Knowing that:

d0 = dot(N, P0) – dot(N, P);
d1 = dot(N, P1) – dot(N, P);

t = -d0 / (d1 – d0)

We can further simplify this by multiplying the numerator and denominator by -1, which doesn’t change the value:

t = d0 / (d0 – d1)

The point of intersection and the new cut vertex is therefore Q = P0 + (d0/(d0 − d1))(P1 − P0).

OpenMesh::Vec3d calculateCutVertex(OpenMesh::Vec3d P0, OpenMesh::Vec3d P1, double d0, double d1) 
{
	double t = d0 / (d0 - d1);
	OpenMesh::Vec3d Q = P0 * (1.0 - t) + P1 * t;
	return Q;
}

Essentially, that’s almost everything. Putting our code together, we get:

std::vector<OpenMesh::VertexHandle> generateCutVertices(
	TriMeshKernel& meshData,
	OpenMesh::FaceHandle fh,
	OpenMesh::Vec3d orig,
	OpenMesh::Vec3d dir)
{
    std::vector<OpenMesh::VertexHandle> vertexHandles;

    for (auto heh : meshData.fh_range(fh)) {
        OpenMesh::VertexHandle vh0 = heh.from();
        OpenMesh::VertexHandle vh1 = heh.to();

        auto P0 = meshData.point(vh0);
	auto P1 = meshData.point(vh1);

	double d0 = OpenMesh::dot(dir, P0) - OpenMesh::dot(dir, orig);
	double d1 = OpenMesh::dot(dir, P1) - OpenMesh::dot(dir, orig);	

	if (d0 <= 0 && d1 <= 0)
	    continue;

	if (d0 > 0 && d1 > 0) {
	    if (!std::ranges::contains(vertexHandles, vh0))
	        vertexHandles.push_back(vh0);
				
	    if (!std::ranges::contains(vertexHandles, vh1))			    
                vertexHandles.push_back(vh1);
	    
            continue;
	}

	auto Q = calculateCutVertex(P0, P1, d0, d1);

	OpenMesh::VertexHandle newVh = meshData.add_vertex(Q);

	if (d0 > 0) {
	    if (!std::ranges::contains(vertexHandles, vh0))
		vertexHandles.push_back(vh0);
				
	    if (!std::ranges::contains(vertexHandles, newVh))
		vertexHandles.push_back(newVh);
	}

	if (d1 > 0) {
	    if (!std::ranges::contains(vertexHandles, newVh))				                
                vertexHandles.push_back(newVh);
				
	    if (!std::ranges::contains(vertexHandles, vh1))
		vertexHandles.push_back(vh1);
	}
    }

    return vertexHandles;
}

But this code contains a certain error related to the fact that one edge of the triangle may be shared between two triangles. Therefore, it’s unnecessary to create a new cut vertex every time; it’s sufficient to reuse the one already created. To do this, it’s enough to establish a simple cache as follows:

std::map<OpenMesh::EdgeHandle, OpenMesh::VertexHandle> cache;

In this way, the code inside the loop changes accordingly:

OpenMesh::VertexHandle newVh;
if (auto it = cache.find(eh); it != std::end(cache))
	newVh = it->second;
else {
	newVh = meshData.add_vertex(Q);
	cache[eh] = newVh;
}

Small improvements

1. If you want to cut the mesh into two halves but don’t want to remove triangles lying in the Negative zone, you just need to handle the Mixed zone and modify the generateCutVertices function so that Positive and Negative zones have their own caches. Additionally, to avoid making significant changes to the code inside the function, you can employ the following trick: add a parameter called sign to determine which zone to consider as Positive or Negative. Then, you can multiply “sign” by the expressions inside each “if” statement within the function. This way, you can treat Positive as Negative zone and vice versa.

if (sign * d0 <= 0 && sign * d1 <= 0)
	continue;

if (sign * d0 > 0 && sign * d1 > 0) {
	if (!std::ranges::contains(vertexHandles, vh0))
		vertexHandles.push_back(vh0);

	if (!std::ranges::contains(vertexHandles, vh1))
		vertexHandles.push_back(vh1);
	continue;
}

auto Q = calculateCutVertex(P0, P1, d0, d1);

OpenMesh::VertexHandle newVh;
if (auto it = cache.find(eh); it != std::end(cache))
	newVh = it->second;
else {
	newVh = meshData.add_vertex(Q);
	cache[eh] = newVh;
}

if (sign * d0 > 0) {
	if (!std::ranges::contains(vertexHandles, vh0))
		vertexHandles.push_back(vh0);

	if (!std::ranges::contains(vertexHandles, newVh))
		vertexHandles.push_back(newVh);
}

if (sign * d1 > 0) {
	if (!std::ranges::contains(vertexHandles, newVh))
		vertexHandles.push_back(newVh);

	if (!std::ranges::contains(vertexHandles, vh1))
		vertexHandles.push_back(vh1);
}

And then calling with two different caches for both sides:

std::vector<OpenMesh::VertexHandle> verticesLeft = generateCutVertices(meshData, fh, cache, P, N, +1.0);

std::vector<OpenMesh::VertexHandle> verticesRight = generateCutVertices(meshData, fh, cache2, P, N, -1.0);

2. Choosing the “absolute zero” of the plane can be tricky. You can adjust with some epsilon and constrain d0, d1 to this epsilon:

if (std::abs(d0) <= eps)
	d0 = 0.0;

if (std::abs(d1) <= eps)
	d1 = 0.0;

This approach helps to avoid issues related to small values and the generation of duplicate vertices.

3. Returning a std::vector of vertices is highly inefficient. Passing through arguments is also not very elegant. In such cases, a data structure like inplace_vector vertices would be more suitable. It allows for efficient addition of vertices without dynamic memory allocation. Since we know that the result of clipping will yield either 3 or 4 vertices, such a structure is ideal. Boost offers such a structure.

Another approach is use pmr. In order to do that just return std::pmr::vector<OpenMesh::VertexHandle> and pass a buffer via std::pmr::memory_resource allocated on the stack into the function. This is also efficient, but for such a simple task, it’s a bit of an overkill.

std::pmr::vector<OpenMesh::VertexHandle> generateCutVertices(/*args as before,*/ std::pmr::memory_resource* resource)
{
	std::pmr::vector<OpenMesh::VertexHandle> vertexHandles(resource);
	
	...
	
	return vertexHandles;
}

// somewhere in the beginning of clip function:
char buffer[1024] = {};
std::pmr::monotonic_buffer_resource pool{ std::data(buffer), std::size(buffer) };

// call 
std::pmr::vector<OpenMesh::VertexHandle> vertices = generateCutVertices(meshData, fh, cache, P, N, 1, &pool);

Using garbage_collection

Don’t forget to use garbage_collection function in OpenMesh. It is a function provided that efficiently removes unused elements from a mesh, such as vertices, edges, and faces, and reclaims memory. It makes sense especially after using Clip operation when many faces could be deleted.

Using update_normal

update_normal is needed whenever there is a change in the mesh geometry or topology that affects the orientation or direction of normals. There is a place where you want to call it as meshData.update_normal(fh). It’s right one after you adding new cut face.

Conclusion

The detailed clipping algorithm effectively handles complex mesh manipulations, using the OpenMesh library to perform precise geometric calculations. This approach ensures efficient processing and provides the capability to handle various mesh configurations in 3D modeling, game development, and more.