# Tesselate triangle

Hi!
Say, I want to increase detalization of a mesh to make a model smoother.

I want to do this by a simple iterative approach: if a certain edge with vertices Va and Vb and normals Na and Nb has an angle between those normals higher than X degrees, I break that edge on two by inserting an extra vertex Vm. It is required that normal Nm at this vertex is equal: Nm=normalize(Na+Nb). So how do I determine the coordinates of that extra vertex Vm?

Should be a trivial task, but I had a hard time searching for a solution…

However you like.

One option is to take the midpoint of the edge of a PN triangle, which would be

Vm = (Va+Vb)/2 - (Wa*Na + Wb*Nb)/8

where

Wa = (Vb - Va) . Na
Wb = (Va - Vb) . Nb

Oh, thank you, GClements, it works!
Care to explain why the displacement divided specifically by 8?

Whoops, I just realized I have incorrectly specified the problem. We can’t say that the normal Nm in the middle point Vm will be the average between normals Na and Nb. Attachment is an example were both normals are pointing the same direction but the Nm shall be perpendicular to them at Vm.
I am getting lost here… I thought I could just divide edge few times and get a smooth curvature as a result. Any suggestions how to do that correctly?

Your title is “Tesselate triangle”. But the figure you included is not a triangle.

It is an edge of the triangle. I tessellate the triangle by subdividing it’s edge.

A PN (point-normal) triangle is a cubic Bezier triangle whose control points are derived from 3 vertices and their associated normals.

A more detailed description can be found e.g. here, but a brief summary follows.

The control points are given by

``````
b300 = P0
b030 = P1
b003 = P2
b210 = (2*P0 + P1 - w01*N0)/3
b120 = (2*P1 + P0 - w10*N1)/3
b021 = (2*P1 + P2 - w12*N1)/3
b012 = (2*P2 + P1 - w21*N2)/3
b102 = (2*P2 + P0 - w20*N2)/3
b201 = (2*P0 + P2 - w02*N0)/3
b111 = (b210 + b120 + b021 + b012 + b102 + b201)/4 - (b300 + b030 + b003)/6

``````

where P0, P1, P2 are the vertices, N0, N1, N2 are the normals, and w[sub]ij[/sub]=dot(P[sub]j[/sub]-P[sub]i[/sub], N[sub]i[/sub])

Each of the 6 edge control points is generated by interpolating the edge in a 2:1 ratio then projecting onto the tangent plane defined by the closest vertex and its associated normal. This ensures that the normals are actually normal to the surface at the vertex. The interior control point b111 is calculated as E+(E-V)/2 where E is the mean of the 6 edge control points and V is the mean of the 3 vertices.

Given the control points and barycentric coordinates (a,b,c), the corresponding point in the surface is calculated by:

``````
P = b300 * c^3
+ b030 * a^3
+ b003 * b^3
+ b210 * 3 * a * c^2
+ b120 * 3 * c * a^2
+ b201 * 3 * b * c^2
+ b021 * 3 * b * a^2
+ b102 * 3 * c * b^2
+ b012 * 3 * a * b^2
+ b111 * 6 * a * b * c

``````

So the midpoint of the edge between P1 and P2 is found by evaluating the above with a=1/2, b=1/2, c=0.

Tangents can be found by substituting c=1-a-b into the above and differentiating with respect to a and b; their cross product gives the surface normal at any point. The expression for each tangent is quadratic, so that for the normal is quartic. But the normal is more commonly approximated by a quadratic Bezier triangle with control points

``````
n200 = N0
n020 = N1
n002 = N2
n110 = N0+N1-v01*(P1-P0)
n011 = N1+N2-v12*(P2-P1)
n101 = N2+N0-v20*(P0-P2)

``````

where v[sub]ij[/sub] = 2*dot(P[sub]j[/sub]-P[sub]i[/sub], N[sub]i[/sub]+N[sub]j[/sub])/dot(P[sub]j[/sub]-P[sub]i[/sub], P[sub]j[/sub]-P[sub]i[/sub])

For given barycentric coordinates, the normal is calculated by:

``````
N = n200 * c^2
+ n020 * a^2
+ n002 * b^2
+ n110 * c * a
+ n011 * a * b
+ n101 * c * b

``````

Note that cubic Bezier triangles typically don’t have C[sub]1[/sub] continuity across edges. To achieve that, you would need a quintic surface (21 control points), which is significantly more expensive both to determine and evaluate.

