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