--- 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=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= 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= 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 } //.............................................................................