Most of the commonly-used subdivision surfaces which do enforce C[sub]1[/sub] continuity across edges (e.g. Catmull-Clark, Loop) only approximate the initial vertices (i.e. each subdivision step moves the original vertices as well as creating additional vertices).

Wow, wow, GClements, you always shoot so much theoretic material for my simple mind that I start feel myself like a black guy in a science lab…
I’ve seen this triangle subdivision with extra 7 points and decided that it is pretty overkill for me. I would like to go with one edge with two corner points subdividing it with 1 extra vertex, and that is it.

So here is what I came up with. We deal with a curve: it starts at one point and exit through another. At both points we know the normals to that curve. So first let’s find a tangent vectors (Ta, Tb) directing the curve at each of those two corner points (Va, Vb). To do this we project the vector VaVb onto planes perpendicular to our normals (Na, Nb):
Ta = VaVb - Na * dot(VaVb, Na);
Tb = VbVa - Nb * dot(VbVa, Nb);
Now let’s add a weight coefficients (wa, wb=1-wa) that show how close a certain point on the curve to a corner point: for a point close to Va the wa will be close to 1 and wb is close to 0, and vice-versa. To make the smooth transition we do a “double” blending: if we close to Va, than we want the curve to follow the direction of Ta and the position should be close to Va, and the same with a proximity with Vb. So here is a blending formula I came up with:

Vm = wa(Va+(1-wa)Ta) + wb(Vb+(1-wb)Tb);

If we simplify this for wa=0.5, we get:

Vm = (Va+Vb)/2 - (dot(Na,(Vb-Va))Na + dot(Nb,(Va-Vb))Nb)/4

Which is the same as the formula at the second post except that the subtracted part is divided by 4, not 8. Shit…

I think I found the solution here. As I understood from the video, the “tension line” is tangent to the curve, and the coordinates of the desired point Vm is a center of that line at 50% “completion”. Therefore I came up with this GLSL-style C++ algorithm:

``````#include "GLSL.hpp"  //GLSL vectors, matrices, operators and functions

//Function breaks an edge onto two by inserting an extra point inbetween the
//specified corner points
void TessellateEdge(const vec3& Va, const vec3& Na,
const vec3& Vb, const vec3& Nb,
vec3* Vm, vec3* Nm){
//Calculate tangents
vec3 Ta = Vb-Va-Na*dot((Vb-Va),Na),
Tb = Va-Vb-Nb*dot((Va-Vb),Nb),
Tm = 0.5*(Vb-Va)+0.25*(Tb-Ta);
//Calculate coordinates of the middle point
*Vm = 0.5*(Va+Vb) + 0.375*(Ta+Tb);
//Calculate normal at the middle point
vec3 perp = cross((Vb-Va),(Na+Nb));
*Nm = normalize(cross(perp,Tm));
}
``````

I double-checked all the formulas, run this in few setups and concluded that it provides the desired result…

One way to look at it:
Wa*Na is the vector from Vb to the closest point on the tangent plane passing through Va with normal Na. Halving that (Wa*Na/2) gives the vector from the midpoint of the edge to that tangent plane.

Similarly, Wb*Nb is the vector from Va to the closest point on the tangent plane passing through Vb with normal Nb. Halving that (Wb*Nb/2) gives the vector from the midpoint of the edge to that tangent plane.

The average of those two is (Wa*Na+Wb*Nb)/4. If we offset the midpoint of the edge, (Va+Vb)/2, by this vector, we get a third point Vc. The points Va, Vc and Vb form a triangle. One edge is the base, the other edges are approximately tangential to the surface at the corresponding endpoint.

If these three points were the control points of a quadratic Bezier curve, the curve would be a parabola which would bisect the line between (Va+Vb)/2 and Vc exactly in half.

To see why this might be the case, consider the curve y=1-x^2. This has roots at x=-1 and x=1 and a maximum of y=1 at x=0. The derivative dy/dx=-2*x which is 2 at x=-1 and -2 at x=1. Thus tangents at the roots intersect at x=0, y=2 to give the interior control point. This gives a control polygon with the vertices (-1,0), (0,2), (1,0) while the curve itself passes through (-1,0), (1,0), (0,0).

Thus to find the midpoint of the curve, we halve the offset vector, which gives us a denominator of 8.

So we halve once because we want the offset of the midpoint rather than the endpoint, once when calculating the average of two midpoints, and once due to the fact that the height of the peak of a parabola is half the height at which its tangents intersect.

Oh yes, now I see, thank you! Funny, but, the more we talk about it the more precise the “However you like.” answer looks like to me… There is clearly a whole bunch of solutions to make a curve bend in different ways and still pass through two PN-handles.

