nurbs1.c

Go to the documentation of this file.

/* --
OpenFX version 1.4 - Modelling, Animation and Rendering Package
Copyright (C) 2000 - 2005 OpenFX Development Team

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

You may contact the OpenFX development team via elecronic mail
at core@openfx.org, or visit our website at http://openfx.org for
further information and support details.
-- */

// File  NURBS1.C  included in NURBS.C

/*
 * Return the current knot the parameter u is less than or equal to.
 * Find this "breakpoint" allows the evaluation routines to concentrate on
 * only those control points actually effecting the curve around u.]
 *
 * m   is the number of points on the curve (or surface direction)
 * k   is the order of the curve (or surface direction)
 * kv  is the knot vector ([0..m+k-1]) to find the break point in.
 */

static long FindBreakPoint(double u, double *kv, long m, long k){
 long i;
 if (u == kv[m+1])   /* Special case for closed interval */
   return m;
 i = m + k;
 while ((u < kv[i]) && (i > 0))i--;
 return( i );
}

/*
 * Compute Bi,k(u), for i = 0..k.
 *  u      is the parameter of the spline to find the basis functions for
 *  brkPoint   is the start of the knot interval ("segment")
 *  kv      is the knot vector
 *  k      is the order of the curve
 *  bvals   is the array of returned basis values.
 *
 * (From Bartels, Beatty & Barsky, p.387)
 */

static void BasisFunctions(double u, long brkPoint, double * kv, long k,
                           double * bvals ){
 register long r, s, i;
 register double omega;
 bvals[0] = 1.0;
 for (r = 2; r <= k; r++){
   i = brkPoint - r + 1;
   bvals[r - 1] = 0.0;
   for (s = r-2; s >= 0; s--){
     i++;
     if (i < 0)
       omega = 0;
     else
       omega = (u - kv[i]) / (kv[i + r - 1] - kv[i]);
     bvals[s + 1] = bvals[s + 1] + (1 - omega) * bvals[s];
     bvals[s] = omega * bvals[s];
   }
 }
}

static void BasisDerivatives(double u, long brkPoint,
                             double * kv, long k, double * dvals ){
/*
 * Compute derivatives of the basis functions Bi,k(u)'
 */
 register long s, i;
 register double omega, knotScale;
 BasisFunctions( u, brkPoint, kv, k - 1, dvals );
 dvals[k-1] = 0.0;       /* BasisFunctions misses this */
 knotScale = kv[brkPoint + 1L] - kv[brkPoint];
 i = brkPoint - k + 1L;
 for (s = k - 2L; s >= 0L; s--){
   i++;
   omega = knotScale * ((double)(k-1L)) / (kv[i+k-1L] - kv[i]);
   dvals[s + 1L] += -omega * dvals[s];
   dvals[s] *= omega;
 }
}
/* eg.
 * Cubic NURBS mimimum requirement 4 control points
 * => num = 4  order = 4   => 8 Index:  knots [0 - 3]  points [0 - 3][0 - 3]
 * in the terms of Watt and Watt spline spec. is n=3 k=4 => knots = 3+4+1
 * Sum i=0,n  (points = n+1)
 * In knot vector first parameter repeated "k" times [0 0 0 0  1 1 1 1]
 */
/*
 * Calculate a point p on NurbSurface n at a specific u, v
 * using the tensor product.
 * If utan and vtan are not NULL, compute the u and v tangents as well.
 *
 * Note the valid parameter range for u and v is
 * (kvU[orderU-1] <= u < kvU[numU]), (kvV[orderV-1] <= v < kvV[numV])
 *
 * eg. Cubic spline with 4 control points: Vector and valid range: O+N = 4+4
 *     [0 0 0 0 1 1 1 1]
 * (kvU[3] <= u < kvU[4]), (kvV[4] <= v < kvV[4])
 * eg. Cubic spline with 5 control points: Vector and valid range: O+N = 4+5
 *     [0 0 0 0 1 2 2 2 2]
 * (kvU[3] <= u < kvU[5]), (kvV[3] <= v < kvV[5])
 * eg. Cubic spline with 8 control points: Vector and valid range: O+N = 4+8
 *     [0 0 0 0 1 2 3 4 5 5 5 5]
 * (kvU[3] <= u < kvU[8]), (kvV[3] <= v < kvV[8])
 * eg. quartic spline with 5 control points: Vector and valid range: O+N = 5+5
 *     [0 0 0 0 0 1 1 1 1 1]
 * (kvU[4] <= u < kvU[5]), (kvV[4] <= v < kvV[5])
 * Note the number of control points must be greater or equal to the order
 * numU >= orderU  and  numV >= orderV
 * Note the maximum value of an element of knot vector is N-O+1 in this case
 * u must lie in range  0 <= u < (N-O+1)
 */

