path: root/chart2/source/view/charttypes/Splines.cxx

/*************************************************************************
 *
 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
 *
 * Copyright 2008 by Sun Microsystems, Inc.
 *
 * OpenOffice.org - a multi-platform office productivity suite
 *
 * $RCSfile: Splines.cxx,v $
 * $Revision: 1.11.44.1 $
 *
 * This file is part of OpenOffice.org.
 *
 * OpenOffice.org is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License version 3
 * only, as published by the Free Software Foundation.
 *
 * OpenOffice.org is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License version 3 for more details
 * (a copy is included in the LICENSE file that accompanied this code).
 *
 * You should have received a copy of the GNU Lesser General Public License
 * version 3 along with OpenOffice.org.  If not, see
 * <http://www.openoffice.org/license.html>
 * for a copy of the LGPLv3 License.
 *
 ************************************************************************/


// MARKER(update_precomp.py): autogen include statement, do not remove
#include "precompiled_chart2.hxx"

#include "Splines.hxx"
#include <rtl/math.hxx>

#include <vector>
#include <algorithm>
#include <functional>

// header for DBG_ASSERT
#include <tools/debug.hxx>

//.............................................................................
namespace chart
{
//.............................................................................
using namespace ::com::sun::star;

namespace
{

template< typename T >
struct lcl_EqualsFirstDoubleOfPair : ::std::binary_function< ::std::pair< double, T >, ::std::pair< double, T >, bool >
{
    inline bool operator() ( const ::std::pair< double, T > & rOne, const ::std::pair< double, T > & rOther )
    {
        return ( ::rtl::math::approxEqual( rOne.first, rOther.first ) );
    }
};

//-----------------------------------------------------------------------------

typedef ::std::pair< double, double >   tPointType;
typedef ::std::vector< tPointType >     tPointVecType;
typedef tPointVecType::size_type        lcl_tSizeType;

class lcl_SplineCalculation
{
public:
    /** @descr creates an object that calculates cublic splines on construction

        @param rSortedPoints  the points for which splines shall be calculated, they need to be sorted in x values
        @param fY1FirstDerivation the resulting spline should have the first
               derivation equal to this value at the x-value of the first point
               of rSortedPoints.  If fY1FirstDerivation is set to infinity, a natural
               spline is calculated.
        @param fYnFirstDerivation the resulting spline should have the first
               derivation equal to this value at the x-value of the last point
               of rSortedPoints
     */
    lcl_SplineCalculation( const tPointVecType & rSortedPoints,
                           double fY1FirstDerivation,
                           double fYnFirstDerivation );

    /** @descr this function corresponds to the function splint in [1].

        [1] Numerical Recipies in C, 2nd edition
            William H. Press, et al.,
            Section 3.3, page 116
    */
    double GetInterpolatedValue( double x );

private:
    /// a copy of the points given in the CTOR
    tPointVecType            m_aPoints;

    /// the result of the Calculate() method
    ::std::vector< double >         m_aSecDerivY;

    double m_fYp1;
    double m_fYpN;

    // these values are cached for performance reasons
    tPointVecType::size_type m_nKLow;
    tPointVecType::size_type m_nKHigh;
    double m_fLastInterpolatedValue;

