NURBS1.C

Go to the documentation of this file.
00001 /* --
00002 OpenFX version 1.0 - Modelling, Animation and Rendering Package
00003 Copyright (C) 2000 OpenFX Development Team
00004 
00005 This program is free software; you can redistribute it and/or
00006 modify it under the terms of the GNU General Public License
00007 as published by the Free Software Foundation; either version 2
00008 of the License, or (at your option) any later version.
00009 
00010 This program is distributed in the hope that it will be useful,
00011 but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013 GNU General Public License for more details.
00014 
00015 You should have received a copy of the GNU General Public License
00016 along with this program; if not, write to the Free Software
00017 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00018 
00019 You may contact the OpenFX development team via elecronic mail
00020 at core@openfx.org, or visit our website at http://openfx.org for
00021 further information and support details.
00022 -- */
00023 
00024 // File  NURBS1.C  included in NURBS.C
00025 
00026 /*
00027  * Return the current knot the parameter u is less than or equal to.
00028  * Find this "breakpoint" allows the evaluation routines to concentrate on
00029  * only those control points actually effecting the curve around u.]
00030  *
00031  * m   is the number of points on the curve (or surface direction)
00032  * k   is the order of the curve (or surface direction)
00033  * kv  is the knot vector ([0..m+k-1]) to find the break point in.
00034  */
00035 
00036 static long FindBreakPoint(double u, double *kv, long m, long k){
00037  long i;
00038  if (u == kv[m+1])   /* Special case for closed interval */
00039    return m;
00040  i = m + k;
00041  while ((u < kv[i]) && (i > 0))i--;
00042  return( i );
00043 }
00044 
00045 /*
00046  * Compute Bi,k(u), for i = 0..k.
00047  *  u      is the parameter of the spline to find the basis functions for
00048  *  brkPoint   is the start of the knot interval ("segment")
00049  *  kv      is the knot vector
00050  *  k      is the order of the curve
00051  *  bvals   is the array of returned basis values.
00052  *
00053  * (From Bartels, Beatty & Barsky, p.387)
00054  */
00055 
00056 static void BasisFunctions(double u, long brkPoint, double * kv, long k,
00057                            double * bvals ){
00058  register long r, s, i;
00059  register double omega;
00060  bvals[0] = 1.0;
00061  for (r = 2; r <= k; r++){
00062    i = brkPoint - r + 1;
00063    bvals[r - 1] = 0.0;
00064    for (s = r-2; s >= 0; s--){
00065      i++;
00066      if (i < 0)
00067        omega = 0;
00068      else
00069        omega = (u - kv[i]) / (kv[i + r - 1] - kv[i]);
00070      bvals[s + 1] = bvals[s + 1] + (1 - omega) * bvals[s];
00071      bvals[s] = omega * bvals[s];
00072    }
00073  }
00074 }
00075 
00076 static void BasisDerivatives(double u, long brkPoint,
00077                              double * kv, long k, double * dvals ){
00078 /*
00079  * Compute derivatives of the basis functions Bi,k(u)'
00080  */
00081  register long s, i;
00082  register double omega, knotScale;
00083  BasisFunctions( u, brkPoint, kv, k - 1, dvals );
00084  dvals[k-1] = 0.0;       /* BasisFunctions misses this */
00085  knotScale = kv[brkPoint + 1L] - kv[brkPoint];
00086  i = brkPoint - k + 1L;
00087  for (s = k - 2L; s >= 0L; s--){
00088    i++;
00089    omega = knotScale * ((double)(k-1L)) / (kv[i+k-1L] - kv[i]);
00090    dvals[s + 1L] += -omega * dvals[s];
00091    dvals[s] *= omega;
00092  }
00093 }
00094 /* eg.
00095  * Cubic NURBS mimimum requirement 4 control points
00096  * => num = 4  order = 4   => 8 Index:  knots [0 - 3]  points [0 - 3][0 - 3]
00097  * in the terms of Watt and Watt spline spec. is n=3 k=4 => knots = 3+4+1
00098  * Sum i=0,n  (points = n+1)
00099  * In knot vector first parameter repeated "k" times [0 0 0 0  1 1 1 1]
00100  */
00101 /*
00102  * Calculate a point p on NurbSurface n at a specific u, v
00103  * using the tensor product.
00104  * If utan and vtan are not NULL, compute the u and v tangents as well.
00105  *
00106  * Note the valid parameter range for u and v is
00107  * (kvU[orderU-1] <= u < kvU[numU]), (kvV[orderV-1] <= v < kvV[numV])
00108  *
00109  * eg. Cubic spline with 4 control points: Vector and valid range: O+N = 4+4
00110  *     [0 0 0 0 1 1 1 1]
00111  * (kvU[3] <= u < kvU[4]), (kvV[4] <= v < kvV[4])
00112  * eg. Cubic spline with 5 control points: Vector and valid range: O+N = 4+5
00113  *     [0 0 0 0 1 2 2 2 2]
00114  * (kvU[3] <= u < kvU[5]), (kvV[3] <= v < kvV[5])
00115  * eg. Cubic spline with 8 control points: Vector and valid range: O+N = 4+8
00116  *     [0 0 0 0 1 2 3 4 5 5 5 5]
00117  * (kvU[3] <= u < kvU[8]), (kvV[3] <= v < kvV[8])
00118  * eg. quartic spline with 5 control points: Vector and valid range: O+N = 5+5
00119  *     [0 0 0 0 0 1 1 1 1 1]
00120  * (kvU[4] <= u < kvU[5]), (kvV[4] <= v < kvV[5])
00121  * Note the number of control points must be greater or equal to the order
00122  * numU >= orderU  and  numV >= orderV
00123  * Note the maximum value of an element of knot vector is N-O+1 in this case
00124  * u must lie in range  0 <= u < (N-O+1)
00125  */
00126 
00127 static void CalcPoint(double u, double v, nurbs *n, vector p,
00128                       vector utan, vector vtan){
00129  long i, j, ri, rj;
00130  vector4 *cp;
00131  double tmp;
00132  double wsqrdiv;
00133  long ubrkPoint, ufirst;
00134  double bu[MAXORDER], buprime[MAXORDER];
00135  long vbrkPoint, vfirst;
00136  double bv[MAXORDER], bvprime[MAXORDER];
00137  vector4 r,rutan,rvtan={0.0,0.0,0.0,0.0};
00138  rutan.x= rvtan.x=r.x=0.0;
00139  rutan.y= rvtan.y=r.y=0.0;
00140  rutan.z= rvtan.z=r.z=0.0;
00141  rutan.w= rvtan.z=r.w=0.0;
00142  /* Evaluate non-uniform basis functions (and derivatives) */
00143  ubrkPoint = FindBreakPoint( u, n->kvU, n->numU-1, n->orderU );
00144  ufirst = ubrkPoint - n->orderU + 1;
00145  BasisFunctions( u, ubrkPoint, n->kvU, n->orderU, bu );
00146  if(utan != NULL)BasisDerivatives( u, ubrkPoint, n->kvU, n->orderU, buprime );
00147  vbrkPoint = FindBreakPoint( v, n->kvV, n->numV-1, n->orderV );
00148  vfirst = vbrkPoint - n->orderV + 1;
00149  BasisFunctions( v, vbrkPoint, n->kvV, n->orderV, bv );
00150  if(vtan != NULL)BasisDerivatives( v, vbrkPoint, n->kvV, n->orderV, bvprime );
00151  /* Weight control points against the basis functions */
00152  for (i = 0; i < n->orderV; i++)
00153  for (j = 0; j < n->orderU; j++){
00154    ri = n->orderV - 1L - i;
00155    rj = n->orderU - 1L - j;
00156    tmp = bu[rj] * bv[ri];
00157    cp = &(n->points[i+vfirst][j+ufirst]);
00158    r.x += cp->x * tmp;
00159    r.y += cp->y * tmp;
00160    r.z += cp->z * tmp;
00161    r.w += cp->w * tmp;
00162    if(utan != NULL){
00163      tmp = buprime[rj] * bv[ri];
00164      rutan.x += cp->x * tmp;
00165      rutan.y += cp->y * tmp;
00166      rutan.z += cp->z * tmp;
00167      rutan.w += cp->w * tmp;
00168    }
00169    if(vtan != NULL){
00170      tmp = bu[rj] * bvprime[ri];
00171      rvtan.x += cp->x * tmp;
00172      rvtan.y += cp->y * tmp;
00173      rvtan.z += cp->z * tmp;
00174      rvtan.w += cp->w * tmp;
00175    }
00176  }
00177  /* Project tangents, using the quotient rule for differentiation */
00178  if(utan != NULL && vtan != NULL)wsqrdiv = 1.0 / (r.w * r.w);
00179  if(utan != NULL){
00180    utan[0] = (r.x * rutan.x - rutan.x * r.x) * wsqrdiv;
00181    utan[1] = (r.y * rutan.y - rutan.y * r.y) * wsqrdiv;
00182    utan[2] = (r.z * rutan.z - rutan.z * r.z) * wsqrdiv;
00183  }
00184  if(vtan != NULL){
00185    vtan[0] = (r.x * rvtan.x - rvtan.x * r.x) * wsqrdiv;
00186    vtan[1] = (r.y * rvtan.y - rvtan.y * r.y) * wsqrdiv;
00187    vtan[2] = (r.z * rvtan.z - rvtan.z * r.z) * wsqrdiv;
00188  }
00189  p[0] = r.x / r.w;
00190  p[1] = r.y / r.w;
00191  p[2] = r.z / r.w;
00192  return;
00193 }
00194 
00195 static BOOL CalcAlpha(double * ukv, double * wkv,
00196                       long m, long n, long k,
00197                       double *** alpha){
00198  register long i, j;
00199  long brkPoint, r, rm1, last, s;
00200  double omega;
00201  double aval[MAXORDER];
00202  if(*alpha == NULL){
00203    if((*alpha = (double **)X__Malloc((long)((k+1) * sizeof(double *))))
00204               == NULL)return FALSE;
00205    for (i = 0; i <= k; i++){
00206      if(((*alpha)[i] = (double *)X__Malloc((long) ((m + n + 1)
00207                      * sizeof(double)))) == NULL)return FALSE;
00208     }
00209  }
00210  for (j = 0; j <= m + n; j++){
00211    brkPoint = FindBreakPoint( wkv[j], ukv, m, k );
00212    aval[0] = 1.0;
00213    for(r = 2; r <= k; r++){
00214      rm1 = r - 1;
00215      last = min(rm1,brkPoint);
00216      i = brkPoint - last;
00217      if(last < rm1)aval[last] = aval[last] * (wkv[j + r - 1] - ukv[i])
00218                               / (ukv[i + r - 1] - ukv[i]);
00219      else  aval[last] = 0.0;
00220      for(s = last - 1; s >= 0; s-- ){
00221        i++;
00222        omega = (wkv[j + r - 1] - ukv[i]) / (ukv[i + r - 1] - ukv[i]);
00223        aval[s + 1] = aval[s+1] + (1 - omega) * aval[s];
00224        aval[s] = omega * aval[s];
00225      }
00226    }
00227    last = min( k - 1, brkPoint );
00228    for (i = 0; i <= k; i++)(*alpha)[i][j] = 0.0;
00229    for (s = 0; s <= last; s++)(*alpha)[last - s][j] = aval[s];
00230  }
00231  return TRUE;
00232 }
00233 
00234 static void RefineSurface(nurbs *src, nurbs *dest, BOOL dirflag){
00235 /*
00236  * Apply the alpha matrix computed above to the rows (or columns)
00237  * of the surface.  If dirflag is true do the U's (row), else do V's (col).
00238  */
00239  register long i, j, out;
00240  register vector4 *dp, *sp;
00241  long i1,brkPoint,maxj,maxout;
00242  register double tmp;
00243  double **alpha = NULL;
00244  /* Compute the alpha matrix and indexing variables for the requested direction */
00245  if(dirflag){
00246    if(!CalcAlpha(src->kvU, dest->kvU, src->numU - 1, dest->numU - src->numU,
00247        src->orderU, &alpha ))return;
00248    maxj = dest->numU;
00249    maxout = src->numV;
00250  }
00251  else{
00252    if(!CalcAlpha(src->kvV, dest->kvV, src->numV - 1, dest->numV - src->numV,
00253          src->orderV, &alpha))return;
00254    maxj = dest->numV;
00255    maxout = dest->numU;
00256  }
00257  /* Apply the alpha matrix to the original control points,
00258     generating new ones */
00259  for(out = 0; out < maxout; out++)
00260  for(j = 0; j < maxj; j++){
00261    if(dirflag){
00262      dp = &(dest->points[out][j]);
00263      brkPoint = FindBreakPoint(dest->kvU[j], src->kvU,
00264                                src->numU-1, src->orderU);
00265      i1 = max( brkPoint - src->orderU + 1, 0 );
00266      sp = &(src->points[out][i1]);
00267    }
00268    else{
00269      dp = &(dest->points[j][out]);
00270      brkPoint = FindBreakPoint(dest->kvV[j], src->kvV,
00271                                src->numV-1, src->orderV );
00272      i1 = max(brkPoint - src->orderV + 1, 0);
00273      sp = &(src->points[i1][out]);
00274    }
00275    dp->x = 0.0;
00276    dp->y = 0.0;
00277    dp->z = 0.0;
00278    dp->w = 0.0;
00279    for(i = i1; i <= brkPoint; i++){
00280      tmp = alpha[i - i1][j];
00281      sp = (dirflag ? &(src->points[out][i]) : &(src->points[i][out]) );
00282      dp->x += tmp * sp->x;
00283      dp->y += tmp * sp->y;
00284      dp->z += tmp * sp->z;
00285      dp->w += tmp * sp->w;
00286    }
00287  }
00288  /* Free up the alpha matrix */
00289  for (i = 0; i <= (dirflag ? src->orderU : src->orderV); i++)
00290    X__Free(alpha[i]);
00291  X__Free(alpha);
00292  return;
00293 }
00294 
00295 #define DIVW(rpt,pt ) \
00296     { (pt[0]) = (rpt)->x / (rpt)->w; \
00297       (pt[1]) = (rpt)->y / (rpt)->w; \
00298       (pt[2]) = (rpt)->z / (rpt)->w; }
00299 #define EPSILON 0.0001   /* Used to determine when things are too small. */
00300 #define DIVPT( p, dn ) { (p[0]) /= (dn); (p[1]) /= (dn); (p[2]) /= (dn); }
00301 #define COPYVEC(b,a)   { b[0] = a[0]; b[1] = a[1]; b[2] = a[2]; }
00302 #define maxV(surf) ((surf)->numV-1L)
00303 #define maxU(surf) ((surf)->numU-1L)
00304 
00305 static void ScreenProject(vector4 *Point4, vector Point3){
00306  DIVW(Point4,Point3);
00307  return;
00308 }
00309 
00310 static long SplitKV(double * srckv,
00311                     double ** destkv,
00312                     long * splitPt,    /* Where the knot interval is split */
00313                     long m, long k ){
00314 /*
00315  * Split a knot vector at the center, by adding multiplicity k knots near
00316  * the middle of the parameter range.  Tries to start with an existing knot,
00317  * but will add a new knot value if there's nothing in "the middle" (e.g.,
00318  * a Bezier curve).
00319  */
00320  long i, last;
00321  long middex, extra, same;   /* "middex" ==> "middle index" */
00322  double midVal;
00323  extra = 0L;
00324  last = (m + k);
00325  middex = last / 2;
00326  midVal = srckv[middex];
00327  /* Search forward and backward to see if multiple knot is already there */
00328  i = middex+1L;
00329  same = 1L;
00330  while ((i < last) && (srckv[i] == midVal)){
00331    i++;
00332    same++;
00333  }
00334  i = middex-1L;
00335  while((i > 0L) && (srckv[i] == midVal)) {
00336    i--;
00337    middex--;       /* middex is start of multiple knot */
00338    same++;
00339  }
00340  if(i <= 0L){       /* No knot in middle, must create it */
00341    midVal = (srckv[0L] + srckv[last]) / 2.0;
00342    middex = last / 2L;
00343    while (srckv[middex + 1L] < midVal)middex++;
00344    same = 0L;
00345  }
00346  extra = k - same;
00347  if((*destkv = (double *)X__Malloc((long)(sizeof(double)
00348                * (m+k+extra+1L)))) == NULL)return 0;
00349  if(same < k){       /* Must add knots */
00350    for(i = 0L; i <= middex; i++) (*destkv)[i] = srckv[i];
00351    for(i = middex+1L; i <= middex+extra; i++)(*destkv)[i] = midVal;
00352    for(i = middex + k - same + 1L; i <= m + k + extra; i++){
00353      (*destkv)[i] = srckv[i - extra];
00354    }
00355  }
00356  else {
00357    for (i = 0L; i <= m + k; i++) (*destkv)[i] = srckv[i];
00358  }
00359  *splitPt = (extra < k) ? middex - 1L : middex;
00360  return extra;
00361 }
00362 
00363 static void ProjectToLine(vector firstPt, vector lastPt, vector midPt){
00364 /*
00365  * Given a line defined by firstPt and lastPt, project midPt onto
00366  * that line.  Used for fixing "cracks".
00367  */
00368  vector base, v0, vm;
00369  double fraction, denom;
00370  VECCOPY(firstPt,base)
00371  VECSUB(lastPt, base, v0)
00372  VECSUB(midPt, base, vm)
00373  denom = (v0[0]*v0[0]+v0[1]*v0[1]+v0[2]*v0[2]);
00374  fraction = (denom == 0.0) ? 0.0 : (DOT(v0,vm)/denom);
00375  midPt[0] = base[0] + fraction * v0[0];
00376  midPt[1] = base[1] + fraction * v0[1];
00377  midPt[2] = base[2] + fraction * v0[2];
00378 }
00379 
00380 static void FixNormals(NurbsSurfSample *s0,
00381                        NurbsSurfSample *s1,
00382                        NurbsSurfSample *s2){
00383 /*
00384  * If a normal has collapsed to zero (normLen == 0.0) then try
00385  * and fix it by looking at its neighbors.  If all the neighbors
00386  * are sick, then re-compute them from the plane they form.
00387  * If that fails too, then we give up...
00388  */
00389  BOOL goodnorm;
00390  long i, j, ok;
00391  double dist;
00392  NurbsSurfSample *V[3];
00393  vector norm={0.0,0.0,0.0};
00394  V[0] = s0; V[1] = s1; V[2] = s2;
00395  /* Find a reasonable normal */
00396  for(ok = 0, goodnorm = FALSE;
00397     (ok < 3L) && !(goodnorm = (V[ok]->normLen > 0.0)); ok++);
00398  if (! goodnorm){  /* All provided normals are zilch, try and invent one */
00399    norm[0] = 0.0; norm[0] = 0.0; norm[0] = 0.0;
00400    for (i = 0; i < 3L; i++){
00401      j = (i + 1L) % 3L;
00402      norm[0] += (V[i]->point[1] - V[j]->point[1]) * (V[i]->point[2] + V[j]->point[2]);
00403      norm[1] += (V[i]->point[2] - V[j]->point[2]) * (V[i]->point[0] + V[j]->point[0]);
00404      norm[2] += (V[i]->point[0] - V[j]->point[0]) * (V[i]->point[1] + V[j]->point[1]);
00405    }
00406    dist = sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]);
00407    if(dist == 0.0)return;      /* This sucker's hopeless... */
00408    DIVPT(norm,dist);
00409    for (i = 0; i < 3; i++){
00410      VECCOPY(norm,V[i]->normal)
00411      V[i]->normLen = dist;
00412    }
00413  }
00414  else{       /* Replace a sick normal with a healthy one nearby */
00415    for (i = 0; i < 3; i++)if ((i != ok) && (V[i]->normLen == 0.0))
00416    VECCOPY(V[ok]->normal,V[i]->normal);
00417  }
00418  return;
00419 }
00420 
00421 static void AdjustNormal(NurbsSurfSample *samp){
00422 /*
00423  * Normalize the normal in a sample.  If it's degenerate,
00424  * flag it as such by setting the normLen to 0.0
00425  */
00426  /* If it's not degenerate, do the normalization now */
00427  samp->normLen = sqrt(samp->normal[0]*samp->normal[0]+
00428                       samp->normal[1]*samp->normal[1]+
00429                       samp->normal[2]*samp->normal[2]);
00430  if (samp->normLen < EPSILON)samp->normLen = 0.0;
00431  else DIVPT(samp->normal,samp->normLen);
00432 }
00433 
00434 static void GetNormal(nurbs *n, long indV, long indU){
00435 /*
00436  * Compute the normal of a corner point of a mesh.  The
00437  * base is the value of the point at the corner, indU and indV
00438  * are the mesh indices of that point (either 0 or numU|numV).
00439  */
00440  vector tmpL,tmpR;   /* "Left" and "Right" of the base point */
00441  NurbsSurfSample *crnr;
00442  if ( (indU && indV) || ((! indU) && (!indV)) ){
00443    if (indU)
00444        crnr = &(n->cnn);
00445    else
00446        crnr = &(n->c00);
00447    DIVW( &(n->points[indV][(indU ? (indU-1L) : 1L)]), tmpL);
00448    DIVW( &(n->points[(indV ? (indV-1L) : 1L)][indU]), tmpR);
00449  }
00450  else {
00451    if(indU) crnr = &(n->c0n);
00452    else     crnr = &(n->cn0);
00453    DIVW( &(n->points[indV][(indU ? (indU-1L) : 1L)]), tmpR);
00454    DIVW( &(n->points[(indV ? (indV-1L) : 1L)][indU]), tmpL);
00455  }
00456  VECSUB(tmpL,crnr->point,tmpL)
00457  VECSUB(tmpR,crnr->point,tmpR)
00458  CROSS(tmpL,tmpR,crnr->normal)
00459  AdjustNormal(crnr);
00460 }
00461 
00462 static void MakeNewCorners(nurbs *parent,
00463       nurbs *kid0,
00464       nurbs *kid1,
00465       BOOL dirflag){
00466 /*
00467  * Build the new corners in the two new surfaces, computing both
00468  * point on the surface along with the normal.   Prevent cracks that may occur.
00469  */
00470  DIVW( &(kid0->points[maxV(kid0)][maxU(kid0)]), kid0->cnn.point);
00471  GetNormal( kid0, maxV(kid0), maxU(kid0) );
00472  if(dirflag){
00473    kid0->strUn = FALSE;   /* Must re-test new edge straightness */
00474    DIVW( &(kid0->points[0L][maxU(kid0)]), kid0->c0n.point);
00475    GetNormal(kid0, 0L, maxU(kid0));
00476    /*
00477     * Normals must be re-calculated for kid1 in case the surface
00478     * was split at a c1 (or c0!) discontinutiy
00479     */
00480    COPYVEC(kid1->c00.point,kid0->c0n.point)
00481    GetNormal( kid1, 0L, 0L );
00482    COPYVEC(kid1->cn0.point,kid0->cnn.point)
00483    GetNormal( kid1, maxV(kid1), 0L );
00484    /*
00485     * Prevent cracks from forming by forcing the points on the seam to
00486     * lie along any straight edges.  (Must do this BEFORE finding normals)
00487     */
00488    if(parent->strV0)
00489      ProjectToLine((parent->c00.point),
00490                    (parent->c0n.point),
00491                    (kid0->c0n.point));
00492    if(parent->strVn)
00493      ProjectToLine((parent->cn0.point),
00494                    (parent->cnn.point),
00495                    (kid0->cnn.point));
00496    COPYVEC(kid1->c00.point,kid0->c0n.point)
00497    COPYVEC(kid1->cn0.point,kid0->cnn.point)
00498    kid1->strU0 = FALSE;
00499  }
00500  else{
00501    kid0->strVn = FALSE;
00502    DIVW( &(kid0->points[maxV(kid0)][0]), kid0->cn0.point);
00503    GetNormal( kid0, maxV(kid0), 0L );
00504    COPYVEC(kid1->c00.point,kid0->cn0.point)
00505    GetNormal( kid1, 0L, 0L );
00506    COPYVEC(kid1->c0n.point,kid0->cnn.point)
00507    GetNormal( kid1, 0L, maxU(kid1) );
00508    if(parent->strU0)
00509      ProjectToLine((parent->c00.point),
00510                    (parent->cn0.point),
00511                    (kid0->cn0.point));
00512    if(parent->strUn)
00513      ProjectToLine((parent->c0n.point),
00514                    (parent->cnn.point),
00515                    (kid0->cnn.point));
00516    COPYVEC(kid1->c00.point,kid0->cn0.point)
00517    COPYVEC(kid1->c0n.point,kid0->cnn.point)
00518    kid1->strV0 = FALSE;
00519  }
00520 }
00521 
00522 static void SplitSurface(nurbs *parent,
00523                          nurbs * kid0, nurbs * kid1,
00524                          BOOL dirflag){ /* If true subdivided U, else V */
00525 /*
00526  * Split a surface into two halves.  First inserts multiplicity k knots
00527  * in the center of the parametric range.  After refinement, the two
00528  * resulting surfaces are copied into separate data structures.    If the
00529  * parent surface had straight edges, the points of the children are
00530  * projected onto those edges.
00531  */
00532  nurbs tmp;
00533  double *newkv;
00534  long i,j,splitPt;
00535  /*
00536   * Add a multiplicty k knot to the knot vector in the direction
00537   * specified by dirflag, and refine the surface.  This creates two
00538   * adjacent surfaces with c0 discontinuity at the seam.
00539   */
00540  tmp = *parent;   /* Copy order, # of points, etc. */
00541  if(dirflag){
00542    tmp.numU = parent->numU + SplitKV( parent->kvU,
00543                    &newkv,
00544                    &splitPt,
00545                    maxU(parent),
00546                    parent->orderU);
00547    AllocNurbs( &tmp, newkv, NULL );
00548    for (i = 0L; i < tmp.numV + tmp.orderV; i++)tmp.kvV[i] = parent->kvV[i];
00549  }
00550  else {
00551    tmp.numV = parent->numV + SplitKV( parent->kvV,
00552                    &newkv,
00553                    &splitPt,
00554                    maxV(parent),
00555                    parent->orderV );
00556    AllocNurbs( &tmp, NULL, newkv );
00557    for (i = 0L; i < tmp.numU + tmp.orderU; i++)tmp.kvU[i] = parent->kvU[i];
00558  }
00559  RefineSurface( parent, &tmp, dirflag );
00560  /*
00561   * Build the two child surfaces, and copy the data from the refined
00562   * version of the parent (tmp) into the two children
00563   */
00564  /* First half */
00565  *kid0 = *parent;   /* copy various edge flags and orders */
00566  kid0->numU = dirflag ? splitPt+1L : parent->numU;
00567  kid0->numV = dirflag ? parent->numV : splitPt+1L;
00568  kid0->kvU = kid0->kvV = NULL;
00569  kid0->points = NULL;
00570  AllocNurbs( kid0, NULL, NULL );
00571  for(i = 0L; i < kid0->numV; i++)   /* Copy the point and kv data */
00572  for(j = 0L; j < kid0->numU; j++)
00573        kid0->points[i][j] = tmp.points[i][j];
00574  for(i = 0L; i < kid0->orderU + kid0->numU; i++)
00575      kid0->kvU[i] = tmp.kvU[i];
00576  for(i = 0L; i < kid0->orderV + kid0->numV; i++)
00577      kid0->kvV[i] = tmp.kvV[i];
00578  /* Second half */
00579  splitPt++;
00580  *kid1 = *parent;
00581  kid1->numU = dirflag ? tmp.numU - splitPt : parent->numU;
00582  kid1->numV = dirflag ? parent->numV : tmp.numV - splitPt;
00583  kid1->kvU = kid1->kvV = NULL;
00584  kid1->points = NULL;
00585  AllocNurbs( kid1, NULL, NULL );
00586  for (i = 0L; i < kid1->numV; i++)   /* Copy the point and kv data */
00587  for (j = 0L; j < kid1->numU; j++)
00588      kid1->points[i][j]
00589      = tmp.points[dirflag ? i: (i + splitPt) ][dirflag ? (j + splitPt) : j];
00590  for (i = 0L; i < kid1->orderU + kid1->numU; i++)
00591    kid1->kvU[i] = tmp.kvU[dirflag ? (i + splitPt) : i];
00592  for (i = 0L; i < kid1->orderV + kid1->numV; i++)
00593    kid1->kvV[i] = tmp.kvV[dirflag ? i : (i + splitPt)];
00594  /* Construct new corners on the boundry between the two kids */
00595  MakeNewCorners(parent,kid0,kid1,dirflag);
00596  ReleaseNurbs(&tmp);       /* Get rid of refined parent */
00597 }
00598 
00599 #define GETPT(i) (( dirflag ? &(n->points[crvInd][i]) : &(n->points[i][crvInd])))
00600 
00601 static BOOL IsCurveStraight(nurbs *n,
00602        double tolerance,
00603        long crvInd,
00604        BOOL dirflag){  /* If true, test in U direction, else test in V */
00605 /*
00606  * Test if a particular row or column of control points in a mesh
00607  * is "straight" with respect to a particular tolerance.  Returns true
00608  * if it is.
00609  */
00610  vector p,vec,prod;
00611  vector cp,e0;
00612  long i,last;
00613  double linelen,dist;
00614  /* Special case: lines are automatically straight. */
00615  if((dirflag ? n->numU : n->numV) == 2L)return( TRUE );
00616  last = (dirflag ? n->numU : n->numV) - 1L;
00617  ScreenProject(GETPT( 0L ),e0);
00618  /* Form an initial line to test the other points against (skiping degen lines) */
00619  linelen = 0.0;
00620  for(i = last; (i > 0L) && (linelen < EPSILON); i--){
00621    ScreenProject(GETPT(i),cp);
00622    VECSUB(cp,e0,vec)
00623    linelen = sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
00624  }
00625  DIVPT(vec,linelen);
00626  if(linelen > EPSILON){  /* If no non-degenerate lines found, it's all degen */
00627    for(i = 1L; i <= last; i++){
00628      /* The cross product of the vector defining the
00629       * initial line with the vector of the current point
00630       * gives the distance to the line. */
00631      ScreenProject(GETPT(i),cp);
00632      VECSUB(cp,e0,p)
00633      CROSS(p,vec,prod)
00634      dist=sqrt(prod[0]*prod[0]+prod[1]*prod[1]+prod[2]*prod[2]);
00635      if (dist > tolerance)return FALSE;
00636    }
00637  }
00638  return  TRUE;
00639 }
00640 
00641 static BOOL TestFlat(nurbs *n, double tolerance){
00642 /*
00643  * Check to see if a surface is flat.  Tests are only performed on edges and
00644  * directions that aren't already straight.  If an edge is flagged as straight
00645  * (from the parent surface) it is assumed it will stay that way.
00646  */
00647  long i;
00648  BOOL straight;
00649  vector cp00,cp0n,cpn0,cpnn,planeEqn;
00650  double dist,d,dd;
00651  /* Check edge straightness */
00652  if(! n->strU0)n->strU0 = IsCurveStraight( n, tolerance, 0L, FALSE);
00653  if(! n->strUn)n->strUn = IsCurveStraight( n, tolerance, maxU(n), FALSE);
00654  if(! n->strV0)n->strV0 = IsCurveStraight( n, tolerance, 0L, TRUE);
00655  if(! n->strVn)n->strVn = IsCurveStraight( n, tolerance, maxV(n), TRUE);
00656  /* Test to make sure control points are straight in U and V */
00657  straight = TRUE;
00658  if( (! n->flatU) && (n->strV0) && (n->strVn) ){
00659    for(i = 1L;
00660       (i < maxV(n)) && (straight = IsCurveStraight( n, tolerance, i, TRUE ));
00661       i++);
00662  }
00663  if (straight && n->strV0 && n->strVn)n->flatU = TRUE;
00664  straight = TRUE;
00665  if( (! n->flatV) && (n->strU0) && (n->strUn) ){
00666    for(i = 1L;
00667       (i < maxU(n)) && (straight = IsCurveStraight( n, tolerance, i, FALSE ));
00668        i++);
00669  }
00670  if (straight && n->strU0 && n->strUn)n->flatV = TRUE;
00671  if( (! n->flatV) || (! n->flatU))return FALSE;
00672  /* The surface can pass the above tests but still be twisted. */
00673  ScreenProject( &(n->points[0L][0L]),            cp00);
00674  ScreenProject( &(n->points[0L][maxU(n)]),       cp0n);
00675  ScreenProject( &(n->points[maxV(n)][0L]),       cpn0);
00676  ScreenProject( &(n->points[maxV(n)][maxU(n)]),  cpnn);
00677  VECSUB(cp0n,cp00,cp0n)    /* Make edges into vectors */
00678  VECSUB(cpn0,cp00,cpn0)
00679  /*
00680   * Compute the plane equation from two adjacent sides, and
00681   * measure the distance from the far point to the plane.  If it's
00682   * larger than tolerance, the surface is twisted.
00683   */
00684  CROSS(cpn0,cp0n,planeEqn)
00685  /* Normalize to keep adds in sync w/ mults */
00686  dd = 1.0/sqrt((planeEqn[0]*planeEqn[0]) +
00687                (planeEqn[1]*planeEqn[1]) +
00688                (planeEqn[2]*planeEqn[2]));
00689  VECSCALE(dd,planeEqn,planeEqn)
00690  d=DOT(planeEqn,cp00);
00691  dist=fabs(DOT(planeEqn,cpnn) - d);
00692  if(dist > tolerance)return FALSE;  /* Surface is twisted */
00693  return TRUE;
00694 }
00695 
00696 static void DoSubdivision(nurbs *n, double tolerance,
00697                           BOOL dirflag, long level ){
00698 /*
00699  * The recursive subdivision algorithm.    Test if the surface is flat.
00700  * If so, split it into triangles.  Otherwise, split it into two halves,
00701  * and invoke the procedure on each half.
00702  */
00703  nurbs left, right;   /* ...or top or bottom. Whatever spins your wheels. */
00704  if(TestFlat(n,tolerance)){
00705    EmitTriangles( n );
00706  }
00707  else {
00708    if ( ((! n->flatV) && (! n->flatU)) || ((n->flatV) && (n->flatU)) )
00709        dirflag = ! dirflag;    /* If twisted or curved in both directions, */
00710    else {  /* then alternate subdivision direction */
00711      if(n->flatU)       /* Only split in directions that aren't flat */
00712        dirflag = FALSE;
00713      else
00714        dirflag = TRUE;
00715    }
00716    SplitSurface(n, &left, &right, dirflag );
00717    DoSubdivision(&left, tolerance, dirflag, level + 1L );
00718    DoSubdivision(&right, tolerance, dirflag, level + 1L );
00719    ReleaseNurbs(&left);
00720    ReleaseNurbs(&right);       /* Deallocate surfaces made by SplitSurface */
00721  }
00722 }
00723 
00724 static void EmitSubdivision(nurbs * surf){ // subdivided nurbs
00725  surf->flatV = FALSE;
00726  surf->flatU = FALSE;
00727  surf->strU0 = FALSE;
00728  surf->strUn = FALSE;
00729  surf->strV0 = FALSE;
00730  surf->strVn = FALSE;
00731  DIVW( &(surf->points[0L][0L]),                       surf->c00.point );
00732  DIVW( &(surf->points[0L][surf->numU-1L]),            surf->c0n.point );
00733  DIVW( &(surf->points[surf->numV-1L][0L]),            surf->cn0.point );
00734  DIVW( &(surf->points[surf->numV-1L][surf->numU-1L]), surf->cnn.point );
00735  GetNormal(surf, 0L, 0L );
00736  GetNormal(surf, 0L, maxU(surf) );
00737  GetNormal(surf, maxV(surf), 0L );
00738  GetNormal(surf, maxV(surf), maxU(surf) );
00739  /* Note surf is deallocated by the subdivision process */
00740  DoSubdivision( surf, SubdivTolerance, TRUE, 0L );
00741  return;
00742 }
00743 

Generated on Sun Apr 27 14:20:11 2014 for OpenFX by  doxygen 1.5.6