qwt_spline_cubic.cpp

/******************************************************************************
 * Qwt Widget Library
 * Copyright (C) 1997   Josef Wilgen
 * Copyright (C) 2002   Uwe Rathmann
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the Qwt License, Version 1.0
 *****************************************************************************/
 
#include "qwt_spline_cubic.h"
#include "qwt_spline_polynomial.h"
 
#include <qpolygon.h>
#include <qpainterpath.h>
 
#define SLOPES_INCREMENTAL 0
#define KAHAN 0
 
namespace QwtSplineCubicP
{
    class KahanSum
    {
      public:
        inline KahanSum( double value = 0.0 ):
            m_sum( value ),
            m_carry( 0.0 )
        {
        }
 
        inline void reset()
        {
            m_sum = m_carry = 0.0;
        }
 
        inline double value() const
        {
            return m_sum;
        }
 
        inline void add( double value )
        {
            const double y = value - m_carry;
            const double t = m_sum + y;
 
            m_carry = ( t - m_sum ) - y;
            m_sum = t;
        }
 
        static inline double sum3( double d1, double d2, double d3 )
        {
            KahanSum sum( d1 );
            sum.add( d2 );
            sum.add( d3 );
 
            return sum.value();
        }
 
        static inline double sum4( double d1, double d2, double d3, double d4 )
        {
            KahanSum sum( d1 );
            sum.add( d2 );
            sum.add( d3 );
            sum.add( d4 );
 
            return sum.value();
        }
 
 
      private:
        double m_sum;
        double m_carry; // The carry from the previous operation
    };
 
    class CurvatureStore
    {
      public:
        inline void setup( int size )
        {
            m_curvatures.resize( size );
            m_cv = m_curvatures.data();
        }
 
        inline void storeFirst( double,
            const QPointF&, const QPointF&, double b1, double )
        {
            m_cv[0] = 2.0 * b1;
        }
 
        inline void storeNext( int index, double,
            const QPointF&, const QPointF&, double, double b2 )
        {
            m_cv[index] = 2.0 * b2;
        }
 
        inline void storeLast( double,
            const QPointF&, const QPointF&, double, double b2 )
        {
            m_cv[m_curvatures.size() - 1] = 2.0 * b2;
        }
 
        inline void storePrevious( int index, double,
            const QPointF&, const QPointF&, double b1, double )
        {
            m_cv[index] = 2.0 * b1;
        }
 
        inline void closeR()
        {
            m_cv[0] = m_cv[m_curvatures.size() - 1];
        }
 
        QVector< double > curvatures() const { return m_curvatures; }
 
      private:
        QVector< double > m_curvatures;
        double* m_cv;
    };
 
    class SlopeStore
    {
      public:
        inline void setup( int size )
        {
            m_slopes.resize( size );
            m_m = m_slopes.data();
        }
 
        inline void storeFirst( double h,
            const QPointF& p1, const QPointF& p2, double b1, double b2 )
        {
            const double s = ( p2.y() - p1.y() ) / h;
            m_m[0] = s - h * ( 2.0 * b1 + b2 ) / 3.0;
#if KAHAN
            m_sum.add( m_m[0] );
#endif
        }
 
        inline void storeNext( int index, double h,
            const QPointF& p1, const QPointF& p2, double b1, double b2 )
        {
#if SLOPES_INCREMENTAL
            Q_UNUSED( p1 )
            Q_UNUSED( p2 )
#if KAHAN
            m_sum.add( ( b1 + b2 ) * h );
            m_m[index] = m_sum.value();
#else
            m_m[index] = m_m[index - 1] + ( b1 + b2 ) * h;
#endif
#else
            const double s = ( p2.y() - p1.y() ) / h;
            m_m[index] = s + h * ( b1 + 2.0 * b2 ) / 3.0;
#endif
        }
 
        inline void storeLast( double h,
            const QPointF& p1, const QPointF& p2, double b1, double b2 )
        {
            const double s = ( p2.y() - p1.y() ) / h;
            m_m[m_slopes.size() - 1] = s + h * ( b1 + 2.0 * b2 ) / 3.0;
#if KAHAN
            m_sum.add( m_m[m_slopes.size() - 1] );
#endif
        }
 
