C++ ON_NurbsSurface::SetCV方法代码示例

static ON_Surface* CreatePlanarSurface( 
                             const ON_3dPoint& SW, const ON_3dPoint& SE,
                             const ON_3dPoint& NE, const ON_3dPoint& NW
                             )
{
  ON_NurbsSurface* pNurbsSurface = new ON_NurbsSurface(
                                        3,     // dimension (>= 1)
                                        FALSE, // not rational
                                        2,     // "u" order (>= 2)
                                        2,     // "v" order (>= 2)
                                        2,     // number of control vertices in "u" dir (>= order)
                                        2      // number of control vertices in "v" dir (>= order)
                                        );
  // corner CVs in counter clockwise order starting in the south west
  pNurbsSurface->SetCV( 0,0, SW );
  pNurbsSurface->SetCV( 1,0, SE );
  pNurbsSurface->SetCV( 1,1, NE );
  pNurbsSurface->SetCV( 0,1, NW );
  // "u" knots
  pNurbsSurface->SetKnot( 0,0, 0.0 );
  pNurbsSurface->SetKnot( 0,1, 1.0 );
  // "v" knots
  pNurbsSurface->SetKnot( 1,0, 0.0 );
  pNurbsSurface->SetKnot( 1,1, 1.0 );

  return pNurbsSurface;
}

static ON_Surface* TwistedCubeSideSurface( 
                             const ON_3dPoint& SW, const ON_3dPoint& SE,
                             const ON_3dPoint& NE, const ON_3dPoint& NW
                             )
{
  ON_NurbsSurface* pNurbsSurface = new ON_NurbsSurface(
                                        3,     // dimension
                                        false, // not rational
                                        2,     // "u" order
                                        2,     // "v" order
                                        2,     // number of control vertices in "u" dir
                                        2      // number of control vertices in "v" dir
                                        );
  // corner CVs in counter clockwise order starting in the south west
  pNurbsSurface->SetCV( 0,0, SW );
  pNurbsSurface->SetCV( 1,0, SE );
  pNurbsSurface->SetCV( 1,1, NE );
  pNurbsSurface->SetCV( 0,1, NW );
  // "u" knots
  pNurbsSurface->SetKnot( 0,0, 0.0 );
  pNurbsSurface->SetKnot( 0,1, 1.0 );
  // "v" knots
  pNurbsSurface->SetKnot( 1,0, 0.0 );
  pNurbsSurface->SetKnot( 1,1, 1.0 );

  return pNurbsSurface;
}

int 
ON_PlaneSurface::GetNurbForm( // returns 0: unable to create NURBS representation
                   //            with desired accuracy.
                   //         1: success - returned NURBS parameterization
                   //            matches the surface's to wthe desired accuracy
                   //         2: success - returned NURBS point locus matches
                   //            the surfaces's to the desired accuracy but, on
                   //            the interior of the surface's domain, the 
                   //            surface's parameterization and the NURBS
                   //            parameterization may not match to the 
                   //            desired accuracy.
        ON_NurbsSurface& nurbs,
        double tolerance
        ) const
{
  ON_BOOL32 rc = IsValid();

  if( !rc )
  {
    if (    m_plane.origin.x != ON_UNSET_VALUE 
         && m_plane.xaxis.x != ON_UNSET_VALUE 
         && m_plane.yaxis.x != ON_UNSET_VALUE
         && m_domain[0].IsIncreasing() && m_domain[1].IsIncreasing()
         && m_extents[0].Length() > 0.0 && m_extents[1].Length() > 0.0
         )
    {
      ON_3dVector N = ON_CrossProduct(m_plane.xaxis,m_plane.yaxis);
      if ( N.Length() <= 1.0e-4 )
      {
        ON_WARNING("ON_PlaneSurface::GetNurbForm - using invalid surface.");
        rc = true;
      }
    }
  }

  if ( rc ) 
  {
    nurbs.m_dim = 3;
    nurbs.m_is_rat = 0;
    nurbs.m_order[0] = nurbs.m_order[1] = 2;
    nurbs.m_cv_count[0] = nurbs.m_cv_count[1] = 2;
    nurbs.m_cv_stride[1] = nurbs.m_dim;
    nurbs.m_cv_stride[0] = nurbs.m_cv_stride[1]*nurbs.m_cv_count[1];
    nurbs.ReserveCVCapacity(12);
    nurbs.ReserveKnotCapacity(0,2);
    nurbs.ReserveKnotCapacity(1,2);
    nurbs.m_knot[0][0] = m_domain[0][0];
    nurbs.m_knot[0][1] = m_domain[0][1];
    nurbs.m_knot[1][0] = m_domain[1][0];
    nurbs.m_knot[1][1] = m_domain[1][1];
    nurbs.SetCV( 0, 0, PointAt( m_domain[0][0], m_domain[1][0] ));
    nurbs.SetCV( 0, 1, PointAt( m_domain[0][0], m_domain[1][1] ));
    nurbs.SetCV( 1, 0, PointAt( m_domain[0][1], m_domain[1][0] ));
    nurbs.SetCV( 1, 1, PointAt( m_domain[0][1], m_domain[1][1] ));
  }

  return rc;
}

