C++ SplineSurface::rational方法代码示例

//===========================================================================
vector<shared_ptr<SplineSurface> >
SurfaceCreators::separateRationalParts(const SplineSurface& sf)
//===========================================================================
{
    bool rat = sf.rational();
    ASSERT(rat);

    int dim= sf.dimension();
    int rdim = dim + 1;
    vector<shared_ptr<SplineSurface> > sep_sfs;
    vector<double> coefs(sf.coefs_begin(), sf.coefs_end());
    int nmb1 = sf.numCoefs_u();
    int nmb2 = sf.numCoefs_v();
    vector<double> rcoefs;
    int num_coefs = nmb1*nmb2;
    vector<double>::const_iterator rcoef_iter = sf.rcoefs_begin();
    for (int ki = 0; ki < num_coefs; ++ki) {
	rcoefs.push_back(rcoef_iter[ki*rdim+1]);
	for (int kj = 0; kj < dim; ++kj) {
	    coefs[ki*dim+kj] /= (rcoefs.back());
	}
    }
    sep_sfs.push_back(shared_ptr<SplineSurface>
		      (new SplineSurface(nmb1, nmb2, sf.order_u(), sf.order_v(),
					 sf.basis_u().begin(), sf.basis_v().begin(),
					 coefs.begin(), dim)));
    sep_sfs.push_back(shared_ptr<SplineSurface>
		      (new SplineSurface(nmb1, nmb2, sf.order_u(), sf.order_v(),
					 sf.basis_u().begin(), sf.basis_v().begin(),
					 rcoefs.begin(), 1)));

    return sep_sfs;
}

//==========================================================================
void cart_to_bary(const SplineSurface& sf, const BaryCoordSystem3D& bc,
		  SplineSurface& sf_bc)
//==========================================================================
{
    ALWAYS_ERROR_IF(sf.dimension() != 3, "Dimension must be 3.");


    int nu = sf.numCoefs_u();
    int nv = sf.numCoefs_v();
    Vector3D cart;
    Vector4D bary;
    vector<double> new_coefs;
    if (!sf.rational()) {
	new_coefs.resize(4 * nu * nv);
	for (int iv = 0; iv < nv; ++iv) {
	    for (int iu = 0; iu < nu; ++iu) {
		int offset = nu * iv + iu;
		cart = Vector3D(sf.coefs_begin() + 3 * offset);
		bary = bc.cartToBary(cart);
		for (int j = 0; j < 4; ++j) {
		    new_coefs[4*offset + j] = bary[j];
		}
	    }
	}
    } else {
	new_coefs.resize(5 * nu * nv);
	for (int iv = 0; iv < nv; ++iv) {
	    for (int iu = 0; iu < nu; ++iu) {
		int offset = nu * iv + iu;
		cart = Vector3D(sf.coefs_begin() + 3 * offset);
		bary = bc.cartToBary(cart);
		double w = sf.rcoefs_begin()[4*offset + 3];
		for (int j = 0; j < 4; ++j) {
		    new_coefs[5*offset + j] = bary[j] * w;
		}
		new_coefs[5*offset + 4] = w;
	    }
	}
    }
    sf_bc = SplineSurface(nu, nv, sf.order_u(), sf.order_v(),
			  sf.basis_u().begin(), sf.basis_v().begin(),
			  new_coefs.begin(), 4, sf.rational());
    return;	
}

shared_ptr<SplineSurface>
GeometryTools::surfaceSum(const SplineSurface& sf1, double fac1,
                          const SplineSurface& sf2, double fac2, double num_tol)