static void CalcPoint(double u, double v, nurbs *n, vector p,
                      vector utan, vector vtan){
 long i, j, ri, rj;
 vector4 *cp;
 double tmp;
 double wsqrdiv;
 long ubrkPoint, ufirst;
 double bu[MAXORDER], buprime[MAXORDER];
 long vbrkPoint, vfirst;
 double bv[MAXORDER], bvprime[MAXORDER];
 vector4 r,rutan,rvtan;
 rutan.x= rvtan.x=r.x=0.0;
 rutan.y= rvtan.y=r.y=0.0;
 rutan.z= rvtan.z=r.z=0.0;
 rutan.w= rvtan.z=r.w=0.0;
 /* Evaluate non-uniform basis functions (and derivatives) */
 ubrkPoint = FindBreakPoint( u, n->kvU, n->numU-1, n->orderU );
 ufirst = ubrkPoint - n->orderU + 1;
 BasisFunctions( u, ubrkPoint, n->kvU, n->orderU, bu );
 if(utan != NULL)BasisDerivatives( u, ubrkPoint, n->kvU, n->orderU, buprime );
 vbrkPoint = FindBreakPoint( v, n->kvV, n->numV-1, n->orderV );
 vfirst = vbrkPoint - n->orderV + 1;
 BasisFunctions( v, vbrkPoint, n->kvV, n->orderV, bv );
 if(vtan != NULL)BasisDerivatives( v, vbrkPoint, n->kvV, n->orderV, bvprime );
 /* Weight control points against the basis functions */
 for (i = 0; i < n->orderV; i++)
 for (j = 0; j < n->orderU; j++){
   ri = n->orderV - 1L - i;
   rj = n->orderU - 1L - j;
   tmp = bu[rj] * bv[ri];
   cp = &(n->points[i+vfirst][j+ufirst]);
   r.x += cp->x * tmp;
   r.y += cp->y * tmp;
   r.z += cp->z * tmp;
   r.w += cp->w * tmp;
   if(utan != NULL){
     tmp = buprime[rj] * bv[ri];
     rutan.x += cp->x * tmp;
     rutan.y += cp->y * tmp;
     rutan.z += cp->z * tmp;
     rutan.w += cp->w * tmp;
   }
   if(vtan != NULL){
     tmp = bu[rj] * bvprime[ri];
     rvtan.x += cp->x * tmp;
     rvtan.y += cp->y * tmp;
     rvtan.z += cp->z * tmp;
     rvtan.w += cp->w * tmp;
   }
 }
 /* Project tangents, using the quotient rule for differentiation */
 if(utan != NULL && vtan != NULL)wsqrdiv = 1.0 / (r.w * r.w);
 if(utan != NULL){
   utan[0] = (r.x * rutan.x - rutan.x * r.x) * wsqrdiv;
   utan[1] = (r.y * rutan.y - rutan.y * r.y) * wsqrdiv;
   utan[2] = (r.z * rutan.z - rutan.z * r.z) * wsqrdiv;
 }
 if(vtan != NULL){
   vtan[0] = (r.x * rvtan.x - rvtan.x * r.x) * wsqrdiv;
   vtan[1] = (r.y * rvtan.y - rvtan.y * r.y) * wsqrdiv;
   vtan[2] = (r.z * rvtan.z - rvtan.z * r.z) * wsqrdiv;
 }
 p[0] = r.x / r.w;
 p[1] = r.y / r.w;
 p[2] = r.z / r.w;
 return;
}