For those who might care - an update:

``````#include "GLSL.hpp"  //GLSL vectors, matrices, operators and functions

//Function breaks an edge onto two by inserting an extra point inbetween the
//specified corner points; false returned if no solution can be found
bool TessellateEdge(const vec3& Va, const vec3& Na,
const vec3& Vb, const vec3& Nb,
vec3* Vm, vec3* Nm){
//Precalculate some values, check for errors
vec3 AB = Vb-Va; //Direction from Va to Vb
if(dot(AB,AB)==0.f || dot(cross(AB,Na),cross(AB,Nb))<=0.f)return false;
//Calculate curve' tangents
vec3 Ta = 0.5f*(AB-Na*dot(AB,Na)),
Tb = 0.5f*(Nb*dot(AB,Nb)-AB),
Tm = AB+0.5f*(Tb-Ta);
//Calculate coordinates of the middle point
if(Vm)*Vm = 0.5*(Va+Vb) + 0.375*(Ta+Tb);
//Calculate normal at the middle point by rotating Na according to Ta->Tm rotation
if(Nm)*Nm = normalize(cross(cross(Ta,Tm),Na)+Na*dot(Ta,Tm));
return true;
}
``````

In the attachment is a test executable for that iterative tessellation algorithm (requires OpenGL4.5). Keys:
Esc - exit;
Enter - fullscreen;
Mouse, Q,E - rotate camera;
R - reset mesh to an initial triangle;
T - tessellate once (break one edge resulting in 1 or 2 additional triangles added).

The algorithm tessellates 1 edge each time T pressed. The edge picked is the one with highest curvature which is defined as
float Curvature = abs(dot( Va-Vb, Na )) + abs(dot( Va-Vb, Nb ));
where Va, Vb are vertices and Na, Nb are normals at Va and Vb.

It does not seem to be perfect: I think the resulting shape slightly dependent on the order at which edges were tessellated, producing wavy shape instead of desired smooth shape. Moreover, after few hundred tessellations (hold T for a few sec) it just starts tessellating a tiny spot over and over again making an ugly kink.

So now I am not sure if the iterative approach can actually produce a smooth shape… Any thoughts?

I suspect that the issue may be related to your calculation for the normal. If you calculate the normal exactly, I suspect that a PN triangle constructed using the interpolated point and normal should actually be one half of the original PN triangle.

However, that requires taking the third point of the triangle into account, which in turn means that the interpolated normal will differ between the two triangles which share the edge. As I’ve mentioned before, you can’t ensure C[sub]1[/sub] continuity with quadratic Bezier triangles (nor with anything less than quintic, for that matter).

Ultimately, if there was a simple subdivision scheme which converged to a smooth surface and which interpolated (rather than approximated) the original vertices, everyone would be using it already. Most of the common approaches merely approximate the original vertices, because it’s much easier to obtain convergence that way.

Say, we have a curve, which shape is defined by 2 PN control points. The shape of that curve actually differ depending on the length we choose for the tangents - we can scale them to any length producing different curves which are still passing through the same PN-handles, right?

The second problem is a subdivision deviation. Each time the curve divided, a new PN-handle created. And subsequent subdivisions use that new handle. What are the chances that all those new PN-handles all lay on the original curve? If we use the original handles to calculate all subdivision locations we explicitly make them located on the original curve.

But what is required to make sure the iterative subdivision produce PN-handles belonging to the original curve? Is it possible at all?
I am sure it is if our curve lays on the sphere or any other predetermined shape. The sphere is not the best example, though, because it can pass through two-PN handles only if their normals lay in the same plane and their average is perpendicular to the line between handles.

But what about ellipse? I think it can be defined by two PN-handles (unless normals are parallel), isn’t it?

Well, there is a way to subdivide a surface using iterative approach and still get a desired smooth shape. The trick is to use the initial control points whenever a new vertex added. But in order to do that each vertex must store an additional attribute (I call it W) showing it’s topological location on the surface.

I deal with triangular patches, therefore I have 3 PN-handles and the input patch is a triangle with vertices A,B,C. The W attribute for those vertices are:
A.W = vec3(1,0,0);
B.W = vec3(0,1,0);
C.W = vec3(0,0,1);
Whenever an edge V[SUB]1[/SUB]-V[SUB]2[/SUB] is subdivided, the new vertex Vm has it’s desired topological location Wm calculated as:
Wm = 0.5*(W[SUB]1[/SUB] + W[SUB]2[/SUB]);
and based on the resulting Wm all other vertex attributes (position, normal) are calculated.
So here is the updated code:

