---

nurbsSub.cpp

/*=============================================================================
        File: nurbsSub.cpp
     Purpose:       
    Revision: $Id: nurbsSub.cpp,v 1.3 2003/01/13 19:42:10 philosophil Exp $
  Created by: Philippe Lavoie          (20 Januray, 1999)
 Modified by: 

 Copyright notice:
          Copyright (C) 1999 Philippe Lavoie
 
          This library is free software; you can redistribute it and/or
          modify it under the terms of the GNU Library General Public
          License as published by the Free Software Foundation; either
          version 2 of the License, or (at your option) any later version.
 
          This library 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
          Library General Public License for more details.
 
          You should have received a copy of the GNU Library General Public
          License along with this library; if not, write to the Free
          Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
=============================================================================*/

#ifndef PLIB_NURBS_NURBSSUB_SOURCE
#define PLIB_NURBS_NURBSSUB_SOURCE

#include "nurbsSub.h"
#include <iostream>
#include <cstdio>
#include <stdio.h>

namespace PLib {

  
  const int FALSE_ = 0 ;
  const int TRUE_ = 1 ; 

  struct BOOL{
    int value ;
    BOOL(): value(FALSE_) {;}
    BOOL(int a) { if(a>0) value=TRUE_; else value=FALSE_;}
    BOOL(const BOOL& b) { value = b.value ; }
    int operator!() { if(value>0) return FALSE_; return TRUE_ ;}
    operator int() const {return value;}
    int& operator=(int a) { value = a ;  return value ; }
  };

  template <class T>
    NurbsSubSurface<T>::NurbsSubSurface(const NurbsSurface<T,3>& s): rsurf(s) {
    surf = 0 ;
    render = 0 ; 
    //initSurf();
    
  }

#define MAXORDER 20         /* Maximum order allowed (for local array sizes) */

  template <class T>
    struct NurbSurface {
      /* Number of Points in the U and V directions, respectivly */
      int numU, numV;
      /* Order of the surface in U and V (must be >= 2, < MAXORDER) */
      int orderU, orderV;
      /* Knot vectors, indexed as [0..numU+orderU-1] and [0..numV+orderV-1] */
      T * kvU, * kvV;
      /* Control points, indexed as points[0..numV-1][0..numU-1] */
      /* Note the w values are *premultiplied* with the x, y and z values */
      Matrix<HPoint_nD<T,3> >& points;
      
      /* These fields are added to support subdivision */
      BOOL strV0, strVn,   /* Edge straightness flags for subdivision */
      strU0, strUn;
      BOOL flatV, flatU;   /* Surface flatness flags for subdivision */
      SurfSample<T> c00, c0n,
      cn0, cnn;    /* Corner data structures for subdivision */
      RenderMesh<T> *render;
      
      static T epsilon ;
      
      NurbSurface(): points(tmpPoints) {   render = 0 ; }
      NurbSurface(Matrix<HPoint_nD<T,3> >& s, int du, int dv);
      
    protected:
      Matrix<HPoint_nD<T,3> > tmpPoints ;
      
    };

#define CHECK( n ) \
  { if (!(n)) { fprintf( stderr, "Ran out of memory\n" ); exit(-1); } }

#define DIVW( rpt, pt ) \
  { (pt)->x() = (rpt)->x() / (rpt)->w(); \
  (pt)->y() = (rpt)->y() / (rpt)->w(); \
  (pt)->z() = (rpt)->z() / (rpt)->w(); }
  
  
  /* Function prototypes */
  
  template <class T> void DrawSubdivision( NurbSurface<T> *, T tolerance );
  template <class T> void DrawEvaluation( NurbSurface<T> * );
  
  template <class T> int FindBreakPoint( T u, T * kv, int m, int k );
  template <class T> void AllocNurb( NurbSurface<T> *, T *, T * );
  template <class T> void CloneNurb( NurbSurface<T> *, NurbSurface<T> * );
  template <class T> void FreeNurb( NurbSurface<T> * );
  template <class T> void RefineSurface( NurbSurface<T> *, NurbSurface<T> *, BOOL );
  
  template <class T> void CalcPoint( T, T, NurbSurface<T> *, Point_nD<T,3> *, Point_nD<T,3> *, Point_nD<T,3> * );
  
  
  template <class T>
    void NurbsSubSurface<T>::drawSubdivision(T tolerance)
    {
      initSurf();
      surf->render = render ; 
      DrawSubdivision(surf,tolerance);
    }
  
  template <class T>
    void NurbsSubSurface<T>::drawSubdivisionPS(const char* f, T tolerance)
    {
      ofstream fout ;
      fout.open(f) ;
      if(fout)
        drawSubdivisionPS(fout,tolerance) ;
      fout.close();
    }

  template <class T>
    void NurbsSubSurface<T>::drawSubdivisionVRML(const char* f, T tolerance, const Color& col)
    {
      ofstream fout ;
      fout.open(f) ;
      if(fout)
        drawSubdivisionVRML(fout,tolerance,col) ;
      fout.close();
    }

  template <class T>
    void NurbsSubSurface<T>::drawSubdivisionVRML97(const char* f, T tolerance, const Color& col)
    {
      ofstream fout ;
      fout.open(f) ;
      if(fout)
        drawSubdivisionVRML97(fout,tolerance,col) ;
      fout.close();
    }

  template <class T>
    void NurbsSubSurface<T>::drawSubdivisionPS(ostream &os, T tolerance)
    {
      if(render)
        delete render ;
      render = new RenderMeshPS<T>(os) ;
      drawSubdivision(tolerance) ;
    }

  template <class T>
    void NurbsSubSurface<T>::drawSubdivisionVRML(ostream &os, T tolerance, const Color& col)
    {
      if(render) 
        delete render ;
      render = new RenderMeshVRML<T>(os,col) ;
      drawSubdivision(tolerance) ;
    }