        inline void storePrevious( int index, double h,
            const QPointF& p1, const QPointF& p2, double b1, double b2 )
        {
#if SLOPES_INCREMENTAL
            Q_UNUSED( p1 )
            Q_UNUSED( p2 )
#if KAHAN
            m_sum.add( -( b1 + b2 ) * h );
            m_m[index] = m_sum.value();
#else
            m_m[index] = m_m[index + 1] - ( b1 + b2 ) * h;
#endif
 
#else
            const double s = ( p2.y() - p1.y() ) / h;
            m_m[index] = s - h * ( 2.0 * b1 + b2 ) / 3.0;
#endif
        }
 
        inline void closeR()
        {
            m_m[0] = m_m[m_slopes.size() - 1];
        }
 
        QVector< double > slopes() const { return m_slopes; }
 
      private:
        QVector< double > m_slopes;
        double* m_m;
#if SLOPES_INCREMENTAL
        KahanSum m_sum;
#endif
    };
}
 
namespace QwtSplineCubicP
{
    class Equation2
    {
      public:
        inline Equation2()
        {
        }
 
        inline Equation2( double p0, double q0, double r0 ):
            p( p0 ),
            q( q0 ),
            r( r0 )
        {
        }
 
        inline void setup( double p0, double q0, double r0 )
        {
            p = p0;
            q = q0;
            r = r0;
        }
 
        inline Equation2 normalized() const
        {
            Equation2 c;
            c.p = 1.0;
            c.q = q / p;
            c.r = r / p;
 
            return c;
        }
 
        inline double resolved1( double x2 ) const
        {
            return ( r - q * x2 ) / p;
        }
 
        inline double resolved2( double x1 ) const
        {
            return ( r - p * x1 ) / q;
        }
 
        inline double resolved1( const Equation2& eq ) const
        {
            // find x1
            double k = q / eq.q;
            return ( r - k * eq.r ) / ( p - k * eq.p );
        }
 
        inline double resolved2( const Equation2& eq ) const
        {
            // find x2
            const double k = p / eq.p;
            return ( r - k * eq.r ) / ( q - k * eq.q );
        }
 
        // p * x1 + q * x2 = r
        double p, q, r;
    };
 
    class Equation3
    {
      public:
        inline Equation3()
        {
        }
 
        inline Equation3( const QPointF& p1, const QPointF& p2, const QPointF& p3 )
        {
            const double h1 = p2.x() - p1.x();
            const double s1 = ( p2.y() - p1.y() ) / h1;
 
            const double h2 = p3.x() - p2.x();
            const double s2 = ( p3.y() - p2.y() ) / h2;
 
            p = h1;
            q = 2 * ( h1 + h2 );
            u = h2;
            r = 3 * ( s2 - s1 );
        }
 
        inline Equation3( double cp, double cq, double du, double dr ):
            p( cp ),
            q( cq ),
            u( du ),
            r( dr )
        {
        }
 
        inline bool operator==( const Equation3& c ) const
        {
            return ( p == c.p ) && ( q == c.q ) &&
                   ( u == c.u ) && ( r == c.r );
        }
 
        inline void setup( double cp, double cq, double du, double dr )
        {
            p = cp;
            q = cq;
            u = du;
            r = dr;
        }
 
        inline Equation3 normalized() const
        {
            Equation3 c;
            c.p = 1.0;
            c.q = q / p;
            c.u = u / p;
            c.r = r / p;
 
            return c;
        }
 
        inline Equation2 substituted1( const Equation3& eq ) const
        {
            // eliminate x1
            const double k = p / eq.p;
            return Equation2( q - k * eq.q, u - k * eq.u, r - k * eq.r );
        }
 
        inline Equation2 substituted2( const Equation3& eq ) const
        {
            // eliminate x2
 
            const double k = q / eq.q;
            return Equation2( p - k * eq.p, u - k * eq.u, r - k * eq.r );
        }
 
        inline Equation2 substituted3( const Equation3& eq ) const
        {
            // eliminate x3
 
            const double k = u / eq.u;
            return Equation2( p - k * eq.p, q - k * eq.q, r - k * eq.r );
        }
 
