Friday, January 15, 2010

Computing Arc-Length of a Spline - Advanced

Here is a snippet of my own software used in another project that computes the arc-length of a spline using the Gauss-Legendre quadrature (numeric method of integration) implemented with recursive template meta-programming in C++:

First-off, you will need a constant 5x2 array of doubles to initialize the abscissae for, which is essentially a pre-computed table of the roots of the Legendre polynomials of order n. In this example I use order of 5, which was accurate enough for my purposes:
 const double ABSCISSAE_5[5][2] = {
-
0.90617984593866399280, 0.23692688505618908751,
-
0.53846931010568309104, 0.47862867049936646804,
0.0
, 0.56888888888888888889,
+
0.53846931010568309104, 0.47862867049936646804,
+
0.90617984593866399280, 0.23692688505618908751
};

Stub out an exception to generate runtime errors when incorrectly implementing the class:
 class CalcEx : public std::exception {};

Below the recursive class template is defined. It uses template-inheritance and since templates are interface-agnostic, you will need to define this method on your class: inline real f (const real t) const. It's essentially providing the implementation for an abstract method declared in the base class.
 template <typename T, typename Derived, int NN>
class
GaussLegendreInt
{

private
:
template
<int N> T abscissa (const unsigned int i, const unsigned int j) const
{

throw
CalcEx ();
}

template
<> T abscissa <5> (const unsigned int i, const unsigned int j) const
{

return
static_cast <T> (ABSCISSAE_5[i][j]);
}


template
<int N> T summate (const T val, const T a, const T b) const
{

return
val + summate <N - 1> (
abscissa <NN> (N - 1, 1) *
static_cast
<const Derived *> (this) -> f(
(
b + a) / 2 +
((
b - a) / 2) * abscissa <NN> (N - 1, 0)
),

a, b
);
}

template
<> T summate <0> (const T val, const T a, const T b) const { return val; }
public
:
inline
T compute (const T a, const T b) const
{

return
((b - a) / 2) * summate <NN> (0, a, b);
}
};


Now it's time to put our algorithm to use and build a class that takes a spline object (usually four points or two points and two tangents) and computes the arc-length using the Gauss-Legendre algorithm we've prepared:
template <typename Spline, typename real>
class
ArcLength : private GaussLegendreInt <real, ArcLength <Spline, real>, 5>
{

private
:
const
Spline & _spline;

public
:
inline
ArcLength(const Spline & spline)
:
_spline(spline) {}

inline
real f (const real t) const
{

return
MAGNITUDE(_spline.computeAt (t));
}


static inline
real calculate (const Spline & spline)
{

return
ArcLength(spline).compute(0, 1);
}
};

The spline object assumes a function called computeAt that returns a vector (2D, 3D or whatever) of which the function MAGNITUDE computes the vector magnitude of. You will define these yourself including your implementation of splines. There are plenty of examples on the Internet on how to implement splines. Some examples include Cubic B-Splines, Catmull-Rom Splines, and NURBS.