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