ON_NurbsSurface
FittingSurface::initNurbsPCA (int order, NurbsDataSurface *m_data, ON_3dVector z)
{
  ON_3dVector mean;
  Eigen::Matrix3d eigenvectors;
  Eigen::Vector3d eigenvalues;

  unsigned s = m_data->interior.size ();

  NurbsTools::pca (m_data->interior, mean, eigenvectors, eigenvalues);

  m_data->mean = mean;
  //m_data->eigenvectors = (*eigenvectors);

  bool flip (false);
  Eigen::Vector3d ez(z[0],z[1],z[2]);
  if (eigenvectors.col (2).dot (ez) < 0.0)
     flip = true;

  eigenvalues = eigenvalues / s; // seems that the eigenvalues are dependent on the number of points (???)

  ON_3dVector sigma(sqrt(eigenvalues[0]), sqrt(eigenvalues[1]), sqrt(eigenvalues[2]));

  ON_NurbsSurface nurbs (3, false, order, order, order, order);
  nurbs.MakeClampedUniformKnotVector (0, 1.0);
  nurbs.MakeClampedUniformKnotVector (1, 1.0);

  // +- 2 sigma -> 95,45 % aller Messwerte
  double dcu = (4.0 * sigma[0]) / (nurbs.Order (0) - 1);
  double dcv = (4.0 * sigma[1]) / (nurbs.Order (1) - 1);

  ON_3dVector cv_t, cv;
  Eigen::Vector3d ecv_t, ecv;
  Eigen::Vector3d emean(mean[0], mean[1], mean[2]);
  for (int i = 0; i < nurbs.Order (0); i++)
  {
    for (int j = 0; j < nurbs.Order (1); j++)
    {
      cv[0] = -2.0 * sigma[0] + dcu * i;
      cv[1] = -2.0 * sigma[1] + dcv * j;
      cv[2] = 0.0;
      ecv (0) = -2.0 * sigma[0] + dcu * i;
      ecv (1) = -2.0 * sigma[1] + dcv * j;
      ecv (2) = 0.0;
      ecv_t = eigenvectors * ecv + emean;
      cv_t[0] = ecv_t (0);
      cv_t[1] = ecv_t (1);
      cv_t[2] = ecv_t (2);
      if (flip)
	nurbs.SetCV (nurbs.Order (0) - 1 - i, j, ON_3dPoint (cv_t[0], cv_t[1], cv_t[2]));
      else
	nurbs.SetCV (i, j, ON_3dPoint (cv_t[0], cv_t[1], cv_t[2]));
    }
  }
  return nurbs;
}

HIDDEN ON_Surface*
sideSurface(const ON_3dPoint& SW, const ON_3dPoint& SE, const ON_3dPoint& NE, const ON_3dPoint& NW)
{
    ON_NurbsSurface *surf = ON_NurbsSurface::New(3, 0, 2, 2, 2, 2);
    surf->SetCV(0, 0, SW);
    surf->SetCV(1, 0, SE);
    surf->SetCV(1, 1, NE);
    surf->SetCV(0, 1, NW);
    surf->SetKnot(0, 0, 0.0);
    surf->SetKnot(0, 1, 1.0);
    surf->SetKnot(1, 0, 0.0);
    surf->SetKnot(1, 1, 1.0);
    return surf;
}

