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;
WSAD,C,Space - slide camera;
Mouse, Q,E - rotate camera;
R - reinitialize mesh;
T - tessellate once (break one edge resulting in 1 or 2 additional triangles added).
Full project code for VS12 can be downloaded from here.