Implementation of a Robust (i.e. with a finite plane thickness) Sutherland–Hodgman algorithm for clipping polygons against an axis-aligned bounding box. I use the following code only for clipping triangles. So the clipped polygon can consist only of at most 9 vertices. The third-party code of Point
, BBox
and Lerp
is not given, since this is rather trivial.
The C++ code works fine (Python version available at) and fast for clipping some triangles. Unfortunately, I typically need to perform 1 billion of these operations. Profiling learns most of the time is spent on these clipping operations (which is not the core of my apllication). So my question: is it possible to make the code more performance efficient (e.g. eliminating most of the branching)?
#pragma once
#define PLANE_THICKNESS_EPSILON 0.00001f
// Classify a given vertex against an axis-aligned plane
//
// @param sign min/max clipping plane
// @param axis axis of the clipping plane
// @param c_v one vertex of the clipping plane
// @param p_v vertex to classify
//
// @return classification of the vertex
template <typename T>
inline int8_t Classify(int8_t sign, uint8_t axis, const Point3<T> &c_v, const Point3d &p_v) {
const double d = sign * (p_v[axis] - c_v[axis]);
if (d > PLANE_THICKNESS_EPSILON) return 1;
else if (d < -PLANE_THICKNESS_EPSILON) return -1;
else return 0;
}
#define POINT_BUFFER_SIZE 9
// Clip the given polygon against an axis-aligned plane
//
// @param p_vs polygon before clipping as a sequence of vertices
// @param nb_p_vs number of vertices before clipping
// @param sign min/max clipping plane
// @param axis axis of the clipping plane
// @param c_v one vertex of the clipping plane
//
// @return p_vs polygon after clipping as a sequence of vertices
// @return nb_p_vs number of vertices after clipping
template <typename T>
inline void Clip3D_plane(Point3d *p_vs, uint8_t *nb_p_vs, int8_t sign, uint8_t axis, const Point3<T> &c_v) {
uint8_t nb = (*nb_p_vs);
if (nb == 0) return;
else if (nb == 1) {
*nb_p_vs = 0;
return;
}
Point3d new_p_vs[POINT_BUFFER_SIZE];
uint8_t k = 0;
bool b = true; // polygon is fully located on clipping plane
Point3d p_v1 = p_vs[nb-1];
int8_t d1 = Classify<T>(sign, axis, c_v, p_v1);
for (uint8_t j = 0; j < nb; ++j) {
const Point3d &p_v2 = p_vs[j];
int8_t d2 = Classify<T>(sign, axis, c_v, p_v2);
if (d2 < 0) {
b = false;
if (d1 > 0) {
const double alpha = (p_v2[axis] - c_v[axis]) / (p_v2[axis] - p_v1[axis]);
new_p_vs[k++] = Lerp(alpha, p_v2, p_v1);
}
else if (d1 == 0 && (k == 0 || new_p_vs[k-1] != p_v1))
new_p_vs[k++] = p_v1;
}
else if (d2 > 0) {
b = false;
if (d1 < 0) {
const double alpha = (p_v2[axis] - c_v[axis]) / (p_v2[axis] - p_v1[axis]);
new_p_vs[k++] = Lerp(alpha, p_v2, p_v1);
}
else if (d1 == 0 && (k == 0 || new_p_vs[k-1] != p_v1))
new_p_vs[k++] = p_v1;
new_p_vs[k++] = p_v2;
}
else {
if (d1 != 0)
new_p_vs[k++] = p_v2;
}
p_v1 = p_v2;
d1 = d2;
}
if (b) return;
*nb_p_vs = k;
for (uint8_t j = 0; j < k; ++j)
p_vs[j] = new_p_vs[j];
}
// Clip the given polygon against an axis-aligned bounding box
//
// @param p_vs polygon before clipping as a sequence of vertices
// @param nb_p_vs number of vertices before clipping
// @param clipper axis-aligned bounding box used for clipping
//
// @return p_vs polygon after clipping as a sequence of vertices
// @return nb_p_vs number of vertices after clipping
inline void Clip3D_AABB(Point3d *p_vs, uint8_t *nb_p_vs, const BBox &clipper) {
for (uint8_t axis = 0; axis < 3; ++axis) {
Clip3D_plane<float>(p_vs, nb_p_vs, 1.0, axis, clipper.pMin);
Clip3D_plane<float>(p_vs, nb_p_vs, -1.0, axis, clipper.pMax);
}
}
All remaining legacy declarations and definitions (stripped as much as possible) to make it compile:
#include "stdint.h"
template <typename T>
class Point3 {
public:
Point3() { x = y = z = 0; }
Point3(T x, T y, T z) : x(x), y(y), z(z) {}
template <typename U>
explicit Point3(const Point3<U> &p) : x((T)p.x), y((T)p.y), z((T)p.z) {}
Point3<T> operator+(const Point3<T> &p) const {
return Point3<T>(x + p.x, y + p.y, z + p.z);
}
template <typename U>
Point3<T> operator*(U f) const {
return Point3<T>(f * x, f * y, f * z);
}
T operator[](int i) const {
if (i == 0) return x;
if (i == 1) return y;
return z;
}
T &operator[](int i) {
if (i == 0) return x;
if (i == 1) return y;
return z;
}
bool operator==(const Point3<T> &p) const {
return x == p.x && y == p.y && z == p.z;
}
bool operator!=(const Point3<T> &p) const {
return x != p.x || y != p.y || z != p.z;
}
T x, y, z;
};
typedef Point3<float> Point;
typedef Point3<double> Point3d;
template <typename T, typename U>
inline Point3<T> operator*(U f, const Point3<T> &p) {
return p * f;
}
template <typename T>
Point3<T> Lerp(double t, const Point3<T> &p0, const Point3<T> &p1) {
return (1 - t) * p0 + t * p1;
}
class BBox {
public:
BBox(const Point &pMin, const Point &pMax) : pMin(pMin), pMax(pMax) {}
Point pMin, pMax;
};