I asked this on SO, but it looks like the question is headed for termination because it doesn't meet their 'rules'. I too am a steward of an open source project -- who seems to be the only one that answers questions. Thanks for all you work on Eigen.

I have a project that spends a lot time evaluating Bezier curves using this De Casteljau implementation in Eigen. As this is the innermost loop for a lot of my work, I'd like it to be as fast as possible.

The result is output in the vector p. It has columns based on the dimensionality of the problem (say x,y,z in 3D).

The Bezier curve control points are input in the matrix cp. Each control point is a row, the columns match p. A Cubic curve would have four rows. Although higher order curves may occur, cubic is typical, so cp is often 4x3.

It starts by creating Q as a copy of cp. The De Casteljau algorithm then operates in-place to calculate q.

Is there a faster way to do this in Eigen?

Code: Select all

template<typename Derived1, typename Derived2>
void de_casteljau(Eigen::MatrixBase<Derived1> &p, const Eigen::MatrixBase<Derived2> &cp, const typename Derived2::Scalar &t)
{
  // do some checks on incoming matrix dimensions
  assert(p.cols()==cp.cols());

  Eigen::Matrix<typename Derived2::Scalar, Eigen::Dynamic, Eigen::Dynamic> Q(cp);
  typename Derived2::Scalar oneminust(1-t);
  typename Derived2::Index i, k;

  for (k=1; k<Q.rows(); ++k)
  {
    for (i=0; i<Q.rows()-k; ++i)
    {
      Q.row(i)=oneminust*Q.row(i)+t*Q.row(i+1);
    }
  }

  p=Q.row(0);
}

I guess I'm hoping for a way to eliminate the inner for loop, or some other expression template magic...

Thanks for any help.

If you only care about the Q.row(0), it's just a linear combination of the rows with binomial coefficients:

Code: Select all

p.setZero();
for (i=0; i<Q.rows(); i++) {
   p += pow(oneminust, i)*pow(t, Q.rows()-1-i)*nchoosek(Q.rows()-1, i)*Q.row(i);
}


where nchoosek(n, k) = factorial(n)/factorial(k)/factorial(n-k) is a function you'll need to write or precompute and store.