Bugzilla – Attachment 47490 Details for
Bug 37889
Adapt smoothing with splines to ODF1.2
Home
|
New
|
Browse
|
Search
|
[?]
|
Reports
|
Help
|
New Account
|
Log In
[x]
|
Forgot Password
Login:
[x]
complete new Splines.cxx for better reading
Splines.cxx (text/plain), 36.66 KB, created by
Regina Henschel
on 2011-06-03 12:16:32 UTC
(
hide
)
Description:
complete new Splines.cxx for better reading
Filename:
MIME Type:
Creator:
Regina Henschel
Created:
2011-06-03 12:16:32 UTC
Size:
36.66 KB
patch
obsolete
>/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ >/************************************************************************* > * > * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. > * > * Copyright 2000, 2010 Oracle and/or its affiliates. > * > * OpenOffice.org - a multi-platform office productivity suite > * > * This file is part of OpenOffice.org. > * > * OpenOffice.org is free software: you can redistribute it and/or modify > * it under the terms of the GNU Lesser General Public License version 3 > * only, as published by the Free Software Foundation. > * > * OpenOffice.org is distributed in the hope that it will be useful, > * but WITHOUT ANY WARRANTY; without even the implied warranty of > * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the > * GNU Lesser General Public License version 3 for more details > * (a copy is included in the LICENSE file that accompanied this code). > * > * You should have received a copy of the GNU Lesser General Public License > * version 3 along with OpenOffice.org. If not, see > * <http://www.openoffice.org/license.html> > * for a copy of the LGPLv3 License. > * > ************************************************************************/ > > >// MARKER(update_precomp.py): autogen include statement, do not remove >#include "precompiled_chart2.hxx" > >#include "Splines.hxx" >#include <rtl/math.hxx> > >#include <vector> >#include <algorithm> >#include <functional> > >// header for DBG_ASSERT >#include <tools/debug.hxx> > >//............................................................................. >namespace chart >{ >//............................................................................. >using namespace ::com::sun::star; > >namespace >{ > >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 ) >{ > DBG_ASSERT( ( 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; > } > DBG_ASSERT( m_nKHigh <= n, "Out of Bounds" ); > } > m_fLastInterpolatedValue = x; > > double h = m_aPoints[ m_nKHigh ].first - m_aPoints[ m_nKLow ].first; > DBG_ASSERT( h != 0, "Bad input to GetInterpolatedValue()" ); > > double a = ( m_aPoints[ m_nKHigh ].first - x ) / h; > double b = ( x - m_aPoints[ m_nKLow ].first ) / h; > > return ( a * m_aPoints[ m_nKLow ].second + > b * m_aPoints[ m_nKHigh ].second + > (( a*a*a - a ) * m_aSecDerivY[ m_nKLow ] + > ( b*b*b - b ) * m_aSecDerivY[ m_nKHigh ] ) * > ( h*h ) / 6.0 ); >} > >//----------------------------------------------------------------------------- > >// 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 ) >{ > DBG_ASSERT( 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 > 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); > > 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 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; > > } // next piece of the series >} > >//............................................................................. >} //namespace chart >//............................................................................. > >/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
<b>You cannot view the attachment while viewing its details because your browser does not support IFRAMEs. <a href="attachment.cgi?id=47490">View the attachment on a separate page</a>.</b>
View Attachment As Raw
Actions:
View
Attachments on
bug 37889
:
47488
|
47489
| 47490 |
47491