        inline Equation2 substituted1( const Equation2& eq ) const
        {
            // eliminate x1
            const double k = p / eq.p;
            return Equation2( q - k * eq.q, u, r - k * eq.r );
        }
 
        inline Equation2 substituted3( const Equation2& eq ) const
        {
            // eliminate x3
 
            const double k = u / eq.q;
            return Equation2( p, q - k * eq.p, r - k * eq.r );
        }
 
 
        inline double resolved1( double x2, double x3 ) const
        {
            return ( r - q * x2 - u * x3 ) / p;
        }
 
        inline double resolved2( double x1, double x3 ) const
        {
            return ( r - u * x3 - p * x1 ) / q;
        }
 
        inline double resolved3( double x1, double x2 ) const
        {
            return ( r - p * x1 - q * x2 ) / u;
        }
 
        // p * x1 + q * x2 + u * x3 = r
        double p, q, u, r;
    };
}
 
#if 0
static QDebug operator<<( QDebug debug, const QwtSplineCubicP::Equation2& eq )
{
    debug.nospace() << "EQ2(" << eq.p << ", " << eq.q << ", " << eq.r << ")";
    return debug.space();
}
 
static QDebug operator<<( QDebug debug, const QwtSplineCubicP::Equation3& eq )
{
    debug.nospace() << "EQ3(" << eq.p << ", "
                    << eq.q << ", " << eq.u << ", " << eq.r << ")";
    return debug.space();
}
#endif
 
namespace QwtSplineCubicP
{
    template< class T >
    class EquationSystem
    {
      public:
        void setStartCondition( double p, double q, double u, double r )
        {
            m_conditionsEQ[0].setup( p, q, u, r );
        }
 
        void setEndCondition( double p, double q, double u, double r )
        {
            m_conditionsEQ[1].setup( p, q, u, r );
        }
 
        const T& store() const
        {
            return m_store;
        }
 
        void resolve( const QPolygonF& p )
        {
            const int n = p.size();
            if ( n < 3 )
                return;
 
            if ( m_conditionsEQ[0].p == 0.0 ||
                ( m_conditionsEQ[0].q == 0.0 && m_conditionsEQ[0].u != 0.0 ) )
            {
                return;
            }
 
            if ( m_conditionsEQ[1].u == 0.0 ||
                ( m_conditionsEQ[1].q == 0.0 && m_conditionsEQ[1].p != 0.0 ) )
            {
                return;
            }
 
            const double h0 = p[1].x() - p[0].x();
            const double h1 = p[2].x() - p[1].x();
            const double hn = p[n - 1].x() - p[n - 2].x();
 
            m_store.setup( n );
 
 
            if ( n == 3 )
            {
                // For certain conditions the first/last point does not
                // necessarily meet the spline equation and we would
                // have many solutions. In this case we resolve using
                // the spline equation - as for all other conditions.
 
                const Equation3 eqSpline0( p[0], p[1], p[2] ); // ???
                const Equation2 eq0 = m_conditionsEQ[0].substituted1( eqSpline0 );
 
                // The equation system can be solved without substitution
                // from the start/end conditions and eqSpline0 ( = eqSplineN ).
 
                double b1;
 
                if ( m_conditionsEQ[0].normalized() == m_conditionsEQ[1].normalized() )
                {
                    // When we have 3 points only and start/end conditions
                    // for 3 points mean the same condition the system
                    // is under-determined and has many solutions.
                    // We chose b1 = 0.0
 
                    b1 = 0.0;
                }
                else
                {
                    const Equation2 eq = m_conditionsEQ[1].substituted1( eqSpline0 );
                    b1 = eq0.resolved1( eq );
                }
 
                const double b2 = eq0.resolved2( b1 );
                const double b0 = eqSpline0.resolved1( b1, b2 );
 
                m_store.storeFirst( h0, p[0], p[1], b0, b1 );
                m_store.storeNext( 1, h0, p[0], p[1], b0, b1 );
                m_store.storeNext( 2, h1, p[1], p[2], b1, b2 );
 
                return;
            }
 
            const Equation3 eqSplineN( p[n - 3], p[n - 2], p[n - 1] );
            const Equation2 eqN = m_conditionsEQ[1].substituted3( eqSplineN );
 
            Equation2 eq = eqN;
            if ( n > 4 )
            {
                const Equation3 eqSplineR( p[n - 4], p[n - 3], p[n - 2] );
                eq = eqSplineR.substituted3( eq );
                eq = substituteSpline( p, eq );
            }
 
            const Equation3 eqSpline0( p[0], p[1], p[2] );
 
            double b0, b1;
            if ( m_conditionsEQ[0].u == 0.0 )
            {
                eq = eqSpline0.substituted3( eq );
 
                const Equation3& eq0 = m_conditionsEQ[0];
                b0 = Equation2( eq0.p, eq0.q, eq0.r ).resolved1( eq );
                b1 = eq.resolved2( b0 );
            }
            else
            {
                const Equation2 eqX = m_conditionsEQ[0].substituted3( eq );
                const Equation2 eqY = eqSpline0.substituted3( eq );
 
                b0 = eqY.resolved1( eqX );
                b1 = eqY.resolved2( b0 );
            }
 
            m_store.storeFirst( h0, p[0], p[1], b0, b1 );
            m_store.storeNext( 1, h0, p[0], p[1], b0, b1 );
 
            const double bn2 = resolveSpline( p, b1 );
 
            const double bn1 = eqN.resolved2( bn2 );
            const double bn0 = m_conditionsEQ[1].resolved3( bn2, bn1 );
 
            const double hx = p[n - 2].x() - p[n - 3].x();
            m_store.storeNext( n - 2, hx, p[n - 3], p[n - 2], bn2, bn1 );
            m_store.storeNext( n - 1, hn, p[n - 2], p[n - 1], bn1, bn0 );
        }
 