void
GlobalOptimizationTDM::updateSurf (double damp)
{
  int ncps (0);

  for (unsigned i = 0; i < m_nurbs.size (); i++)
  {
    ON_NurbsSurface* nurbs = m_nurbs[i];

    int ncp = nurbs->CVCount ();

    for (int A = 0; A < ncp; A++)
    {

      int I = gl2gr (*nurbs, A);
      int J = gl2gc (*nurbs, A);

      ON_3dPoint cp_prev;
      nurbs->GetCV (I, J, cp_prev);

      ON_3dPoint cp;
      cp.x = cp_prev.x + damp * (m_solver.x (3 * (ncps + A) + 0, 0) - cp_prev.x);
      cp.y = cp_prev.y + damp * (m_solver.x (3 * (ncps + A) + 1, 0) - cp_prev.y);
      cp.z = cp_prev.z + damp * (m_solver.x (3 * (ncps + A) + 2, 0) - cp_prev.z);

      nurbs->SetCV (I, J, cp);
    }

    ncps += ncp;
  }
}

ON_NurbsSurface
FittingSurface::initNurbsPCA (int order, NurbsDataSurface *m_data, Eigen::Vector3d z)
{
  Eigen::Vector3d mean;
  Eigen::Matrix3d eigenvectors;
  Eigen::Vector3d eigenvalues;

  unsigned s = m_data->interior.size ();

  NurbsTools::pca (m_data->interior, mean, eigenvectors, eigenvalues);

  m_data->mean = mean;
  m_data->eigenvectors = eigenvectors;

  bool flip (false);
  if (eigenvectors.col (2).dot (z) < 0.0)
    flip = true;

  eigenvalues = eigenvalues / s; // seems that the eigenvalues are dependent on the number of points (???)

  Eigen::Vector3d sigma (sqrt (eigenvalues (0)), sqrt (eigenvalues (1)), sqrt (eigenvalues (2)));

  ON_NurbsSurface nurbs (3, false, order, order, order, order);
  nurbs.MakeClampedUniformKnotVector (0, 1.0);
  nurbs.MakeClampedUniformKnotVector (1, 1.0);

  // +- 2 sigma -> 95,45 % aller Messwerte
  double dcu = (4.0 * sigma (0)) / (nurbs.Order (0) - 1);
  double dcv = (4.0 * sigma (1)) / (nurbs.Order (1) - 1);

  Eigen::Vector3d cv_t, cv;
  for (int i = 0; i < nurbs.Order (0); i++)
  {
    for (int j = 0; j < nurbs.Order (1); j++)
    {
      cv (0) = -2.0 * sigma (0) + dcu * i;
      cv (1) = -2.0 * sigma (1) + dcv * j;
      cv (2) = 0.0;
      cv_t = eigenvectors * cv + mean;
      if (flip)
        nurbs.SetCV (nurbs.Order (0) - 1 - i, j, ON_3dPoint (cv_t (0), cv_t (1), cv_t (2)));
      else
        nurbs.SetCV (i, j, ON_3dPoint (cv_t (0), cv_t (1), cv_t (2)));
    }
  }
  return nurbs;
}

