[Libreoffice-commits] .: chart2/source
Fridrich Strba
fridrich at kemper.freedesktop.org
Mon Jun 6 10:32:28 PDT 2011
chart2/source/view/charttypes/Splines.cxx | 759 +++++++++++++++++++++++-------
chart2/source/view/charttypes/Splines.hxx | 6
2 files changed, 606 insertions(+), 159 deletions(-)
New commits:
commit f18c69f3f40bf443b59bd006db220267ba5174d0
Author: Regina Henschel <rb.henschel at t-online.de>
Date: Mon Jun 6 19:31:56 2011 +0200
Adapt smoothing with splines to ODF1.2
diff --git a/chart2/source/view/charttypes/Splines.cxx b/chart2/source/view/charttypes/Splines.cxx
index 6376eb1..52f2020 100644
--- a/chart2/source/view/charttypes/Splines.cxx
+++ b/chart2/source/view/charttypes/Splines.cxx
@@ -46,17 +46,6 @@ using namespace ::com::sun::star;
namespace
{
-template< typename T >
-struct lcl_EqualsFirstDoubleOfPair : ::std::binary_function< ::std::pair< double, T >, ::std::pair< double, T >, bool >
-{
- inline bool operator() ( const ::std::pair< double, T > & rOne, const ::std::pair< double, T > & rOther )
- {
- return ( ::rtl::math::approxEqual( rOne.first, rOther.first ) );
- }
-};
-
-//-----------------------------------------------------------------------------
-
typedef ::std::pair< double, double > tPointType;
typedef ::std::vector< tPointType > tPointVecType;
typedef tPointVecType::size_type lcl_tSizeType;
@@ -79,6 +68,16 @@ public:
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
@@ -98,8 +97,8 @@ private:
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].
@@ -109,6 +108,22 @@ private:
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();
};
//-----------------------------------------------------------------------------
@@ -124,21 +139,32 @@ lcl_SplineCalculation::lcl_SplineCalculation(
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 );
@@ -156,9 +182,9 @@ void lcl_SplineCalculation::Calculate()
((( 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 ];
@@ -196,19 +222,164 @@ void lcl_SplineCalculation::Calculate()
// 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 )
{
OSL_PRECOND( ( 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;
@@ -218,7 +389,7 @@ double lcl_SplineCalculation::GetInterpolatedValue( double x )
// 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
@@ -252,63 +423,142 @@ double lcl_SplineCalculation::GetInterpolatedValue( double x )
//-----------------------------------------------------------------------------
-//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;
+
+ // initialize with indicator function degree 0
+ rowN[p] = 1.0; // all others are zero
- sal_Int32 i0 = (sal_Int32)floor(x) + k - 1;
- b [i0] = 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];
+ }
- 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];
+ // 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 )
{
OSL_PRECOND( nGranularity > 0, "Granularity is invalid" );
@@ -316,7 +566,7 @@ void SplineCalculater::CalculateCubicSplines(
rResult.SequenceY.realloc(0);
rResult.SequenceZ.realloc(0);
- sal_Int32 nOuterCount = rInput.SequenceX.getLength();
+ sal_uInt32 nOuterCount = rInput.SequenceX.getLength();
if( !nOuterCount )
return;
@@ -324,30 +574,23 @@ void SplineCalculater::CalculateCubicSplines(
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);
@@ -355,7 +598,7 @@ void SplineCalculater::CalculateCubicSplines(
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 ];
@@ -370,16 +613,20 @@ void SplineCalculater::CalculateCubicSplines(
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"
{
@@ -388,10 +635,10 @@ void SplineCalculater::CalculateCubicSplines(
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);
@@ -402,13 +649,13 @@ void SplineCalculater::CalculateCubicSplines(
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];
@@ -422,9 +669,10 @@ void SplineCalculater::CalculateCubicSplines(
{
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++;
}
}
@@ -432,23 +680,39 @@ void SplineCalculater::CalculateCubicSplines(
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
+ 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
@@ -460,74 +724,257 @@ void SplineCalculater::CalculateBSplines(
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 *b = new double [n + k + 1]; // values of blending functions
+ 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
+ }
- const double* t = createTVector(n, k); // knot vector
+ 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;
+ }
- 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
-
- 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++)
+ // 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)
{
- fX +=pOldX[nPointIndex]*b[nPointIndex];
- fY +=pOldY[nPointIndex]*b[nPointIndex];
- fZ +=pOldZ[nPointIndex]*b[nPointIndex];
+ // Matrix N is singular, although this is mathematically impossible
+ bIsSuccessful = false;
}
- pNewX[nNewSector] = fX;
- pNewY[nNewSector] = fY;
- pNewZ[nNewSector] = fZ;
-
- fCurveparam += fCurveStep; //for next looping
+ 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;
}
- // add last control point to BSpline curve
- pNewX[nNewSectorCount] = pOldX[n];
- pNewY[nNewSectorCount] = pOldY[n];
- pNewZ[nNewSectorCount] = pOldZ[n];
-
+ 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
}
//.............................................................................
diff --git a/chart2/source/view/charttypes/Splines.hxx b/chart2/source/view/charttypes/Splines.hxx
index 27806f2..0aad52a 100644
--- a/chart2/source/view/charttypes/Splines.hxx
+++ b/chart2/source/view/charttypes/Splines.hxx
@@ -46,13 +46,13 @@ public:
static void CalculateCubicSplines(
const ::com::sun::star::drawing::PolyPolygonShape3D& rPoints
, ::com::sun::star::drawing::PolyPolygonShape3D& rResult
- , sal_Int32 nGranularity );
+ , sal_uInt32 nGranularity );
static void CalculateBSplines(
const ::com::sun::star::drawing::PolyPolygonShape3D& rPoints
, ::com::sun::star::drawing::PolyPolygonShape3D& rResult
- , sal_Int32 nGranularity
- , sal_Int32 nSplineDepth );
+ , sal_uInt32 nGranularity
+ , sal_uInt32 nSplineDepth );
};
More information about the Libreoffice-commits
mailing list