      private:
        Equation2 substituteSpline( const QPolygonF& points, const Equation2& eq )
        {
            const int n = points.size();
 
            m_eq.resize( n - 2 );
            m_eq[n - 3] = eq;
 
            // eq[i].resolved2( b[i-1] ) => b[i]
 
            double slope2 = ( points[n - 3].y() - points[n - 4].y() ) / eq.p;
 
            for ( int i = n - 4; i > 1; i-- )
            {
                const Equation2& eq2 = m_eq[i + 1];
                Equation2& eq1 = m_eq[i];
 
                eq1.p = points[i].x() - points[i - 1].x();
                const double slope1 = ( points[i].y() - points[i - 1].y() ) / eq1.p;
 
                const double v = eq2.p / eq2.q;
 
                eq1.q = 2.0 * ( eq1.p + eq2.p ) - v * eq2.p;
                eq1.r = 3.0 * ( slope2 - slope1 ) - v * eq2.r;
 
                slope2 = slope1;
            }
 
            return m_eq[2];
        }
 
        double resolveSpline( const QPolygonF& points, double b1 )
        {
            const int n = points.size();
            const QPointF* p = points.constData();
 
            for ( int i = 2; i < n - 2; i++ )
            {
                // eq[i].resolved2( b[i-1] ) => b[i]
                const double b2 = m_eq[i].resolved2( b1 );
                m_store.storeNext( i, m_eq[i].p, p[i - 1], p[i], b1, b2 );
 
                b1 = b2;
            }
 
            return b1;
        }
 
      private:
        Equation3 m_conditionsEQ[2];
        QVector< Equation2 > m_eq;
        T m_store;
    };
 
    template< class T >
    class EquationSystem2
    {
      public:
        const T& store() const
        {
            return m_store;
        }
 
        void resolve( const QPolygonF& p )
        {
            const int n = p.size();
 
            if ( p[n - 1].y() != p[0].y() )
            {
                // TODO ???
            }
 
            const double h0 = p[1].x() - p[0].x();
            const double s0 = ( p[1].y() - p[0].y() ) / h0;
 
            if ( n == 3 )
            {
                const double h1 = p[2].x() - p[1].x();
                const double s1 = ( p[2].y() - p[1].y() ) / h1;
 
                const double b = 3.0 * ( s0 - s1 ) / ( h0 + h1 );
 
                m_store.setup( 3 );
                m_store.storeLast( h1, p[1], p[2], -b, b );
                m_store.storePrevious( 1, h1, p[1], p[2], -b, b );
                m_store.closeR();
 
                return;
            }
 
            const double hn = p[n - 1].x() - p[n - 2].x();
 
            Equation2 eqn, eqX;
            substitute( p, eqn, eqX );
 
            const double b0 = eqn.resolved2( eqX );
            const double bn = eqn.resolved1( b0 );
 
            m_store.setup( n );
            m_store.storeLast( hn, p[n - 2], p[n - 1], bn, b0 );
            m_store.storePrevious( n - 2, hn, p[n - 2], p[n - 1], bn, b0 );
 
            resolveSpline( p, b0, bn );
 
            m_store.closeR();
        }
 
