Attachment 47488 Details for Bug 37889 – patch for Splines.cxx against LO 3.4.0.2

[patch] patch for Splines.cxx against LO 3.4.0.2

Splines_cxx.patch (text/plain), 38.47 KB, created by Regina Henschel on 2011-06-03 12:12:41 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Regina Henschel

Created: 2011-06-03 12:12:41 UTC

Size: 38.47 KB

patch

obsolete

>--- Splines.cxx.original	2011-05-19 11:09:16.000000000 +0200
>+++ Splines.cxx	2011-06-03 20:49:57.187500000 +0200
>@@ -49,17 +49,6 @@
> 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;
>@@ -82,6 +71,16 @@
>                            double fY1FirstDerivation,
>                            double fYnFirstDerivation );
> 
>+
>+    /** @descr creates an object that calculates cublic splines on construction
>+               for the special case of periodic cubic spline
>+
>+        @param rSortedPoints  the points for which splines shall be calculated,
>+               they need to be sorted in x values. First and last y value must be equal
>+     */
>+    lcl_SplineCalculation( const tPointVecType & rSortedPoints);
>+
>+
>     /** @descr this function corresponds to the function splint in [1].
> 
>         [1] Numerical Recipies in C, 2nd edition
>@@ -101,8 +100,8 @@
>     double m_fYpN;
> 
>     // these values are cached for performance reasons
>-    tPointVecType::size_type m_nKLow;
>-    tPointVecType::size_type m_nKHigh;
>+    lcl_tSizeType m_nKLow;
>+    lcl_tSizeType m_nKHigh;
>     double m_fLastInterpolatedValue;
> 
>     /** @descr this function corresponds to the function spline in [1].
>@@ -112,6 +111,22 @@
>             Section 3.3, page 115
>     */
>     void Calculate();
>+    
>+    /** @descr this function corresponds to the algoritm 4.76 in [2] and
>+        theorem 5.3.7 in [3]
>+        
>+        [2] Engeln-MÃ¼llges, Gisela: Numerik-Algorithmen: Verfahren, Beispiele, Anwendungen
>+            Springer, Berlin; Auflage: 9., Ã¼berarb. und erw. A. (8. Dezember 2004)
>+            Section 4.10.2, page 175
>+
>+        [3] Hanrath, Wilhelm: Mathematik III / Numerik, Vorlesungsskript zur
>+            Veranstaltung im WS 2007/2008
>+            Fachhochschule Aachen, 2009-09-19
>+            Numerik_01.pdf, downloaded 2011-04-19 via
>+            http://www.fh-aachen.de/index.php?id=11424&no_cache=1&file=5016&uid=44191
>+            Section 5.3, page 129
>+    */
>+    void CalculatePeriodic();
> };
> 
> //-----------------------------------------------------------------------------
>@@ -127,21 +142,32 @@
>           m_nKHigh( rSortedPoints.size() - 1 )
> {
>     ::rtl::math::setInf( &m_fLastInterpolatedValue, sal_False );
>-
>-    // 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();
> }
> 
>+//-----------------------------------------------------------------------------
>+
>+lcl_SplineCalculation::lcl_SplineCalculation(
>+    const tPointVecType & rSortedPoints)
>+        : m_aPoints( rSortedPoints ),
>+          m_fYp1( 0.0 ),  /*dummy*/
>+          m_fYpN( 0.0 ),  /*dummy*/
>+          m_nKLow( 0 ),
>+          m_nKHigh( rSortedPoints.size() - 1 )
>+{
>+    ::rtl::math::setInf( &m_fLastInterpolatedValue, sal_False );
>+    CalculatePeriodic();
>+}
>+//-----------------------------------------------------------------------------
>+
>+
> 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;
>+    const lcl_tSizeType n = m_aPoints.size() - 1;
>     ::std::vector< double > u( n );
>     m_aSecDerivY.resize( n + 1, 0.0 );
> 
>@@ -159,9 +185,9 @@
>             ((( m_aPoints[ 1 ].second - m_aPoints[ 0 ].second ) / xDiff ) - m_fYp1 );
>     }
> 
>-    for( tPointVecType::size_type i = 1; i < n; ++i )
>+    for( lcl_tSizeType i = 1; i < n; ++i )
>     {
>-        ::std::pair< double, double >
>+        tPointType
>             p_i = m_aPoints[ i ],
>             p_im1 = m_aPoints[ i - 1 ],
>             p_ip1 = m_aPoints[ i + 1 ];
>@@ -199,19 +225,164 @@
>     // 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 )
>+    for( lcl_tSizeType k = n; k > 0; --k )
>     {
>         ( m_aSecDerivY[ k - 1 ] *= m_aSecDerivY[ k ] ) += u[ k - 1 ];
>     }
> }
> 
>+void lcl_SplineCalculation::CalculatePeriodic()
>+{
>+    if( m_aPoints.size() <= 1 )
>+        return;
>+
>+    // n is the last valid index to m_aPoints
>+    const lcl_tSizeType n = m_aPoints.size() - 1;
>+    
>+    // u is used for vector f in A*c=f in [3], vector a in  Ax=a in [2],
>+    // vector z in Rtranspose z = a and Dr=z in [2]
>+    ::std::vector< double > u( n + 1, 0.0 );
>+    
>+    // used for vector c in A*c=f and vector x in Ax=a in [2]
>+    m_aSecDerivY.resize( n + 1, 0.0 );
>+    
>+    // diagonal of matrix A, used index 1 to n
>+    ::std::vector< double > Adiag( n + 1, 0.0 );
>+    
>+    // secondary diagonal of matrix A with index 1 to n-1 and upper right element in A[n]
>+    ::std::vector< double > Aupper( n + 1, 0.0 );
>+    
>+    // diagonal of matrix D in A=(R transpose)*D*R in [2], used index 1 to n
>+    ::std::vector< double > Ddiag( n+1, 0.0 );
>+    
>+    // right column of matrix R, used index 1 to n-2
>+    ::std::vector< double > Rright( n-1, 0.0 );
>+    
>+    // secondary diagonal of matrix R, used index 1 to n-1
>+    ::std::vector< double > Rupper( n, 0.0 );
>+
>+    if (n<4)
>+    {
>+        if (n==3)
>+        {   // special handling of three polynomials, that are four points
>+            double xDiff0 = m_aPoints[ 1 ].first - m_aPoints[ 0 ].first ;
>+            double xDiff1 = m_aPoints[ 2 ].first - m_aPoints[ 1 ].first ;
>+            double xDiff2 = m_aPoints[ 3 ].first - m_aPoints[ 2 ].first ;
>+            double xDiff2p1 = xDiff2 + xDiff1;
>+            double xDiff0p2 = xDiff0 + xDiff2;
>+            double xDiff1p0 = xDiff1 + xDiff0;
>+            double fFaktor = 1.5 / (xDiff0*xDiff1 + xDiff1*xDiff2 + xDiff2*xDiff0);
>+            double yDiff0 = (m_aPoints[ 1 ].second - m_aPoints[ 0 ].second) / xDiff0;
>+            double yDiff1 = (m_aPoints[ 2 ].second - m_aPoints[ 1 ].second) / xDiff1;
>+            double yDiff2 = (m_aPoints[ 0 ].second - m_aPoints[ 2 ].second) / xDiff2;
>+            m_aSecDerivY[ 1 ] = fFaktor * (yDiff1*xDiff2p1 - yDiff0*xDiff0p2);
>+            m_aSecDerivY[ 2 ] = fFaktor * (yDiff2*xDiff0p2 - yDiff1*xDiff1p0);
>+            m_aSecDerivY[ 3 ] = fFaktor * (yDiff0*xDiff1p0 - yDiff2*xDiff2p1);
>+            m_aSecDerivY[ 0 ] = m_aSecDerivY[ 3 ];
>+        }
>+        else if (n==2)
>+        {
>+        // special handling of two polynomials, that are three points
>+            double xDiff0 = m_aPoints[ 1 ].first - m_aPoints[ 0 ].first;
>+            double xDiff1 = m_aPoints[ 2 ].first - m_aPoints[ 1 ].first;
>+            double fHelp = 3.0 * (m_aPoints[ 0 ].second - m_aPoints[ 1 ].second) / (xDiff0*xDiff1);
>+            m_aSecDerivY[ 1 ] = fHelp ;
>+            m_aSecDerivY[ 2 ] = -fHelp ;
>+            m_aSecDerivY[ 0 ] = m_aSecDerivY[ 2 ] ;
>+        }
>+        else
>+        {
>+            // should be handled with natural spline, periodic not possible.
>+        }
>+    }
>+    else
>+    {
>+        double xDiff_i =1.0; // values are dummy;
>+        double xDiff_im1 =1.0;
>+        double yDiff_i = 1.0;
>+        double yDiff_im1 = 1.0;
>+        // fill matrix A and fill right side vector u
>+        for( lcl_tSizeType i=1; i<n; ++i )
>+        {
>+            xDiff_im1 = m_aPoints[ i ].first - m_aPoints[ i-1 ].first;
>+            xDiff_i = m_aPoints[ i+1 ].first - m_aPoints[ i ].first;
>+            yDiff_im1 = (m_aPoints[ i ].second - m_aPoints[ i-1 ].second) / xDiff_im1;
>+            yDiff_i = (m_aPoints[ i+1 ].second - m_aPoints[ i ].second) / xDiff_i;
>+            Adiag[ i ] = 2 * (xDiff_im1 + xDiff_i);
>+            Aupper[ i ] = xDiff_i;
>+            u [ i ] = 3 * (yDiff_i - yDiff_im1);
>+        }
>+        xDiff_im1 = m_aPoints[ n ].first - m_aPoints[ n-1 ].first;
>+        xDiff_i = m_aPoints[ 1 ].first - m_aPoints[ 0 ].first;
>+        yDiff_im1 = (m_aPoints[ n ].second - m_aPoints[ n-1 ].second) / xDiff_im1;
>+        yDiff_i = (m_aPoints[ 1 ].second - m_aPoints[ 0 ].second) / xDiff_i;
>+        Adiag[ n ] = 2 * (xDiff_im1 + xDiff_i);
>+        Aupper[ n ] = xDiff_i;
>+        u [ n ] = 3 * (yDiff_i - yDiff_im1);
>+        
>+        // decomposite A=(R transpose)*D*R
>+        Ddiag[1] = Adiag[1];
>+        Rupper[1] = Aupper[1] / Ddiag[1];
>+        Rright[1] = Aupper[n] / Ddiag[1];
>+        for( lcl_tSizeType i=2; i<=n-2; ++i )
>+        {
>+            Ddiag[i] = Adiag[i] - Aupper[ i-1 ] * Rupper[ i-1 ];
>+            Rupper[ i ] = Aupper[ i ] / Ddiag[ i ];
>+            Rright[ i ] = - Rright[ i-1 ] * Aupper[ i-1 ] / Ddiag[ i ];
>+        }
>+        Ddiag[ n-1 ] = Adiag[ n-1 ] - Aupper[ n-2 ] * Rupper[ n-2 ];
>+        Rupper[ n-1 ] = ( Aupper[ n-1 ] - Aupper[ n-2 ] * Rright[ n-2] ) / Ddiag[ n-1 ];
>+        double fSum = 0.0;
>+        for ( lcl_tSizeType i=1; i<=n-2; ++i )
>+        {
>+            fSum += Ddiag[ i ] * Rright[ i ] * Rright[ i ];
>+        }
>+        Ddiag[ n ] = Adiag[ n ] - fSum - Ddiag[ n-1 ] * Rupper[ n-1 ] * Rupper[ n-1 ]; // bug in [2]!
>+        
>+        // solve forward (R transpose)*z=u, overwrite u with z
>+        for ( lcl_tSizeType i=2; i<=n-1; ++i )
>+        {
>+            u[ i ] -= u[ i-1 ]* Rupper[ i-1 ];
>+        }
>+        fSum = 0.0;
>+        for ( lcl_tSizeType i=1; i<=n-2; ++i )
>+        {
>+            fSum += Rright[ i ] * u[ i ];
>+        }
>+        u[ n ] = u[ n ] - fSum - Rupper[ n - 1] * u[ n-1 ];
>+        
>+        // solve forward D*r=z, z is in u, overwrite u with r
>+        for ( lcl_tSizeType i=1; i<=n; ++i )
>+        {
>+            u[ i ] = u[i] / Ddiag[ i ];
>+        }
>+        
>+        // solve backward R*x= r, r is in u
>+        m_aSecDerivY[ n ] = u[ n ];
>+        m_aSecDerivY[ n-1 ] = u[ n-1 ] - Rupper[ n-1 ] * m_aSecDerivY[ n ];
>+        for ( lcl_tSizeType i=n-2; i>=1; --i)
>+        {
>+            m_aSecDerivY[ i ] = u[ i ] - Rupper[ i ] * m_aSecDerivY[ i+1 ] - Rright[ i ] * m_aSecDerivY[ n ];
>+        }
>+        // periodic
>+        m_aSecDerivY[ 0 ] = m_aSecDerivY[ n ];
>+    }
>+    
>+    // adapt m_aSecDerivY for usage in GetInterpolatedValue()
>+    for( lcl_tSizeType i = 0; i <= n ; ++i )
>+    {
>+        m_aSecDerivY[ i ] *= 2.0;
>+    }
>+
>+}
>+
> 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;
>+    const lcl_tSizeType n = m_aPoints.size() - 1;
>     if( x < m_fLastInterpolatedValue )
>     {
>         m_nKLow = 0;
>@@ -221,7 +392,7 @@
>         // first initialization is done in CTOR
>         while( m_nKHigh - m_nKLow > 1 )
>         {
>-            tPointVecType::size_type k = ( m_nKHigh + m_nKLow ) / 2;
>+            lcl_tSizeType k = ( m_nKHigh + m_nKLow ) / 2;
>             if( m_aPoints[ k ].first > x )
>                 m_nKHigh = k;
>             else
>@@ -255,63 +426,142 @@
> 
> //-----------------------------------------------------------------------------
> 
>-//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;
>+// helper methods for B-spline
>+
>+// Create parameter t_0 to t_n using the centripetal method with a power of 0.5
>+bool createParameterT(const tPointVecType aUniquePoints, double* t)
>+{   // precondition: no adjacent identical points
>+    // postcondition: 0 = t_0 < t_1 < ... < t_n = 1
>+    bool bIsSuccessful = true;
>+    const lcl_tSizeType n = aUniquePoints.size() - 1;
>+    t[0]=0.0;
>+    double dx = 0.0;
>+    double dy = 0.0;
>+    double fDiffMax = 1.0; //dummy values
>+    double fDenominator = 0.0; // initialized for summing up
>+    for (lcl_tSizeType i=1; i<=n ; ++i)
>+    {   // 4th root(dx^2+dy^2)
>+        dx = aUniquePoints[i].first - aUniquePoints[i-1].first;
>+        dy = aUniquePoints[i].second - aUniquePoints[i-1].second;
>+        // scaling to avoid underflow or overflow
>+        fDiffMax = (fabs(dx)>fabs(dy)) ? fabs(dx) : fabs(dy);
>+        if (fDiffMax == 0.0)
>+        {
>+            bIsSuccessful = false;
>+            break;
>+        }
>         else
>-            t[i] = n-k+2;
>+        {
>+            dx /= fDiffMax;
>+            dy /= fDiffMax;
>+            fDenominator += sqrt(sqrt(dx * dx + dy * dy)) * sqrt(fDiffMax);
>+        }
>     }
>-    return t;
>-}
>+    if (fDenominator == 0.0)
>+    {
>+        bIsSuccessful = false;
>+    }
>+    if (bIsSuccessful)
>+    {
>+        for (lcl_tSizeType j=1; j<=n ; ++j)
>+        {
>+            double fNumerator = 0.0;
>+            for (lcl_tSizeType i=1; i<=j ; ++i)
>+            {
>+                dx = aUniquePoints[i].first - aUniquePoints[i-1].first;
>+                dy = aUniquePoints[i].second - aUniquePoints[i-1].second;
>+                fDiffMax = (abs(dx)>abs(dy)) ? abs(dx) : abs(dy);
>+                // same as above, so should not be zero
>+                dx /= fDiffMax;
>+                dy /= fDiffMax;
>+                fNumerator += sqrt(sqrt(dx * dx + dy * dy)) * sqrt(fDiffMax);
>+            }
>+            t[j] = fNumerator / fDenominator;
> 
>-//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;
>+        }
>+        // postcondition check
>+        t[n] = 1.0;
>+        double fPrevious = 0.0;
>+        for (lcl_tSizeType i=1; i <= n && bIsSuccessful ; ++i)
>+        {
>+            if (fPrevious >= t[i])
>+            {
>+                bIsSuccessful = false;
>+            }
>+            else
>+            {
>+                fPrevious = t[i];
>+            }
>+        }
>+    }
>+    return bIsSuccessful;
> }
> 
>-//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;
>+void createKnotVector(const lcl_tSizeType n, const sal_uInt32 p, double* t, double* u)
>+{  // precondition: 0 = t_0 < t_1 < ... < t_n = 1
>+        for (lcl_tSizeType j = 0; j <= p; ++j)
>+        {
>+            u[j] = 0.0;
>+        }
>+        double fSum = 0.0;
>+        for (lcl_tSizeType j = 1; j <= n-p; ++j )
>+        {
>+            fSum = 0.0;
>+            for (lcl_tSizeType i = j; i <= j+p-1; ++i)
>+            {
>+                fSum += t[i];
>+            }
>+            u[j+p] = fSum / p ;
>+        }
>+        for (lcl_tSizeType j = n+1; j <= n+1+p; ++j)
>+        {
>+            u[j] = 1.0;
>+        }
> }
> 
>-//calculate weight vector
>-void BVector(double x, sal_Int32 n, sal_Int32 k, double *b, const double *t)
>+void applyNtoParameterT(const lcl_tSizeType i,const double tk,const sal_uInt32 p,const double* u, double* rowN)
> {
>-    sal_Int32 i = 0;
>-    for( i=0; i<=n+k; i++ )
>-        b[i]=0;
>+    // get N_p(t_k) recursively, only N_(i-p) till N_(i) are relevant, all other N_# are zero
>+    double fRightFactor = 0.0;
>+    double fLeftFactor = 0.0;
> 
>-    sal_Int32 i0 = (sal_Int32)floor(x) + k - 1;
>-    b [i0] = 1;
>+    // initialize with indicator function degree 0
>+    rowN[p] = 1.0; // all others are zero
> 
>-    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];
>+    // calculate up to degree p
>+    for (sal_uInt32 s = 1; s <= p; ++s)
>+    {
>+        // first element
>+        fRightFactor = ( u[i+1] - tk ) / ( u[i+1]- u[i-s+1] );
>+        // i-s "true index" - (i-p)"shift" = p-s
>+        rowN[p-s] = fRightFactor * rowN[p-s+1];
>+
>+        // middle elements
>+        for (sal_uInt32 j = s-1; j>=1 ; --j)
>+        {
>+            fLeftFactor = ( tk - u[i-j] ) / ( u[i-j+s] - u[i-j] ) ;
>+            fRightFactor = ( u[i-j+s+1] - tk ) / ( u[i-j+s+1] - u[i-j+1] );
>+            // i-j "true index" - (i-p)"shift" = p-j
>+            rowN[p-j] = fLeftFactor * rowN[p-j] + fRightFactor *  rowN[p-j+1];
>+        }
>+
>+        // last element
>+        fLeftFactor = ( tk - u[i] ) / ( u[i+s] - u[i] );
>+        // i "true index" - (i-p)"shift" = p
>+        rowN[p] = fLeftFactor * rowN[p];
>+    }
> }
> 
> } //  anonymous namespace
> 
> //-----------------------------------------------------------------------------
> 
>+// Calculates uniform parametric splines with subinterval length 1,
>+// according ODF1.2 part 1, chapter 'chart interpolation'.
> void SplineCalculater::CalculateCubicSplines(
>     const drawing::PolyPolygonShape3D& rInput
>     , drawing::PolyPolygonShape3D& rResult
>-    , sal_Int32 nGranularity )
>+    , sal_uInt32 nGranularity )
> {
>     DBG_ASSERT( nGranularity > 0, "Granularity is invalid" );
> 
>@@ -319,7 +569,7 @@
>     rResult.SequenceY.realloc(0);
>     rResult.SequenceZ.realloc(0);
> 
>-    sal_Int32 nOuterCount = rInput.SequenceX.getLength();
>+    sal_uInt32 nOuterCount = rInput.SequenceX.getLength();
>     if( !nOuterCount )
>         return;
> 
>@@ -327,30 +577,23 @@
>     rResult.SequenceY.realloc(nOuterCount);
>     rResult.SequenceZ.realloc(nOuterCount);
> 
>-    for( sal_Int32 nOuter = 0; nOuter < nOuterCount; ++nOuter )
>+    for( sal_uInt32 nOuter = 0; nOuter < nOuterCount; ++nOuter )
>     {
>         if( rInput.SequenceX[nOuter].getLength() <= 1 )
>             continue; //we need at least two points
> 
>-        sal_Int32 nMaxIndexPoints = rInput.SequenceX[nOuter].getLength()-1; // is >=1
>+        sal_uInt32 nMaxIndexPoints = rInput.SequenceX[nOuter].getLength()-1; // is >=1
>         const double* pOldX = rInput.SequenceX[nOuter].getConstArray();
>         const double* pOldY = rInput.SequenceY[nOuter].getConstArray();
>         const double* pOldZ = rInput.SequenceZ[nOuter].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++ )
>+        for( sal_uInt32 nIndex=1; nIndex<=nMaxIndexPoints; nIndex++ )
>         {
>-            // The euclidian distance leads to curve loops for functions having single extreme points
>-            // 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);
>@@ -358,7 +601,7 @@
>         aInputY.resize(nMaxIndexPoints+1);
>         tPointVecType aInputZ;
>         aInputZ.resize(nMaxIndexPoints+1);
>-        for (sal_Int32 nN=0;nN<=nMaxIndexPoints; nN++ )
>+        for (sal_uInt32 nN=0;nN<=nMaxIndexPoints; nN++ )
>         {
>           aInputX[ nN ].first=aParameter[nN];
>           aInputX[ nN ].second=pOldX[ nN ];
>@@ -373,16 +616,20 @@
>         double fXDerivation;
>         double fYDerivation;
>         double fZDerivation;
>+        lcl_SplineCalculation* aSplineX;
>+        lcl_SplineCalculation* aSplineY;
>+        // lcl_SplineCalculation* aSplineZ; the z-coordinates of all points in
>+        // a data series are equal. No spline calculation needed, but copy
>+        // coordinate to output
>+        
>         if( pOldX[ 0 ] == pOldX[nMaxIndexPoints] &&
>             pOldY[ 0 ] == pOldY[nMaxIndexPoints] &&
>-            pOldZ[ 0 ] == pOldZ[nMaxIndexPoints] )
>-        {
>-            // #i101050# avoid a corner in closed lines, which are smoothed by spline
>-            // This derivation are special for parameter of kind 0,1,2,3... If you
>-            // change generating parameters (see above), then adapt derivations too.)
>-            fXDerivation = 0.5 * (pOldX[1]-pOldX[nMaxIndexPoints-1]);
>-            fYDerivation = 0.5 * (pOldY[1]-pOldY[nMaxIndexPoints-1]);
>-            fZDerivation = 0.5 * (pOldZ[1]-pOldZ[nMaxIndexPoints-1]);
>+            pOldZ[ 0 ] == pOldZ[nMaxIndexPoints] &&
>+            nMaxIndexPoints >=2 )
>+        {   // periodic spline
>+            aSplineX = new lcl_SplineCalculation( aInputX) ;
>+            aSplineY = new lcl_SplineCalculation( aInputY) ;
>+            // aSplineZ = new lcl_SplineCalculation( aInputZ) ;
>         }
>         else // generate the kind "natural spline"
>         {
>@@ -391,10 +638,10 @@
>             fXDerivation = fInfty;
>             fYDerivation = fInfty;
>             fZDerivation = fInfty;
>+            aSplineX = new lcl_SplineCalculation( aInputX, fXDerivation, fXDerivation );
>+            aSplineY = new lcl_SplineCalculation( aInputY, fYDerivation, fYDerivation );
>+            // aSplineZ = new lcl_SplineCalculation( aInputZ, fZDerivation, fZDerivation );
>         }
>-        lcl_SplineCalculation aSplineX( aInputX, fXDerivation, fXDerivation );
>-        lcl_SplineCalculation aSplineY( aInputY, fYDerivation, fYDerivation );
>-        lcl_SplineCalculation aSplineZ( aInputZ, fZDerivation, fZDerivation );
> 
>         // fill result polygon with calculated values
>         rResult.SequenceX[nOuter].realloc( nMaxIndexPoints*nGranularity + 1);
>@@ -405,13 +652,13 @@
>         double* pNewY = rResult.SequenceY[nOuter].getArray();
>         double* pNewZ = rResult.SequenceZ[nOuter].getArray();
> 
>-        sal_Int32 nNewPointIndex = 0; // Index in result points
>+        sal_uInt32 nNewPointIndex = 0; // Index in result points
>         // needed for inner loop
>         double    fInc;   // step for intermediate points
>-        sal_Int32 nj;     // for loop
>+        sal_uInt32 nj;     // for loop
>         double    fParam; // a intermediate parameter value
> 
>-        for( sal_Int32 ni = 0; ni < nMaxIndexPoints; ni++ )
>+        for( sal_uInt32 ni = 0; ni < nMaxIndexPoints; ni++ )
>         {
>             // given point is surely a curve point
>             pNewX[nNewPointIndex] = pOldX[ni];
>@@ -425,9 +672,10 @@
>             {
>                 fParam = aParameter[ni] + ( fInc * static_cast< double >( nj ) );
> 
>-                pNewX[nNewPointIndex]=aSplineX.GetInterpolatedValue( fParam );
>-                pNewY[nNewPointIndex]=aSplineY.GetInterpolatedValue( fParam );
>-                pNewZ[nNewPointIndex]=aSplineZ.GetInterpolatedValue( fParam );
>+                pNewX[nNewPointIndex]=aSplineX->GetInterpolatedValue( fParam );
>+                pNewY[nNewPointIndex]=aSplineY->GetInterpolatedValue( fParam );
>+                // pNewZ[nNewPointIndex]=aSplineZ->GetInterpolatedValue( fParam );
>+                pNewZ[nNewPointIndex] = pOldZ[ni];
>                 nNewPointIndex++;
>             }
>         }
>@@ -435,23 +683,39 @@
>         pNewX[nNewPointIndex] = pOldX[nMaxIndexPoints];
>         pNewY[nNewPointIndex] = pOldY[nMaxIndexPoints];
>         pNewZ[nNewPointIndex] = pOldZ[nMaxIndexPoints];
>+        delete aSplineX;
>+        delete aSplineY;
>+        // delete aSplineZ;
>     }
> }
> 
>+
>+// The implementation follows closely ODF1.2 spec, chapter chart:interpolation
>+// using the same names as in spec as far as possible, without prefix.
>+// More details can be found on
>+// Dr. C.-K. Shene: CS3621 Introduction to Computing with Geometry Notes
>+// Unit 9: Interpolation and Approximation/Curve Global Interpolation
>+// Department of Computer Science, Michigan Technological University
>+// http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/
>+// [last called 2011-05-20]
> void SplineCalculater::CalculateBSplines(
>             const ::com::sun::star::drawing::PolyPolygonShape3D& rInput
>             , ::com::sun::star::drawing::PolyPolygonShape3D& rResult
>-            , sal_Int32 nGranularity
>-            , sal_Int32 nDegree )
>+            , sal_uInt32 nResolution
>+			, sal_uInt32 nDegree )
> {
>-    // #issue 72216#
>-    // k is the order of the BSpline, nDegree is the degree of its polynoms
>-    sal_Int32 k = nDegree + 1;
>+    // nResolution is ODF1.2 file format attribut chart:spline-resolution and
>+    // ODF1.2 spec variable k. Causion, k is used as index in the spec in addition.
>+    // nDegree is ODF1.2 file format attribut chart:spline-order and
>+    // ODF1.2 spec variable p
>+    DBG_ASSERT( nResolution > 1, "resolution 1 is equivalent to attribut none" );
>+    DBG_ASSERT( nDegree >= 1, "degree of polynomials must be positive");
>+    sal_uInt32 p = nDegree;
>+
>+	rResult.SequenceX.realloc(0);
>+	rResult.SequenceY.realloc(0);
>+	rResult.SequenceZ.realloc(0);
> 
>-    rResult.SequenceX.realloc(0);
>-    rResult.SequenceY.realloc(0);
>-    rResult.SequenceZ.realloc(0);
>-    
>     sal_Int32 nOuterCount = rInput.SequenceX.getLength();
>     if( !nOuterCount )
>         return; // no input
>@@ -463,74 +727,257 @@
>     for( sal_Int32 nOuter = 0; nOuter < nOuterCount; ++nOuter )
>     {
>         if( rInput.SequenceX[nOuter].getLength() <= 1 )
>-            continue; // need at least 2 control points
>-        
>-        sal_Int32 n = rInput.SequenceX[nOuter].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
>+            continue; // need at least 2 points, next piece of the series
> 
>+        // Copy input to vector of points and remove adjacent double points. The
>+        // Z-coordinate is equal for all points in a series and holds the depth
>+        // in 3D mode, simple copying is enough.
>+        lcl_tSizeType nMaxIndexPoints = rInput.SequenceX[nOuter].getLength()-1; // is >=1
>         const double* pOldX = rInput.SequenceX[nOuter].getConstArray();
>         const double* pOldY = rInput.SequenceY[nOuter].getConstArray();
>         const double* pOldZ = rInput.SequenceZ[nOuter].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 fZCoordinate = pOldZ[0];
>+        tPointVecType aPointsIn;
>+        aPointsIn.resize(nMaxIndexPoints+1);
>+        for (lcl_tSizeType i = 0; i <= nMaxIndexPoints; ++i )
>+        {
>+          aPointsIn[ i ].first = pOldX[i];
>+          aPointsIn[ i ].second = pOldY[i];
>+        }
>+        aPointsIn.erase( ::std::unique( aPointsIn.begin(), aPointsIn.end()),
>+                     aPointsIn.end() );
> 
>-        double *b       = new double [n + k + 1]; // values of blending functions
>+        // n is the last valid index to the reduced aPointsIn
>+        // There are n+1 valid data points.
>+        const lcl_tSizeType n = aPointsIn.size() - 1;
>+        if (n < 1 || p > n)
>+            continue; // need at least 2 points, degree p needs at least n+1 points
>+                      // next piece of series
> 
>-        const double* t = createTVector(n, k); // knot vector
>+        double* t = new double [n+1];
>+        if (!createParameterT(aPointsIn, t))
>+        {
>+            delete[] t;
>+            continue; // next piece of series
>+        }
> 
>-        rResult.SequenceX[nOuter].realloc(nNewSectorCount+1);
>-        rResult.SequenceY[nOuter].realloc(nNewSectorCount+1);
>-        rResult.SequenceZ[nOuter].realloc(nNewSectorCount+1);
>-        double* pNewX = rResult.SequenceX[nOuter].getArray();
>-        double* pNewY = rResult.SequenceY[nOuter].getArray();
>-        double* pNewZ = rResult.SequenceZ[nOuter].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
>+        lcl_tSizeType m = n + p + 1;
>+        double* u = new double [m+1];
>+        createKnotVector(n, p, t, u);
> 
>-        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++)
>+        // The matrix N contains the B-spline basis functions applied to parameters.
>+        // In each row only p+1 adjacent elements are non-zero. The starting
>+        // column in a higher row is equal or greater than in the lower row.
>+        // To store this matrix the non-zero elements are shifted to column 0
>+        // and the amount of shifting is remembered in an array.
>+        double** aMatN = new double*[n+1];
>+        for (lcl_tSizeType row = 0; row <=n; ++row)
>+        {
>+            aMatN[row] = new double[p+1];
>+            for (sal_uInt32 col = 0; col <= p; ++col)
>+            aMatN[row][col] = 0.0;
>+        }
>+        lcl_tSizeType* aShift = new lcl_tSizeType[n+1];
>+        aMatN[0][0] = 1.0; //all others are zero
>+        aShift[0] = 0;
>+        aMatN[n][0] = 1.0;
>+        aShift[n] = n;
>+        for (lcl_tSizeType k = 1; k<=n-1; ++k)
>+        { // all basis functions are applied to t_k,
>+            // results are elements in row k in matrix N
>+
>+            // find the one interval with u_i <= t_k < u_(i+1)
>+            // remember u_0 = ... = u_p = 0.0 and u_(m-p) = ... u_m = 1.0 and 0<t_k<1
>+            lcl_tSizeType i = p;
>+            while (!(u[i] <= t[k] && t[k] < u[i+1]))
>             {
>-                fX +=pOldX[nPointIndex]*b[nPointIndex];
>-                fY +=pOldY[nPointIndex]*b[nPointIndex];
>-                fZ +=pOldZ[nPointIndex]*b[nPointIndex];
>+                ++i;
>             }
>-            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];
> 
>+            // index in reduced matrix aMatN = (index in full matrix N) - (i-p)
>+            aShift[k] = i - p;
>+
>+            applyNtoParameterT(i, t[k], p, u, aMatN[k]);
>+        } // next row k
>+
>+        // Get matrix C of control points from the matrix equation aMatN * C = aPointsIn
>+        // aPointsIn is overwritten with C.
>+        // Gaussian elimination is possible without pivoting, see reference
>+        lcl_tSizeType r = 0; // true row index
>+        lcl_tSizeType c = 0; // true column index
>+        double fDivisor = 1.0; // used for diagonal element
>+        double fEliminate = 1.0; // used for the element, that will become zero
>+        double fHelp;
>+        tPointType aHelp;
>+        lcl_tSizeType nHelp; // used in triangle change
>+        bool bIsSuccessful = true;
>+        for (c = 0 ; c <= n && bIsSuccessful; ++c)
>+        {
>+            // search for first non-zero downwards
>+            r = c;
>+            while ( aMatN[r][c-aShift[r]] == 0 &&  r < n)
>+            {
>+                ++r;
>+            }
>+            if (aMatN[r][c-aShift[r]] == 0.0)
>+            {
>+                // Matrix N is singular, although this is mathematically impossible
>+                bIsSuccessful = false;
>+            }
>+            else
>+            {
>+                // exchange total row r with total row c if necessary
>+                if (r != c)
>+                {
>+                    for ( sal_uInt32 i = 0; i <= p ; ++i)
>+                    {
>+                        fHelp = aMatN[r][i];
>+                        aMatN[r][i] = aMatN[c][i];
>+                        aMatN[c][i] = fHelp;
>+                    }
>+                    aHelp = aPointsIn[r];
>+                    aPointsIn[r] = aPointsIn[c];
>+                    aPointsIn[c] = aHelp;
>+                    nHelp = aShift[r];
>+                    aShift[r] = aShift[c];
>+                    aShift[c] = nHelp;
>+                }
>+
>+                // divide row c, so that element(c,c) becomes 1
>+                fDivisor = aMatN[c][c-aShift[c]]; // not zero, see above
>+                for (sal_uInt32 i = 0; i <= p; ++i)
>+                {
>+                    aMatN[c][i] /= fDivisor;
>+                }
>+                aPointsIn[c].first /= fDivisor;
>+                aPointsIn[c].second /= fDivisor;
>+
>+                // eliminate forward, examine row c+1 to n-1 (worst case)
>+                // stop if first non-zero element in row has an higher column as c
>+                // look at nShift for that, elements in nShift are equal or increasing
>+                for ( r = c+1; aShift[r]<=c && r < n; ++r)
>+                {
>+                    fEliminate = aMatN[r][0];
>+                    if (fEliminate != 0.0) // else accidentally zero, nothing to do
>+                    {
>+                        for (sal_uInt32 i = 1; i <= p; ++i)
>+                        {
>+                            aMatN[r][i-1] = aMatN[r][i] - fEliminate * aMatN[c][i];
>+                        }
>+                        aMatN[r][p]=0;
>+                        aPointsIn[r].first -= fEliminate * aPointsIn[c].first;
>+                        aPointsIn[r].second -= fEliminate * aPointsIn[c].second;
>+                        ++aShift[r];
>+                    }
>+                }
>+            }
>+        }// upper triangle form is reached
>+        if( bIsSuccessful)
>+        {
>+            // eliminate backwards, begin with last column
>+            for (lcl_tSizeType c = n; c >= 1; --c )
>+            {
>+                // In row c the diagonal element(c,c) == 1 and all elements left from
>+                // diagonal are zero and do not influence other rows.
>+                // Full matrix N has semibandwidth < p, therefore element(r,c) is
>+                // zero, if abs(r-c)>=p.  abs(r-c)=c-r, because r<c.
>+                r = c - 1;
>+                while ( r !=0 && c-r < p )
>+                {
>+                    fEliminate = aMatN[r][ c - aShift[r] ];
>+                    if ( fEliminate != 0.0) // else element is accidentically zero, no action needed
>+                    {
>+                        // row r -= fEliminate * row c only relevant for right side
>+                        aMatN[r][c - aShift[r]] = 0.0;
>+                        aPointsIn[r].first -= fEliminate * aPointsIn[c].first;
>+                        aPointsIn[r].second -= fEliminate * aPointsIn[c].second;
>+                    }
>+                    --r;
>+                }
>+            }
>+        }   // aPointsIn contains the control points now.
>+        if (bIsSuccessful)
>+        {
>+            // calculate the intermediate points according given resolution
>+            // using deBoor-Cox algorithm
>+            lcl_tSizeType nNewSize = nResolution * n + 1;
>+            rResult.SequenceX[nOuter].realloc(nNewSize);
>+	        rResult.SequenceY[nOuter].realloc(nNewSize);
>+	        rResult.SequenceZ[nOuter].realloc(nNewSize);
>+	        double* pNewX = rResult.SequenceX[nOuter].getArray();
>+	        double* pNewY = rResult.SequenceY[nOuter].getArray();
>+	        double* pNewZ = rResult.SequenceZ[nOuter].getArray();
>+	        pNewX[0] = aPointsIn[0].first;
>+	        pNewY[0] = aPointsIn[0].second;
>+	        pNewZ[0] = fZCoordinate; // Precondition: z-coordinates of all points of a series are equal
>+	        pNewX[nNewSize -1 ] = aPointsIn[n].first;
>+	        pNewY[nNewSize -1 ] = aPointsIn[n].second;
>+	        pNewZ[nNewSize -1 ] = fZCoordinate;
>+            double* aP = new double[m+1];
>+            lcl_tSizeType nLow = 0;
>+	        for ( lcl_tSizeType nTIndex = 0; nTIndex <= n-1; ++nTIndex)
>+	        {
>+                for (sal_uInt32 nResolutionStep = 1;
>+                     nResolutionStep <= nResolution && !( nTIndex == n-1 && nResolutionStep == nResolution);
>+                     ++nResolutionStep)
>+                {
>+                    lcl_tSizeType nNewIndex = nTIndex * nResolution + nResolutionStep;
>+                    double ux = t[nTIndex] + nResolutionStep * ( t[nTIndex+1] - t[nTIndex]) /nResolution;
>+
>+                    // get index nLow, so that u[nLow]<= ux < u[nLow +1]
>+                    // continue from previous nLow
>+                    while ( u[nLow] <= ux)
>+                    {
>+                        ++nLow;
>+                    }
>+                    --nLow;
>+
>+                    // x-coordinate
>+                    for (lcl_tSizeType i = nLow-p; i <= nLow; ++i)
>+                    {
>+                        aP[i] = aPointsIn[i].first;
>+                    }
>+                    for (sal_uInt32 lcl_Degree = 1; lcl_Degree <= p; ++lcl_Degree)
>+                    {
>+                        double fFactor = 0.0;
>+                        for (lcl_tSizeType i = nLow; i >= nLow + lcl_Degree - p; --i)
>+                        {
>+                            fFactor = ( ux - u[i] ) / ( u[i+p+1-lcl_Degree] - u[i]);
>+                            aP[i] = (1 - fFactor)* aP[i-1] + fFactor * aP[i];
>+                        }
>+                    }
>+                    pNewX[nNewIndex] = aP[nLow];
>+
>+                    // y-coordinate
>+                    for (lcl_tSizeType i = nLow - p; i <= nLow; ++i)
>+                    {
>+                        aP[i] = aPointsIn[i].second;
>+                    }
>+                    for (sal_uInt32 lcl_Degree = 1; lcl_Degree <= p; ++lcl_Degree)
>+                    {
>+                        double fFactor = 0.0;
>+                        for (lcl_tSizeType i = nLow; i >= nLow +lcl_Degree - p; --i)
>+                        {
>+                            fFactor = ( ux - u[i] ) / ( u[i+p+1-lcl_Degree] - u[i]);
>+                            aP[i] = (1 - fFactor)* aP[i-1] + fFactor * aP[i];
>+                        }
>+                    }
>+                    pNewY[nNewIndex] = aP[nLow];
>+                    pNewZ[nNewIndex] = fZCoordinate;
>+                }
>+            }
>+            delete[] aP;
>+        }
>+        delete[] aShift;
>+        for (lcl_tSizeType row = 0; row <=n; ++row)
>+        {
>+            delete[] aMatN[row];
>+        }
>+        delete[] aMatN;
>+        delete[] u;
>         delete[] t;
>-        delete[] b;
>-    }
>+
>+    } // next piece of the series
> }
> 
> //.............................................................................

Actions: View | Splinter Review

Attachments on bug 37889: 47488 | 47489 | 47490 | 47491

Use of freedesktop.org services, including Bugzilla, is subject to our Code of Conduct. How we collect and use information is described in our Privacy Policy.