``````
//Tesselate curve AB at k percent calculating the position (Vm) and
//normal (Nm) of M;
//k is an interpolation percent (k=0 -> M=A, k=1 -> M=B);
//the Curvature controls the curve' extrusion: 0 for flat interpolation
//and >0 for a smooth extrusion (0.5 is a recommended value);
//false returned on error
bool TessellateCurve(const vec3& Va, const vec3& Na,
const vec3& Vb, const vec3& Nb,
vec3* Vm, vec3* Nm,
float k, float Curvature=0.5f){
//Check for input errors
if(k<0.f||k>1.f||Curvature<=0.f||Curvature>=1.f)return false;
//Precalculate some values, check for errors
vec3 AB = Vb-Va; //Direction from Va to Vb
if(dot(AB,AB)==0.f || dot(cross(AB,Na),cross(AB,Nb))<=0.f)return false;
//Calculate curve' tangents
vec3 Ta = Curvature*(AB-Na*dot(AB,Na)),
Tb = Curvature*(Nb*dot(AB,Nb)-AB),
DF = AB+Tb-Ta,
J = k*((2.f-k)*DF - k*Tb) + Va + Ta,
H = k*(k*DF + (2.f-k)*Ta) + Va,
Tm = J-H;
//Calculate coordinates of the middle point
if(Vm)*Vm = H*(1.f-k)+J*k;
//Calculate normal at the middle point
if(Nm)*Nm = normalize(cross(cross(Ta,Tm),Na)+Na*dot(Ta,Tm));
return true;
}
``````
``````
//Evaluate point on the bezier surface defined by 3 vertices;
//Wm is a topological location of the vertex M:
//Wm={1,0,0} -> M=A,
//Wm={0,1,0} -> M=B,
//Wm={0,0,1} -> M=C,
//the Curvature controls the curve' extrusion: 0 for flat interpolation
//and >0 for a smooth extrusion (0.5 is a recommended value);
//false returned on error
bool Evaluate(const vec3& Va, const vec3& Na,
const vec3& Vb, const vec3& Nb,
const vec3& Vc, const vec3& Nc,
const vec3& Wm,
vec3* Vm, vec3* Nm,
float Curvature=0.5f){
//Check for input errors
if(Wm.x<0.f||Wm.y<0.f||Wm.z<0.f)return false;
float Wabc = Wm.x+Wm.y+Wm.z; //Must be 1.f
if(Wabc==0.f)return false;
Wm /= Wabc; //Ensure the sum is equal to 1.f
//Calculate weight sum pairs
float Wab = Wm.x+Wm.y,
Wbc = Wm.y+Wm.z,
Wac = Wm.x+Wm.z;
//Prepare the containers for the intermediate points
vec3 Vab = Vc, Nab = Nc,
Vbc = Va, Nbc = Na,
Vac = Vb, Nac = Nb;
//Calculate the intermediate points
if(Wab>0.f)if(!TessellateCurve(Va,Na, Vb,Nb, &Vab,&Nab, Wm.y/Wab, Curvature))return false;
if(Wbc>0.f)if(!TessellateCurve(Vb,Nb, Vc,Nc, &Vbc,&Nbc, Wm.z/Wbc, Curvature))return false;
if(Wac>0.f)if(!TessellateCurve(Vc,Nc, Va,Na, &Vac,&Nac, Wm.x/Wac, Curvature))return false;
//Calculate the evaluated point 3 different ways and approximate the result
vec3 Vm1, Nm1,
Vm2, Nm2,
Vm3, Nm3;
if(!TessellateCurve(Va,Na, Vbc,Nbc, &Vm1,&Nm1, Wbc, Curvature))return false;
if(!TessellateCurve(Vb,Nb, Vac,Nac, &Vm2,&Nm2, Wac, Curvature))return false;
if(!TessellateCurve(Vc,Nc, Vab,Nab, &Vm3,&Nm3, Wab, Curvature))return false;
Vm = (Vm1+Vm2+Vm3)/3.f;
Vn = normalize(Nm1+Nm2+Nm3);
return true;
}

``````

Here is the result in the attachment showing the resulting iterative subdivision process step-by-step (requires OpenGL4.5). Keys:
Esc - exit;
Enter - fullscreen;
Mouse, Q,E - rotate camera;
R - reinitialize mesh;
T - tessellate once (break one edge resulting in 1 or 2 additional triangles added).