      private:
 
        void substitute( const QPolygonF& points, Equation2& eqn, Equation2& eqX )
        {
            const int n = points.size();
 
            const double hn = points[n - 1].x() - points[n - 2].x();
 
            const Equation3 eqSpline0( points[0], points[1], points[2] );
            const Equation3 eqSplineN(
                QPointF( points[0].x() - hn, points[n - 2].y() ), points[0], points[1] );
 
            m_eq.resize( n - 1 );
 
            double dq = 0;
            double dr = 0;
 
            m_eq[1] = eqSpline0;
 
            double slope1 = ( points[2].y() - points[1].y() ) / m_eq[1].u;
 
            // a) p1 * b[0] + q1 * b[1] + u1 * b[2] = r1
            // b) p2 * b[n-2] + q2 * b[0] + u2 * b[1] = r2
            // c) pi * b[i-1] + qi * b[i] + ui * b[i+1] = ri
            //
            // Using c) we can substitute b[i] ( starting from 2 ) by b[i+1]
            // until we reach n-1. As we know, that b[0] == b[n-1] we found
            // an equation where only 2 coefficients ( for b[n-2], b[0] ) are left unknown.
            // Each step we have an equation that depends on b[0], b[i] and b[i+1]
            // that can also be used to substitute b[i] in b). Ding so we end up with another
            // equation depending on b[n-2], b[0] only.
            // Finally 2 equations with 2 coefficients can be solved.
 
            for ( int i = 2; i < n - 1; i++ )
            {
                const Equation3& eq1 = m_eq[i - 1];
                Equation3& eq2 = m_eq[i];
 
                dq += eq1.p * eq1.p / eq1.q;
                dr += eq1.p * eq1.r / eq1.q;
 
                eq2.u = points[i + 1].x() - points[i].x();
                const double slope2 = ( points[i + 1].y() - points[i].y() ) / eq2.u;
 
                const double k = eq1.u / eq1.q;
 
                eq2.p = -eq1.p * k;
                eq2.q = 2.0 * ( eq1.u + eq2.u ) - eq1.u * k;
                eq2.r = 3.0 * ( slope2 - slope1 ) - eq1.r * k;
 
                slope1 = slope2;
            }
 
 
            // b[0] * m_p[n-2] + b[n-2] * m_q[n-2] + b[n-1] * pN = m_r[n-2]
            eqn.setup( m_eq[n - 2].q, m_eq[n - 2].p + eqSplineN.p, m_eq[n - 2].r );
 
            // b[n-2] * pN + b[0] * ( qN - dq ) + b[n-2] * m_p[n-2] = rN - dr
            eqX.setup( m_eq[n - 2].p + eqSplineN.p, eqSplineN.q - dq, eqSplineN.r - dr );
        }
 
        void resolveSpline( const QPolygonF& points, double b0, double bi )
        {
            const int n = points.size();
 
            for ( int i = n - 3; i >= 1; i-- )
            {
                const Equation3& eq = m_eq[i];
 
                const double b = eq.resolved2( b0, bi );
                m_store.storePrevious( i, eq.u, points[i], points[i + 1], b, bi );
 
                bi = b;
            }
        }
 
        void resolveSpline2( const QPolygonF& points,
            double b0, double bi, QVector< double >& m )
        {
            const int n = points.size();
 
            bi = m_eq[0].resolved3( b0, bi );
 
            for ( int i = 1; i < n - 2; i++ )
            {
                const Equation3& eq = m_eq[i];
 
                const double b = eq.resolved3( b0, bi );
                m[i + 1] = m[i] + ( b + bi ) * m_eq[i].u;
 
                bi = b;
            }
        }
 