static BOOL CalcAlpha(double * ukv, double * wkv,
                      long m, long n, long k,
                      double *** alpha){
 register long i, j;
 long brkPoint, r, rm1, last, s;
 double omega;
 double aval[MAXORDER];
 if(*alpha == NULL){
   if((*alpha = (double **)X__Malloc((long)((k+1) * sizeof(double *))))
              == NULL)return FALSE;
   for (i = 0; i <= k; i++){
     if(((*alpha)[i] = (double *)X__Malloc((long) ((m + n + 1)
                     * sizeof(double)))) == NULL)return FALSE;
    }
 }
 for (j = 0; j <= m + n; j++){
   brkPoint = FindBreakPoint( wkv[j], ukv, m, k );
   aval[0] = 1.0;
   for(r = 2; r <= k; r++){
     rm1 = r - 1;
     last = min(rm1,brkPoint);
     i = brkPoint - last;
     if(last < rm1)aval[last] = aval[last] * (wkv[j + r - 1] - ukv[i])
                              / (ukv[i + r - 1] - ukv[i]);
     else  aval[last] = 0.0;
     for(s = last - 1; s >= 0; s-- ){
       i++;
       omega = (wkv[j + r - 1] - ukv[i]) / (ukv[i + r - 1] - ukv[i]);
       aval[s + 1] = aval[s+1] + (1 - omega) * aval[s];
       aval[s] = omega * aval[s];
     }
   }
   last = min( k - 1, brkPoint );
   for (i = 0; i <= k; i++)(*alpha)[i][j] = 0.0;
   for (s = 0; s <= last; s++)(*alpha)[last - s][j] = aval[s];
 }
 return TRUE;
}

static void RefineSurface(nurbs *src, nurbs *dest, BOOL dirflag){
/*
 * Apply the alpha matrix computed above to the rows (or columns)
 * of the surface.  If dirflag is true do the U's (row), else do V's (col).
 */
 register long i, j, out;
 register vector4 *dp, *sp;
 long i1,brkPoint,maxj,maxout;
 register double tmp;
 double **alpha = NULL;
 /* Compute the alpha matrix and indexing variables for the requested direction */
 if(dirflag){
   if(!CalcAlpha(src->kvU, dest->kvU, src->numU - 1, dest->numU - src->numU,
       src->orderU, &alpha ))return;
   maxj = dest->numU;
   maxout = src->numV;
 }
 else{
   if(!CalcAlpha(src->kvV, dest->kvV, src->numV - 1, dest->numV - src->numV,
         src->orderV, &alpha))return;
   maxj = dest->numV;
   maxout = dest->numU;
 }
 /* Apply the alpha matrix to the original control points,
    generating new ones */
 for(out = 0; out < maxout; out++)
 for(j = 0; j < maxj; j++){
   if(dirflag){
     dp = &(dest->points[out][j]);
     brkPoint = FindBreakPoint(dest->kvU[j], src->kvU,
                               src->numU-1, src->orderU);
     i1 = max( brkPoint - src->orderU + 1, 0 );
     sp = &(src->points[out][i1]);
   }
   else{
     dp = &(dest->points[j][out]);
     brkPoint = FindBreakPoint(dest->kvV[j], src->kvV,
                               src->numV-1, src->orderV );
     i1 = max(brkPoint - src->orderV + 1, 0);
     sp = &(src->points[i1][out]);
   }
   dp->x = 0.0;
   dp->y = 0.0;
   dp->z = 0.0;
   dp->w = 0.0;
   for(i = i1; i <= brkPoint; i++){
     tmp = alpha[i - i1][j];
     sp = (dirflag ? &(src->points[out][i]) : &(src->points[i][out]) );
     dp->x += tmp * sp->x;
     dp->y += tmp * sp->y;
     dp->z += tmp * sp->z;
     dp->w += tmp * sp->w;
   }
 }
 /* Free up the alpha matrix */
 for (i = 0; i <= (dirflag ? src->orderU : src->orderV); i++)
   X__Free(alpha[i]);
 X__Free(alpha);
 return;
}

