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