        void resolveSpline3( const QPolygonF& points,
            double b0, double b1, QVector< double >& m )
        {
            const int n = points.size();
 
            double h0 = ( points[1].x() - points[0].x() );
            double s0 = ( points[1].y() - points[0].y() ) / h0;
 
            m[1] = m[0] + ( b0 + b1 ) * h0;
 
            for ( int i = 1; i < n - 1; i++ )
            {
                const double h1 = ( points[i + 1].x() - points[i].x() );
                const double s1 = ( points[i + 1].y() - points[i].y() ) / h1;
 
                const double r = 3.0 * ( s1 - s0 );
 
                const double b2 = ( r - h0 * b0 - 2.0 * ( h0 + h1 ) * b1 ) / h1;
                m[i + 1] = m[i] + ( b1 + b2 ) * h1;
 
                h0 = h1;
                s0 = s1;
                b0 = b1;
                b1 = b2;
            }
        }
 
        void resolveSpline4( const QPolygonF& points,
            double b2, double b1, QVector< double >& m )
        {
            const int n = points.size();
 
            double h2 = ( points[n - 1].x() - points[n - 2].x() );
            double s2 = ( points[n - 1].y() - points[n - 2].y() ) / h2;
 
            for ( int i = n - 2; i > 1; i-- )
            {
                const double h1 = ( points[i].x() - points[i - 1].x() );
                const double s1 = ( points[i].y() - points[i - 1].y() ) / h1;
 
                const double r = 3.0 * ( s2 - s1 );
                const double k = 2.0 * ( h1 + h2 );
 
                const double b0 = ( r - h2 * b2 - k * b1 ) / h1;
 
                m[i - 1] = m[i] - ( b0 + b1 ) * h1;
 
                h2 = h1;
                s2 = s1;
                b2 = b1;
                b1 = b0;
            }
        }
 
      public:
        QVector< Equation3 > m_eq;
        T m_store;
    };
}
 