//********************************************************************
// Addition of two signed SplineSurfaces, i.e. this function can
// also be used for subtraction. The surfaces is assumed to live on
// the same parameter domain, but may have different knot vectors.
//********************************************************************
{
    // Check input
    ALWAYS_ERROR_IF(fabs(sf1.startparam_u() - sf2.startparam_u()) > num_tol ||
                    fabs(sf1.endparam_u() - sf2.endparam_u()) > num_tol ||
                    fabs(sf1.startparam_v() - sf2.startparam_v()) > num_tol ||
                    fabs(sf1.endparam_v() - sf2.endparam_v()) > num_tol,
                    "Inconsistent parameter domain.");

    // For the time being
    if (sf1.rational() || sf2.rational()) {
        THROW("Sum of rational surfaces is not implemented");
    }

    // Make copy of surfaces
    vector<shared_ptr<SplineSurface> > surfaces;
    surfaces.reserve(2);
    shared_ptr<SplineSurface> sf;
// #ifdef _MSC_VER
//     sf = shared_ptr<SplineSurface>(dynamic_cast<SplineSurface*>(sf1.clone()));
// #else
    sf = shared_ptr<SplineSurface>(sf1.clone());
// #endif
    surfaces.push_back(sf);
// #ifdef _MSC_VER
//     sf = shared_ptr<SplineSurface>(dynamic_cast<SplineSurface*>(sf2.clone()));
// #else
    sf = shared_ptr<SplineSurface>(sf2.clone());
// #endif
    surfaces.push_back(sf);

    // Make sure that the surfaces live on the same knot vector
    GeometryTools::unifySurfaceSplineSpace(surfaces, num_tol);

    // Add signed coefficients
    vector<double> coefs;
    int nmb_coefs_u = surfaces[0]->numCoefs_u();
    int nmb_coefs_v = surfaces[0]->numCoefs_v();
    int dim = surfaces[0]->dimension();
    coefs.resize(dim*nmb_coefs_u*nmb_coefs_v);
    int ki;
    std::vector<double>::iterator s1 = surfaces[0]->coefs_begin();
    std::vector<double>::iterator s2 = surfaces[1]->coefs_begin();
    for (ki=0; ki<dim*nmb_coefs_u*nmb_coefs_v; ki++)
        coefs[ki] = fac1*s1[ki] + fac2*s2[ki];

    // Create output curve
    shared_ptr<SplineSurface>
    surfacesum(new SplineSurface(nmb_coefs_u, nmb_coefs_v,
                                 surfaces[0]->order_u(),
                                 surfaces[0]->order_v(),
                                 surfaces[0]->basis_u().begin(),
                                 surfaces[0]->basis_v().begin(),
                                 &coefs[0], dim, false));

    return surfacesum;
}

//==========================================================================
void make_matrix(const SplineSurface& surf, int deg,
		 vector<vector<double> >& mat)
//==========================================================================
{
    // Create BernsteinMulti. In the rational case the weights are
    // included in an "extra" coordinate.
    int dim = surf.dimension();
    bool rational = surf.rational();
    vector<BernsteinMulti> beta;
    spline_to_bernstein(surf, beta);

    // Make vector of basis functions (with the surface plugged in) by
    // using recursion
    int num = (deg+1) * (deg+2) * (deg+3) / 6;
    vector<BernsteinMulti> basis(num);
    vector<BernsteinMulti> tmp(num);
    basis[0] = BernsteinMulti(1.0);
    BernsteinMulti zero_multi = BernsteinMulti(0.0);
    for (int r = 1; r <= deg; ++r) {
	int m = -1;
	int tmp_num = (r + 1) * (r + 2) * (r + 3) / 6;
	fill(tmp.begin(), tmp.begin() + tmp_num, zero_multi);
	for (int i = 0; i < r; ++i) {
	    int k = (i + 1) * (i + 2) / 2;
	    for (int j = 0; j <= i; ++j) {
		for (int l = 0; l <= j; ++l) {
		    ++m;
		    tmp[m] += beta[0] * basis[m];
		    tmp[m + k] += beta[1] * basis[m];
		    tmp[m + 1 + j + k] += beta[2] * basis[m];
		    tmp[m + 2 + j + k] += beta[3] * basis[m];
		}
	    }
	}
	basis.swap(tmp);
    }

    // Fill up the matrix mat
    int deg_u = surf.order_u() - 1;
    int deg_v = surf.order_v() - 1;
    int numbas = (deg * deg_u + 1) * (deg * deg_v + 1);
    mat.resize(numbas);
    for (int row = 0; row < numbas; ++row) {
	mat[row].resize(num);
	for (int col = 0; col < num; ++col) {
	    mat[row][col] = basis[col][row];
	}
    }

    // If rational, include diagonal scaling matrix. Dividing the
    // D-matrix by the weights has the same effect as multiplying the
    // basis with the same weights. (Included for numerical reasons only -
    // it makes the basis a partition of unity.)
    if (rational) {
        BernsteinMulti weights = BernsteinMulti(1.0);
	for (int i = 1; i <= deg; ++i)
	    weights *= beta[dim];
	for (int row = 0; row < numbas; ++row) {
	    double scaling = 1.0 / weights[row];
	    for (int col = 0; col < num; ++col) {
		mat[row][col] *= scaling;
	    }
	}
    }

//     // Check Frobenius norm
//     double norm = 0.0;
//     for (int irow = 0; irow < numbas; ++irow) {
//  	for (int icol = 0; icol < num; ++icol) {
//  	    norm += mat[irow][icol] * mat[irow][icol];
//  	}
//     }
//     norm = sqrt(norm);
//     cout << "Frobenius norm = " << norm << endl;

    return;
}

