Best unique subdivision algorithm

Hello OpenGL forum

This is my first topic and message here, hope won’t disappoint you :slight_smile:

I’m working with native OpenGL/WebGL context ( by this I meant not using any 3rd party libraries ).
And I have achieved the step where I need some useful subdivision algorithm.

Because of fact that I’m really new to this subject, I don’t know what algorithm is best. I’ve searched some information, there are different topics like:

  • loop subdivision surface scheme ( which was developed by Charles Loop, 1987 )
  • Delaunay triangulation
  • Voronoi diagram
  • Tessellation algorithm

and etc…

I’m kinda surprised and don’t even know what subdivision algorithm is better to learn previously. So, I’m asking you to give me such a piece of advice.

Does exist some simple and unique subdivision algorithm?

I’m preparing my own OpenGL framework and there are some primitives like cube, triangle, pyramid, sphere, cone, torus and some kind of custom geometry
( but with no extrusion ).

I suppose that for custom geometry the best option will be some variant of Delaunay triangulation for 3D.

But what about primitives?
Maybe there is a clean, great, useful solution for this?


Delaunay triangulation and Voronoi diagrams aren’t subdivision algorithms. They’re ways of partitioning a plane based upon a set of points.

For subdividing arbitrary geometry, the main factors distinguishing subdivision schemes are the types of geometry to which they can be applied (e.g. arbitrary polygons or just triangles), the degrees of continuity provided, whether they offer additional parameters (e.g. variable edge sharpness), and whether they interpolate or approximate the original vertices.

For primitives (e.g. sphere, cylinder, torus, etc), you’d typically just choose a topology then calculate the vertex coordinates using either a parametric equation or an implicit equation.

OpenGL 4 tessellation shaders favour schemes which sample a surface defined by an equation (e.g. Bezier patch, Bezier triangle, NURBS surface, etc) over recursive subdivision (Loop, Catmull-Clark, Doo-Sabin, etc).

The recursive subdivision itself is fairly straightforward, and you have parts of it already. Taking one triangle of the octahedron, we seed the recursive calls by passing the coordinates of the full triangle to the recursive function, together with the desired subdivision level (something in the range of 3 to 5 will give reasonable spheres):

subdivide(1.0f, 0.0f, 0.0f,
0.0f, 1.0f, 0.0f,
0.0f, 0.0f, 1.0f,

Then, in the recursive function, we first check if we reached level 0, in which case we have a final triangle. Otherwise we split the triangle into 4, and make the recursive call for each of them. The vertices of the 4 sub-triangles are formed from the original triangle vertices and the middle of their edges, based on the diagrams you already found:

 /  \
/    \

/ \ /
/ /

The recursive function then looks like this:

void subdivide(float v1x, float v1y, float v1z,
float v2x, float v2y, float v2z,
float v3x, float v3y, float v3z,
int level) {
if (level == 0) {
// Reached desired tessellation level, emit triangle.
drawTriangle(v1x, v1y, v1z,
v2x, v2y, v2z,
v3x, v3y, v3z);
} else {
// Calculate middle of first edge…
float v12x = 0.5f * (v1x + v2x);
float v12y = 0.5f * (v1y + v2y);
float v12z = 0.5f * (v1z + v2z);
// … and renormalize it to get a point on the sphere.
float s = 1.0f / sqrt(v12x * v12x + v12y * v12y + v12z * v12z);
v12x *= s;
v12y *= s;
v12z *= s;

    // Same thing for the middle of the other two edges.
    float v13x = 0.5f * (v1x + v3x);
    float v13y = 0.5f * (v1y + v3y);
    float v13z = 0.5f * (v1z + v3z);
    float s = 1.0f / sqrt(v13x * v13x + v13y * v13y + v13z * v13z);
    v13x *= s;
    v13y *= s;
    v13z *= s;

    float v23x = 0.5f * (v2x + v3x);
    float v23y = 0.5f * (v2y + v3y);
    float v23z = 0.5f * (v2z + v3z);
    float s = 1.0f / sqrt(v23x * v23x + v23y * v23y + v23z * v23z);
    v23x *= s;
    v23y *= s;
    v23z *= s;

    // Make the recursive calls.
    subdivide(v1x, v1y, v1z,
              v12x, v12y, v12z,
              v13x, v13y, v13z,
              level - 1);
    subdivide(v12x, v12y, v12z,
              v2x, v2y, v2z,
              v23x, v23y, v23z,
              level - 1);
    subdivide(v13x, v13y, v13z,
              v23x, v23y, v23z,
              v3x, v3y, v3z,
              level - 1);
    subdivide(v12x, v12y, v12z,
              v23x, v23y, v23z,
              v13x, v13y, v13z,
              level - 1);


This gives you all the vertices you need. But it calculates each vertex of the tessellation multiple times. Where it gets slightly painful is if you want to calculate each vertex only once, and store them in a vertex buffer. To do that, you’ll have to also pass indices to the recursive function, and do some math on those indices so that you can place each resulting vertex at a unique index. It gets even more cumbersome if you want nice triangle strips from all the resulting unique vertices.

Unfortunately all those details are somewhat beyond the scope of an answer here. But I hope that the above will at least clearly explain the underlying math and logic needed to do the actual subdivision.