  template <class T>
    void NurbsSubSurface<T>::drawSubdivisionVRML97(ostream &os, T tolerance, const Color& col)
    {
      if(render) 
        delete render ;
      render = new RenderMeshVRML97<T>(os,col) ;
      drawSubdivision(tolerance) ;
    }

  template <class T>
  //void NurbsSubSurface<T>::drawSubdivisionPoints(vector<Point_nD<T,3> > &pnts, T tolerance)
  void NurbsSubSurface<T>::drawSubdivisionPoints(BasicArray<Point_nD<T,3> >&pnts, T tolerance)
    {
      if(render) 
        delete render ;
      render = new RenderMeshPoints<T>(pnts) ;
      drawSubdivision(tolerance) ;
    }

  template <class T>
    void NurbsSubSurface<T>::initSurf()
    {
      if(surf)
        delete surf ;
      surf = new NurbSurface<T> ;
      surf->numU = rsurf.ctrlPnts().rows() ;
      surf->numV = rsurf.ctrlPnts().cols() ; 
      surf->orderU = rsurf.degreeU()+1 ;
      surf->orderV = rsurf.degreeV()+1 ;
      surf->kvU = new T[surf->numU + surf->orderU];
      surf->kvV = new T[surf->numV + surf->orderV];
      surf->points.resize(surf->numV,surf->numU);
      
      surf->flatV = FALSE_;
      surf->flatU = FALSE_;
      surf->strU0 = FALSE_;
      surf->strUn = FALSE_;
      surf->strV0 = FALSE_;
      surf->strVn = FALSE_;

      surf->render = render ; 
     
      int i;
      for(i=rsurf.knotU().n()-1;i>=0;--i){
        surf->kvU[i] = rsurf.knotU()[i] ;
      }
      for(i=rsurf.knotV().n()-1;i>=0;--i){
        surf->kvV[i] = rsurf.knotV()[i] ;
      }
      for(i=rsurf.ctrlPnts().rows()-1;i>=0;--i)
        for(int j=rsurf.ctrlPnts().cols()-1;j>=0;--j)
          surf->points(j,i) = rsurf.ctrlPnts()(i,j) ; 
    }
  
  template <class T>
    NurbSurface<T>::NurbSurface(Matrix<HPoint_nD<T,3> >& s, int du, int dv): points(s), orderU(du+1), orderV(dv+1){
    numU = points.rows() ;
    numV = points.cols() ;
    kvU = new T[numU + orderU] ;
    kvV = new T[numV + orderV] ; 
    render = 0 ; 
  }
  
  template <class T>
    NurbsSubSurface<T>::~NurbsSubSurface()
    {      
      if(surf){
        delete []surf->kvU ;
        delete []surf->kvV ;
      }
    }
  
  template <class T>
    void RenderMeshPS<T>::screenProject(const HPoint_nD<T,3> &worldPt, Point_nD<T,3> &screenPt )
    {
      screenPt.x() = worldPt.x() / worldPt.w() * 100 + 200;
      screenPt.y() = worldPt.y() / worldPt.w() * 100 + 200;
      screenPt.z() = worldPt.z() / worldPt.w() * 100 + 200;
    }

  template <class T>
    void RenderMeshVRML<T>::screenProject(const HPoint_nD<T,3> &worldPt, Point_nD<T,3> &screenPt )
    {
      screenPt = project(worldPt) ;
    }


  template <class T>
    void RenderMeshVRML97<T>::screenProject(const HPoint_nD<T,3> &worldPt, Point_nD<T,3> &screenPt )
    {      
      screenPt = project(worldPt) ;
      if(init){
        p_min = p_max = screenPt ;
        init = 0 ; 
      }
      if(screenPt.x()<p_min.x()) p_min.x() = screenPt.x();
      if(screenPt.y()<p_min.y()) p_min.y() = screenPt.y();
      if(screenPt.z()<p_min.z()) p_min.z() = screenPt.z();
      if(screenPt.x()>p_max.x()) p_max.x() = screenPt.x();
      if(screenPt.y()>p_max.y()) p_max.y() = screenPt.y();
      if(screenPt.z()>p_max.z()) p_max.z() = screenPt.z();
    }


  template <class T>
    void RenderMeshPoints<T>::screenProject(const HPoint_nD<T,3> &worldPt, Point_nD<T,3> &screenPt )
    {
      screenPt = project(worldPt) ;
    }


  inline void LineTo( ostream& out, short x, short y )
    {
      out << (int)x << " " << (int)y << " lineto\n" ; 
    }
  
  inline void MoveTo( ostream& out, short x, short y )
    {
      out << (int)x << " " << (int)y << " moveto\n" ; 
    }
  
  template <class T>
  void RenderMeshPS<T>::drawHeader(){
    out << "%!PS-Adobe-2.1\n";
    out << "%%Title: code5_0.ps (20)\n" ;
    out << "%%Creator: color_grid_generator\n" ; 
    out << "%%BoundingBox: 0 0 500 500\n" ;
    out << "%%Pages: 0\n" ;
    out << "%%EndComments\n" ;
    out << "0 setlinewidth\n" ;
    out << "0 0 0 setrgbcolor\n" ;
  }

  template <class T>
    void RenderMeshPS<T>::drawFooter(){
    out << "showpage\n%%EOF\n" ; 
  }


  template <class T>
    void RenderMeshPS<T>::drawTriangle( const SurfSample<T> &v0, const SurfSample<T> &v1, const SurfSample<T> & v2 )
    {
      out << "newpath\n" ;
      MoveTo( out, (short) (v0.point.x() * 100 + 200), (short) (v0.point.y() * 100 + 200) );
      LineTo( out, (short) (v1.point.x() * 100 + 200), (short) (v1.point.y() * 100 + 200) );
      LineTo( out, (short) (v2.point.x() * 100 + 200), (short) (v2.point.y() * 100 + 200) );
      LineTo( out, (short) (v0.point.x() * 100 + 200), (short) (v0.point.y() * 100 + 200) );
      out << "stroke\n" ;
    }