ON_NurbsSurface
FittingCylinder::initNurbsCylinderWithAxes (int order, NurbsDataSurface *data, Eigen::Matrix3d &axes)
{
  Eigen::Vector3d mean;

  unsigned s = unsigned (data->interior.size ());
  mean = NurbsTools::computeMean (data->interior);

  data->mean = mean;

  Eigen::Vector3d v_max (0.0, 0.0, 0.0);
  Eigen::Vector3d v_min (DBL_MAX, DBL_MAX, DBL_MAX);
  for (unsigned i = 0; i < s; i++)
  {
    Eigen::Vector3d p (axes.inverse () * (data->interior[i] - mean));

    if (p (0) > v_max (0))
      v_max (0) = p (0);
    if (p (1) > v_max (1))
      v_max (1) = p (1);
    if (p (2) > v_max (2))
      v_max (2) = p (2);

    if (p (0) < v_min (0))
      v_min (0) = p (0);
    if (p (1) < v_min (1))
      v_min (1) = p (1);
    if (p (2) < v_min (2))
      v_min (2) = p (2);
  }

  int ncpsU (order);
  int ncpsV (2 * order + 4);
  ON_NurbsSurface nurbs (3, false, order, order, ncpsU, ncpsV);
  nurbs.MakeClampedUniformKnotVector (0, 1.0);
  nurbs.MakePeriodicUniformKnotVector (1, 1.0 / (ncpsV - order + 1));

  double dcu = (v_max (0) - v_min (0)) / (ncpsU - 1);
  double dcv = (2.0 * M_PI) / (ncpsV - order + 1);

  double ry = std::max<double> (std::fabs (v_min (1)), std::fabs (v_max (1)));
  double rz = std::max<double> (std::fabs (v_min (2)), std::fabs (v_max (2)));

  Eigen::Vector3d cv_t, cv;
  for (int i = 0; i < ncpsU; i++)
  {
    for (int j = 0; j < ncpsV; j++)
    {
      cv (0) = v_min (0) + dcu * i;
      cv (1) = ry * sin (dcv * j);
      cv (2) = rz * cos (dcv * j);
      cv_t = axes * cv + mean;
      nurbs.SetCV (i, j, ON_3dPoint (cv_t (0), cv_t (1), cv_t (2)));
    }
  }

  return nurbs;
}

extern "C" void
rt_nurb_brep(ON_Brep **b, const struct rt_db_internal *ip, const struct bn_tol *)
{
    int i, j, k;
    struct rt_nurb_internal *nip;

    RT_CK_DB_INTERNAL(ip);
    nip = (struct rt_nurb_internal *)ip->idb_ptr;
    RT_NURB_CK_MAGIC(nip);

    ON_TextLog log(stderr);

    for (i = 0; i < nip->nsrf; i++) {
	struct face_g_snurb *surface = nip->srfs[i];
	NMG_CK_SNURB(surface);

	ON_NurbsSurface *nurb = ON_NurbsSurface::New(3, true, surface->order[0], surface->order[1], surface->s_size[0], surface->s_size[1]);

	/* set 'u' knots */
	/* skip first and last (duplicates?) */
	for (j = 1; j < surface->u.k_size - 1; j++) {
	    nurb->SetKnot(0, j-1, surface->u.knots[j]);
	    /* bu_log("u knot %d is %f\n", j-1, surface->u.knots[j]); */
	}
	/* set 'v' knots */
	/* skip first and last (duplicates?) */
	for (j = 1; j < surface->v.k_size - 1; j++) {
	    nurb->SetKnot(1, j-1, surface->v.knots[j]);
	    /* bu_log("v knot %d is %f\n", j-1, surface->u.knots[j]); */
	}

	/* set control points */
	for (j = 0; j < surface->s_size[0]; j++) {
	    for (k = 0; k < surface->s_size[1]; k++) {
		ON_3dPoint point = &RT_NURB_GET_CONTROL_POINT(surface, j, k);
		nurb->SetCV(k, j, point);
	    }
	}

	/* nurb->Dump(log); */
	bu_log("NURBS surface %d %s valid\n", i, nurb->IsValid(&log) ? "is" : "is not");

	(*b)->m_S.Append(nurb);
	int sindex = (*b)->m_S.Count();
	(*b)->NewFace(sindex - 1);
	int findex = (*b)->m_F.Count();
	(*b)->NewOuterLoop(findex - 1);
    }

    bu_log("BREP object %s a single surface\n", (*b)->IsSurface() ? "is" : "is not");
    bu_log("BREP object %s valid\n", (*b)->IsValid(&log) ? "is" : "is not");
    bu_log("BREP object %s valid topology\n", (*b)->IsValidTopology(&log) ? "is" : "is not");
    bu_log("BREP object %s valid geometry\n", (*b)->IsValidGeometry(&log) ? "is" : "is not");
    bu_log("BREP object %s solid\n", (*b)->IsSolid() ? "is" : "is not");
    bu_log("BREP object %s manifold\n", (*b)->IsManifold() ? "is" : "is not");
}

