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

void
NurbsTools::computeBoundingBox (const ON_NurbsSurface &nurbs, Eigen::Vector3d &_min, Eigen::Vector3d &_max)
{
  _min = Eigen::Vector3d (DBL_MAX, DBL_MAX, DBL_MAX);
  _max = Eigen::Vector3d (-DBL_MAX, -DBL_MAX, -DBL_MAX);
  for (int i = 0; i < nurbs.CVCount (0); i++)
  {
    for (int j = 0; j < nurbs.CVCount (1); j++)
    {
      ON_3dPoint p;
      nurbs.GetCV (i, j, p);

      if (p.x < _min (0))
        _min (0) = p.x;
      if (p.y < _min (1))
        _min (1) = p.y;
      if (p.z < _min (2))
        _min (2) = p.z;

      if (p.x > _max (0))
        _max (0) = p.x;
      if (p.y > _max (1))
        _max (1) = p.y;
      if (p.z > _max (2))
        _max (2) = p.z;
    }
  }
}

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;
  }
}

void ON_GL( const ON_NurbsSurface& s,
              GLUnurbsObj* nobj, // created with gluNewNurbsRenderer )
              GLenum type,       // = 0 (and type is automatically set)
              int bPermitKnotScaling,
              double* knot_scale0,
              double* knot_scale1
             )
{
  int i, j, k;

  // The "bPermitScaling" parameters to the ON_GL() call that
  // fills in the knot vectors is set to false because any
  // rescaling that is applied to a surface domain must also
  // be applied to parameter space trimming curve geometry.

  // GL "s" knots
  GLint sknot_count = s.KnotCount(0) + 2;
  GLfloat* sknot = (GLfloat*)onmalloc( sknot_count*sizeof(*sknot) );
  ON_GL( s.Order(0), s.CVCount(0), s.Knot(0), sknot, 
           bPermitKnotScaling, knot_scale0 );

  // GL "t" knots
  GLint tknot_count = s.KnotCount(1) + 2;
  GLfloat* tknot = (GLfloat*)onmalloc( tknot_count*sizeof(*tknot) );
  ON_GL( s.Order(1), s.CVCount(1), s.Knot(1), tknot,
           bPermitKnotScaling, knot_scale1 );

  // control vertices
  const int cv_size= s.CVSize();
  const int cv_count[2] = {s.CVCount(0), s.CVCount(1)};
  GLint s_stride = cv_size*cv_count[1];
  GLint t_stride = cv_size;
  GLfloat* ctlarray = (GLfloat*)onmalloc( s_stride*cv_count[0]*sizeof(*ctlarray) );
  for ( i = 0; i < cv_count[0]; i++ ) {
    for ( j = 0; j < cv_count[1]; j++ ) {
      const double*  cv = s.CV(i,j);
      GLfloat* gl_cv = ctlarray + s_stride*i + t_stride*j;
      for ( k = 0; k < cv_size; k++ ) {
        gl_cv[k] = (GLfloat)cv[k];
      }
    }
  }
  
  GLint sorder = s.Order(0);
  GLint torder = s.Order(1);

  if ( type == 0 ) {
    // set GL surface type for 3d CVs in homogeneous/euclidean form.
    type = ( s.IsRational() ) ? GL_MAP2_VERTEX_4 : GL_MAP2_VERTEX_3;
  }

  gluNurbsSurface (
    nobj,
    sknot_count,
    sknot,
    tknot_count,
    tknot,
    s_stride, 	
    t_stride, 	
    ctlarray, 	
    sorder, 	
    torder, 	
    type	
  );	

  onfree( ctlarray );
  onfree( tknot );
  onfree( sknot );
}

void
GlobalOptimizationTDM::addCageCornerRegularisation (unsigned id, int ncps, double weight, unsigned &row)
{
  ON_NurbsSurface *nurbs = m_nurbs[id];

  { // NORTH-WEST
    int i = 0;
    int j = 0;

    //      m_solver.f(row+0, 0, 0.0);
    //      m_solver.f(row+1, 0, 0.0);
    //      m_solver.f(row+2, 0, 0.0);

    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 0, -2.0 * weight);
    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i + 1, j + 0)) + 0, 1.0 * weight);
    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 1)) + 0, 1.0 * weight);

    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 1, -2.0 * weight);
    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i + 1, j + 0)) + 1, 1.0 * weight);
    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 1)) + 1, 1.0 * weight);

    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 2, -2.0 * weight);
    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i + 1, j + 0)) + 2, 1.0 * weight);
    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 1)) + 2, 1.0 * weight);

    row += 3;
  }

  { // NORTH-EAST
    int i = nurbs->CVCount (0) - 1;
    int j = 0;

    //      m_solver.f(row+0, 0, 0.0);
    //      m_solver.f(row+1, 0, 0.0);
    //      m_solver.f(row+2, 0, 0.0);

    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 0, -2.0 * weight);
    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i - 1, j + 0)) + 0, 1.0 * weight);
    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 1)) + 0, 1.0 * weight);

    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 1, -2.0 * weight);
    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i - 1, j + 0)) + 1, 1.0 * weight);
    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 1)) + 1, 1.0 * weight);

    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 2, -2.0 * weight);
    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i - 1, j + 0)) + 2, 1.0 * weight);
    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 1)) + 2, 1.0 * weight);

    row += 3;
  }

  { // SOUTH-EAST
    int i = nurbs->CVCount (0) - 1;
    int j = nurbs->CVCount (1) - 1;

    //      m_solver.f(row+0, 0, 0.0);
    //      m_solver.f(row+1, 0, 0.0);
    //      m_solver.f(row+2, 0, 0.0);

    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 0, -2.0 * weight);
    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i - 1, j + 0)) + 0, 1.0 * weight);
    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i + 0, j - 1)) + 0, 1.0 * weight);

    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 1, -2.0 * weight);
    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i - 1, j + 0)) + 1, 1.0 * weight);
    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i + 0, j - 1)) + 1, 1.0 * weight);

    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 2, -2.0 * weight);
    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i - 1, j + 0)) + 2, 1.0 * weight);
    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i + 0, j - 1)) + 2, 1.0 * weight);

    row += 3;
  }

  { // SOUTH-WEST
    int i = 0;
    int j = nurbs->CVCount (1) - 1;

    //      m_solver.f(row+0, 0, 0.0);
    //      m_solver.f(row+1, 0, 0.0);
    //      m_solver.f(row+2, 0, 0.0);

    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 0, -2.0 * weight);
    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i + 1, j + 0)) + 0, 1.0 * weight);
    m_solver.K (row + 0, 3 * (ncps + grc2gl (*nurbs, i + 0, j - 1)) + 0, 1.0 * weight);

    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 1, -2.0 * weight);
    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i + 1, j + 0)) + 1, 1.0 * weight);
    m_solver.K (row + 1, 3 * (ncps + grc2gl (*nurbs, i + 0, j - 1)) + 1, 1.0 * weight);

    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i + 0, j + 0)) + 2, -2.0 * weight);
    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i + 1, j + 0)) + 2, 1.0 * weight);
    m_solver.K (row + 2, 3 * (ncps + grc2gl (*nurbs, i + 0, j - 1)) + 2, 1.0 * weight);

    row += 3;
  }

  //  if (!m_quiet && !(row % 100))
  //    printf("[GlobalOptimizationTDM::addCageCornerRegularisation] row: %d / %d\n", row, m_nrows);
}

//.........这里部分代码省略.........
    (*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);

}

//.........这里部分代码省略.........
    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);
}