int main(int argc, char** argv)
{

  if (argc != 6)
    {
      cout << "Usage: " << argv[0] << " surfaceinfile surface3doutfile points3doutfile num_u num_v" << endl;
      exit(-1);
    }

  ifstream filein(argv[1]);
  ALWAYS_ERROR_IF(filein.bad(), "Bad or no curvee input filename");
  ObjectHeader head;
  filein >> head;
  if (head.classType() != SplineSurface::classType()) {
    THROW("Not a spline surface");
  }

  SplineSurface sf;
  filein >> sf;

  ofstream fileoutsurf(argv[2]);
  ALWAYS_ERROR_IF(fileoutsurf.bad(), "Bad surface output filename");

  ofstream fileoutpts(argv[3]);
  ALWAYS_ERROR_IF(fileoutpts.bad(), "Bad points output filename");

  int num_u = atoi(argv[4]);
  int num_v = atoi(argv[5]);

  vector<double> pts, param_u, param_v;

  sf.gridEvaluator(num_u, num_v, pts, param_u, param_v);

  vector<double> coefs3d;
  vector<Point> pts3d;
  int dim = sf.dimension();
  bool rational = sf.rational();

  int ctrl_pts = sf.numCoefs_u() * sf.numCoefs_v();
  vector<double>::const_iterator it = sf.ctrl_begin();

  for (int i = 0; i < ctrl_pts; ++i)
    {
      if (dim <= 3)
	for (int j = 0; j < 3; ++j)
	  {
	    if (j>=dim)
	      coefs3d.push_back(0.0);
	    else
	      {
		coefs3d.push_back(*it);
		++it;
	      }
	  }
      else
	{
	  for (int j = 0; j < 3; ++j, ++it)
	    coefs3d.push_back(*it);
	  it += (dim-3);
	}
      if (rational)
	{
	  coefs3d.push_back(*it);
	  ++it;
	}
    }

  int pts_pos = 0;
  for (int i = 0; i < num_u*num_v; ++i)
    {
      double x, y, z;
      if (dim == 0)
	x = 0.0;
      else
	x = pts[pts_pos];
      if (dim <= 1)
	y = 0.0;
      else
	y = pts[pts_pos+1];
      if (dim <= 2)
	z = 0.0;
      else
	z = pts[pts_pos+2];
      pts_pos += dim;
      pts3d.push_back(Point(x, y, z));
    }

  SplineSurface sf3d(sf.basis_u(), sf.basis_v(), coefs3d.begin(), 3, rational);

  sf3d.writeStandardHeader(fileoutsurf);
  sf3d.write(fileoutsurf);

  fileoutpts << "400 1 0 4 255 255 0 255" << endl;
  fileoutpts << pts3d.size() << endl;
  for (int i = 0; i < (int)pts3d.size(); ++i)
    fileoutpts << pts3d[i] << endl;
}