int ON_Cone::GetNurbForm( ON_NurbsSurface& s ) const
{
  int rc = 0;
  if ( IsValid() ) {
    ON_Circle c = CircleAt(height);
    ON_NurbsCurve n;
    c.GetNurbForm(n);
    ON_3dPoint apex = ApexPoint();
    ON_4dPoint cv;
    int i, j0, j1;

    s.Create(3,TRUE,3,2,9,2);
    for ( i = 0; i < 10; i++ )
      s.m_knot[0][i] = n.m_knot[i];

    if ( height >= 0.0 ) {
      s.m_knot[1][0] = 0.0;
      s.m_knot[1][1] = height;
      j0 = 0;
      j1 = 1;
    }
    else {
      s.m_knot[1][0] = height;
      s.m_knot[1][1] = 0.0;
      j0 = 1;
      j1 = 0;
    }

    for ( i = 0; i < 9; i++ ) {
      cv = n.CV(i);
      s.SetCV(i, j1, ON::homogeneous_rational, &cv.x );
      cv.x = apex.x*cv.w;
      cv.y = apex.y*cv.w;
      cv.z = apex.z*cv.w;
      s.SetCV(i, j0, cv);
    }
    rc = 2;
  }
  return rc;
}

int ON_Cylinder::GetNurbForm( ON_NurbsSurface& s ) const
{
  int rc = 0;
  if ( IsValid() && height[0] != height[1] ) {
    ON_NurbsCurve n0, n1;
    int i;
    ON_Circle c0 = CircleAt(height[0]);
    ON_Circle c1 = CircleAt(height[1]);


    if ( height[0] <= height[1] ) {
      c0.GetNurbForm(n0);
      c1.GetNurbForm(n1);
    }
    else {
      c0.GetNurbForm(n1);
      c1.GetNurbForm(n0);
    }

    if (    n0.m_dim != n1.m_dim
         || n0.m_is_rat != n1.m_is_rat
         || n0.m_order != n1.m_order
         || n0.m_cv_count != n1.m_cv_count )
      return 0;

    s.Create(3,TRUE, n0.m_order, 2, n0.m_cv_count, 2 );
    if ( height[0] <= height[1] ) {
      s.m_knot[1][0] = height[0];
      s.m_knot[1][1] = height[1];
    }
    else {
      s.m_knot[1][0] = height[1];
      s.m_knot[1][1] = height[0];
    }

    for ( i = 0; i < n0.KnotCount(); i++ )
      s.m_knot[0][i] = n0.m_knot[i];

    for ( i = 0; i < n0.m_cv_count; i++ ) {
      s.SetCV(i,0,ON::homogeneous_rational,n0.CV(i));
      s.SetCV(i,1,ON::homogeneous_rational,n1.CV(i));
    }
    rc = 2;
  }
  return rc;
}

ON_NurbsSurface
FittingSurface::initNurbs4Corners (int order, ON_3dPoint ll, ON_3dPoint lr, ON_3dPoint ur, ON_3dPoint ul)
{
  std::map<int, std::map<int, ON_3dPoint> > cv_map;

  double dc = 1.0 / (order - 1);

  for (int i = 0; i < order; i++)
  {
    double di = dc * i;
    cv_map[i][0] = ll + (lr - ll) * di;
    cv_map[0][i] = ll + (ul - ll) * di;
    cv_map[i][order - 1] = ul + (ur - ul) * di;
    cv_map[order - 1][i] = lr + (ur - lr) * di;
  }

  for (int i = 1; i < order - 1; i++)
  {
    for (int j = 1; j < order - 1; j++)
    {
      ON_3dPoint du = cv_map[0][j] + (cv_map[order - 1][j] - cv_map[0][j]) * dc * i;
      ON_3dPoint dv = cv_map[i][0] + (cv_map[i][order - 1] - cv_map[i][0]) * dc * j;
      cv_map[i][j] = du * 0.5 + dv * 0.5;
    }
  }

  ON_NurbsSurface nurbs (3, false, order, order, order, order);
  nurbs.MakeClampedUniformKnotVector (0, 1.0);
  nurbs.MakeClampedUniformKnotVector (1, 1.0);

  for (int i = 0; i < order; i++)
  {
    for (int j = 0; j < order; j++)
    {
      nurbs.SetCV (i, j, cv_map[i][j]);
    }
  }
  return nurbs;
}

