// Triangular patch interpolator
Point3 Patch::interp(PatchMesh *pMesh,float pu, float pv, float pw)
{
   // It had better be a triangular patch!
   assert(type == PATCH_TRI);
   Point3 p;
   PatchVert *vp = pMesh->verts;
   PatchVec *vecp = pMesh->vecs;
   floatpu2 = pu * pu;
   floatpu3 = pu2 * pu;
   floatpu4 = pu3 * pu;
   floatpv2 = pv * pv;
   floatpv3 = pv2 * pv;
   floatpv4 = pv3 * pv;
   floatpw2 = pw * pw;
   floatpw3 = pw2 * pw;
   floatpw4 = pw3 * pw;
  
   // Hold on to your hats -- Here it comes!
   p = vp[v[0]].p * pw4 +
   aux[0] * 4.0f * pu * pw3 +
   aux[1] * 6.0f * pu2 * pw2 +
   aux[2] * 4.0f * pu3 * pw +
   vp[v[1]].p * pu4 +
   aux[3] * 4.0f * pv * pu3 +
   aux[4] * 6.0f * pu2 * pv2 +
   aux[5] * 4.0f * pv3 * pu +
   vp[v[2]].p * pv4 +
   aux[6] * 4.0f * pw * pv3 +
   aux[7] * 6.0f * pv2 * pw2 +
   aux[8] * 4.0f * pw3 * pv +
   vecp[interior[0]].p * 12.0f * pu * pv * pw2 +
   vecp[interior[1]].p * 12.0f * pu2 * pv * pw +
   vecp[interior[2]].p * 12.0f * pu * pv2 * pw;
   returnp;
}
 
// Quadrilateral patch interpolator
Point3 Patch::interp(PatchMesh *pMesh,float pu, float pv)
{
   // It had better be a quad patch!
   assert(type == PATCH_QUAD);
   Point3 p;
   PatchVert *vp = pMesh->verts;
   PatchVec *vecp = pMesh->vecs;
   floatpu2 = pu * pu;
   floatpu1 = 1.0f - pu;
   floatpu12 = pu1 * pu1;
   floatu0 = pu12 * pu1;
   floatu1 = 3.0f * pu * pu12;
   floatu2 = 3.0f * pu2 * pu1;
   floatu3 = pu2 * pu;
   floatpv2 = pv * pv;
   floatpv1 = 1.0f - pv;
   floatpv12 = pv1 * pv1;
   floatv0 = pv12 * pv1;
   floatv1 = 3.0f * pv * pv12;
   floatv2 = 3.0f * pv2 * pv1;
   floatv3 = pv2 * pv;
  
   // Hold on to your hats -- Here it comes!
   p = vp[v[0]].p * u0 * v0 +
   vecp[vec[7]].p * u1 * v0 +
   vecp[vec[6]].p * u2 * v0 +
   vp[v[3]].p * u3 * v0 +
   vecp[vec[0]].p * u0 * v1 +
   vecp[interior[0]].p * u1 * v1 +
   vecp[interior[3]].p * u2 * v1 +
   vecp[vec[5]].p * u3 * v1 +
   vecp[vec[1]].p * u0 * v2 +
   vecp[interior[1]].p * u1 * v2 +
   vecp[interior[2]].p * u2 * v2 +
   vecp[vec[4]].p * u3 * v2 +
   vp[v[1]].p * u0 * v3 +
   vecp[vec[2]].p * u1 * v3 +
   vecp[vec[3]].p * u2 * v3 +
   vp[v[2]].p * u3 * v3;
   returnp;
}

float fpd = (float)patchDivs;
for(px = 0; px < numPatches; ++px)
{
   Patch &p = patches[px];
   switch(p.type)
   {
     case PATCH_TRI:
     {
      for(intax = patchDivs - 1; ax > 0; --ax)
       {
         for(int bx = 1; bx < ax; ++bx)
         {
          float u = (float)bx / fpd;
          float v = (float)(patchDivs - ax) / fpd;
          float w = 1.0f - u - v; // Barycentric validity guaranteed!
          mesh.setVert(vert++, p.interp(this, u, v, w));
         }
       }
     }
     break;
     case PATCH_QUAD:
     {
      for(intu = 1; u < patchDivs; ++u)
       {
         float fu = (float)u / fpd;
         for(int v = 1; v < patchDivs; ++v)
         {
          float fv = (float)v / fpd;
          mesh.setVert(vert++, p.interp(this, fu, fv));
         }
       }
     }
     break;
   }
}