//===========================================================================
void SplineUtils::refinedBezierCoefsCubic(SplineSurface& spline_sf,
					  int ind_u_min, int ind_v_min,
					  vector<double>& bez_coefs)
//===========================================================================
{
    assert(!spline_sf.rational());

    if (bez_coefs.size() != 48)
	bez_coefs.resize(48);
    std::fill(bez_coefs.begin(), bez_coefs.end(), 0.0);

    // Values for inpute spline surface.
    int dim = spline_sf.dimension();
    int order_u = spline_sf.order_u();
    int order_v = spline_sf.order_u();
    int num_coefs_u = spline_sf.numCoefs_u();
    int num_coefs_v = spline_sf.numCoefs_v();

    // Checking that input index is within range.
    assert(ind_u_min >= order_u - 1 && ind_u_min < num_coefs_u);
    assert(ind_v_min >= order_v - 1 && ind_v_min < num_coefs_v);

    BsplineBasis& basis_u = spline_sf.basis_u();
    BsplineBasis& basis_v = spline_sf.basis_v();
    double* knot_u = &basis_u.begin()[0];
    double* knot_v = &basis_v.begin()[0];

    // We expect the knot index to refer to the last occurence.
    assert(knot_u[ind_u_min] != knot_u[ind_u_min+1]);
    assert(knot_v[ind_v_min] != knot_v[ind_v_min+1]);

    // We expect knot mult to be 1 or 4.
    int knot_mult_umin = (knot_u[ind_u_min-1] == knot_u[ind_u_min]) ? 4 : 1;
    int knot_mult_umax = (knot_u[ind_u_min+1] == knot_u[ind_u_min+2]) ? 4 : 1;
    int knot_mult_vmin = (knot_v[ind_v_min-1] == knot_v[ind_v_min]) ? 4 : 1;
    int knot_mult_vmax = (knot_v[ind_v_min+1] == knot_v[ind_v_min+2]) ? 4 : 1;

    bool kreg_at_ustart = (knot_mult_umin == 4);
    bool kreg_at_uend = (knot_mult_umax == 4);
    vector<double> transf_mat_u(16, 0.0);
    // if (!kreg_at_ustart && !kreg_at_uend)
    splineToBezierTransfMat(knot_u + ind_u_min - 3,
			    transf_mat_u);

#ifndef NDEBUG
    std::cout << "\ntransf_mat_u=" << std::endl;
    for (size_t kj = 0; kj < 4; ++kj)
    {
	for (size_t ki = 0; ki < 4; ++ki)
	    std::cout << transf_mat_u[kj*4+ki] << " ";
	std::cout << std::endl;
    }
    std::cout << std::endl;
#endif // NDEBUG

    // else
    // 	cubicTransfMat(knot_u + ind_u_min - 3,
    // 		       kreg_at_ustart, kreg_at_uend,
    // 		       transf_mat_u);

    bool kreg_at_vstart = (knot_mult_vmin == 4);
    bool kreg_at_vend = (knot_mult_vmax == 4);
    vector<double> transf_mat_v(16, 0.0);
    // if (!kreg_at_ustart && !kreg_at_uend)
    splineToBezierTransfMat(knot_v + ind_v_min - 3,
			    transf_mat_v);

#ifndef NDEBUG
    std::cout << "\ntransf_mat_v=" << std::endl;
    for (size_t kj = 0; kj < 4; ++kj)
    {
	for (size_t ki = 0; ki < 4; ++ki)
	    std::cout << transf_mat_v[kj*4+ki] << " ";
	std::cout << std::endl;
    }
    std::cout << std::endl;
#endif // NDEBUG

    extractBezierCoefs(&spline_sf.coefs_begin()[0],
		       num_coefs_u, num_coefs_v,
		       ind_u_min, ind_v_min,
		       transf_mat_u, transf_mat_v,
		       bez_coefs);

    return;
}