static void qwtSetupEndEquations(
    int conditionBegin, double valueBegin, int conditionEnd, double valueEnd,
    const QPolygonF& points, QwtSplineCubicP::Equation3 eq[2] )
{
    const int n = points.size();
 
    const double h0 = points[1].x() - points[0].x();
    const double s0 = ( points[1].y() - points[0].y() ) / h0;
 
    const double hn = ( points[n - 1].x() - points[n - 2].x() );
    const double sn = ( points[n - 1].y() - points[n - 2].y() ) / hn;
 
    switch( conditionBegin )
    {
        case QwtSpline::Clamped1:
        {
            // first derivative at end points given
 
            // 3 * a1 * h + b1 = b2
            // a1 * h * h + b1 * h + c1 = s
 
            // c1 = slopeBegin
            // => b1 * ( 2 * h / 3.0 ) + b2 * ( h / 3.0 ) = s - slopeBegin
 
            // c2 = slopeEnd
            // => b1 * ( 1.0 / 3.0 ) + b2 * ( 2.0 / 3.0 ) = ( slopeEnd - s ) / h;
 
            eq[0].setup( 2 * h0 / 3.0, h0 / 3.0, 0.0, s0 - valueBegin );
            break;
        }
        case QwtSpline::Clamped2:
        {
            // second derivative at end points given
 
            // b0 = 0.5 * cvStart
            // => b0 * 1.0 + b1 * 0.0 = 0.5 * cvStart
 
            // b1 = 0.5 * cvEnd
            // => b0 * 0.0 + b1 * 1.0 = 0.5 * cvEnd
 
            eq[0].setup( 1.0, 0.0, 0.0, 0.5 * valueBegin );
            break;
        }
        case QwtSpline::Clamped3:
        {
            // third derivative at end point given
 
            // 3 * a * h0 + b[0] = b[1]
 
            // a = marg_0 / 6.0
            // => b[0] * 1.0 + b[1] * ( -1.0 ) = -0.5 * v0 * h0
 
            // a = marg_n / 6.0
            // => b[n-2] * 1.0 + b[n-1] * ( -1.0 ) = -0.5 * v1 * h5
 
            eq[0].setup( 1.0, -1.0, 0.0, -0.5 * valueBegin * h0 );
 
            break;
        }
        case QwtSpline::LinearRunout:
        {
            const double r0 = qBound( 0.0, valueBegin, 1.0 );
            if ( r0 == 0.0 )
            {
                // clamping s0
                eq[0].setup( 2 * h0 / 3.0, h0 / 3.0, 0.0, 0.0 );
            }
            else
            {
                eq[0].setup( 1.0 + 2.0 / r0, 2.0 + 1.0 / r0, 0.0, 0.0 );
            }
            break;
        }
        case QwtSplineC2::NotAKnot:
        case QwtSplineC2::CubicRunout:
        {
            // building one cubic curve from 3 points
 
            double v0;
 
            if ( conditionBegin == QwtSplineC2::CubicRunout )
            {
                // first/last point are the endpoints of the curve
 
                // b0 = 2 * b1 - b2
                // => 1.0 * b0 - 2 * b1 + 1.0 * b2 = 0.0
 
                v0 = 1.0;
            }
            else
            {
                // first/last points are on the curve,
                // the imaginary endpoints have the same distance as h0/hn
 
                v0 = h0 / ( points[2].x() - points[1].x() );
            }
 
            eq[0].setup( 1.0, -( 1.0 + v0 ), v0, 0.0 );
            break;
        }
        default:
        {
            // a natural spline, where the
            // second derivative at end points set to 0.0
            eq[0].setup( 1.0, 0.0, 0.0, 0.0 );
            break;
        }
    }
 
    switch( conditionEnd )
    {
        case QwtSpline::Clamped1:
        {
            // first derivative at end points given
            eq[1].setup( 0.0, 1.0 / 3.0 * hn, 2.0 / 3.0 * hn, valueEnd - sn );
            break;
        }
        case QwtSpline::Clamped2:
        {
            // second derivative at end points given
            eq[1].setup( 0.0, 0.0, 1.0, 0.5 * valueEnd );
            break;
        }
        case QwtSpline::Clamped3:
        {
            // third derivative at end point given
            eq[1].setup( 0.0, 1.0, -1.0, -0.5 * valueEnd * hn );
            break;
        }
        case QwtSpline::LinearRunout:
        {
            const double rn = qBound( 0.0, valueEnd, 1.0 );
            if ( rn == 0.0 )
            {
                // clamping sn
                eq[1].setup( 0.0, 1.0 / 3.0 * hn, 2.0 / 3.0 * hn, 0.0 );
            }
            else
            {
                eq[1].setup( 0.0, 2.0 + 1.0 / rn, 1.0 + 2.0 / rn, 0.0 );
            }
 
            break;
        }
        case QwtSplineC2::NotAKnot:
        case QwtSplineC2::CubicRunout:
        {
            // building one cubic curve from 3 points
 
            double vn;
 
            if ( conditionEnd == QwtSplineC2::CubicRunout )
            {
                // last point is the endpoints of the curve
                vn = 1.0;
            }
            else
            {
                // last points on the curve,
                // the imaginary endpoints have the same distance as hn
 
                vn = hn / ( points[n - 2].x() - points[n - 3].x() );
            }
 
            eq[1].setup( vn, -( 1.0 + vn ), 1.0, 0.0 );
            break;
        }
        default:
        {
            // a natural spline, where the
            // second derivative at end points set to 0.0
            eq[1].setup( 0.0, 0.0, 1.0, 0.0 );
            break;
        }
    }
}
 
class QwtSplineCubic::PrivateData
{
};
 
QwtSplineCubic::QwtSplineCubic()
    : m_data( NULL )
{
    // a natural spline
 
    setBoundaryCondition( QwtSpline::AtBeginning, QwtSpline::Clamped2 );
    setBoundaryValue( QwtSpline::AtBeginning, 0.0 );
 
    setBoundaryCondition( QwtSpline::AtEnd, QwtSpline::Clamped2 );
    setBoundaryValue( QwtSpline::AtEnd, 0.0 );
}
 
QwtSplineCubic::~QwtSplineCubic()
{
}
 
uint QwtSplineCubic::locality() const
{
    return 0;
}
 