#define DIVW(rpt,pt ) \
    { (pt[0]) = (rpt)->x / (rpt)->w; \
      (pt[1]) = (rpt)->y / (rpt)->w; \
      (pt[2]) = (rpt)->z / (rpt)->w; }
#define EPSILON 0.0001   /* Used to determine when things are too small. */
#define DIVPT( p, dn ) { (p[0]) /= (dn); (p[1]) /= (dn); (p[2]) /= (dn); }
#define COPYVEC(b,a)   { b[0] = a[0]; b[1] = a[1]; b[2] = a[2]; }
#define maxV(surf) ((surf)->numV-1L)
#define maxU(surf) ((surf)->numU-1L)

static void ScreenProject(vector4 *Point4, vector Point3){
 DIVW(Point4,Point3);
 return;
}

static long SplitKV(double * srckv,
                    double ** destkv,
                    long * splitPt,    /* Where the knot interval is split */
                    long m, long k ){
/*
 * Split a knot vector at the center, by adding multiplicity k knots near
 * the middle of the parameter range.  Tries to start with an existing knot,
 * but will add a new knot value if there's nothing in "the middle" (e.g.,
 * a Bezier curve).
 */
 long i, last;
 long middex, extra, same;   /* "middex" ==> "middle index" */
 double midVal;
 extra = 0L;
 last = (m + k);
 middex = last / 2;
 midVal = srckv[middex];
 /* Search forward and backward to see if multiple knot is already there */
 i = middex+1L;
 same = 1L;
 while ((i < last) && (srckv[i] == midVal)){
   i++;
   same++;
 }
 i = middex-1L;
 while((i > 0L) && (srckv[i] == midVal)) {
   i--;
   middex--;       /* middex is start of multiple knot */
   same++;
 }
 if(i <= 0L){       /* No knot in middle, must create it */
   midVal = (srckv[0L] + srckv[last]) / 2.0;
   middex = last / 2L;
   while (srckv[middex + 1L] < midVal)middex++;
   same = 0L;
 }
 extra = k - same;
 if((*destkv = (double *)X__Malloc((long)(sizeof(double)
               * (m+k+extra+1L)))) == NULL)return 0;
 if(same < k){       /* Must add knots */
   for(i = 0L; i <= middex; i++) (*destkv)[i] = srckv[i];
   for(i = middex+1L; i <= middex+extra; i++)(*destkv)[i] = midVal;
   for(i = middex + k - same + 1L; i <= m + k + extra; i++){
     (*destkv)[i] = srckv[i - extra];
   }
 }
 else {
   for (i = 0L; i <= m + k; i++) (*destkv)[i] = srckv[i];
 }
 *splitPt = (extra < k) ? middex - 1L : middex;
 return extra;
}

static void ProjectToLine(vector firstPt, vector lastPt, vector midPt){
/*
 * Given a line defined by firstPt and lastPt, project midPt onto
 * that line.  Used for fixing "cracks".
 */
 vector base, v0, vm;
 double fraction, denom;
 VECCOPY(firstPt,base)
 VECSUB(lastPt, base, v0)
 VECSUB(midPt, base, vm)
 denom = (v0[0]*v0[0]+v0[1]*v0[1]+v0[2]*v0[2]);
 fraction = (denom == 0.0) ? 0.0 : (DOT(v0,vm)/denom);
 midPt[0] = base[0] + fraction * v0[0];
 midPt[1] = base[1] + fraction * v0[1];
 midPt[2] = base[2] + fraction * v0[2];
}

