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.

1 comment:

  1. hi :)

    I looked everywhere on the web to calculate the arc length of spline.

    but I found only complicated example: (
    Could you give a little piece of code to calculate it ?

    thank you in advance

    ReplyDelete