    /** @descr this function corresponds to the function spline in [1].

        [1] Numerical Recipies in C, 2nd edition
            William H. Press, et al.,
            Section 3.3, page 115
    */
    void Calculate();
};

//-----------------------------------------------------------------------------

lcl_SplineCalculation::lcl_SplineCalculation(
    const tPointVecType & rSortedPoints,
    double fY1FirstDerivation,
    double fYnFirstDerivation )
        : m_aPoints( rSortedPoints ),
          m_fYp1( fY1FirstDerivation ),
          m_fYpN( fYnFirstDerivation ),
          m_nKLow( 0 ),
          m_nKHigh( rSortedPoints.size() - 1 )
{
    ::rtl::math::setInf( &m_fLastInterpolatedValue, sal_False );

    // #108301# remove points that have equal x-values
    m_aPoints.erase( ::std::unique( m_aPoints.begin(), m_aPoints.end(),
                             lcl_EqualsFirstDoubleOfPair< double >() ),
                     m_aPoints.end() );
    Calculate();
}

void lcl_SplineCalculation::Calculate()
{
    if( m_aPoints.size() <= 1 )
        return;

    // n is the last valid index to m_aPoints
    const tPointVecType::size_type n = m_aPoints.size() - 1;
    ::std::vector< double > u( n );
    m_aSecDerivY.resize( n + 1, 0.0 );

    if( ::rtl::math::isInf( m_fYp1 ) )
    {
        // natural spline
        m_aSecDerivY[ 0 ] = 0.0;
        u[ 0 ] = 0.0;
    }
    else
    {
        m_aSecDerivY[ 0 ] = -0.5;
        double xDiff = ( m_aPoints[ 1 ].first - m_aPoints[ 0 ].first );
        u[ 0 ] = ( 3.0 / xDiff ) *
            ((( m_aPoints[ 1 ].second - m_aPoints[ 0 ].second ) / xDiff ) - m_fYp1 );
    }

    for( tPointVecType::size_type i = 1; i < n; ++i )
    {
        ::std::pair< double, double >
            p_i = m_aPoints[ i ],
            p_im1 = m_aPoints[ i - 1 ],
            p_ip1 = m_aPoints[ i + 1 ];

        double sig = ( p_i.first - p_im1.first ) /
            ( p_ip1.first - p_im1.first );
        double p = sig * m_aSecDerivY[ i - 1 ] + 2.0;

        m_aSecDerivY[ i ] = ( sig - 1.0 ) / p;
        u[ i ] =
            ( ( p_ip1.second - p_i.second ) /
              ( p_ip1.first - p_i.first ) ) -
            ( ( p_i.second - p_im1.second ) /
              ( p_i.first - p_im1.first ) );
        u[ i ] =
            ( 6.0 * u[ i ] / ( p_ip1.first - p_im1.first )
              - sig * u[ i - 1 ] ) / p;
    }

    // initialize to values for natural splines (used for m_fYpN equal to
    // infinity)
    double qn = 0.0;
    double un = 0.0;

    if( ! ::rtl::math::isInf( m_fYpN ) )
    {
        qn = 0.5;
        double xDiff = ( m_aPoints[ n ].first - m_aPoints[ n - 1 ].first );
        un = ( 3.0 / xDiff ) *
            ( m_fYpN - ( m_aPoints[ n ].second - m_aPoints[ n - 1 ].second ) / xDiff );
    }

    m_aSecDerivY[ n ] = ( un - qn * u[ n - 1 ] ) * ( qn * m_aSecDerivY[ n - 1 ] + 1.0 );

    // note: the algorithm in [1] iterates from n-1 to 0, but as size_type
    // may be (usuall is) an unsigned type, we can not write k >= 0, as this
    // is always true.
    for( tPointVecType::size_type k = n; k > 0; --k )
    {
        ( m_aSecDerivY[ k - 1 ] *= m_aSecDerivY[ k ] ) += u[ k - 1 ];
    }
}

double lcl_SplineCalculation::GetInterpolatedValue( double x )
{
    DBG_ASSERT( ( m_aPoints[ 0 ].first <= x ) &&
                ( x <= m_aPoints[ m_aPoints.size() - 1 ].first ),
                "Trying to extrapolate" );

    const tPointVecType::size_type n = m_aPoints.size() - 1;
    if( x < m_fLastInterpolatedValue )
    {
        m_nKLow = 0;
        m_nKHigh = n;

        // calculate m_nKLow and m_nKHigh
        // first initialization is done in CTOR
        while( m_nKHigh - m_nKLow > 1 )
        {
            tPointVecType::size_type k = ( m_nKHigh + m_nKLow ) / 2;
            if( m_aPoints[ k ].first > x )
                m_nKHigh = k;
            else
                m_nKLow = k;
        }
    }
    else
    {
        while( ( m_aPoints[ m_nKHigh ].first < x ) &&
               ( m_nKHigh <= n ) )
        {
            ++m_nKHigh;
            ++m_nKLow;
        }
        DBG_ASSERT( m_nKHigh <= n, "Out of Bounds" );
    }
    m_fLastInterpolatedValue = x;

    double h = m_aPoints[ m_nKHigh ].first - m_aPoints[ m_nKLow ].first;
    DBG_ASSERT( h != 0, "Bad input to GetInterpolatedValue()" );

    double a = ( m_aPoints[ m_nKHigh ].first - x ) / h;
    double b = ( x - m_aPoints[ m_nKLow ].first  ) / h;

    return ( a * m_aPoints[ m_nKLow ].second +
             b * m_aPoints[ m_nKHigh ].second +
             (( a*a*a - a ) * m_aSecDerivY[ m_nKLow ] +
              ( b*b*b - b ) * m_aSecDerivY[ m_nKHigh ] ) *
             ( h*h ) / 6.0 );
}

//-----------------------------------------------------------------------------

//create knot vector for B-spline
double* createTVector( sal_Int32 n, sal_Int32 k )
{
    double* t = new double [n + k + 1];
    for (sal_Int32 i=0; i<=n+k; i++ )
    {
        if(i < k)
            t[i] = 0;
        else if(i <= n)
            t[i] = i-k+1;
        else
            t[i] = n-k+2;
    }
    return t;
}

//calculate left knot vector
double TLeft (double x, sal_Int32 i, sal_Int32 k, const double *t )
{
    double deltaT = t[i + k - 1] - t[i];
    return (deltaT == 0.0)
               ? 0.0
               : (x - t[i]) / deltaT;
}

//calculate right knot vector
double TRight(double x, sal_Int32 i, sal_Int32 k, const double *t )
{
    double deltaT = t[i + k] - t[i + 1];
    return (deltaT == 0.0)
               ? 0.0
               : (t[i + k] - x) / deltaT;
}

//calculate weight vector
void BVector(double x, sal_Int32 n, sal_Int32 k, double *b, const double *t)
{
    sal_Int32 i = 0;
    for( i=0; i<=n+k; i++ )
        b[i]=0;

    sal_Int32 i0 = (sal_Int32)floor(x) + k - 1;
    b [i0] = 1;

    for( sal_Int32 j=2; j<=k; j++ )
        for( i=0; i<=i0; i++ )
            b[i] = TLeft(x, i, j, t) * b[i] + TRight(x, i, j, t) * b [i + 1];
}

} //  anonymous namespace

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------

void SplineCalculater::CalculateCubicSplines(
    const drawing::PolyPolygonShape3D& rInput
    , drawing::PolyPolygonShape3D& rResult
    , sal_Int32 nGranularity )
{
    DBG_ASSERT( nGranularity > 0, "Granularity is invalid" );
    rResult.SequenceX.realloc(0);
    rResult.SequenceY.realloc(0);
    rResult.SequenceZ.realloc(0);

    if( !rInput.SequenceX.getLength() )
        return;
    if( rInput.SequenceX[0].getLength() <= 1 )
        return; //we need at least two points

    sal_Int32 nMaxIndexPoints = rInput.SequenceX[0].getLength()-1; // is >=1
    const double* pOldX = rInput.SequenceX[0].getConstArray();
    const double* pOldY = rInput.SequenceY[0].getConstArray();
    const double* pOldZ = rInput.SequenceZ[0].getConstArray();

    // #i13699# The curve gets a parameter and then for each coordinate a
    // separate spline will be calculated using the parameter as first argument
    // and the point coordinate as second argument. Therefore the points need
    // not to be sorted in its x-coordinates. The parameter is sorted by
    // construction.

    ::std::vector < double > aParameter(nMaxIndexPoints+1);
    aParameter[0]=0.0;
    for( sal_Int32 nIndex=1; nIndex<=nMaxIndexPoints; nIndex++ )
    {
        // The euclidian distance leads to curve loops for functions having single extreme points
//         aParameter[nIndex]=aParameter[nIndex-1]+
//             sqrt( (pOldX[nIndex]-pOldX[nIndex-1])*(pOldX[nIndex]-pOldX[nIndex-1])+
//                   (pOldY[nIndex]-pOldY[nIndex-1])*(pOldY[nIndex]-pOldY[nIndex-1])+
//                   (pOldZ[nIndex]-pOldZ[nIndex-1])*(pOldZ[nIndex]-pOldZ[nIndex-1]));

        // use increment of 1 instead
        aParameter[nIndex]=aParameter[nIndex-1]+1;
    }
    // Split the calculation to X, Y and Z coordinate
    tPointVecType aInputX;
    aInputX.resize(nMaxIndexPoints+1);
    tPointVecType aInputY;
    aInputY.resize(nMaxIndexPoints+1);
    tPointVecType aInputZ;
    aInputZ.resize(nMaxIndexPoints+1);
    for (sal_Int32 nN=0;nN<=nMaxIndexPoints; nN++ )
    {
      aInputX[ nN ].first=aParameter[nN];
      aInputX[ nN ].second=pOldX[ nN ];
      aInputY[ nN ].first=aParameter[nN];
      aInputY[ nN ].second=pOldY[ nN ];
      aInputZ[ nN ].first=aParameter[nN];
      aInputZ[ nN ].second=pOldZ[ nN ];
    }

    // generate a spline for each coordinate. It holds the complete
    // information to calculate each point of the curve

    // generate the kind "natural spline"
    double fInfty;
    ::rtl::math::setInf( &fInfty, sal_False );
    lcl_SplineCalculation aSplineX( aInputX, fInfty, fInfty );
    lcl_SplineCalculation aSplineY( aInputY, fInfty, fInfty );
    lcl_SplineCalculation aSplineZ( aInputZ, fInfty, fInfty );

    // fill result polygon with calculated values
    rResult.SequenceX.realloc(1);
    rResult.SequenceY.realloc(1);
    rResult.SequenceZ.realloc(1);
    rResult.SequenceX[0].realloc( nMaxIndexPoints*nGranularity + 1);
    rResult.SequenceY[0].realloc( nMaxIndexPoints*nGranularity + 1);
    rResult.SequenceZ[0].realloc( nMaxIndexPoints*nGranularity + 1);

    double* pNewX = rResult.SequenceX[0].getArray();
    double* pNewY = rResult.SequenceY[0].getArray();
    double* pNewZ = rResult.SequenceZ[0].getArray();

    sal_Int32 nNewPointIndex = 0; // Index in result points
    // needed for inner loop
    double    fInc;   // step for intermediate points
    sal_Int32 nj;     // for loop
    double    fParam; // a intermediate parameter value

    for( sal_Int32 ni = 0; ni < nMaxIndexPoints; ni++ )
    {
        // given point is surely a curve point
        pNewX[nNewPointIndex] = pOldX[ni];
        pNewY[nNewPointIndex] = pOldY[ni];
        pNewZ[nNewPointIndex] = pOldZ[ni];
        nNewPointIndex++;

        // calculate intermediate points
        fInc = ( aParameter[ ni+1 ] - aParameter[ni] ) / static_cast< double >( nGranularity );
        for(nj = 1; nj < nGranularity; nj++)
        {
            fParam = aParameter[ni] + ( fInc * static_cast< double >( nj ) );

            pNewX[nNewPointIndex]=aSplineX.GetInterpolatedValue( fParam );
            pNewY[nNewPointIndex]=aSplineY.GetInterpolatedValue( fParam );
            pNewZ[nNewPointIndex]=aSplineZ.GetInterpolatedValue( fParam );
            nNewPointIndex++;
        }
    }
    // add last point
    pNewX[nNewPointIndex] = pOldX[nMaxIndexPoints];
    pNewY[nNewPointIndex] = pOldY[nMaxIndexPoints];
    pNewZ[nNewPointIndex] = pOldZ[nMaxIndexPoints];
}

void SplineCalculater::CalculateBSplines(
            const ::com::sun::star::drawing::PolyPolygonShape3D& rInput
            , ::com::sun::star::drawing::PolyPolygonShape3D& rResult
            , sal_Int32 nGranularity
            , sal_Int32 nDegree )
{
    // #issue 72216#
    // k is the order of the BSpline, nDegree is the degree of its polynoms
    sal_Int32 k = nDegree + 1;

    rResult.SequenceX.realloc(0);
    rResult.SequenceY.realloc(0);
    rResult.SequenceZ.realloc(0);

    if( !rInput.SequenceX.getLength() )
        return; // no input

    if( rInput.SequenceX[0].getLength() <= 1 )
        return; // need at least 2 control points

    sal_Int32 n = rInput.SequenceX[0].getLength()-1; // maximum index of control points

    double fCurveparam =0.0; // parameter for the curve
    // 0<= fCurveparam < fMaxCurveparam
    double fMaxCurveparam = 2.0+ n - k;
    if (fMaxCurveparam <= 0.0)
        return; // not enough control points for desired spline order

    if (nGranularity < 1)
        return; //need at least 1 line for each part beween the control points

    const double* pOldX = rInput.SequenceX[0].getConstArray();
    const double* pOldY = rInput.SequenceY[0].getConstArray();
    const double* pOldZ = rInput.SequenceZ[0].getConstArray();

    // keep this amount of steps to go well with old version
    sal_Int32 nNewSectorCount = nGranularity * n;
    double fCurveStep = fMaxCurveparam/static_cast< double >(nNewSectorCount);

    double *b       = new double [n + k + 1]; // values of blending functions

    const double* t = createTVector(n, k); // knot vector

    rResult.SequenceX.realloc(1);
    rResult.SequenceY.realloc(1);
    rResult.SequenceZ.realloc(1);
    rResult.SequenceX[0].realloc(nNewSectorCount+1);
    rResult.SequenceY[0].realloc(nNewSectorCount+1);
    rResult.SequenceZ[0].realloc(nNewSectorCount+1);
    double* pNewX = rResult.SequenceX[0].getArray();
    double* pNewY = rResult.SequenceY[0].getArray();
    double* pNewZ = rResult.SequenceZ[0].getArray();

    // variables needed inside loop, when calculating one point of output
    sal_Int32 nPointIndex =0; //index of given contol points
    double fX=0.0;
    double fY=0.0;
    double fZ=0.0; //coordinates of a new BSpline point

    for(sal_Int32 nNewSector=0; nNewSector<nNewSectorCount; nNewSector++)
    { // in first looping fCurveparam has value 0.0

        // Calculate the values of the blending functions for actual curve parameter
        BVector(fCurveparam, n, k, b, t);

        // output point(fCurveparam) = sum over {input point * value of blending function}
        fX = 0.0;
        fY = 0.0;
        fZ = 0.0;
        for (nPointIndex=0;nPointIndex<=n;nPointIndex++)
        {
            fX +=pOldX[nPointIndex]*b[nPointIndex];
            fY +=pOldY[nPointIndex]*b[nPointIndex];
            fZ +=pOldZ[nPointIndex]*b[nPointIndex];
        }
        pNewX[nNewSector] = fX;
        pNewY[nNewSector] = fY;
        pNewZ[nNewSector] = fZ;

        fCurveparam += fCurveStep; //for next looping
    }
    // add last control point to BSpline curve
    pNewX[nNewSectorCount] = pOldX[n];
    pNewY[nNewSectorCount] = pOldY[n];
    pNewZ[nNewSectorCount] = pOldZ[n];

    delete[] t;
    delete[] b;
}

//.............................................................................
} //namespace chart
//.............................................................................