static void FixNormals(NurbsSurfSample *s0,
                       NurbsSurfSample *s1,
                       NurbsSurfSample *s2){
/*
 * If a normal has collapsed to zero (normLen == 0.0) then try
 * and fix it by looking at its neighbors.  If all the neighbors
 * are sick, then re-compute them from the plane they form.
 * If that fails too, then we give up...
 */
 BOOL goodnorm;
 long i, j, ok;
 double dist;
 NurbsSurfSample *V[3];
 vector norm;
 norm[0]=norm[1]=norm[2]=0.0;
 V[0] = s0; V[1] = s1; V[2] = s2;
 /* Find a reasonable normal */
 for(ok = 0, goodnorm = FALSE;
    (ok < 3L) && !(goodnorm = (V[ok]->normLen > 0.0)); ok++);
 if (! goodnorm){  /* All provided normals are zilch, try and invent one */
   norm[0] = 0.0; norm[0] = 0.0; norm[0] = 0.0;
   for (i = 0; i < 3L; i++){
     j = (i + 1L) % 3L;
     norm[0] += (V[i]->point[1] - V[j]->point[1]) * (V[i]->point[2] + V[j]->point[2]);
     norm[1] += (V[i]->point[2] - V[j]->point[2]) * (V[i]->point[0] + V[j]->point[0]);
     norm[2] += (V[i]->point[0] - V[j]->point[0]) * (V[i]->point[1] + V[j]->point[1]);
   }
   dist = sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]);
   if(dist == 0.0)return;      /* This sucker's hopeless... */
   DIVPT(norm,dist);
   for (i = 0; i < 3; i++){
     VECCOPY(norm,V[i]->normal)
     V[i]->normLen = dist;
   }
 }
 else{       /* Replace a sick normal with a healthy one nearby */
   for (i = 0; i < 3; i++)if ((i != ok) && (V[i]->normLen == 0.0))
   VECCOPY(V[ok]->normal,V[i]->normal);
 }
 return;
}

static void AdjustNormal(NurbsSurfSample *samp){
/*
 * Normalize the normal in a sample.  If it's degenerate,
 * flag it as such by setting the normLen to 0.0
 */
 /* If it's not degenerate, do the normalization now */
 samp->normLen = sqrt(samp->normal[0]*samp->normal[0]+
                      samp->normal[1]*samp->normal[1]+
                      samp->normal[2]*samp->normal[2]);
 if (samp->normLen < EPSILON)samp->normLen = 0.0;
 else DIVPT(samp->normal,samp->normLen);
}

static void GetNormal(nurbs *n, long indV, long indU){
/*
 * Compute the normal of a corner point of a mesh.  The
 * base is the value of the point at the corner, indU and indV
 * are the mesh indices of that point (either 0 or numU|numV).
 */
 vector tmpL,tmpR;   /* "Left" and "Right" of the base point */
 NurbsSurfSample *crnr;
 if ( (indU && indV) || ((! indU) && (!indV)) ){
   if (indU)
       crnr = &(n->cnn);
   else
       crnr = &(n->c00);
   DIVW( &(n->points[indV][(indU ? (indU-1L) : 1L)]), tmpL);
   DIVW( &(n->points[(indV ? (indV-1L) : 1L)][indU]), tmpR);
 }
 else {
   if(indU) crnr = &(n->c0n);
   else     crnr = &(n->cn0);
   DIVW( &(n->points[indV][(indU ? (indU-1L) : 1L)]), tmpR);
   DIVW( &(n->points[(indV ? (indV-1L) : 1L)][indU]), tmpL);
 }
 VECSUB(tmpL,crnr->point,tmpL)
 VECSUB(tmpR,crnr->point,tmpR)
 CROSS(tmpL,tmpR,crnr->normal)
 AdjustNormal(crnr);
}