//.........这里部分代码省略.........
    bp->SetExtents(1, bp->Domain(1));
    (*b)->SetTrimIsoFlags(bface);
    delete ellcurve1;

    //  Now, the hard part.  Need an elliptical hyperbolic NURBS surface
    //  First step is to create a nurbs curve.

    double intercept_calc = (eip->ehy_c)*(eip->ehy_c)/(MAGNITUDE(eip->ehy_H) + eip->ehy_c);
    double intercept_dist = MAGNITUDE(eip->ehy_H) + eip->ehy_c - intercept_calc;
    double intercept_length = intercept_dist - MAGNITUDE(eip->ehy_H);
    double MX = MAGNITUDE(eip->ehy_H);
    double MP = MX + intercept_length;
    double w = (MX/MP)/(1-MX/MP);

    point_t ep1, ep2, ep3;
    VSET(ep1, -eip->ehy_r1, 0, 0);
    VSET(ep2, 0, 0, w*intercept_dist);
    VSET(ep3, eip->ehy_r1, 0, 0);
    ON_3dPoint onp1 = ON_3dPoint(ep1);
    ON_3dPoint onp2 = ON_3dPoint(ep2);
    ON_3dPoint onp3 = ON_3dPoint(ep3);

    ON_3dPointArray cpts(3);
    cpts.Append(onp1);
    cpts.Append(onp2);
    cpts.Append(onp3);
    ON_BezierCurve *bcurve = new ON_BezierCurve(cpts);
    bcurve->MakeRational();
    bcurve->SetWeight(1, w);

    ON_NurbsCurve* tnurbscurve = ON_NurbsCurve::New();
    bcurve->GetNurbForm(*tnurbscurve);
    ON_NurbsCurve* hypbnurbscurve = ON_NurbsCurve::New();
    const ON_Interval subinterval = ON_Interval(0, 0.5);
    tnurbscurve->GetNurbForm(*hypbnurbscurve, 0.0, &subinterval);

    // Next, rotate that curve around the height vector.

    point_t revpoint1, revpoint2;
    VSET(revpoint1, 0, 0, 0);
    VSET(revpoint2, 0, 0, MX);
    ON_3dPoint rpnt1 = ON_3dPoint(revpoint1);
    ON_3dPoint rpnt2 = ON_3dPoint(revpoint2);

    ON_Line revaxis = ON_Line(rpnt1, rpnt2);
    ON_RevSurface* hyp_surf = ON_RevSurface::New();
    hyp_surf->m_curve = hypbnurbscurve;
    hyp_surf->m_axis = revaxis;
    hyp_surf->m_angle = ON_Interval(0, 2*ON_PI);

    // Get the NURBS form of the surface
    ON_NurbsSurface *ehycurvedsurf = ON_NurbsSurface::New();
    hyp_surf->GetNurbForm(*ehycurvedsurf, 0.0);

    delete hyp_surf;
    delete tnurbscurve;
    delete bcurve;

    // Transformations

    for (int i = 0; i < ehycurvedsurf->CVCount(0); i++) {
	for (int j = 0; j < ehycurvedsurf->CVCount(1); j++) {
	    point_t cvpt;
	    ON_4dPoint ctrlpt;
	    ehycurvedsurf->GetCV(i, j, ctrlpt);

	    // Scale the control points of the
	    // resulting surface to map to the shorter axis.
	    VSET(cvpt, ctrlpt.x, ctrlpt.y * eip->ehy_r2/eip->ehy_r1, ctrlpt.z);

	    // Rotate according to the directions of Au and H
	    vect_t Hu;
	    mat_t R;
	    point_t new_cvpt;

	    VSCALE(Hu, eip->ehy_H, 1/MAGNITUDE(eip->ehy_H));
	    MAT_IDN(R);
	    VMOVE(&R[0], eip->ehy_Au);
	    VMOVE(&R[4], y_dir);
	    VMOVE(&R[8], Hu);
	    VEC3X3MAT(new_cvpt, cvpt, R);
	    VMOVE(cvpt, new_cvpt);

	    // Translate according to V
	    vect_t scale_v;
	    VSCALE(scale_v, eip->ehy_V, ctrlpt.w);
	    VADD2(cvpt, cvpt, scale_v);

	    ON_4dPoint newpt = ON_4dPoint(cvpt[0], cvpt[1], cvpt[2], ctrlpt.w);
	    ehycurvedsurf->SetCV(i, j, newpt);
	}
    }

    (*b)->m_S.Append(ehycurvedsurf);
    int surfindex = (*b)->m_S.Count();
    ON_BrepFace& face = (*b)->NewFace(surfindex - 1);
    (*b)->FlipFace(face);
    int faceindex = (*b)->m_F.Count();
    (*b)->NewOuterLoop(faceindex-1);
}

