1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 /* 3 * This file is part of the LibreOffice project. 4 * 5 * This Source Code Form is subject to the terms of the Mozilla Public 6 * License, v. 2.0. If a copy of the MPL was not distributed with this 7 * file, You can obtain one at http://mozilla.org/MPL/2.0/. 8 * 9 * This file incorporates work covered by the following license notice: 10 * 11 * Licensed to the Apache Software Foundation (ASF) under one or more 12 * contributor license agreements. See the NOTICE file distributed 13 * with this work for additional information regarding copyright 14 * ownership. The ASF licenses this file to you under the Apache 15 * License, Version 2.0 (the "License"); you may not use this file 16 * except in compliance with the License. You may obtain a copy of 17 * the License at http://www.apache.org/licenses/LICENSE-2.0 . 18 */ 19 20 #include "Splines.hxx" 21 #include <rtl/math.hxx> 22 #include <osl/diagnose.h> 23 #include <com/sun/star/drawing/PolyPolygonShape3D.hpp> 24 25 #include <vector> 26 #include <algorithm> 27 #include <memory> 28 29 namespace chart 30 { 31 using namespace ::com::sun::star; 32 33 namespace 34 { 35 36 typedef std::pair< double, double > tPointType; 37 typedef std::vector< tPointType > tPointVecType; 38 typedef tPointVecType::size_type lcl_tSizeType; 39 40 class lcl_SplineCalculation 41 { 42 public: 43 /** @descr creates an object that calculates cubic splines on construction 44 45 @param rSortedPoints the points for which splines shall be calculated, they need to be sorted in x values 46 @param fY1FirstDerivation the resulting spline should have the first 47 derivation equal to this value at the x-value of the first point 48 of rSortedPoints. If fY1FirstDerivation is set to infinity, a natural 49 spline is calculated. 50 @param fYnFirstDerivation the resulting spline should have the first 51 derivation equal to this value at the x-value of the last point 52 of rSortedPoints 53 */ 54 lcl_SplineCalculation( const tPointVecType & rSortedPoints, 55 double fY1FirstDerivation, 56 double fYnFirstDerivation ); 57 58 /** @descr creates an object that calculates cubic splines on construction 59 for the special case of periodic cubic spline 60 61 @param rSortedPoints the points for which splines shall be calculated, 62 they need to be sorted in x values. First and last y value must be equal 63 */ 64 explicit lcl_SplineCalculation( const tPointVecType & rSortedPoints); 65 66 /** @descr this function corresponds to the function splint in [1]. 67 68 [1] Numerical Recipes in C, 2nd edition 69 William H. Press, et al., 70 Section 3.3, page 116 71 */ 72 double GetInterpolatedValue( double x ); 73 74 private: 75 /// a copy of the points given in the CTOR 76 tPointVecType m_aPoints; 77 78 /// the result of the Calculate() method 79 std::vector< double > m_aSecDerivY; 80 81 double m_fYp1; 82 double m_fYpN; 83 84 // these values are cached for performance reasons 85 lcl_tSizeType m_nKLow; 86 lcl_tSizeType m_nKHigh; 87 double m_fLastInterpolatedValue; 88 89 /** @descr this function corresponds to the function spline in [1]. 90 91 [1] Numerical Recipes in C, 2nd edition 92 William H. Press, et al., 93 Section 3.3, page 115 94 */ 95 void Calculate(); 96 97 /** @descr this function corresponds to the algorithm 4.76 in [2] and 98 theorem 5.3.7 in [3] 99 100 [2] Engeln-Müllges, Gisela: Numerik-Algorithmen: Verfahren, Beispiele, Anwendungen 101 Springer, Berlin; Auflage: 9., überarb. und erw. A. (8. Dezember 2004) 102 Section 4.10.2, page 175 103 104 [3] Hanrath, Wilhelm: Mathematik III / Numerik, Vorlesungsskript zur 105 Veranstaltung im WS 2007/2008 106 Fachhochschule Aachen, 2009-09-19 107 Numerik_01.pdf, downloaded 2011-04-19 via 108 http://www.fh-aachen.de/index.php?id=11424&no_cache=1&file=5016&uid=44191 109 Section 5.3, page 129 110 */ 111 void CalculatePeriodic(); 112 }; 113 114 lcl_SplineCalculation::lcl_SplineCalculation( 115 const tPointVecType & rSortedPoints, 116 double fY1FirstDerivation, 117 double fYnFirstDerivation ) 118 : m_aPoints( rSortedPoints ), 119 m_fYp1( fY1FirstDerivation ), 120 m_fYpN( fYnFirstDerivation ), 121 m_nKLow( 0 ), 122 m_nKHigh( rSortedPoints.size() - 1 ), 123 m_fLastInterpolatedValue(0.0) 124 { 125 ::rtl::math::setInf( &m_fLastInterpolatedValue, false ); 126 Calculate(); 127 } 128 129 lcl_SplineCalculation::lcl_SplineCalculation( 130 const tPointVecType & rSortedPoints) 131 : m_aPoints( rSortedPoints ), 132 m_fYp1( 0.0 ), /*dummy*/ 133 m_fYpN( 0.0 ), /*dummy*/ 134 m_nKLow( 0 ), 135 m_nKHigh( rSortedPoints.size() - 1 ), 136 m_fLastInterpolatedValue(0.0) 137 { 138 ::rtl::math::setInf( &m_fLastInterpolatedValue, false ); 139 CalculatePeriodic(); 140 } 141 142 void lcl_SplineCalculation::Calculate() 143 { 144 if( m_aPoints.size() <= 1 ) 145 return; 146 147 // n is the last valid index to m_aPoints 148 const lcl_tSizeType n = m_aPoints.size() - 1; 149 std::vector< double > u( n ); 150 m_aSecDerivY.resize( n + 1, 0.0 ); 151 152 if( std::isinf( m_fYp1 ) ) 153 { 154 // natural spline 155 m_aSecDerivY[ 0 ] = 0.0; 156 u[ 0 ] = 0.0; 157 } 158 else 159 { 160 m_aSecDerivY[ 0 ] = -0.5; 161 double xDiff = m_aPoints[ 1 ].first - m_aPoints[ 0 ].first; 162 u[ 0 ] = ( 3.0 / xDiff ) * 163 ((( m_aPoints[ 1 ].second - m_aPoints[ 0 ].second ) / xDiff ) - m_fYp1 ); 164 } 165 166 for( lcl_tSizeType i = 1; i < n; ++i ) 167 { 168 tPointType 169 p_i = m_aPoints[ i ], 170 p_im1 = m_aPoints[ i - 1 ], 171 p_ip1 = m_aPoints[ i + 1 ]; 172 173 double sig = ( p_i.first - p_im1.first ) / 174 ( p_ip1.first - p_im1.first ); 175 double p = sig * m_aSecDerivY[ i - 1 ] + 2.0; 176 177 m_aSecDerivY[ i ] = ( sig - 1.0 ) / p; 178 u[ i ] = 179 ( ( p_ip1.second - p_i.second ) / 180 ( p_ip1.first - p_i.first ) ) - 181 ( ( p_i.second - p_im1.second ) / 182 ( p_i.first - p_im1.first ) ); 183 u[ i ] = 184 ( 6.0 * u[ i ] / ( p_ip1.first - p_im1.first ) 185 - sig * u[ i - 1 ] ) / p; 186 } 187 188 // initialize to values for natural splines (used for m_fYpN equal to 189 // infinity) 190 double qn = 0.0; 191 double un = 0.0; 192 193 if( ! std::isinf( m_fYpN ) ) 194 { 195 qn = 0.5; 196 double xDiff = m_aPoints[ n ].first - m_aPoints[ n - 1 ].first; 197 un = ( 3.0 / xDiff ) * 198 ( m_fYpN - ( m_aPoints[ n ].second - m_aPoints[ n - 1 ].second ) / xDiff ); 199 } 200 201 m_aSecDerivY[ n ] = ( un - qn * u[ n - 1 ] ) * ( qn * m_aSecDerivY[ n - 1 ] + 1.0 ); 202 203 // note: the algorithm in [1] iterates from n-1 to 0, but as size_type 204 // may be (usually is) an unsigned type, we can not write k >= 0, as this 205 // is always true. 206 for( lcl_tSizeType k = n; k > 0; --k ) 207 { 208 m_aSecDerivY[ k - 1 ] = (m_aSecDerivY[ k - 1 ] * m_aSecDerivY[ k ] ) + u[ k - 1 ]; 209 } 210 } 211 212 void lcl_SplineCalculation::CalculatePeriodic() 213 { 214 if( m_aPoints.size() <= 1 ) 215 return; 216 217 // n is the last valid index to m_aPoints 218 const lcl_tSizeType n = m_aPoints.size() - 1; 219 220 // u is used for vector f in A*c=f in [3], vector a in Ax=a in [2], 221 // vector z in Rtranspose z = a and Dr=z in [2] 222 std::vector< double > u( n + 1, 0.0 ); 223 224 // used for vector c in A*c=f and vector x in Ax=a in [2] 225 m_aSecDerivY.resize( n + 1, 0.0 ); 226 227 // diagonal of matrix A, used index 1 to n 228 std::vector< double > Adiag( n + 1, 0.0 ); 229 230 // secondary diagonal of matrix A with index 1 to n-1 and upper right element in A[n] 231 std::vector< double > Aupper( n + 1, 0.0 ); 232 233 // diagonal of matrix D in A=(R transpose)*D*R in [2], used index 1 to n 234 std::vector< double > Ddiag( n+1, 0.0 ); 235 236 // right column of matrix R, used index 1 to n-2 237 std::vector< double > Rright( n-1, 0.0 ); 238 239 // secondary diagonal of matrix R, used index 1 to n-1 240 std::vector< double > Rupper( n, 0.0 ); 241 242 if (n<4) 243 { 244 if (n==3) 245 { // special handling of three polynomials, that are four points 246 double xDiff0 = m_aPoints[ 1 ].first - m_aPoints[ 0 ].first ; 247 double xDiff1 = m_aPoints[ 2 ].first - m_aPoints[ 1 ].first ; 248 double xDiff2 = m_aPoints[ 3 ].first - m_aPoints[ 2 ].first ; 249 double xDiff2p1 = xDiff2 + xDiff1; 250 double xDiff0p2 = xDiff0 + xDiff2; 251 double xDiff1p0 = xDiff1 + xDiff0; 252 double fFactor = 1.5 / (xDiff0*xDiff1 + xDiff1*xDiff2 + xDiff2*xDiff0); 253 double yDiff0 = (m_aPoints[ 1 ].second - m_aPoints[ 0 ].second) / xDiff0; 254 double yDiff1 = (m_aPoints[ 2 ].second - m_aPoints[ 1 ].second) / xDiff1; 255 double yDiff2 = (m_aPoints[ 0 ].second - m_aPoints[ 2 ].second) / xDiff2; 256 m_aSecDerivY[ 1 ] = fFactor * (yDiff1*xDiff2p1 - yDiff0*xDiff0p2); 257 m_aSecDerivY[ 2 ] = fFactor * (yDiff2*xDiff0p2 - yDiff1*xDiff1p0); 258 m_aSecDerivY[ 3 ] = fFactor * (yDiff0*xDiff1p0 - yDiff2*xDiff2p1); 259 m_aSecDerivY[ 0 ] = m_aSecDerivY[ 3 ]; 260 } 261 else if (n==2) 262 { 263 // special handling of two polynomials, that are three points 264 double xDiff0 = m_aPoints[ 1 ].first - m_aPoints[ 0 ].first; 265 double xDiff1 = m_aPoints[ 2 ].first - m_aPoints[ 1 ].first; 266 double fHelp = 3.0 * (m_aPoints[ 0 ].second - m_aPoints[ 1 ].second) / (xDiff0*xDiff1); 267 m_aSecDerivY[ 1 ] = fHelp ; 268 m_aSecDerivY[ 2 ] = -fHelp ; 269 m_aSecDerivY[ 0 ] = m_aSecDerivY[ 2 ] ; 270 } 271 else 272 { 273 // should be handled with natural spline, periodic not possible. 274 } 275 } 276 else 277 { 278 double xDiff_i =1.0; // values are dummy; 279 double xDiff_im1 =1.0; 280 double yDiff_i = 1.0; 281 double yDiff_im1 = 1.0; 282 // fill matrix A and fill right side vector u 283 for( lcl_tSizeType i=1; i<n; ++i ) 284 { 285 xDiff_im1 = m_aPoints[ i ].first - m_aPoints[ i-1 ].first; 286 xDiff_i = m_aPoints[ i+1 ].first - m_aPoints[ i ].first; 287 yDiff_im1 = (m_aPoints[ i ].second - m_aPoints[ i-1 ].second) / xDiff_im1; 288 yDiff_i = (m_aPoints[ i+1 ].second - m_aPoints[ i ].second) / xDiff_i; 289 Adiag[ i ] = 2 * (xDiff_im1 + xDiff_i); 290 Aupper[ i ] = xDiff_i; 291 u [ i ] = 3 * (yDiff_i - yDiff_im1); 292 } 293 xDiff_im1 = m_aPoints[ n ].first - m_aPoints[ n-1 ].first; 294 xDiff_i = m_aPoints[ 1 ].first - m_aPoints[ 0 ].first; 295 yDiff_im1 = (m_aPoints[ n ].second - m_aPoints[ n-1 ].second) / xDiff_im1; 296 yDiff_i = (m_aPoints[ 1 ].second - m_aPoints[ 0 ].second) / xDiff_i; 297 Adiag[ n ] = 2 * (xDiff_im1 + xDiff_i); 298 Aupper[ n ] = xDiff_i; 299 u [ n ] = 3 * (yDiff_i - yDiff_im1); 300 301 // decomposite A=(R transpose)*D*R 302 Ddiag[1] = Adiag[1]; 303 Rupper[1] = Aupper[1] / Ddiag[1]; 304 Rright[1] = Aupper[n] / Ddiag[1]; 305 for( lcl_tSizeType i=2; i<=n-2; ++i ) 306 { 307 Ddiag[i] = Adiag[i] - Aupper[ i-1 ] * Rupper[ i-1 ]; 308 Rupper[ i ] = Aupper[ i ] / Ddiag[ i ]; 309 Rright[ i ] = - Rright[ i-1 ] * Aupper[ i-1 ] / Ddiag[ i ]; 310 } 311 Ddiag[ n-1 ] = Adiag[ n-1 ] - Aupper[ n-2 ] * Rupper[ n-2 ]; 312 Rupper[ n-1 ] = ( Aupper[ n-1 ] - Aupper[ n-2 ] * Rright[ n-2] ) / Ddiag[ n-1 ]; 313 double fSum = 0.0; 314 for ( lcl_tSizeType i=1; i<=n-2; ++i ) 315 { 316 fSum += Ddiag[ i ] * Rright[ i ] * Rright[ i ]; 317 } 318 Ddiag[ n ] = Adiag[ n ] - fSum - Ddiag[ n-1 ] * Rupper[ n-1 ] * Rupper[ n-1 ]; // bug in [2]! 319 320 // solve forward (R transpose)*z=u, overwrite u with z 321 for ( lcl_tSizeType i=2; i<=n-1; ++i ) 322 { 323 u[ i ] -= u[ i-1 ]* Rupper[ i-1 ]; 324 } 325 fSum = 0.0; 326 for ( lcl_tSizeType i=1; i<=n-2; ++i ) 327 { 328 fSum += Rright[ i ] * u[ i ]; 329 } 330 u[ n ] = u[ n ] - fSum - Rupper[ n - 1] * u[ n-1 ]; 331 332 // solve forward D*r=z, z is in u, overwrite u with r 333 for ( lcl_tSizeType i=1; i<=n; ++i ) 334 { 335 u[ i ] = u[i] / Ddiag[ i ]; 336 } 337 338 // solve backward R*x= r, r is in u 339 m_aSecDerivY[ n ] = u[ n ]; 340 m_aSecDerivY[ n-1 ] = u[ n-1 ] - Rupper[ n-1 ] * m_aSecDerivY[ n ]; 341 for ( lcl_tSizeType i=n-2; i>=1; --i) 342 { 343 m_aSecDerivY[ i ] = u[ i ] - Rupper[ i ] * m_aSecDerivY[ i+1 ] - Rright[ i ] * m_aSecDerivY[ n ]; 344 } 345 // periodic 346 m_aSecDerivY[ 0 ] = m_aSecDerivY[ n ]; 347 } 348 349 // adapt m_aSecDerivY for usage in GetInterpolatedValue() 350 for( lcl_tSizeType i = 0; i <= n ; ++i ) 351 { 352 m_aSecDerivY[ i ] *= 2.0; 353 } 354 355 } 356 357 double lcl_SplineCalculation::GetInterpolatedValue( double x ) 358 { 359 OSL_PRECOND( ( m_aPoints[ 0 ].first <= x ) && 360 ( x <= m_aPoints[ m_aPoints.size() - 1 ].first ), 361 "Trying to extrapolate" ); 362 363 const lcl_tSizeType n = m_aPoints.size() - 1; 364 if( x < m_fLastInterpolatedValue ) 365 { 366 m_nKLow = 0; 367 m_nKHigh = n; 368 369 // calculate m_nKLow and m_nKHigh 370 // first initialization is done in CTOR 371 while( m_nKHigh - m_nKLow > 1 ) 372 { 373 lcl_tSizeType k = ( m_nKHigh + m_nKLow ) / 2; 374 if( m_aPoints[ k ].first > x ) 375 m_nKHigh = k; 376 else 377 m_nKLow = k; 378 } 379 } 380 else 381 { 382 while( ( m_nKHigh <= n ) && 383 ( m_aPoints[ m_nKHigh ].first < x ) ) 384 { 385 ++m_nKHigh; 386 ++m_nKLow; 387 } 388 OSL_ENSURE( m_nKHigh <= n, "Out of Bounds" ); 389 } 390 m_fLastInterpolatedValue = x; 391 392 double h = m_aPoints[ m_nKHigh ].first - m_aPoints[ m_nKLow ].first; 393 OSL_ENSURE( h != 0, "Bad input to GetInterpolatedValue()" ); 394 395 double a = ( m_aPoints[ m_nKHigh ].first - x ) / h; 396 double b = ( x - m_aPoints[ m_nKLow ].first ) / h; 397 398 return ( a * m_aPoints[ m_nKLow ].second + 399 b * m_aPoints[ m_nKHigh ].second + 400 (( a*a*a - a ) * m_aSecDerivY[ m_nKLow ] + 401 ( b*b*b - b ) * m_aSecDerivY[ m_nKHigh ] ) * 402 ( h*h ) / 6.0 ); 403 } 404 405 // helper methods for B-spline 406 407 // Create parameter t_0 to t_n using the centripetal method with a power of 0.5 408 bool createParameterT(const tPointVecType& rUniquePoints, double* t) 409 { // precondition: no adjacent identical points 410 // postcondition: 0 = t_0 < t_1 < ... < t_n = 1 411 bool bIsSuccessful = true; 412 const lcl_tSizeType n = rUniquePoints.size() - 1; 413 t[0]=0.0; 414 double dx = 0.0; 415 double dy = 0.0; 416 double fDiffMax = 1.0; //dummy values 417 double fDenominator = 0.0; // initialized for summing up 418 for (lcl_tSizeType i=1; i<=n ; ++i) 419 { // 4th root(dx^2+dy^2) 420 dx = rUniquePoints[i].first - rUniquePoints[i-1].first; 421 dy = rUniquePoints[i].second - rUniquePoints[i-1].second; 422 // scaling to avoid underflow or overflow 423 fDiffMax = std::max(fabs(dx), fabs(dy)); 424 if (fDiffMax == 0.0) 425 { 426 bIsSuccessful = false; 427 break; 428 } 429 else 430 { 431 dx /= fDiffMax; 432 dy /= fDiffMax; 433 fDenominator += sqrt(sqrt(dx * dx + dy * dy)) * sqrt(fDiffMax); 434 } 435 } 436 if (fDenominator == 0.0) 437 { 438 bIsSuccessful = false; 439 } 440 if (bIsSuccessful) 441 { 442 for (lcl_tSizeType j=1; j<=n ; ++j) 443 { 444 double fNumerator = 0.0; 445 for (lcl_tSizeType i=1; i<=j ; ++i) 446 { 447 dx = rUniquePoints[i].first - rUniquePoints[i-1].first; 448 dy = rUniquePoints[i].second - rUniquePoints[i-1].second; 449 fDiffMax = std::max(fabs(dx), fabs(dy)); 450 // same as above, so should not be zero 451 dx /= fDiffMax; 452 dy /= fDiffMax; 453 fNumerator += sqrt(sqrt(dx * dx + dy * dy)) * sqrt(fDiffMax); 454 } 455 t[j] = fNumerator / fDenominator; 456 457 } 458 // postcondition check 459 t[n] = 1.0; 460 double fPrevious = 0.0; 461 for (lcl_tSizeType i=1; i <= n && bIsSuccessful ; ++i) 462 { 463 if (fPrevious >= t[i]) 464 { 465 bIsSuccessful = false; 466 } 467 else 468 { 469 fPrevious = t[i]; 470 } 471 } 472 } 473 return bIsSuccessful; 474 } 475 476 void createKnotVector(const lcl_tSizeType n, const sal_uInt32 p, const double* t, double* u) 477 { // precondition: 0 = t_0 < t_1 < ... < t_n = 1 478 for (lcl_tSizeType j = 0; j <= p; ++j) 479 { 480 u[j] = 0.0; 481 } 482 for (lcl_tSizeType j = 1; j <= n-p; ++j ) 483 { 484 double fSum = 0.0; 485 for (lcl_tSizeType i = j; i <= j+p-1; ++i) 486 { 487 fSum += t[i]; 488 } 489 assert(p != 0); 490 u[j+p] = fSum / p ; 491 } 492 for (lcl_tSizeType j = n+1; j <= n+1+p; ++j) 493 { 494 u[j] = 1.0; 495 } 496 } 497 498 void applyNtoParameterT(const lcl_tSizeType i,const double tk,const sal_uInt32 p,const double* u, double* rowN) 499 { 500 // get N_p(t_k) recursively, only N_(i-p) till N_(i) are relevant, all other N_# are zero 501 502 // initialize with indicator function degree 0 503 rowN[p] = 1.0; // all others are zero 504 505 // calculate up to degree p 506 for (sal_uInt32 s = 1; s <= p; ++s) 507 { 508 // first element 509 double fLeftFactor = 0.0; 510 double fRightFactor = ( u[i+1] - tk ) / ( u[i+1]- u[i-s+1] ); 511 // i-s "true index" - (i-p)"shift" = p-s 512 rowN[p-s] = fRightFactor * rowN[p-s+1]; 513 514 // middle elements 515 for (sal_uInt32 j = s-1; j>=1 ; --j) 516 { 517 fLeftFactor = ( tk - u[i-j] ) / ( u[i-j+s] - u[i-j] ) ; 518 fRightFactor = ( u[i-j+s+1] - tk ) / ( u[i-j+s+1] - u[i-j+1] ); 519 // i-j "true index" - (i-p)"shift" = p-j 520 rowN[p-j] = fLeftFactor * rowN[p-j] + fRightFactor * rowN[p-j+1]; 521 } 522 523 // last element 524 fLeftFactor = ( tk - u[i] ) / ( u[i+s] - u[i] ); 525 // i "true index" - (i-p)"shift" = p 526 rowN[p] = fLeftFactor * rowN[p]; 527 } 528 } 529 530 } // anonymous namespace 531 532 // Calculates uniform parametric splines with subinterval length 1, 533 // according ODF1.2 part 1, chapter 'chart interpolation'. 534 void SplineCalculater::CalculateCubicSplines( 535 const drawing::PolyPolygonShape3D& rInput 536 , drawing::PolyPolygonShape3D& rResult 537 , sal_uInt32 nGranularity ) 538 { 539 OSL_PRECOND( nGranularity > 0, "Granularity is invalid" ); 540 541 sal_uInt32 nOuterCount = rInput.SequenceX.getLength(); 542 543 rResult.SequenceX.realloc(nOuterCount); 544 rResult.SequenceY.realloc(nOuterCount); 545 rResult.SequenceZ.realloc(nOuterCount); 546 547 if( !nOuterCount ) 548 return; 549 550 for( sal_uInt32 nOuter = 0; nOuter < nOuterCount; ++nOuter ) 551 { 552 if( rInput.SequenceX[nOuter].getLength() <= 1 ) 553 continue; //we need at least two points 554 555 sal_uInt32 nMaxIndexPoints = rInput.SequenceX[nOuter].getLength()-1; // is >=1 556 const double* pOldX = rInput.SequenceX[nOuter].getConstArray(); 557 const double* pOldY = rInput.SequenceY[nOuter].getConstArray(); 558 const double* pOldZ = rInput.SequenceZ[nOuter].getConstArray(); 559 560 std::vector < double > aParameter(nMaxIndexPoints+1); 561 aParameter[0]=0.0; 562 for( sal_uInt32 nIndex=1; nIndex<=nMaxIndexPoints; nIndex++ ) 563 { 564 aParameter[nIndex]=aParameter[nIndex-1]+1; 565 } 566 567 // Split the calculation to X, Y and Z coordinate 568 tPointVecType aInputX; 569 aInputX.resize(nMaxIndexPoints+1); 570 tPointVecType aInputY; 571 aInputY.resize(nMaxIndexPoints+1); 572 tPointVecType aInputZ; 573 aInputZ.resize(nMaxIndexPoints+1); 574 for (sal_uInt32 nN=0;nN<=nMaxIndexPoints; nN++ ) 575 { 576 aInputX[ nN ].first=aParameter[nN]; 577 aInputX[ nN ].second=pOldX[ nN ]; 578 aInputY[ nN ].first=aParameter[nN]; 579 aInputY[ nN ].second=pOldY[ nN ]; 580 aInputZ[ nN ].first=aParameter[nN]; 581 aInputZ[ nN ].second=pOldZ[ nN ]; 582 } 583 584 // generate a spline for each coordinate. It holds the complete 585 // information to calculate each point of the curve 586 std::unique_ptr<lcl_SplineCalculation> aSplineX; 587 std::unique_ptr<lcl_SplineCalculation> aSplineY; 588 // lcl_SplineCalculation* aSplineZ; the z-coordinates of all points in 589 // a data series are equal. No spline calculation needed, but copy 590 // coordinate to output 591 592 if( pOldX[ 0 ] == pOldX[nMaxIndexPoints] && 593 pOldY[ 0 ] == pOldY[nMaxIndexPoints] && 594 pOldZ[ 0 ] == pOldZ[nMaxIndexPoints] && 595 nMaxIndexPoints >=2 ) 596 { // periodic spline 597 aSplineX.reset(new lcl_SplineCalculation( aInputX)); 598 aSplineY.reset(new lcl_SplineCalculation( aInputY)); 599 // aSplineZ = new lcl_SplineCalculation( aInputZ) ; 600 } 601 else // generate the kind "natural spline" 602 { 603 double fInfty; 604 ::rtl::math::setInf( &fInfty, false ); 605 double fXDerivation = fInfty; 606 double fYDerivation = fInfty; 607 aSplineX.reset(new lcl_SplineCalculation( aInputX, fXDerivation, fXDerivation )); 608 aSplineY.reset(new lcl_SplineCalculation( aInputY, fYDerivation, fYDerivation )); 609 } 610 611 // fill result polygon with calculated values 612 rResult.SequenceX[nOuter].realloc( nMaxIndexPoints*nGranularity + 1); 613 rResult.SequenceY[nOuter].realloc( nMaxIndexPoints*nGranularity + 1); 614 rResult.SequenceZ[nOuter].realloc( nMaxIndexPoints*nGranularity + 1); 615 616 double* pNewX = rResult.SequenceX[nOuter].getArray(); 617 double* pNewY = rResult.SequenceY[nOuter].getArray(); 618 double* pNewZ = rResult.SequenceZ[nOuter].getArray(); 619 620 sal_uInt32 nNewPointIndex = 0; // Index in result points 621 622 for( sal_uInt32 ni = 0; ni < nMaxIndexPoints; ni++ ) 623 { 624 // given point is surely a curve point 625 pNewX[nNewPointIndex] = pOldX[ni]; 626 pNewY[nNewPointIndex] = pOldY[ni]; 627 pNewZ[nNewPointIndex] = pOldZ[ni]; 628 nNewPointIndex++; 629 630 // calculate intermediate points 631 double fInc = ( aParameter[ ni+1 ] - aParameter[ni] ) / static_cast< double >( nGranularity ); 632 for(sal_uInt32 nj = 1; nj < nGranularity; nj++) 633 { 634 double fParam = aParameter[ni] + ( fInc * static_cast< double >( nj ) ); 635 636 pNewX[nNewPointIndex]=aSplineX->GetInterpolatedValue( fParam ); 637 pNewY[nNewPointIndex]=aSplineY->GetInterpolatedValue( fParam ); 638 // pNewZ[nNewPointIndex]=aSplineZ->GetInterpolatedValue( fParam ); 639 pNewZ[nNewPointIndex] = pOldZ[ni]; 640 nNewPointIndex++; 641 } 642 } 643 // add last point 644 pNewX[nNewPointIndex] = pOldX[nMaxIndexPoints]; 645 pNewY[nNewPointIndex] = pOldY[nMaxIndexPoints]; 646 pNewZ[nNewPointIndex] = pOldZ[nMaxIndexPoints]; 647 } 648 } 649 650 // The implementation follows closely ODF1.2 spec, chapter chart:interpolation 651 // using the same names as in spec as far as possible, without prefix. 652 // More details can be found on 653 // Dr. C.-K. Shene: CS3621 Introduction to Computing with Geometry Notes 654 // Unit 9: Interpolation and Approximation/Curve Global Interpolation 655 // Department of Computer Science, Michigan Technological University 656 // http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/ 657 // [last called 2011-05-20] 658 void SplineCalculater::CalculateBSplines( 659 const css::drawing::PolyPolygonShape3D& rInput 660 , css::drawing::PolyPolygonShape3D& rResult 661 , sal_uInt32 nResolution 662 , sal_uInt32 nDegree ) 663 { 664 // nResolution is ODF1.2 file format attribute chart:spline-resolution and 665 // ODF1.2 spec variable k. Causion, k is used as index in the spec in addition. 666 // nDegree is ODF1.2 file format attribute chart:spline-order and 667 // ODF1.2 spec variable p 668 OSL_ASSERT( nResolution > 1 ); 669 OSL_ASSERT( nDegree >= 1 ); 670 671 // limit the b-spline degree at 15 to prevent insanely large sets of points 672 sal_uInt32 p = std::min<sal_uInt32>(nDegree, 15); 673 674 sal_Int32 nOuterCount = rInput.SequenceX.getLength(); 675 676 rResult.SequenceX.realloc(nOuterCount); 677 rResult.SequenceY.realloc(nOuterCount); 678 rResult.SequenceZ.realloc(nOuterCount); 679 680 if( !nOuterCount ) 681 return; // no input 682 683 for( sal_Int32 nOuter = 0; nOuter < nOuterCount; ++nOuter ) 684 { 685 if( rInput.SequenceX[nOuter].getLength() <= 1 ) 686 continue; // need at least 2 points, next piece of the series 687 688 // Copy input to vector of points and remove adjacent double points. The 689 // Z-coordinate is equal for all points in a series and holds the depth 690 // in 3D mode, simple copying is enough. 691 lcl_tSizeType nMaxIndexPoints = rInput.SequenceX[nOuter].getLength()-1; // is >=1 692 const double* pOldX = rInput.SequenceX[nOuter].getConstArray(); 693 const double* pOldY = rInput.SequenceY[nOuter].getConstArray(); 694 const double* pOldZ = rInput.SequenceZ[nOuter].getConstArray(); 695 double fZCoordinate = pOldZ[0]; 696 tPointVecType aPointsIn; 697 aPointsIn.resize(nMaxIndexPoints+1); 698 for (lcl_tSizeType i = 0; i <= nMaxIndexPoints; ++i ) 699 { 700 aPointsIn[ i ].first = pOldX[i]; 701 aPointsIn[ i ].second = pOldY[i]; 702 } 703 aPointsIn.erase( std::unique( aPointsIn.begin(), aPointsIn.end()), 704 aPointsIn.end() ); 705 706 // n is the last valid index to the reduced aPointsIn 707 // There are n+1 valid data points. 708 const lcl_tSizeType n = aPointsIn.size() - 1; 709 if (n < 1 || p > n) 710 continue; // need at least 2 points, degree p needs at least n+1 points 711 // next piece of series 712 713 std::unique_ptr<double[]> t(new double [n+1]); 714 if (!createParameterT(aPointsIn, t.get())) 715 { 716 continue; // next piece of series 717 } 718 719 lcl_tSizeType m = n + p + 1; 720 std::unique_ptr<double[]> u(new double [m+1]); 721 createKnotVector(n, p, t.get(), u.get()); 722 723 // The matrix N contains the B-spline basis functions applied to parameters. 724 // In each row only p+1 adjacent elements are non-zero. The starting 725 // column in a higher row is equal or greater than in the lower row. 726 // To store this matrix the non-zero elements are shifted to column 0 727 // and the amount of shifting is remembered in an array. 728 std::unique_ptr<double*[]> aMatN(new double*[n+1]); 729 for (lcl_tSizeType row = 0; row <=n; ++row) 730 { 731 aMatN[row] = new double[p+1]; 732 for (sal_uInt32 col = 0; col <= p; ++col) 733 aMatN[row][col] = 0.0; 734 } 735 std::unique_ptr<lcl_tSizeType[]> aShift(new lcl_tSizeType[n+1]); 736 aMatN[0][0] = 1.0; //all others are zero 737 aShift[0] = 0; 738 aMatN[n][0] = 1.0; 739 aShift[n] = n; 740 for (lcl_tSizeType k = 1; k<=n-1; ++k) 741 { // all basis functions are applied to t_k, 742 // results are elements in row k in matrix N 743 744 // find the one interval with u_i <= t_k < u_(i+1) 745 // remember u_0 = ... = u_p = 0.0 and u_(m-p) = ... u_m = 1.0 and 0<t_k<1 746 lcl_tSizeType i = p; 747 while (u[i] > t[k] || t[k] >= u[i+1]) 748 { 749 ++i; 750 } 751 752 // index in reduced matrix aMatN = (index in full matrix N) - (i-p) 753 aShift[k] = i - p; 754 755 applyNtoParameterT(i, t[k], p, u.get(), aMatN[k]); 756 } // next row k 757 758 // Get matrix C of control points from the matrix equation aMatN * C = aPointsIn 759 // aPointsIn is overwritten with C. 760 // Gaussian elimination is possible without pivoting, see reference 761 lcl_tSizeType r = 0; // true row index 762 lcl_tSizeType c = 0; // true column index 763 double fDivisor = 1.0; // used for diagonal element 764 double fEliminate = 1.0; // used for the element, that will become zero 765 double fHelp; 766 tPointType aHelp; 767 lcl_tSizeType nHelp; // used in triangle change 768 bool bIsSuccessful = true; 769 for (c = 0 ; c <= n && bIsSuccessful; ++c) 770 { 771 // search for first non-zero downwards 772 r = c; 773 while ( r < n && aMatN[r][c-aShift[r]] == 0 ) 774 { 775 ++r; 776 } 777 if (aMatN[r][c-aShift[r]] == 0.0) 778 { 779 // Matrix N is singular, although this is mathematically impossible 780 bIsSuccessful = false; 781 } 782 else 783 { 784 // exchange total row r with total row c if necessary 785 if (r != c) 786 { 787 for ( sal_uInt32 i = 0; i <= p ; ++i) 788 { 789 fHelp = aMatN[r][i]; 790 aMatN[r][i] = aMatN[c][i]; 791 aMatN[c][i] = fHelp; 792 } 793 aHelp = aPointsIn[r]; 794 aPointsIn[r] = aPointsIn[c]; 795 aPointsIn[c] = aHelp; 796 nHelp = aShift[r]; 797 aShift[r] = aShift[c]; 798 aShift[c] = nHelp; 799 } 800 801 // divide row c, so that element(c,c) becomes 1 802 fDivisor = aMatN[c][c-aShift[c]]; // not zero, see above 803 for (sal_uInt32 i = 0; i <= p; ++i) 804 { 805 aMatN[c][i] /= fDivisor; 806 } 807 aPointsIn[c].first /= fDivisor; 808 aPointsIn[c].second /= fDivisor; 809 810 // eliminate forward, examine row c+1 to n-1 (worst case) 811 // stop if first non-zero element in row has an higher column as c 812 // look at nShift for that, elements in nShift are equal or increasing 813 for ( r = c+1; r < n && aShift[r]<=c ; ++r) 814 { 815 fEliminate = aMatN[r][0]; 816 if (fEliminate != 0.0) // else accidentally zero, nothing to do 817 { 818 for (sal_uInt32 i = 1; i <= p; ++i) 819 { 820 aMatN[r][i-1] = aMatN[r][i] - fEliminate * aMatN[c][i]; 821 } 822 aMatN[r][p]=0; 823 aPointsIn[r].first -= fEliminate * aPointsIn[c].first; 824 aPointsIn[r].second -= fEliminate * aPointsIn[c].second; 825 ++aShift[r]; 826 } 827 } 828 } 829 }// upper triangle form is reached 830 if( bIsSuccessful) 831 { 832 // eliminate backwards, begin with last column 833 for (lcl_tSizeType cc = n; cc >= 1; --cc ) 834 { 835 // In row cc the diagonal element(cc,cc) == 1 and all elements left from 836 // diagonal are zero and do not influence other rows. 837 // Full matrix N has semibandwidth < p, therefore element(r,c) is 838 // zero, if abs(r-cc)>=p. abs(r-cc)=cc-r, because r<cc. 839 r = cc - 1; 840 while ( r !=0 && cc-r < p ) 841 { 842 fEliminate = aMatN[r][ cc - aShift[r] ]; 843 if ( fEliminate != 0.0) // else element is accidentally zero, no action needed 844 { 845 // row r -= fEliminate * row cc only relevant for right side 846 aMatN[r][cc - aShift[r]] = 0.0; 847 aPointsIn[r].first -= fEliminate * aPointsIn[cc].first; 848 aPointsIn[r].second -= fEliminate * aPointsIn[cc].second; 849 } 850 --r; 851 } 852 } 853 // aPointsIn contains the control points now. 854 855 // calculate the intermediate points according given resolution 856 // using deBoor-Cox algorithm 857 lcl_tSizeType nNewSize = nResolution * n + 1; 858 rResult.SequenceX[nOuter].realloc(nNewSize); 859 rResult.SequenceY[nOuter].realloc(nNewSize); 860 rResult.SequenceZ[nOuter].realloc(nNewSize); 861 double* pNewX = rResult.SequenceX[nOuter].getArray(); 862 double* pNewY = rResult.SequenceY[nOuter].getArray(); 863 double* pNewZ = rResult.SequenceZ[nOuter].getArray(); 864 pNewX[0] = aPointsIn[0].first; 865 pNewY[0] = aPointsIn[0].second; 866 pNewZ[0] = fZCoordinate; // Precondition: z-coordinates of all points of a series are equal 867 pNewX[nNewSize -1 ] = aPointsIn[n].first; 868 pNewY[nNewSize -1 ] = aPointsIn[n].second; 869 pNewZ[nNewSize -1 ] = fZCoordinate; 870 std::unique_ptr<double[]> aP(new double[m+1]); 871 lcl_tSizeType nLow = 0; 872 for ( lcl_tSizeType nTIndex = 0; nTIndex <= n-1; ++nTIndex) 873 { 874 for (sal_uInt32 nResolutionStep = 1; 875 nResolutionStep <= nResolution && ( nTIndex != n-1 || nResolutionStep != nResolution); 876 ++nResolutionStep) 877 { 878 lcl_tSizeType nNewIndex = nTIndex * nResolution + nResolutionStep; 879 double ux = t[nTIndex] + nResolutionStep * ( t[nTIndex+1] - t[nTIndex]) /nResolution; 880 881 // get index nLow, so that u[nLow]<= ux < u[nLow +1] 882 // continue from previous nLow 883 while ( u[nLow] <= ux) 884 { 885 ++nLow; 886 } 887 --nLow; 888 889 // x-coordinate 890 for (lcl_tSizeType i = nLow-p; i <= nLow; ++i) 891 { 892 aP[i] = aPointsIn[i].first; 893 } 894 for (sal_uInt32 lcl_Degree = 1; lcl_Degree <= p; ++lcl_Degree) 895 { 896 for (lcl_tSizeType i = nLow; i >= nLow + lcl_Degree - p; --i) 897 { 898 double fFactor = ( ux - u[i] ) / ( u[i+p+1-lcl_Degree] - u[i]); 899 aP[i] = (1 - fFactor)* aP[i-1] + fFactor * aP[i]; 900 } 901 } 902 pNewX[nNewIndex] = aP[nLow]; 903 904 // y-coordinate 905 for (lcl_tSizeType i = nLow - p; i <= nLow; ++i) 906 { 907 aP[i] = aPointsIn[i].second; 908 } 909 for (sal_uInt32 lcl_Degree = 1; lcl_Degree <= p; ++lcl_Degree) 910 { 911 for (lcl_tSizeType i = nLow; i >= nLow +lcl_Degree - p; --i) 912 { 913 double fFactor = ( ux - u[i] ) / ( u[i+p+1-lcl_Degree] - u[i]); 914 aP[i] = (1 - fFactor)* aP[i-1] + fFactor * aP[i]; 915 } 916 } 917 pNewY[nNewIndex] = aP[nLow]; 918 pNewZ[nNewIndex] = fZCoordinate; 919 } 920 } 921 } 922 for (lcl_tSizeType row = 0; row <=n; ++row) 923 { 924 delete[] aMatN[row]; 925 } 926 } // next piece of the series 927 } 928 929 } //namespace chart 930 931 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */ 932
