| /************************************************************** |
| * |
| * Licensed to the Apache Software Foundation (ASF) under one |
| * or more contributor license agreements. See the NOTICE file |
| * distributed with this work for additional information |
| * regarding copyright ownership. The ASF licenses this file |
| * to you under the Apache License, Version 2.0 (the |
| * "License"); you may not use this file except in compliance |
| * with the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, |
| * software distributed under the License is distributed on an |
| * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY |
| * KIND, either express or implied. See the License for the |
| * specific language governing permissions and limitations |
| * under the 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 |
| { |
| |
| 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 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 |
| 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 |
| lcl_tSizeType m_nKLow; |
| lcl_tSizeType 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(); |
| |
| /** @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(); |
| }; |
| |
| //----------------------------------------------------------------------------- |
| |
| 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 ); |
| 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 lcl_tSizeType 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( lcl_tSizeType i = 1; i < n; ++i ) |
| { |
| tPointType |
| 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( 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 ) |
| { |
| OSL_PRECOND( ( m_aPoints[ 0 ].first <= x ) && |
| ( x <= m_aPoints[ m_aPoints.size() - 1 ].first ), |
| "Trying to extrapolate" ); |
| |
| const lcl_tSizeType 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 ) |
| { |
| lcl_tSizeType 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; |
| } |
| OSL_ENSURE( m_nKHigh <= n, "Out of Bounds" ); |
| } |
| m_fLastInterpolatedValue = x; |
| |
| double h = m_aPoints[ m_nKHigh ].first - m_aPoints[ m_nKLow ].first; |
| OSL_ENSURE( 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 ); |
| } |
| |
| //----------------------------------------------------------------------------- |
| |
| // 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 |
| { |
| dx /= fDiffMax; |
| dy /= fDiffMax; |
| fDenominator += sqrt(sqrt(dx * dx + dy * dy)) * sqrt(fDiffMax); |
| } |
| } |
| 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; |
| |
| } |
| // 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; |
| } |
| |
| 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; |
| } |
| } |
| |
| void applyNtoParameterT(const lcl_tSizeType i,const double tk,const sal_uInt32 p,const double* u, double* rowN) |
| { |
| // 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; |
| |
| // initialize with indicator function degree 0 |
| rowN[p] = 1.0; // all others are zero |
| |
| // 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_uInt32 nGranularity ) |
| { |
| OSL_PRECOND( nGranularity > 0, "Granularity is invalid" ); |
| |
| rResult.SequenceX.realloc(0); |
| rResult.SequenceY.realloc(0); |
| rResult.SequenceZ.realloc(0); |
| |
| sal_uInt32 nOuterCount = rInput.SequenceX.getLength(); |
| if( !nOuterCount ) |
| return; |
| |
| rResult.SequenceX.realloc(nOuterCount); |
| rResult.SequenceY.realloc(nOuterCount); |
| rResult.SequenceZ.realloc(nOuterCount); |
| |
| for( sal_uInt32 nOuter = 0; nOuter < nOuterCount; ++nOuter ) |
| { |
| if( rInput.SequenceX[nOuter].getLength() <= 1 ) |
| continue; //we need at least two points |
| |
| 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(); |
| |
| ::std::vector < double > aParameter(nMaxIndexPoints+1); |
| aParameter[0]=0.0; |
| for( sal_uInt32 nIndex=1; nIndex<=nMaxIndexPoints; nIndex++ ) |
| { |
| 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_uInt32 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 |
| 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] && |
| 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" |
| { |
| double fInfty; |
| ::rtl::math::setInf( &fInfty, sal_False ); |
| 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 ); |
| } |
| |
| // fill result polygon with calculated values |
| rResult.SequenceX[nOuter].realloc( nMaxIndexPoints*nGranularity + 1); |
| rResult.SequenceY[nOuter].realloc( nMaxIndexPoints*nGranularity + 1); |
| rResult.SequenceZ[nOuter].realloc( nMaxIndexPoints*nGranularity + 1); |
| |
| double* pNewX = rResult.SequenceX[nOuter].getArray(); |
| double* pNewY = rResult.SequenceY[nOuter].getArray(); |
| double* pNewZ = rResult.SequenceZ[nOuter].getArray(); |
| |
| sal_uInt32 nNewPointIndex = 0; // Index in result points |
| // needed for inner loop |
| double fInc; // step for intermediate points |
| sal_uInt32 nj; // for loop |
| double fParam; // a intermediate parameter value |
| |
| for( sal_uInt32 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 ); |
| pNewZ[nNewPointIndex] = pOldZ[ni]; |
| nNewPointIndex++; |
| } |
| } |
| // add last point |
| 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_uInt32 nResolution |
| , sal_uInt32 nDegree ) |
| { |
| // 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 |
| OSL_ASSERT( nResolution > 1 ); |
| OSL_ASSERT( nDegree >= 1 ); |
| sal_uInt32 p = nDegree; |
| |
| rResult.SequenceX.realloc(0); |
| rResult.SequenceY.realloc(0); |
| rResult.SequenceZ.realloc(0); |
| |
| sal_Int32 nOuterCount = rInput.SequenceX.getLength(); |
| if( !nOuterCount ) |
| return; // no input |
| |
| rResult.SequenceX.realloc(nOuterCount); |
| rResult.SequenceY.realloc(nOuterCount); |
| rResult.SequenceZ.realloc(nOuterCount); |
| |
| for( sal_Int32 nOuter = 0; nOuter < nOuterCount; ++nOuter ) |
| { |
| if( rInput.SequenceX[nOuter].getLength() <= 1 ) |
| 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(); |
| 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() ); |
| |
| // 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 |
| |
| double* t = new double [n+1]; |
| if (!createParameterT(aPointsIn, t)) |
| { |
| delete[] t; |
| continue; // next piece of series |
| } |
| |
| lcl_tSizeType m = n + p + 1; |
| double* u = new double [m+1]; |
| createKnotVector(n, p, t, u); |
| |
| // 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])) |
| { |
| ++i; |
| } |
| |
| // 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 cc = n; cc >= 1; --cc ) |
| { |
| // In row cc the diagonal element(cc,cc) == 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-cc)>=p. abs(r-cc)=cc-r, because r<cc. |
| r = cc - 1; |
| while ( r !=0 && cc-r < p ) |
| { |
| fEliminate = aMatN[r][ cc - aShift[r] ]; |
| if ( fEliminate != 0.0) // else element is accidentically zero, no action needed |
| { |
| // row r -= fEliminate * row cc only relevant for right side |
| aMatN[r][cc - aShift[r]] = 0.0; |
| aPointsIn[r].first -= fEliminate * aPointsIn[cc].first; |
| aPointsIn[r].second -= fEliminate * aPointsIn[cc].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; |
| |
| } // next piece of the series |
| } |
| |
| //............................................................................. |
| } //namespace chart |
| //............................................................................. |
| |