//.........这里部分代码省略.........
    (*b)->m_S.Append(tp);
    int tsi = (*b)->m_S.Count() - 1;
    ON_BrepFace& tface = (*b)->NewFace(tsi);
    (*b)->NewPlanarFaceLoop(tface.m_face_index, ON_BrepLoop::outer, boundary, true);
    ON_BrepLoop* tloop = (*b)->m_L.Last();
    tp->SetDomain(0, tloop->m_pbox.m_min.x, tloop->m_pbox.m_max.x);
    tp->SetDomain(1, tloop->m_pbox.m_min.y, tloop->m_pbox.m_max.y);
    tp->SetExtents(0, bp->Domain(0));
    tp->SetExtents(1, bp->Domain(1));
    (*b)->SetTrimIsoFlags(tface);
    delete tcurve;

    //  Now, the hard part.  Need an elliptical hyperbolic NURBS surface.
    //  First step is to create a nurbs curve.

    double MX = eip->hyp_b * eip->hyp_bnr;
    point_t ep1, ep2, ep3;
    VSET(ep1, -eip->hyp_b, 0, 0.5*MAGNITUDE(eip->hyp_Hi));
    VSET(ep2, -MX*eip->hyp_bnr, 0, 0);
    VSET(ep3, -eip->hyp_b, 0, -0.5*MAGNITUDE(eip->hyp_Hi));

    ON_3dPoint onp1 = ON_3dPoint(ep1);
    ON_3dPoint onp2 = ON_3dPoint(ep2);
    ON_3dPoint onp3 = ON_3dPoint(ep3);

    ON_3dPointArray cpts(3);
    cpts.Append(onp1);
    cpts.Append(onp2);
    cpts.Append(onp3);
    ON_BezierCurve *bezcurve = new ON_BezierCurve(cpts);
    bezcurve->MakeRational();
    bezcurve->SetWeight(1, bezcurve->Weight(0)/eip->hyp_bnr);

    ON_NurbsCurve* tnurbscurve = ON_NurbsCurve::New();
    bezcurve->GetNurbForm(*tnurbscurve);
    delete bezcurve;

    ON_3dPoint revpnt1 = ON_3dPoint(0, 0, -0.5*MAGNITUDE(eip->hyp_Hi));
    ON_3dPoint revpnt2 = ON_3dPoint(0, 0, 0.5*MAGNITUDE(eip->hyp_Hi));

    ON_Line revaxis = ON_Line(revpnt1, revpnt2);
    ON_RevSurface* hyp_surf = ON_RevSurface::New();
    hyp_surf->m_curve = tnurbscurve;
    hyp_surf->m_axis = revaxis;
    hyp_surf->m_angle = ON_Interval(0, 2*ON_PI);

    // Get the NURBS form of the surface
    ON_NurbsSurface *hypcurvedsurf = ON_NurbsSurface::New();
    hyp_surf->GetNurbForm(*hypcurvedsurf, 0.0);
    delete hyp_surf;

    for (int i = 0; i < hypcurvedsurf->CVCount(0); i++) {
	for (int j = 0; j < hypcurvedsurf->CVCount(1); j++) {
	    point_t cvpt;
	    ON_4dPoint ctrlpt;
	    hypcurvedsurf->GetCV(i, j, ctrlpt);

	    // Scale and shear
	    vect_t proj_ah;
	    vect_t proj_ax;
	    fastf_t factor;

	    VPROJECT(eip->hyp_A, eip->hyp_Hi, proj_ah, proj_ax);
	    VSET(cvpt, ctrlpt.x * MAGNITUDE(proj_ax)/eip->hyp_b, ctrlpt.y, ctrlpt.z);
	    factor = VDOT(eip->hyp_A, eip->hyp_Hi)>0 ? 1.0 : -1.0;
	    cvpt[2] += factor*cvpt[0]/MAGNITUDE(proj_ax)*MAGNITUDE(proj_ah) + 0.5*MAGNITUDE(eip->hyp_Hi)*ctrlpt.w;

	    // Rotate
	    vect_t Au, Bu, Hu;
	    mat_t R;
	    point_t new_cvpt;

	    VSCALE(Bu, y_dir, 1/MAGNITUDE(y_dir));
	    VSCALE(Hu, eip->hyp_Hi, 1/MAGNITUDE(eip->hyp_Hi));
	    VCROSS(Au, Bu, Hu);
	    VUNITIZE(Au);
	    MAT_IDN(R);
	    VMOVE(&R[0], Au);
	    VMOVE(&R[4], Bu);
	    VMOVE(&R[8], Hu);
	    VEC3X3MAT(new_cvpt, cvpt, R);
	    VMOVE(cvpt, new_cvpt);

	    // Translate
	    vect_t scale_v;
	    VSCALE(scale_v, eip->hyp_Vi, ctrlpt.w);
	    VADD2(cvpt, cvpt, scale_v);
	    ON_4dPoint newpt = ON_4dPoint(cvpt[0], cvpt[1], cvpt[2], ctrlpt.w);
	    hypcurvedsurf->SetCV(i, j, newpt);
	}
    }

    (*b)->m_S.Append(hypcurvedsurf);
    int surfindex = (*b)->m_S.Count();
    ON_BrepFace& face = (*b)->NewFace(surfindex - 1);
    (*b)->FlipFace(face);
    int faceindex = (*b)->m_F.Count();
    (*b)->NewOuterLoop(faceindex-1);

}