  template <class T>
    void RenderMeshPS<T>::drawLine( const SurfSample<T> &v0, const SurfSample<T> &v1)
    {
      out << "newpath\n" ; 
      MoveTo(out , (short) (v0.point.x() * 100 + 200), (short) (v0.point.y() * 100 + 200) );
      LineTo( out, (short) (v1.point.x() * 100 + 200), (short) (v1.point.y() * 100 + 200) );
      out << "stroke\n" ; 
    }


  template <class T>
    void RenderMeshVRML97<T>::drawHeader()
    {
      size = 0 ; 
      out << "#VRML V2.0 utf8\n" ;
      out << "#  Generated by Phil's NURBS library\n" ;
      out << "\nGroup {\n" ;
      out << "\n  children [\n" ; 
      out << "\tDEF T Transform {\n"; 
      out << "\t  children [\n" ;
      out << "\t\tShape {\n" ;
      out << "\t\t appearance Appearance {\n" ; 
      out << "\t\t\tmaterial Material{ diffuseColor " << float(color.r/255.0) << ' ' << float(color.g/255.0) << ' ' << float(color.b/255.0) << " }\n" ;
      out << "\t\t }\n" ; 
      out << "\t\t geometry IndexedFaceSet {\n" ;
      out << "\t\t\tsolid FALSE\n" ; 
      out << "\t\t\tcoord Coordinate {\n" ;
      out << "\t\t\t point [\n" ;
    }
 
  template <class T>
    void RenderMeshVRML97<T>::drawFooter(){

    out << "\t\t\t ]\n" ; // point [
    out << "\t\t\t}\n" ; // coord
    
    out << "\t\t\t coordIndex [\n" ;
  
    for(int i=0;i<size;++i){
      out << "\t\t\t\t" << 3*i << ", " << 3*i+1 << ", " << 3*i+2 << ", -1,\n" ;
    }
    out << "\t\t\t ]\n" ;
    out << "\t\t\t}\n" ; // IndexedFaceSet
    out << "\t\t}\n" ;
    out << "\t ]\n" ; 
    out << "\t}\n" ;
    out << "  ]\n" ; 
    out << "}\n" ; 
    
    Point_nD<T,3> p_mid((p_max.x()+p_min.x())/T(2),
                        (p_max.y()+p_min.y())/T(2),
                        (p_max.z()+p_min.z())/T(2));
    
    T x_axis = p_max.x() - p_min.x() ;
    T y_axis = p_max.y() - p_min.y() ; 
    
    T axis = (x_axis< y_axis) ? y_axis : x_axis ;
    axis *= T(2) ;
    
    out << "Viewpoint {\n\t position " << p_mid.x() << ' ' << p_mid.y() << ' ' << p_max.z()+axis << "\n\t description \"top\"\n}\n" ; 
    
    
    out << "NavigationInfo { type \"EXAMINE\" }\n" ; 

  }

  template <class T>
    void RenderMeshVRML<T>::drawHeader()
    {
      size = 0 ; 
      out << "#VRML V1.0 ascii\n" ;
      out << "#  Generated by Phil's NURBS library\n" ;
      out << "\nSeparator {\n" << "\tMaterialBinding { value PER_VERTEX_INDEXED }\n"  ;
      out << "\tMaterial{ ambientColor 0.25 0.25 0.25\n\t\tdiffuseColor " << float(color.r/255.0) << ' ' << float(color.g/255.0) << ' ' << float(color.b/255.0) << " }\n" ;
      out << "\tCoordinate3 {\n" ;
      out << "\t\tpoint [\n" ;
    }
 
  template <class T>
    void RenderMeshVRML<T>::drawFooter(){

    out << "\t\t]\n" ; // point [
    out << "\t}\n" ; // cordinate3
    
    out << "\tIndexedFaceSet{\n" ;
    out << "\t\tcoordIndex[\n" ;
    
    for(int i=0;i<size;++i){
      out << "\t\t\t" << 3*i << ", " << 3*i+1 << ", " << 3*i+2 << ", -1,\n" ;
    }
    out << "\t\t]\n" ;
    out << "\t}\n" ; // IndexedFaceSet
    
    out << "}\n" ;
  }


  template <class T>
    void RenderMeshVRML<T>::drawTriangle(const SurfSample<T> &v0, const SurfSample<T> &v1, const SurfSample<T> & v2 )
    {
      ++size ;
      out << "\t\t\t" << v0.point.x() << " " << v0.point.y() << " " << v0.point.z() << ",\n" ;
      out << "\t\t\t" << v1.point.x() << " " << v1.point.y() << " " << v1.point.z() << ",\n" ;
      out << "\t\t\t" << v2.point.x() << " " << v2.point.y() << " " << v2.point.z() << ",\n" ;
      
    }
  
  template <class T>
    void RenderMeshVRML97<T>::drawTriangle(const SurfSample<T> &v0, const SurfSample<T> &v1, const SurfSample<T> & v2 )
    {
      ++size ;
      out << "\t\t\t\t" << v0.point.x() << " " << v0.point.y() << " " << v0.point.z() << ",\n" ;
      out << "\t\t\t\t" << v1.point.x() << " " << v1.point.y() << " " << v1.point.z() << ",\n" ;
      out << "\t\t\t\t" << v2.point.x() << " " << v2.point.y() << " " << v2.point.z() << ",\n" ;
      
    }
  
  template <class T>
    void RenderMeshPoints<T>::drawHeader()
    {
      points.clear();
    }
 
  template <class T>
    void RenderMeshPoints<T>::drawFooter(){
  }