QVector< double > QwtSplineCubic::slopes( const QPolygonF& points ) const
{
    using namespace QwtSplineCubicP;
 
    if ( points.size() <= 2 )
        return QVector< double >();
 
    if ( ( boundaryType() == QwtSpline::PeriodicPolygon )
        || ( boundaryType() == QwtSpline::ClosedPolygon ) )
    {
        EquationSystem2< SlopeStore > eqs;
        eqs.resolve( points );
 
        return eqs.store().slopes();
    }
 
    if ( points.size() == 3 )
    {
        if ( boundaryCondition( QwtSpline::AtBeginning ) == QwtSplineCubic::NotAKnot
            || boundaryCondition( QwtSpline::AtEnd ) == QwtSplineCubic::NotAKnot )
        {
#if 0
            const double h0 = points[1].x() - points[0].x();
            const double h1 = points[2].x() - points[1].x();
 
            const double s0 = ( points[1].y() - points[0].y() ) / h0;
            const double s1 = ( points[2].y() - points[1].y() ) / h1;
 
            /*
               the system is under-determined and we only
               compute a quadratic spline.
             */
 
            const double b = ( s1 - s0 ) / ( h0 + h1 );
 
            QVector< double > m( 3 );
            m[0] = s0 - h0 * b;
            m[1] = s1 - h1 * b;
            m[2] = s1 + h1 * b;
 
            return m;
#else
            return QVector< double >();
#endif
        }
    }
 
    Equation3 eq[2];
    qwtSetupEndEquations(
        boundaryCondition( QwtSpline::AtBeginning ),
        boundaryValue( QwtSpline::AtBeginning ),
        boundaryCondition( QwtSpline::AtEnd ),
        boundaryValue( QwtSpline::AtEnd ),
        points, eq );
 
    EquationSystem< SlopeStore > eqs;
    eqs.setStartCondition( eq[0].p, eq[0].q, eq[0].u, eq[0].r );
    eqs.setEndCondition( eq[1].p, eq[1].q, eq[1].u, eq[1].r );
    eqs.resolve( points );
 
    return eqs.store().slopes();
}
 
QVector< double > QwtSplineCubic::curvatures( const QPolygonF& points ) const
{
    using namespace QwtSplineCubicP;
 
    if ( points.size() <= 2 )
        return QVector< double >();
 
    if ( ( boundaryType() == QwtSpline::PeriodicPolygon )
        || ( boundaryType() == QwtSpline::ClosedPolygon ) )
    {
        EquationSystem2< CurvatureStore > eqs;
        eqs.resolve( points );
 
        return eqs.store().curvatures();
    }
 
    if ( points.size() == 3 )
    {
        if ( boundaryCondition( QwtSpline::AtBeginning ) == QwtSplineC2::NotAKnot
            || boundaryCondition( QwtSpline::AtEnd ) == QwtSplineC2::NotAKnot )
        {
            return QVector< double >();
        }
    }
 
    Equation3 eq[2];
    qwtSetupEndEquations(
        boundaryCondition( QwtSpline::AtBeginning ),
        boundaryValue( QwtSpline::AtBeginning ),
        boundaryCondition( QwtSpline::AtEnd ),
        boundaryValue( QwtSpline::AtEnd ),
        points, eq );
 
    EquationSystem< CurvatureStore > eqs;
    eqs.setStartCondition( eq[0].p, eq[0].q, eq[0].u, eq[0].r );
    eqs.setEndCondition( eq[1].p, eq[1].q, eq[1].u, eq[1].r );
    eqs.resolve( points );
 
    return eqs.store().curvatures();
}
 
QPainterPath QwtSplineCubic::painterPath( const QPolygonF& points ) const
{
    // as QwtSplineCubic can calculate slopes directly we can
    // use the implementation of QwtSplineC1 without any performance loss.
 
    return QwtSplineC1::painterPath( points );
}
 
QVector< QLineF > QwtSplineCubic::bezierControlLines( const QPolygonF& points ) const
{
    // as QwtSplineCubic can calculate slopes directly we can
    // use the implementation of QwtSplineC1 without any performance loss.
 
    return QwtSplineC1::bezierControlLines( points );
}
 
QVector< QwtSplinePolynomial > QwtSplineCubic::polynomials( const QPolygonF& points ) const
{
    return QwtSplineC2::polynomials( points );
}