I just confirmed that the function TessellateCurve incorrectly computes the normal Nm for the subdivision point Vm. But I have no idea how to do that the right way… If I simply interpolate it using k, then there is no guarantee that the resulting normal will be perpendicular to tangent Tm at the point Vm. I can force it to be so, but not always. Any suggestions?

The surface is a cubic function in barycentric coordinates (a,b,c). Barycentric coordinates have the constraint a+b+c=1 meaning c=1-a-b, so you only actually have two independent variables, a and b. Calculating partial derivatives with respect to a and b gives you the tangents as cubic functions in a and b (more precisely, the derivative w.r.t. a will have terms up to a2 and b3 while the derivative w.r.t. b will have terms up to a3 and b2). The normal is the cross product of the tangents, which is a quinitic function in a and b.

But PN triangles don’t normally use the correct normal. A quadratic approximation is both cheaper to compute and (more importantly) can be used to hide the fact that you don’t have C[sub]1[/sub] continuity at the edges (i.e. there’s a sharp corner at the edge between triangles). If you derive the equations for the tangents, you’ll note that the tangent perpendicular to the edge between P0,N0 and P1,N1 depends not only on those points but also on the opposing point P2,N2. Similarly for the other two edges. Which means that when you have two triangles sharing a common edge, the tangents on either side won’t be co-planar and thus the normals won’t be co-linear.

Using an approximation where the tangents along an edge depend only upon the two endpoints of that edge (with their associated normals) ensures that the tangents are the same on either side, meaning that the normal is the same on either side, meaning that there’s no sudden discontinuity in the lighting at the edge, meaning that the surface appears smooth.

You can find the usual approximation for the normal in one of my earlier posts in this thread.

Well, I do need continuity across edges so two patches sharing two PN-handles can be stacked together without a visible seem. And my approach makes it so - I just checked it with two patches - ideally smooth transition both in shape and normal. The reason for that is that the edge points are independent from the opposite handle. The only problem is the calculation of Nm inside “TessellateCurve” function - I had to change it to linear approximation:

if(Nm)*Nm = normalize((1-k)Na+kNb);

Only because in the example code above it is calculated clearly the wrong way - even for k=1 the resulting Nm does not correspond to Nb while it should.

Therefore the last problem I need to solve, again, is the very first one: subdividing a PN-defined 3D cubic bezier curve calculating a new PN handle at the percentage k.

Success! In order to calculate the normal I use the second service curve created by shifting the control points of the original one. The resulting normal is calculated as a cross-product of a tangent at Vm and the vector (U-Vm), where U is a point on a service curve for the same k value. Here is the modified code of TessellateCurve function:

``````
//Tesselate curve AB at k percent calculating the position (Vm) and
//normal (Nm) of M;
//k is an interpolation percent (k=0 -> M=A, k=1 -> M=B);
//the Curvature controls the curve' extrusion: 0 for flat interpolation
//and >0 for a smooth extrusion (0.5 is a recommended value);
//false returned on error
bool TessellateCurve(const vec3& Va, const vec3& Na,
const vec3& Vb, const vec3& Nb,
vec3* Vm, vec3* Nm,
float k, float Curvature=0.5f){
//Check for input errors
if(k<0.f||k>1.f||Curvature<=0.f||Curvature>=1.f)return false;
//Precalculate some values, check for errors
vec3 AB = Vb-Va; //Direction from Va to Vb
if(dot(AB,AB)==0.f || dot(cross(AB,Na),cross(AB,Nb))<=0.f)return false;
//Calculate curve' tangents
vec3 Ta = Curvature*(AB-Na*dot(AB,Na)),
Tb = Curvature*(Nb*dot(AB,Nb)-AB),
DF = AB+Tb-Ta,
J = k*((2.f-k)*DF - k*Tb) + Va + Ta,
H = k*(k*DF + (2.f-k)*Ta) + Va,
Tm = J-H;
//Calculate coordinates of the middle point
if(Vm)*Vm = H*(1.f-k)+J*k;
//Calculate normal at the middle point
vec3 P = Va+normalize(cross(Na,AB)),
Q = Vb+normalize(cross(Nb,AB)),
PQ = Q-P,
RS = PQ+Tb-Ta,
L = k*((2.f-k)*RS - k*Tb) + P + Ta,
W = k*(k*RS + (2.f-k)*Ta) + P,
U = W*(1.f-k)+L*k;
if(Nm)*Nm = normalize(cross( Tm, U - (*Vm) ));
return true;
}
``````

And here is proof in the attachment with two triangular patches sharing 2 handles. Both patches are tessellated (press T) independently but form a continuous surface anyway, as desired!