  template <class T>
    void RenderMeshPoints<T>::drawTriangle(const SurfSample<T> &v0, const SurfSample<T> &v1, const SurfSample<T> & v2 )
    {
      // naive method
      points.push_back(v0.point);
      points.push_back(v1.point);
      points.push_back(v2.point);
    }
  
template <class T>
int
FindBreakPoint( T u, T * kv, int m, int k )
{
  int 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)
 */

template <class T> void
BasisFunctions( T u, int brkPoint, T * kv, int k, T * bvals )
{
  int r, s, i;
  T 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];
        }
    }
}

/*
 * Compute derivatives of the basis functions Bi,k(u)'
 */
template <class T> void
BasisDerivatives( T u, int brkPoint, T * kv, int k, T * dvals )
{
  register int s, i;
  T omega, knotScale;

  BasisFunctions( u, brkPoint, kv, k - 1, dvals );
  
  dvals[k-1] = 0.0;         /* BasisFunctions misses this */
  
  knotScale = kv[brkPoint + 1] - kv[brkPoint];
  
  i = brkPoint - k + 1;
  for (s = k - 2; s >= 0; s--)
    {
      i++;
      omega = knotScale * ((T)(k-1)) / (kv[i+k-1] - kv[i]);
      dvals[s + 1] += -omega * dvals[s];
      dvals[s] *= omega;
    }
}

/*
 * Calculate a point p on NurbSurface<T> n at a specific u, v using the tensor product.
 * If utan and vtan are not nil, compute the u and v tangents as well.
 *
 * Note the valid parameter range for u and v is
 * (kvU[orderU] <= u < kvU[numU), (kvV[orderV] <= v < kvV[numV])
 */

template <class T>
void
CalcPoint(T u, T v, NurbSurface<T> * n,
           Point_nD<T,3> * p, Point_nD<T,3> * utan, Point_nD<T,3> * vtan)
{
  int i, j, ri, rj;
  HPoint_nD<T,3> * cp;
  T tmp;
  T wsqrdiv;
  int ubrkPoint, ufirst;
  T bu[MAXORDER], buprime[MAXORDER];
  int vbrkPoint, vfirst;
  T bv[MAXORDER], bvprime[MAXORDER];
  HPoint_nD<T,3> r, rutan, rvtan;
  
  r.x() = 0.0;
  r.y() = 0.0;
  r.z() = 0.0;
  r.w() = 0.0;
  
  rutan = r;
  rvtan = r;
  
  /* 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)
    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)
    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 - 1 - i;
          rj = n->orderU - 1 - 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)
            {
              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)
            {
              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 */
  
  wsqrdiv = 1.0 / (r.w() * r.w());
  if (utan)
    {
      utan->x() = (r.w() * rutan.x() - rutan.w() * r.x()) * wsqrdiv;
      utan->y() = (r.w() * rutan.y() - rutan.w() * r.y()) * wsqrdiv;
      utan->z() = (r.w() * rutan.z() - rutan.w() * r.z()) * wsqrdiv;
    }
  if (vtan)
    {
      vtan->x() = (r.w() * rvtan.x() - rvtan.w() * r.x()) * wsqrdiv;
      vtan->y() = (r.w() * rvtan.y() - rvtan.w() * r.y()) * wsqrdiv;
      vtan->z() = (r.w() * rvtan.z() - rvtan.w() * r.z()) * wsqrdiv;
    }
  
  p->x() = r.x() / r.w();
  p->y() = r.y() / r.w();
  p->z() = r.z() / r.w();
}

/*
 * Draw a mesh of points by evaluating the surface at evenly spaced
 * points.
 */
template <class T>
void
DrawEvaluation( NurbSurface<T> * n )
{
  Point_nD<T,3> p, utan, vtan;
  register int i, j;
  register T u, v, d;
  SurfSample<T> ** pts ;
  
  int Granularity = 10;  /* Controls the number of steps in u and v */
  
  /* Allocate storage for the grid of points generated */
  
  CHECK( pts = new (SurfSample<T>*)[Granularity+1]);
  CHECK( pts[0] = new SurfSample<T>[(Granularity+1)*(Granularity+1)]);
  
  for (i = 1; i <= Granularity; i++)
    pts[i] = &(pts[0][(Granularity+1) * i]);
  
  /* Compute points on curve */
  
  for (i = 0; i <= Granularity; i++)
    {
      v = ((T) i / (T) Granularity)
            * (n->kvV[n->numV] - n->kvV[n->orderV-1])
        + n->kvV[n->orderV-1];
      
      for (j = 0; j <= Granularity; j++)
        {
          u = ((T) j / (T) Granularity)
            * (n->kvU[n->numU] - n->kvU[n->orderU-1])
            + n->kvU[n->orderU-1];
          
          CalcPoint( u, v, n, &(pts[i][j].point), &utan, &vtan );
          p = crossProduct(utan,vtan) ; //(void) V3Cross( &utan, &vtan, &p );
          d = norm(p) ; // d = V3Length( &p );
          if (d != 0.0)
            {
              p.x() /= d;
              p.y() /= d;
              p.z() /= d;
            }
          else
            {
              p.x() = 0;
              p.y() = 0;
              p.z() = 0;
            }
          pts[i][j].normLen = d;
          pts[i][j].normal = p;
          pts[i][j].u = u;
          pts[i][j].v = v;
        }
    }
  
  /* Draw the grid */
  
  for (i = 0; i < Granularity; i++)
    for (j = 0; j < Granularity; j++)
      {
        n->render->drawTriangle( pts[i][j], pts[i+1][j+1], pts[i+1][j] );
        n->render->drawTriangle( pts[i][j], pts[i][j+1],   pts[i+1][j+1] );
      }
  
  delete [] pts[0];
  delete [] pts ;
}

/*
 * NurbRefine.c - Given a refined knot vector, add control points to a surface.
 *
 * John Peterson
 */


/*
 * Given the original knot vector ukv, and a new knotvector vkv, compute
 * the "alpha matrix" used to generate the corresponding new control points.
 * This routines allocates the alpha matrix if it isn't allocated already.
 *
 * This is from Bartels, Beatty & Barsky, p. 407
 */