static void MakeNewCorners(nurbs *parent,
      nurbs *kid0,
      nurbs *kid1,
      BOOL dirflag){
/*
 * Build the new corners in the two new surfaces, computing both
 * point on the surface along with the normal.   Prevent cracks that may occur.
 */
 DIVW( &(kid0->points[maxV(kid0)][maxU(kid0)]), kid0->cnn.point);
 GetNormal( kid0, maxV(kid0), maxU(kid0) );
 if(dirflag){
   kid0->strUn = FALSE;   /* Must re-test new edge straightness */
   DIVW( &(kid0->points[0L][maxU(kid0)]), kid0->c0n.point);
   GetNormal(kid0, 0L, maxU(kid0));
   /*
    * Normals must be re-calculated for kid1 in case the surface
    * was split at a c1 (or c0!) discontinutiy
    */
   COPYVEC(kid1->c00.point,kid0->c0n.point)
   GetNormal( kid1, 0L, 0L );
   COPYVEC(kid1->cn0.point,kid0->cnn.point)
   GetNormal( kid1, maxV(kid1), 0L );
   /*
    * Prevent cracks from forming by forcing the points on the seam to
    * lie along any straight edges.  (Must do this BEFORE finding normals)
    */
   if(parent->strV0)
     ProjectToLine((parent->c00.point),
                   (parent->c0n.point),
                   (kid0->c0n.point));
   if(parent->strVn)
     ProjectToLine((parent->cn0.point),
                   (parent->cnn.point),
                   (kid0->cnn.point));
   COPYVEC(kid1->c00.point,kid0->c0n.point)
   COPYVEC(kid1->cn0.point,kid0->cnn.point)
   kid1->strU0 = FALSE;
 }
 else{
   kid0->strVn = FALSE;
   DIVW( &(kid0->points[maxV(kid0)][0]), kid0->cn0.point);
   GetNormal( kid0, maxV(kid0), 0L );
   COPYVEC(kid1->c00.point,kid0->cn0.point)
   GetNormal( kid1, 0L, 0L );
   COPYVEC(kid1->c0n.point,kid0->cnn.point)
   GetNormal( kid1, 0L, maxU(kid1) );
   if(parent->strU0)
     ProjectToLine((parent->c00.point),
                   (parent->cn0.point),
                   (kid0->cn0.point));
   if(parent->strUn)
     ProjectToLine((parent->c0n.point),
                   (parent->cnn.point),
                   (kid0->cnn.point));
   COPYVEC(kid1->c00.point,kid0->cn0.point)
   COPYVEC(kid1->c0n.point,kid0->cnn.point)
   kid1->strV0 = FALSE;
 }
}

static void SplitSurface(nurbs *parent,
                         nurbs * kid0, nurbs * kid1,
                         BOOL dirflag){ /* If true subdivided U, else V */
/*
 * Split a surface into two halves.  First inserts multiplicity k knots
 * in the center of the parametric range.  After refinement, the two
 * resulting surfaces are copied into separate data structures.    If the
 * parent surface had straight edges, the points of the children are
 * projected onto those edges.
 */
 nurbs tmp;
 double *newkv;
 long i,j,splitPt;
 /*
  * Add a multiplicty k knot to the knot vector in the direction
  * specified by dirflag, and refine the surface.  This creates two
  * adjacent surfaces with c0 discontinuity at the seam.
  */
 tmp = *parent;   /* Copy order, # of points, etc. */
 if(dirflag){
   tmp.numU = parent->numU + SplitKV( parent->kvU,
                   &newkv,
                   &splitPt,
                   maxU(parent),
                   parent->orderU);
   AllocNurbs( &tmp, newkv, NULL );
   for (i = 0L; i < tmp.numV + tmp.orderV; i++)tmp.kvV[i] = parent->kvV[i];
 }
 else {
   tmp.numV = parent->numV + SplitKV( parent->kvV,
                   &newkv,
                   &splitPt,
                   maxV(parent),
                   parent->orderV );
   AllocNurbs( &tmp, NULL, newkv );
   for (i = 0L; i < tmp.numU + tmp.orderU; i++)tmp.kvU[i] = parent->kvU[i];
 }
 RefineSurface( parent, &tmp, dirflag );
 /*
  * Build the two child surfaces, and copy the data from the refined
  * version of the parent (tmp) into the two children
  */
 /* First half */
 *kid0 = *parent;   /* copy various edge flags and orders */
 kid0->numU = dirflag ? splitPt+1L : parent->numU;
 kid0->numV = dirflag ? parent->numV : splitPt+1L;
 kid0->kvU = kid0->kvV = NULL;
 kid0->points = NULL;
 AllocNurbs( kid0, NULL, NULL );
 for(i = 0L; i < kid0->numV; i++)   /* Copy the point and kv data */
 for(j = 0L; j < kid0->numU; j++)
       kid0->points[i][j] = tmp.points[i][j];
 for(i = 0L; i < kid0->orderU + kid0->numU; i++)
     kid0->kvU[i] = tmp.kvU[i];
 for(i = 0L; i < kid0->orderV + kid0->numV; i++)
     kid0->kvV[i] = tmp.kvV[i];
 /* Second half */
 splitPt++;
 *kid1 = *parent;
 kid1->numU = dirflag ? tmp.numU - splitPt : parent->numU;
 kid1->numV = dirflag ? parent->numV : tmp.numV - splitPt;
 kid1->kvU = kid1->kvV = NULL;
 kid1->points = NULL;
 AllocNurbs( kid1, NULL, NULL );
 for (i = 0L; i < kid1->numV; i++)   /* Copy the point and kv data */
 for (j = 0L; j < kid1->numU; j++)
     kid1->points[i][j]
     = tmp.points[dirflag ? i: (i + splitPt) ][dirflag ? (j + splitPt) : j];
 for (i = 0L; i < kid1->orderU + kid1->numU; i++)
   kid1->kvU[i] = tmp.kvU[dirflag ? (i + splitPt) : i];
 for (i = 0L; i < kid1->orderV + kid1->numV; i++)
   kid1->kvV[i] = tmp.kvV[dirflag ? i : (i + splitPt)];
 /* Construct new corners on the boundry between the two kids */
 MakeNewCorners(parent,kid0,kid1,dirflag);
 ReleaseNurbs(&tmp);       /* Get rid of refined parent */
}