ON_NurbsSurface
FittingCylinder::initNurbsPCACylinder (int order, NurbsDataSurface *data)
{
  Eigen::Vector3d mean;
  Eigen::Matrix3d eigenvectors;
  Eigen::Vector3d eigenvalues;

  unsigned s = unsigned (data->interior.size ());

  NurbsTools::pca (data->interior, mean, eigenvectors, eigenvalues);

  data->mean = mean;
  data->eigenvectors = eigenvectors;

  eigenvalues = eigenvalues / s; // seems that the eigenvalues are dependent on the number of points (???)

  Eigen::Vector3d v_max (0.0, 0.0, 0.0);
  Eigen::Vector3d v_min (DBL_MAX, DBL_MAX, DBL_MAX);
  for (unsigned i = 0; i < s; i++)
  {
    Eigen::Vector3d p (eigenvectors.inverse () * (data->interior[i] - mean));

    if (p (0) > v_max (0))
      v_max (0) = p (0);
    if (p (1) > v_max (1))
      v_max (1) = p (1);
    if (p (2) > v_max (2))
      v_max (2) = p (2);

    if (p (0) < v_min (0))
      v_min (0) = p (0);
    if (p (1) < v_min (1))
      v_min (1) = p (1);
    if (p (2) < v_min (2))
      v_min (2) = p (2);
  }

  int ncpsU (order);
  int ncpsV (2 * order + 4);
  ON_NurbsSurface nurbs (3, false, order, order, ncpsU, ncpsV);
  nurbs.MakeClampedUniformKnotVector (0, 1.0);
  nurbs.MakePeriodicUniformKnotVector (1, 1.0 / (ncpsV - order + 1));

  double dcu = (v_max (0) - v_min (0)) / (ncpsU - 1);
  double dcv = (2.0 * M_PI) / (ncpsV - order + 1);

  double ry = std::max<double> (std::fabs (v_min (1)), std::fabs (v_max (1)));
  double rz = std::max<double> (std::fabs (v_min (2)), std::fabs (v_max (2)));

  Eigen::Vector3d cv_t, cv;
  for (int i = 0; i < ncpsU; i++)
  {
    for (int j = 0; j < ncpsV; j++)
    {
      cv (0) = v_min (0) + dcu * i;
      cv (1) = ry * sin (dcv * j);
      cv (2) = rz * cos (dcv * j);
      cv_t = eigenvectors * cv + mean;
      nurbs.SetCV (i, j, ON_3dPoint (cv_t (0), cv_t (1), cv_t (2)));
    }
  }

  return nurbs;
}