template <class T> void
CalcAlpha( T * ukv, T * wkv, int m, int n, int k, T *** alpha )
{
  int i, j;
  int brkPoint, r, rm1, last, s;
  T omega;
  T aval[MAXORDER];
  
  if (! *alpha) /* Must allocate alpha */
    {
      CHECK( *alpha = new (T*)[k+1]);
      for (i = 0; i <= k; i++)
        CHECK( (*alpha)[i] = new T[(m + n + 1)]);
    }
  
  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 = minimum( 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 = minimum( 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];
    }
}

/*
 * 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).
 */
template <class T>
void
RefineSurface( NurbSurface<T> * src, NurbSurface<T> * dest, BOOL dirflag )
{
  int i, j, out;
  HPoint_nD<T,3> * dp, * sp;
  int i1, brkPoint, maxj, maxout;
  register T tmp;
  T ** alpha = 0;
  
  // Compute the alpha matrix and indexing variables for the requested direction 
  
  if (dirflag)
    {
      CalcAlpha( src->kvU, dest->kvU, src->numU - 1, dest->numU - src->numU,
                 src->orderU, &alpha );
      maxj = dest->numU;
      maxout = src->numV;
    }
  else
    {
      CalcAlpha( src->kvV, dest->kvV, src->numV - 1, dest->numV - src->numV,
                 src->orderV, &alpha );
      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 = maximum( 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 = maximum( 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++)
    delete [] alpha[i] ;
  delete [] alpha ;
}

/*
 * NurbUtils.c - Code for Allocating, freeing, & copying NURB surfaces.
 *
 * John Peterson
 */


/*
 * Allocate memory for a NURB (assumes numU, numV, orderU
 * and orderV have been set).  If ukv or vkv are not NIL, they
 * are assumed to be pointers to valid knot vectors.
 */

template <class T>
void
AllocNurb( NurbSurface<T> * n, T * ukv, T * vkv )
{
  int i;
  
  if (! ukv)
    n->kvU = new T[n->numU + n->orderU] ;
  else
    n->kvU = ukv;
  if (! vkv)
    n->kvV = new T[n->numV + n->orderV];
  else
    n->kvV = vkv;
  
  n->points.resize(n->numV,n->numU) ;
}

/*
 * Release storage for a patch
 */

template <class T>
void
FreeNurb( NurbSurface<T> * n )
{
  int i;
  
  if (n->kvU) delete [] n->kvU ;
  n->kvU = 0;
  if (n->kvV) delete [] n->kvV ;
  n->kvV = 0;
  delete n ; 
  n = 0 ; 
  // Don't touch render, it might still be used.
}

/*
 * Clone a nurb (deep copy)
 */

template <class T>
void
CloneNurb( NurbSurface<T> * src, NurbSurface<T> * dst )
{
  int i, j;
  T * srcp, *dstp;
  
  //*dst = *src;        /* Copy fields that don't change */
  dst->numU = src->numU ;
  dst->numV = src->numV ;
  dst->orderU = src->orderU ;
  dst->orderV = src->orderV ;
  
  dst->strU0 = src->strU0 ; 
  dst->strUn = src->strUn ; 
  dst->strV0 = src->strV0 ; 
  dst->strVn = src->strVn ; 
  
  dst->kvU = 0;
  dst->kvV = 0; /* So they get allocated */
  dst->points = 0;
  
  AllocNurb( dst, (T*)0, (T*)0 );
  
  /* Copy kv's */
  srcp = src->kvU;
  dstp = dst->kvU;
  for (i = 0; i < src->numU + src->orderU; i++)
    *dstp++ = *srcp++;
  
  srcp = src->kvV;
  dstp = dst->kvV;
  for (i = 0; i < src->numV + src->orderV; i++)
    *dstp++ = *srcp++;

  /* Copy control points */
  for (i = 0; i < src->numV; i++)
    for (j = 0; j < src->numU; j++)
      dst->points(i,j) = src->points(i,j);
}

/*
 * NurbSubdiv.c - Perform adaptive subdivision on a NURB surface.
 *
 * John Peterson
 *
 */

#define DIVPT( p, dn ) { ((p).x()) /= (dn); ((p).y()) /= (dn); ((p).z()) /= (dn); }

#define maxV(surf) ((surf)->numV-1)
#define maxU(surf) ((surf)->numU-1)

/*
 * 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).
 */
template <class T> int
SplitKV( T * srckv,
         T ** destkv,
         int * splitPt,    /* Where the knot interval is split */
         int m, int k )
{
  int i, last;
  int middex, extra, same;      /* "middex" ==> "middle index" */
  T midVal;
  
  extra = 0;
  last = (m + k);
  
  middex = last / 2;
  midVal = srckv[middex];
  
  /* Search forward and backward to see if multiple knot is already there */
  
  i = middex+1;
  same = 1;
  while ((i < last) && (srckv[i] == midVal)) {
    i++;
    same++;
  }
  
  i = middex-1;
  while ((i > 0) && (srckv[i] == midVal)) {
    i--;
    middex--;       /* middex is start of multiple knot */
    same++;
  }
  
  if (i <= 0)       /* No knot in middle, must create it */
    {
      midVal = (srckv[0] + srckv[last]) / 2.0;
      middex = last / 2;
      while (srckv[middex + 1] < midVal)
        middex++;
      same = 0;
    }
  
  extra = k - same;
  *destkv = new T[m+k+extra+1];
  
  if (same < k)     /* Must add knots */
    {
      for (i = 0; i <= middex; i++)
        (*destkv)[i] = srckv[i];
      
      for (i = middex+1; i <= middex+extra; i++)
        (*destkv)[i] = midVal;
      
      for (i = middex + k - same + 1; i <= m + k + extra; i++)
        (*destkv)[i] = srckv[i - extra];
    }
  else
    {
      for (i = 0; i <= m + k; i++)
        (*destkv)[i] = srckv[i];
    }
  
  *splitPt = (extra < k) ? middex - 1 : middex;
  return( extra );
}

/*
 * Given a line defined by firstPt and lastPt, project midPt onto
 * that line.  Used for fixing "cracks".
 */
template <class T> void
ProjectToLine( Point_nD<T,3> * firstPt, Point_nD<T,3> * lastPt, Point_nD<T,3> * midPt )
{
  Point_nD<T,3> base, v0, vm;
  T fraction, denom;
  
  base = *firstPt;
  
  v0 = *lastPt - base ; // (void) V3Sub( lastPt, &base, &v0 );
  vm = *midPt - base ; // (void) V3Sub( midPt, &base, &vm );
  
  denom = norm2(v0) ; // V3SquaredLength( &v0 );
  //fraction = (denom == 0.0) ? 0.0 : (V3Dot( &v0, &vm ) / denom);
  fraction = (denom == 0.0) ? 0.0 : (v0*vm ) / denom;
  
  midPt->x() = base.x() + fraction * v0.x();
  midPt->y() = base.y() + fraction * v0.y();
  midPt->z() = base.z() + fraction * v0.z();
}

/*
 * 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...
 */
template <class T> void
FixNormals( SurfSample<T> * s0, SurfSample<T> * s1, SurfSample<T> * s2 )
{
  BOOL goodnorm;
  int i, j, ok;
  T dist;
  SurfSample<T> * V[3];
  Point_nD<T,3> normal;
  
  V[0] = s0; V[1] = s1; V[2] = s2;
  
  /* Find a reasonable normalal */
  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 */
    {
      normal.x() = 0.0; normal.y() = 0.0; normal.z() = 0.0;
      
      for (i = 0; i < 3L; i++)
        {
          j = (i + 1) % 3L;
          normal.x() += (V[i]->point.y() - V[j]->point.y()) * (V[i]->point.z() + V[j]->point.z());
          normal.y() += (V[i]->point.z() - V[j]->point.z()) * (V[i]->point.x() + V[j]->point.x());
          normal.z() += (V[i]->point.x() - V[j]->point.x()) * (V[i]->point.y() + V[j]->point.y());
        }
      //dist = V3Length( &norm );
      dist = norm(normal) ; 
      if (dist == 0.0)
        return;         /* This sucker's hopeless... */
      
      DIVPT( normal, dist );
      
      for (i = 0; i < 3; i++)
        {
          V[i]->normal = 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))
          V[i]->normal = V[ok]->normal;
    }
  return;
}

/*
 * Normalize the normal in a sample.  If it's degenerate,
 * flag it as such by setting the normLen to 0.0
 */
template <class T> void
AdjustNormal( SurfSample<T> * samp )
{
  // If it's not degenerate, do the normalization now */
  samp->normLen = norm(samp->normal) ; // V3Length( &(samp->normal) );

  if (samp->normLen < samp->epsilon )
    samp->normLen = 0.0;
  else
    DIVPT( (samp->normal), samp->normLen );
}

/*
 * 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).
 */
template <class T> void
GetNormal( NurbSurface<T> * n, int indV, int indU )
{
  Point_nD<T,3> tmpL, tmpR;     /* "Left" and "Right" of the base point */
  SurfSample<T> * crnr;
  
  if ( (indU && indV) || ((! indU) && (!indV)) )
    {
      if (indU)
        crnr = &(n->cnn);
      else
        crnr = &(n->c00);
      DIVW( &(n->points(indV,(indU ? (indU-1) : 1))), &tmpL );
      DIVW( &(n->points((indV ? (indV-1) : 1),indU)), &tmpR );
    }
  else
    {
      if (indU)
        crnr = &(n->c0n);
      else
        crnr = &(n->cn0);
      DIVW( &(n->points(indV,(indU ? (indU-1) : 1))), &tmpR );
      DIVW( &(n->points((indV ? (indV-1) : 1),indU)), &tmpL );
    }
  
  tmpL -= crnr->point ; //(void) V3Sub( &tmpL, &(crnr->point), &tmpL );
  tmpR -= crnr->point ;//(void) V3Sub( &tmpR, &(crnr->point), &tmpR );
  crnr->normal = crossProduct(tmpL,tmpR); //(void) V3Cross( &tmpL, &tmpR, &(crnr->normal) );
  AdjustNormal( crnr );
}

/*
 * Build the new corners in the two new surfaces, computing both
 * point on the surface aint with the normal.   Prevent cracks that may occur.
 */
template <class T> void
MakeNewCorners( NurbSurface<T> * parent,
                NurbSurface<T> * kid0,
                NurbSurface<T> * kid1,
                BOOL dirflag )
{
  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(0,maxU(kid0))), &(kid0->c0n.point) );
      GetNormal( kid0, 0, maxU(kid0) );
      /*
       * Normals must be re-calculated for kid1 in case the surface
       * was split at a c1 (or c0!) discontinutiy
       */
      kid1->c00.point = kid0->c0n.point;
      GetNormal( kid1, 0, 0 );
      kid1->cn0.point = kid0->cnn.point;
      GetNormal( kid1, maxV(kid1), 0 );
      
      /*
       * Prevent cracks from forming by forcing the points on the seam to
       * lie aint 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) );
      
      kid1->c00.point = kid0->c0n.point;
      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), 0 );
      kid1->c00.point = kid0->cn0.point;
      GetNormal( kid1, 0, 0 );
      kid1->c0n.point = kid0->cnn.point;
      GetNormal( kid1, 0, 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) );
      
      kid1->c00.point = kid0->cn0.point;
      kid1->c0n.point = kid0->cnn.point;
      kid1->strV0 = FALSE_;
    }
}

/*
 * 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.
 */