#define GETPT(i) (( dirflag ? &(n->points[crvInd][i]) : &(n->points[i][crvInd])))

static BOOL IsCurveStraight(nurbs *n,
       double tolerance,
       long crvInd,
       BOOL dirflag){  /* If true, test in U direction, else test in V */
/*
 * Test if a particular row or column of control points in a mesh
 * is "straight" with respect to a particular tolerance.  Returns true
 * if it is.
 */
 vector p,vec,prod;
 vector cp,e0;
 long i,last;
 double linelen,dist;
 /* Special case: lines are automatically straight. */
 if((dirflag ? n->numU : n->numV) == 2L)return( TRUE );
 last = (dirflag ? n->numU : n->numV) - 1L;
 ScreenProject(GETPT( 0L ),e0);
 /* Form an initial line to test the other points against (skiping degen lines) */
 linelen = 0.0;
 for(i = last; (i > 0L) && (linelen < EPSILON); i--){
   ScreenProject(GETPT(i),cp);
   VECSUB(cp,e0,vec)
   linelen = sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
 }
 DIVPT(vec,linelen);
 if(linelen > EPSILON){  /* If no non-degenerate lines found, it's all degen */
   for(i = 1L; i <= last; i++){
     /* The cross product of the vector defining the
      * initial line with the vector of the current point
      * gives the distance to the line. */
     ScreenProject(GETPT(i),cp);
     VECSUB(cp,e0,p)
     CROSS(p,vec,prod)
     dist=sqrt(prod[0]*prod[0]+prod[1]*prod[1]+prod[2]*prod[2]);
     if (dist > tolerance)return FALSE;
   }
 }
 return  TRUE;
}