template <class T> void
SplitSurface( NurbSurface<T> * parent,
              NurbSurface<T> * kid0, NurbSurface<T> * kid1,
              BOOL dirflag )        /* If TRUE_ subdivided in U, else in V */
{
  NurbSurface<T> *tmp;
  T * newkv;
  int i, j, splitPt;
  
  tmp = new NurbSurface<T> ;

  //
  // 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. 
  tmp->numU = parent->numU ;
  tmp->numV = parent->numV ;
  tmp->orderU = parent->orderU ;
  tmp->orderV = parent->orderV ;
  
  tmp->strU0 = parent->strU0 ; 
  tmp->strUn = parent->strUn ; 
  tmp->strV0 = parent->strV0 ; 
  tmp->strVn = parent->strVn ; 
  
  tmp->render = parent->render ; 
  
  if (dirflag)
    {
      tmp->numU = parent->numU + SplitKV( parent->kvU,
                                         &newkv,
                                         &splitPt,
                                         maxU(parent),
                                         parent->orderU );
      AllocNurb( tmp, newkv, (T*)0 );
      for (i = 0; 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 );
      AllocNurb( tmp, (T*)0, newkv );
      for (i = 0; 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->orderU = parent->orderU ;
  kid0->orderV = parent->orderV ;
  
  kid0->strU0 = parent->strU0 ; 
  kid0->strUn = parent->strUn ; 
  kid0->strV0 = parent->strV0 ; 
  kid0->strVn = parent->strVn ; 
  
  kid0->c00 = parent->c00 ;
  kid0->c0n = parent->c0n ;
  kid0->cn0 = parent->cn0 ;
  kid0->cnn = parent->cnn ;
  
  kid0->render = parent->render ; 
  
  kid0->numU = dirflag ? splitPt+1 : parent->numU;
  kid0->numV = dirflag ? parent->numV : splitPt+1;
  kid0->kvU = kid0->kvV = 0;
  AllocNurb( kid0, (T*)0, (T*)0 );
  
  for (i = 0; i < kid0->numV; i++)      // Copy the point and kv data 
    for (j = 0; j < kid0->numU; j++)
      kid0->points(i,j) = tmp->points(i,j) ;
  for (i = 0; i < kid0->orderU + kid0->numU; i++)
    kid0->kvU[i] = tmp->kvU[i];
  for (i = 0; i < kid0->orderV + kid0->numV; i++)
    kid0->kvV[i] = tmp->kvV[i];
  
  // Second half 
  
  splitPt++;
  //*kid1 = *parent;
  kid1->orderU = parent->orderU ;
  kid1->orderV = parent->orderV ;
  
  kid1->strU0 = parent->strU0 ; 
  kid1->strUn = parent->strUn ; 
  kid1->strV0 = parent->strV0 ; 
  kid1->strVn = parent->strVn ; 

  kid1->c00 = parent->c00 ;
  kid1->c0n = parent->c0n ;
  kid1->cn0 = parent->cn0 ;
  kid1->cnn = parent->cnn ;
  
  kid1->render = parent->render ; 
  
  kid1->numU = dirflag ? tmp->numU - splitPt : parent->numU;
  kid1->numV = dirflag ? parent->numV : tmp->numV - splitPt;
  kid1->kvU = kid1->kvV = 0;
  AllocNurb( kid1, (T*)0, (T*)0 );
  
  for (i = 0; i < kid1->numV; i++)      // Copy the point and kv data 
    for (j = 0; j < kid1->numU; j++)
      kid1->points(i,j)
        = tmp->points(dirflag ? i: (i + splitPt) ,dirflag ? (j + splitPt) : j);
  for (i = 0; i < kid1->orderU + kid1->numU; i++)
    kid1->kvU[i] = tmp->kvU[dirflag ? (i + splitPt) : i];
  for (i = 0; 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 );
  
  FreeNurb( tmp );          // Get rid of refined parent 
    
}

/*
 * 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.
 */

#define GETPT( i )  (( dirflag ? (n->points(crvInd,i)) : (n->points(i,crvInd)) ))

#define EPSILON n->epsilon

template <class T> BOOL
IsCurveStraight( NurbSurface<T> * n,
                 T tolerance,
                 int crvInd,
                 BOOL dirflag )  /* If TRUE_, test in U direction, else test in V */
{
  Point_nD<T,3> p, vec, prod;
  Point_nD<T,3> cp, e0;
  int i, last;
  T linelen, dist;
  
  /* Special case: lines are automatically straight. */
  if ((dirflag ? n->numU : n->numV) == 2)
    return( TRUE_ );
  
  last = (dirflag ? n->numU : n->numV) - 1;
  n->render->screenProject( GETPT( 0 ), e0 );
  
  /* Form an initial line to test the other points against (skiping degen lines) */
  
  linelen = 0.0;
  for (i = last; (i > 0) && (linelen < EPSILON); i--)
    {
      n->render->screenProject( GETPT( i ), cp );
      vec = cp - e0 ; 
      linelen = norm(vec) ;
    }
  
  DIVPT( vec, linelen );

  if (linelen > EPSILON)        /* If no non-degenerate lines found, it's all degen */
    for (i = 1; 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. */
        n->render->screenProject( GETPT( i ), cp );
        p = cp - e0 ; 
        
        prod = crossProduct(p,vec) ; 
        dist = norm(prod) ; 
        
        if (dist > tolerance)
          return( FALSE_ );
      }
  
  return( TRUE_ );
}

/*
 * 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.
 */
template <class T> BOOL
TestFlat( NurbSurface<T> * n, T tolerance )
{
  int i;
  BOOL straight;
  Point_nD<T,3> cp00, cp0n, cpn0, cpnn, planeEqn;
  T dist,d ;
  
  /* Check edge straightness */
  
  if (! n->strU0)
    n->strU0 = IsCurveStraight( n, tolerance, 0, FALSE_ );
  if (! n->strUn)
    n->strUn = IsCurveStraight( n, tolerance, maxU(n), FALSE_ );
  if (! n->strV0)
    n->strV0 = IsCurveStraight( n, tolerance, 0, 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 = 1;
         (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 = 1;
         (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. 
  
  n->render->screenProject( (n->points(0,0)),       cp00 );
  n->render->screenProject( (n->points(0,maxU(n))),         cp0n );
  n->render->screenProject( (n->points(maxV(n),0)),         cpn0 );
  n->render->screenProject( (n->points(maxV(n),maxU(n))),  cpnn );
  
  cp0n -= cp00 ; // Make edges into vectors
  cpn0 -= cp00 ; 
  
  
  // 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.
  
  planeEqn = crossProduct(cpn0,cp0n) ; 
  planeEqn = planeEqn.unitLength() ; // Normalize to keep adds in sync w/ mults 
  
  d = planeEqn * cp00 ; 
  dist = fabs( ( planeEqn * cpnn ) - d );
  
  if ( dist > tolerance ) // Surface is twisted 
    return( FALSE_ );
  else
    return( TRUE_ );
}

/*
 * Turn a sufficiently flat surface into triangles.
 */
template <class T> 
void EmitTriangles( NurbSurface<T> * n )
{
  Point_nD<T,3> vecnn, vec0n;   // Diagonal vectors 
  T len2nn, len20n;             // Diagonal lengths squared 
  T u0, un, v0, vn;             // Texture coords;
                                   
  //
  // Measure the distance aint the two diagonals to decide the best
  // way to cut the rectangle into triangles.
  //
  
  vecnn = n->c00.point - n->cnn.point ; 
  vec0n = n->c0n.point - n->cn0.point ; 
  
  len2nn = norm2(vecnn) ; 
  len20n = norm2(vec0n) ; 
  
  if (maximum(len2nn, len20n) < n->epsilon)
    return;                             // Triangles are too small to render 
  
  //
  // Assign the texture coordinates
  //
  u0 = n->kvU[n->orderU-1];
  un = n->kvU[n->numU];
  v0 = n->kvV[n->orderV-1];
  vn = n->kvV[n->numV];
  n->c00.u = u0; n->c00.v = v0;
  n->c0n.u = un; n->c0n.v = v0;
  n->cn0.u = u0; n->cn0.v = vn;
  n->cnn.u = un; n->cnn.v = vn;
  
  //
  // If any normals are sick, fix them now.
  //
  if ((n->c00.normLen == 0.0) || (n->cnn.normLen == 0.0) || (n->cn0.normLen == 0.0))
    FixNormals( &(n->c00), &(n->cnn), &(n->cn0) );
  if (n->c0n.normLen == 0.0)
    FixNormals( &(n->c00), &(n->c0n), &(n->cnn) );
  
  if ( len2nn < len20n )
    {
      n->render->drawTriangle( n->c00, n->cnn, n->cn0 );
      n->render->drawTriangle( n->c00, n->c0n, n->cnn );
    }                                            
  else                                   
    {                                            
      n->render->drawTriangle( n->c0n, n->cnn, n->cn0 );
      n->render->drawTriangle( n->c0n, n->cn0, n->c00 );
    }
}

/*
 * 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.
 */
template <class T> void
DoSubdivision( NurbSurface<T> * n, T tolerance, BOOL dirflag, int level )
{
  NurbSurface<T> *left, *right;
  
  left = new NurbSurface<T>;
  right = new NurbSurface<T>;

  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 + 1 );
      DoSubdivision( right, tolerance, dirflag, level + 1 );
      FreeNurb( left );
      FreeNurb( right );
    }
}

/*
 * Main entry point for subdivision */
template <class T> void
DrawSubdivision( NurbSurface<T> * surf, T tolerance )
{
  surf->flatV = FALSE_;
  surf->flatU = FALSE_;
  surf->strU0 = FALSE_;
  surf->strUn = FALSE_;
  surf->strV0 = FALSE_;
  surf->strVn = FALSE_;
  
  //
  // Initialize the projected corners of the surface
  // and the normals.
  //
  DIVW( &(surf->points(0,0)), &surf->c00.point );
  DIVW( &(surf->points(0,surf->numU-1)), &surf->c0n.point );
  DIVW( &(surf->points(surf->numV-1,0)), &surf->cn0.point );
  DIVW( &(surf->points(surf->numV-1,surf->numU-1)), &surf->cnn.point );

  GetNormal( surf, 0, 0 );
  GetNormal( surf, 0, maxU(surf) );
  GetNormal( surf, maxV(surf), 0 );
  GetNormal( surf, maxV(surf), maxU(surf) );
  
  RenderMesh<T> *render ;

  render = surf->render ;
  render->drawHeader();
  DoSubdivision( surf, tolerance, TRUE_, 0 );
  // Note surf is deallocated by the subdivision process 
  render->drawFooter();
}

template <class T>
SurfSample<T>& SurfSample<T>::operator=(const SurfSample<T>& s) {
  point = s.point ;
  normal = s.normal ;
  normLen = s.normLen ;
  u = s.u ;
  v = s.v ;
  return *this ;
}

/*
template <class T>
class NurbsCurveTess : public NurbsCurve<T,2> {
public:
  NurbsCurveTess(const NurbsSurface<T,3>& rs) ;

  void tesselate(T tol) ; 

protected:
  int isStraight() ; 

  const NurbsSurface<T,3>& rsurf ; 
  T tolerance ; 
}

template <class T>
NurbsCurveTess::NurbsCurveTess(const NurbsSurface<T,3>& rs) : rsurf(rs) {
  tolerance = 0.01 ; 
}

template <class T>
void NurbsCurveTess::tesselate(T tol){
  tolerance = tol ; 
  if(knot.n() != 2*(deg_+1)){
    
  }
}
*/

}

#endif // #define PLIB_NURBS_NURBSSUB_SOURCE

---