static BOOL TestFlat(nurbs *n, double tolerance){
/*
 * Check to see if a surface is flat.  Tests are only performed on edges and
 * directions that aren't already straight.  If an edge is flagged as straight
 * (from the parent surface) it is assumed it will stay that way.
 */
 long i;
 BOOL straight;
 vector cp00,cp0n,cpn0,cpnn,planeEqn;
 double dist,d,dd;
 /* Check edge straightness */
 if(! n->strU0)n->strU0 = IsCurveStraight( n, tolerance, 0L, FALSE);
 if(! n->strUn)n->strUn = IsCurveStraight( n, tolerance, maxU(n), FALSE);
 if(! n->strV0)n->strV0 = IsCurveStraight( n, tolerance, 0L, TRUE);
 if(! n->strVn)n->strVn = IsCurveStraight( n, tolerance, maxV(n), TRUE);
 /* Test to make sure control points are straight in U and V */
 straight = TRUE;
 if( (! n->flatU) && (n->strV0) && (n->strVn) ){
   for(i = 1L;
      (i < maxV(n)) && (straight = IsCurveStraight( n, tolerance, i, TRUE ));
      i++);
 }
 if (straight && n->strV0 && n->strVn)n->flatU = TRUE;
 straight = TRUE;
 if( (! n->flatV) && (n->strU0) && (n->strUn) ){
   for(i = 1L;
      (i < maxU(n)) && (straight = IsCurveStraight( n, tolerance, i, FALSE ));
       i++);
 }
 if (straight && n->strU0 && n->strUn)n->flatV = TRUE;
 if( (! n->flatV) || (! n->flatU))return FALSE;
 /* The surface can pass the above tests but still be twisted. */
 ScreenProject( &(n->points[0L][0L]),            cp00);
 ScreenProject( &(n->points[0L][maxU(n)]),       cp0n);
 ScreenProject( &(n->points[maxV(n)][0L]),       cpn0);
 ScreenProject( &(n->points[maxV(n)][maxU(n)]),  cpnn);
 VECSUB(cp0n,cp00,cp0n)    /* Make edges into vectors */
 VECSUB(cpn0,cp00,cpn0)
 /*
  * Compute the plane equation from two adjacent sides, and
  * measure the distance from the far point to the plane.  If it's
  * larger than tolerance, the surface is twisted.
  */
 CROSS(cpn0,cp0n,planeEqn)
 /* Normalize to keep adds in sync w/ mults */
 dd = 1.0/sqrt((planeEqn[0]*planeEqn[0]) +
               (planeEqn[1]*planeEqn[1]) +
               (planeEqn[2]*planeEqn[2]));
 VECSCALE(dd,planeEqn,planeEqn)
 d=DOT(planeEqn,cp00);
 dist=fabs(DOT(planeEqn,cpnn) - d);
 if(dist > tolerance)return FALSE;  /* Surface is twisted */
 return TRUE;
}

static void DoSubdivision(nurbs *n, double tolerance,
                          BOOL dirflag, long level ){
/*
 * The recursive subdivision algorithm.    Test if the surface is flat.
 * If so, split it into triangles.  Otherwise, split it into two halves,
 * and invoke the procedure on each half.
 */
 nurbs left, right;   /* ...or top or bottom. Whatever spins your wheels. */
 if(TestFlat(n,tolerance)){
   EmitTriangles( n );
 }
 else {
   if ( ((! n->flatV) && (! n->flatU)) || ((n->flatV) && (n->flatU)) )
       dirflag = ! dirflag;    /* If twisted or curved in both directions, */
   else {  /* then alternate subdivision direction */
     if(n->flatU)       /* Only split in directions that aren't flat */
       dirflag = FALSE;
     else
       dirflag = TRUE;
   }
   SplitSurface(n, &left, &right, dirflag );
   DoSubdivision(&left, tolerance, dirflag, level + 1L );
   DoSubdivision(&right, tolerance, dirflag, level + 1L );
   ReleaseNurbs(&left);
   ReleaseNurbs(&right);       /* Deallocate surfaces made by SplitSurface */
 }
}

static void EmitSubdivision(nurbs * surf){ // subdivided nurbs
 surf->flatV = FALSE;
 surf->flatU = FALSE;
 surf->strU0 = FALSE;
 surf->strUn = FALSE;
 surf->strV0 = FALSE;
 surf->strVn = FALSE;
 DIVW( &(surf->points[0L][0L]),                       surf->c00.point );
 DIVW( &(surf->points[0L][surf->numU-1L]),            surf->c0n.point );
 DIVW( &(surf->points[surf->numV-1L][0L]),            surf->cn0.point );
 DIVW( &(surf->points[surf->numV-1L][surf->numU-1L]), surf->cnn.point );
 GetNormal(surf, 0L, 0L );
 GetNormal(surf, 0L, maxU(surf) );
 GetNormal(surf, maxV(surf), 0L );
 GetNormal(surf, maxV(surf), maxU(surf) );
 /* Note surf is deallocated by the subdivision process */
 DoSubdivision( surf, SubdivTolerance, TRUE, 0L );
 return;
}