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 <math.h>
21 :
22 : #include <tools/poly.hxx>
23 : #include <boost/scoped_array.hpp>
24 :
25 : #include <sgvspln.hxx>
26 :
27 : extern "C" {
28 :
29 : /*.pn 277 */
30 : /*.hlAppendix: C - programs*/
31 : /*.hrConstants- and macrodefinitions*/
32 : /*.fe The include file u_const.h should be stored in the directory, */
33 : /*.fe where the compiler searches for include files. */
34 :
35 : /*----------------------- FILE u_const.h ---------------------------*/
36 :
37 : #define IEEE
38 :
39 : /* IEEE - standard for representation of floating-point numbers:
40 :
41 : 8 byte long floating point numbers with
42 :
43 : 53 bit mantissa ==> mantissa range: 2^52 different numbers
44 : with 0.1 <= number < 1.0,
45 : 1 sign-bit
46 : 11 bit exponent ==> exponent range: -1024...+1023
47 :
48 : The first line (#define IEEE) should be deleted if the machine
49 : or the compiler does not use floating-point numbers according
50 : to the IEEE standard. In which case also MAXEXPON, MINEXPON (see
51 : below) should be adapted.
52 : */
53 :
54 : #ifdef IEEE /*-------------- if IEEE norm --------------------*/
55 :
56 : #define MACH_EPS 2.220446049250313e-016 /* machine precision */
57 : /* IBM-AT: = 2^-52 */
58 : /* MACH_EPS is the smallest positive, by the machine representable
59 : number x, which fulfills the equation: 1.0 + x > 1.0 */
60 :
61 : #define MAXROOT 9.48075190810918e+153
62 :
63 : #else /*------------------ otherwise--------------------*/
64 :
65 : double exp (double);
66 : double atan (double);
67 : double pow (double,double);
68 : double sqrt (double);
69 :
70 : double masch() /* calculate MACH_EPS machine independence */
71 : {
72 : double eps = 1.0, x = 2.0, y = 1.0;
73 : while ( y < x )
74 : { eps *= 0.5;
75 : x = 1.0 + eps;
76 : }
77 : eps *= 2.0;
78 : return (eps);
79 : }
80 :
81 : short basis() /* calculate BASE machine independence */
82 : {
83 : double x = 1.0, one = 1.0, b = 1.0;
84 :
85 : while ( (x + one) - x == one ) x *= 2.0;
86 : while ( (x + b) == x ) b *= 2.0;
87 :
88 : return ( (short) ((x + b) - x) );
89 : }
90 :
91 : #define BASIS basis() /* base of number representation */
92 :
93 : /* If the machine (the compiler) does not use the IEEE-representation
94 : for floating-point numbers, the next 2 constants should be adapted.
95 : */
96 :
97 : #define MAXEXPON 1023.0 /* largest exponent */
98 : #define MINEXPON -1024.0 /* smallest exponent */
99 :
100 : #define MACH_EPS masch()
101 :
102 : #define POSMAX pow ((double) BASIS, MAXEXPON)
103 : #define MAXROOT sqrt(POSMAX)
104 :
105 : #endif /*-------------- END of ifdef --------------------*/
106 :
107 : /* defines for function macros: */
108 :
109 : #define abs(X) ((X) >= 0 ? (X) : -(X)) /* absolute number X */
110 : #define sign(X, Y) (Y < 0 ? -abs(X) : abs(X)) /* sign of Y times */
111 : /* abs(X) */
112 :
113 : /*-------------------- END of FILE u_const.h -----------------------*/
114 :
115 : /*.HL appendix: C - programs*/
116 : /*.HR Systems of equations for tridiagonal matrices*/
117 :
118 : /*.FE P 3.7 tridiagonal systems of equations */
119 :
120 : /*---------------------- MODULE tridiagonal -----------------------*/
121 :
122 0 : sal_uInt16 TriDiagGS(bool rep, sal_uInt16 n, double* lower,
123 : double* diag, double* upper, double* b)
124 : /*************************/
125 : /* Gaussian methods for */
126 : /* tridiagonal matrices */
127 : /*************************/
128 :
129 : /*====================================================================*/
130 : /* */
131 : /* trdiag determines solution x of the system of linear equations */
132 : /* A * x = b with tridiagonal n x n coefficient matrix A, which is */
133 : /* stored in 3 vectors lower, upper and diag as per below: */
134 : /* */
135 : /* ( diag[0] upper[0] 0 0 . . . 0 ) */
136 : /* ( lower[1] diag[1] upper[1] 0 . . . ) */
137 : /* ( 0 lower[2] diag[2] upper[2] 0 . ) */
138 : /* A = ( . 0 lower[3] . . . ) */
139 : /* ( . . . . . 0 ) */
140 : /* ( . . . . . ) */
141 : /* ( . . . upper[n-2] ) */
142 : /* ( 0 . . . 0 lower[n-1] diag[n-1] ) */
143 : /* */
144 : /*====================================================================*/
145 : /* */
146 : /* Usage: */
147 : /* ====== */
148 : /* Mainly for diagonal determinant triangular matrix, as they */
149 : /* occur in Spline-interpolations. */
150 : /* For diagonal dominant matrices always a left-upper row */
151 : /* reduction exists; for non diagonal dominant triangular */
152 : /* matrices we should pull forward the function band, as this */
153 : /* works with row pivot searches, which is numerical more stable.*/
154 : /* */
155 : /*====================================================================*/
156 : /* */
157 : /* Input parameters: */
158 : /* ================ */
159 : /* n dimension of the matrix ( > 1 ) sal_uInt16 n */
160 : /* */
161 : /* lower lower antidiagonal double lower[n] */
162 : /* diag main diagonal double diag[n] */
163 : /* upper upper antidiagonal double upper[n] */
164 : /* */
165 : /* for rep = true lower, diag and upper contain the */
166 : /* triangulation of the start matrix. */
167 : /* */
168 : /* b right side of equation double b[n] */
169 : /* rep = false first call bool rep */
170 : /* = true next call */
171 : /* for the same matrix, */
172 : /* but different b. */
173 : /* */
174 : /* Output parameters: */
175 : /* ================= */
176 : /* b solution vector of the system; double b[n] */
177 : /* the original right side is overwritten */
178 : /* */
179 : /* lower ) contain for rep = false the decomposition of the */
180 : /* diag ) matrix; the original values of the lower and */
181 : /* upper ) diagonals are overwritten */
182 : /* */
183 : /* The determinant of the matrix is for rep = false defined by */
184 : /* determinant A = diag[0] * ... * diag[n-1] */
185 : /* */
186 : /* Return value: */
187 : /* ============= */
188 : /* = 0 all ok */
189 : /* = 1 n < 2 chosen */
190 : /* = 2 triangular decomposition of matrix does not exist */
191 : /* */
192 : /*====================================================================*/
193 : /* */
194 : /* Functions used: */
195 : /* =============== */
196 : /* */
197 : /* From the C library: fabs() */
198 : /* */
199 : /*====================================================================*/
200 :
201 : /*.cp 5 */
202 : {
203 : sal_uInt16 i;
204 : short j;
205 :
206 : // double fabs(double);
207 :
208 0 : if ( n < 2 ) return(1); /* n at least 2 */
209 :
210 : /* if rep = false, */
211 : /* determine the */
212 : /* triangular */
213 : /* decomposition of */
214 0 : if (!rep) /* matrix and determinant*/
215 : {
216 0 : for (i = 1; i < n; i++)
217 0 : { if ( fabs(diag[i-1]) < MACH_EPS ) /* do not decompose */
218 0 : return(2); /* if one diag[i] = 0 */
219 0 : lower[i] /= diag[i-1];
220 0 : diag[i] -= lower[i] * upper[i-1];
221 : }
222 : }
223 :
224 0 : if ( fabs(diag[n-1]) < MACH_EPS ) return(2);
225 :
226 0 : for (i = 1; i < n; i++) /* forward elimination */
227 0 : b[i] -= lower[i] * b[i-1];
228 :
229 0 : b[n-1] /= diag[n-1]; /* reverse elimination */
230 0 : for (j = n-2; j >= 0; j--) {
231 0 : i=j;
232 0 : b[i] = ( b[i] - upper[i] * b[i+1] ) / diag[i];
233 : }
234 0 : return(0);
235 : }
236 :
237 : /*----------------------- END OF TRIDIAGONAL ------------------------*/
238 :
239 : /*.HL Appendix: C - Programs*/
240 : /*.HRSystems of equations with cyclic tridiagonal matrices*/
241 :
242 : /*.FE P 3.8 Systems with cyclic tridiagonal matrices */
243 :
244 : /*---------------- Module cyclic tridiagonal -----------------------*/
245 :
246 12 : sal_uInt16 ZyklTriDiagGS(bool rep, sal_uInt16 n, double* lower, double* diag,
247 : double* upper, double* lowrow, double* ricol, double* b)
248 : /******************************/
249 : /* Systems with cyclic */
250 : /* tridiagonal matrices */
251 : /******************************/
252 :
253 : /*====================================================================*/
254 : /* */
255 : /* tzdiag determines the solution x of the linear equation system */
256 : /* A * x = b with cyclic tridiagonal n x n coefficient- */
257 : /* matrix A, which is stored in the 5 vectors: lower, upper, diag, */
258 : /* lowrow and ricol as per below: */
259 : /* */
260 : /* ( diag[0] upper[0] 0 0 . . 0 ricol[0] ) */
261 : /* ( lower[1] diag[1] upper[1] 0 . . 0 ) */
262 : /* ( 0 lower[2] diag[2] upper[2] 0 . ) */
263 : /* A = ( . 0 lower[3] . . . . ) */
264 : /* ( . . . . . 0 ) */
265 : /* ( . . . . . ) */
266 : /* ( 0 . . . upper[n-2] ) */
267 : /* ( lowrow[0] 0 . . 0 lower[n-1] diag[n-1] ) */
268 : /* */
269 : /* Memory for lowrow[1],..,lowrow[n-3] und ricol[1],...,ricol[n-3] */
270 : /* should be provided separately, as this should be available to */
271 : /* store the decomposition matrix, which is overwritting */
272 : /* the 5 vectors mentioned. */
273 : /* */
274 : /*====================================================================*/
275 : /* */
276 : /* Usage: */
277 : /* ====== */
278 : /* Predominantly for diagonal dominant cyclic tridiagonal- */
279 : /* matrices as they occur in spline-interpolations. */
280 : /* For diagonal dominant matices only a LU-decomposition exists. */
281 : /* */
282 : /*====================================================================*/
283 : /* */
284 : /* Input parameters: */
285 : /* ================= */
286 : /* n Dimension of the matrix ( > 2 ) sal_uInt16 n */
287 : /* lower lower antidiagonal double lower[n] */
288 : /* diag main diagonal double diag[n] */
289 : /* upper upper antidiagonal double upper[n] */
290 : /* b right side of the system double b[n] */
291 : /* rep = FALSE first call bool rep */
292 : /* = TRUE repeated call */
293 : /* for equal matrix, */
294 : /* but different b. */
295 : /* */
296 : /* Output parameters: */
297 : /* ================== */
298 : /* b solution vector of the system, double b[n] */
299 : /* the original right side is overwritten */
300 : /* */
301 : /* lower ) contain for rep = false the solution of the matrix;*/
302 : /* diag ) the original values of lower and diagonal will be */
303 : /* upper ) overwritten */
304 : /* lowrow ) double lowrow[n-2] */
305 : /* ricol ) double ricol[n-2] */
306 : /* */
307 : /* The determinant of the matrix is for rep = false */
308 : /* det A = diag[0] * ... * diag[n-1] defined . */
309 : /* */
310 : /* Return value: */
311 : /* ============= */
312 : /* = 0 all ok */
313 : /* = 1 n < 3 chosen */
314 : /* = 2 Decomposition matrix does not exist */
315 : /* */
316 : /*====================================================================*/
317 : /* */
318 : /* Used functions: */
319 : /* =============== */
320 : /* */
321 : /* from the C library: fabs() */
322 : /* */
323 : /*====================================================================*/
324 :
325 : /*.cp 5 */
326 : {
327 : double temp; // fabs(double);
328 : sal_uInt16 i;
329 : short j;
330 :
331 12 : if ( n < 3 ) return(1);
332 :
333 12 : if (!rep) /* If rep = false, */
334 : { /* calculate decomposition */
335 12 : lower[0] = upper[n-1] = 0.0; /* of the matrix. */
336 :
337 12 : if ( fabs (diag[0]) < MACH_EPS ) return(2);
338 : /* Do not decompose if the */
339 12 : temp = 1.0 / diag[0]; /* value of a diagonal */
340 12 : upper[0] *= temp; /* element is smaller then */
341 12 : ricol[0] *= temp; /* MACH_EPS */
342 :
343 132 : for (i = 1; i < n-2; i++)
344 120 : { diag[i] -= lower[i] * upper[i-1];
345 120 : if ( fabs(diag[i]) < MACH_EPS ) return(2);
346 120 : temp = 1.0 / diag[i];
347 120 : upper[i] *= temp;
348 120 : ricol[i] = -lower[i] * ricol[i-1] * temp;
349 : }
350 :
351 12 : diag[n-2] -= lower[n-2] * upper[n-3];
352 12 : if ( fabs(diag[n-2]) < MACH_EPS ) return(2);
353 :
354 132 : for (i = 1; i < n-2; i++)
355 120 : lowrow[i] = -lowrow[i-1] * upper[i-1];
356 :
357 12 : lower[n-1] -= lowrow[n-3] * upper[n-3];
358 12 : upper[n-2] = ( upper[n-2] - lower[n-2] * ricol[n-3] ) / diag[n-2];
359 :
360 144 : for (temp = 0.0, i = 0; i < n-2; i++)
361 132 : temp -= lowrow[i] * ricol[i];
362 12 : diag[n-1] += temp - lower[n-1] * upper[n-2];
363 :
364 12 : if ( fabs(diag[n-1]) < MACH_EPS ) return(2);
365 : }
366 :
367 12 : b[0] /= diag[0]; /* forward elimination */
368 144 : for (i = 1; i < n-1; i++)
369 132 : b[i] = ( b[i] - b[i-1] * lower[i] ) / diag[i];
370 :
371 144 : for (temp = 0.0, i = 0; i < n-2; i++)
372 132 : temp -= lowrow[i] * b[i];
373 :
374 12 : b[n-1] = ( b[n-1] + temp - lower[n-1] * b[n-2] ) / diag[n-1];
375 :
376 12 : b[n-2] -= b[n-1] * upper[n-2]; /* backward elimination */
377 144 : for (j = n-3; j >= 0; j--) {
378 132 : i=j;
379 132 : b[i] -= upper[i] * b[i+1] + ricol[i] * b[n-1];
380 : }
381 12 : return(0);
382 : }
383 :
384 : /*------------------ END of CYCLIC TRIDIAGONAL ---------------------*/
385 :
386 : } // extern "C"
387 :
388 : // Calculates the coefficients of natural cubic splines with n intervals.
389 0 : sal_uInt16 NaturalSpline(sal_uInt16 n, double* x, double* y,
390 : double Marg0, double MargN,
391 : sal_uInt8 MargCond,
392 : double* b, double* c, double* d)
393 : {
394 : sal_uInt16 i;
395 0 : boost::scoped_array<double> a;
396 0 : boost::scoped_array<double> h;
397 : sal_uInt16 error;
398 :
399 0 : if (n<2) return 1;
400 0 : if ( (MargCond & ~3) ) return 2;
401 0 : a.reset(new double[n+1]);
402 0 : h.reset(new double[n+1]);
403 0 : for (i=0;i<n;i++) {
404 0 : h[i]=x[i+1]-x[i];
405 0 : if (h[i]<=0.0) return 1;
406 : }
407 0 : for (i=0;i<n-1;i++) {
408 0 : a[i]=3.0*((y[i+2]-y[i+1])/h[i+1]-(y[i+1]-y[i])/h[i]);
409 0 : b[i]=h[i];
410 0 : c[i]=h[i+1];
411 0 : d[i]=2.0*(h[i]+h[i+1]);
412 : }
413 0 : switch (MargCond) {
414 : case 0: {
415 0 : if (n==2) {
416 0 : a[0]=a[0]/3.0;
417 0 : d[0]=d[0]*0.5;
418 : } else {
419 0 : a[0] =a[0]*h[1]/(h[0]+h[1]);
420 0 : a[n-2]=a[n-2]*h[n-2]/(h[n-1]+h[n-2]);
421 0 : d[0] =d[0]-h[0];
422 0 : d[n-2]=d[n-2]-h[n-1];
423 0 : c[0] =c[0]-h[0];
424 0 : b[n-2]=b[n-2]-h[n-1];
425 : }
426 : }
427 : case 1: {
428 0 : a[0] =a[0]-1.5*((y[1]-y[0])/h[0]-Marg0);
429 0 : a[n-2]=a[n-2]-1.5*(MargN-(y[n]-y[n-1])/h[n-1]);
430 0 : d[0] =d[0]-h[0]*0.5;
431 0 : d[n-2]=d[n-2]-h[n-1]*0.5;
432 : }
433 : case 2: {
434 0 : a[0] =a[0]-h[0]*Marg0*0.5;
435 0 : a[n-2]=a[n-2]-h[n-1]*MargN*0.5;
436 : }
437 : case 3: {
438 0 : a[0] =a[0]+Marg0*h[0]*h[0]*0.5;
439 0 : a[n-2]=a[n-2]-MargN*h[n-1]*h[n-1]*0.5;
440 0 : d[0] =d[0]+h[0];
441 0 : d[n-2]=d[n-2]+h[n-1];
442 : }
443 : } // switch MargCond
444 0 : if (n==2) {
445 0 : c[1]=a[0]/d[0];
446 : } else {
447 0 : error=TriDiagGS(false,n-1,b,d,c,a.get());
448 0 : if (error!=0) return error+2;
449 0 : for (i=0;i<n-1;i++) c[i+1]=a[i];
450 : }
451 0 : switch (MargCond) {
452 : case 0: {
453 0 : if (n==2) {
454 0 : c[2]=c[1];
455 0 : c[0]=c[1];
456 : } else {
457 0 : c[0]=c[1]+h[0]*(c[1]-c[2])/h[1];
458 0 : c[n]=c[n-1]+h[n-1]*(c[n-1]-c[n-2])/h[n-2];
459 : }
460 : }
461 : case 1: {
462 0 : c[0]=1.5*((y[1]-y[0])/h[0]-Marg0);
463 0 : c[0]=(c[0]-c[1]*h[0]*0.5)/h[0];
464 0 : c[n]=1.5*((y[n]-y[n-1])/h[n-1]-MargN);
465 0 : c[n]=(c[n]-c[n-1]*h[n-1]*0.5)/h[n-1];
466 : }
467 : case 2: {
468 0 : c[0]=Marg0*0.5;
469 0 : c[n]=MargN*0.5;
470 : }
471 : case 3: {
472 0 : c[0]=c[1]-Marg0*h[0]*0.5;
473 0 : c[n]=c[n-1]+MargN*h[n-1]*0.5;
474 : }
475 : } // switch MargCond
476 0 : for (i=0;i<n;i++) {
477 0 : b[i]=(y[i+1]-y[i])/h[i]-h[i]*(c[i+1]+2.0*c[i])/3.0;
478 0 : d[i]=(c[i+1]-c[i])/(3.0*h[i]);
479 : }
480 0 : return 0;
481 : }
482 :
483 : // calculates the coefficients of periodical cubic splines with n intervals.
484 12 : sal_uInt16 PeriodicSpline(sal_uInt16 n, double* x, double* y,
485 : double* b, double* c, double* d)
486 : { // array dimensions should range from [0..n]!
487 : sal_uInt16 Error;
488 : sal_uInt16 i,im1,nm1; //integer
489 : double hr,hl;
490 12 : boost::scoped_array<double> a;
491 24 : boost::scoped_array<double> lowrow;
492 24 : boost::scoped_array<double> ricol;
493 :
494 12 : if (n<2) return 4;
495 12 : nm1=n-1;
496 12 : for (i=0;i<=nm1;i++) if (x[i+1]<=x[i]) return 2; // should be strictly monotonically decreasing!
497 12 : if (y[n]!=y[0]) return 3; // begin and end should be equal!
498 :
499 12 : a.reset(new double[n+1]);
500 12 : lowrow.reset(new double[n+1]);
501 12 : ricol.reset(new double[n+1]);
502 :
503 12 : if (n==2) {
504 0 : c[1]=3.0*((y[2]-y[1])/(x[2]-x[1]));
505 0 : c[1]=c[1]-3.0*((y[i]-y[0])/(x[1]-x[0]));
506 0 : c[1]=c[1]/(x[2]-x[0]);
507 0 : c[2]=-c[1];
508 : } else {
509 156 : for (i=1;i<=nm1;i++) {
510 144 : im1=i-1;
511 144 : hl=x[i]-x[im1];
512 144 : hr=x[i+1]-x[i];
513 144 : b[im1]=hl;
514 144 : d[im1]=2.0*(hl+hr);
515 144 : c[im1]=hr;
516 144 : a[im1]=3.0*((y[i+1]-y[i])/hr-(y[i]-y[im1])/hl);
517 : }
518 12 : hl=x[n]-x[nm1];
519 12 : hr=x[1]-x[0];
520 12 : b[nm1]=hl;
521 12 : d[nm1]=2.0*(hl+hr);
522 12 : lowrow[0]=hr;
523 12 : ricol[0]=hr;
524 12 : a[nm1]=3.0*((y[1]-y[0])/hr-(y[n]-y[nm1])/hl);
525 12 : Error=ZyklTriDiagGS(false,n,b,d,c,lowrow.get(),ricol.get(),a.get());
526 12 : if ( Error != 0 )
527 : {
528 0 : return(Error+4);
529 : }
530 12 : for (i=0;i<=nm1;i++) c[i+1]=a[i];
531 : }
532 12 : c[0]=c[n];
533 168 : for (i=0;i<=nm1;i++) {
534 156 : hl=x[i+1]-x[i];
535 156 : b[i]=(y[i+1]-y[i])/hl;
536 156 : b[i]=b[i]-hl*(c[i+1]+2.0*c[i])/3.0;
537 156 : d[i]=(c[i+1]-c[i])/hl/3.0;
538 : }
539 24 : return 0;
540 : }
541 :
542 : // calculate the coefficients of parametric natural of periodical cubic splines
543 : // with n intervals
544 6 : sal_uInt16 ParaSpline(sal_uInt16 n, double* x, double* y, sal_uInt8 MargCond,
545 : double Marg01, double Marg02,
546 : double MargN1, double MargN2,
547 : bool CondT, double* T,
548 : double* bx, double* cx, double* dx,
549 : double* by, double* cy, double* dy)
550 : {
551 : sal_uInt16 Error;
552 : sal_uInt16 i;
553 : double deltX,deltY,delt,
554 6 : alphX = 0,alphY = 0,
555 6 : betX = 0,betY = 0;
556 :
557 6 : if (n<2) return 1;
558 6 : if ((MargCond & ~3) && (MargCond != 4)) return 2; // invalid boundary condition
559 6 : if (!CondT) {
560 6 : T[0]=0.0;
561 84 : for (i=0;i<n;i++) {
562 78 : deltX=x[i+1]-x[i]; deltY=y[i+1]-y[i];
563 78 : delt =deltX*deltX+deltY*deltY;
564 78 : if (delt<=0.0) return 3; // two identical adjacent points!
565 78 : T[i+1]=T[i]+sqrt(delt);
566 : }
567 : }
568 6 : switch (MargCond) {
569 0 : case 0: break;
570 : case 1: case 2: {
571 0 : alphX=Marg01; betX=MargN1;
572 0 : alphY=Marg02; betY=MargN2;
573 0 : } break;
574 : case 3: {
575 6 : if (x[n]!=x[0]) return 3;
576 6 : if (y[n]!=y[0]) return 4;
577 6 : } break;
578 : case 4: {
579 0 : if (abs(Marg01)>=MAXROOT) {
580 0 : alphX=0.0;
581 0 : alphY=sign(1.0,y[1]-y[0]);
582 : } else {
583 0 : alphX=sign(sqrt(1.0/(1.0+Marg01*Marg01)),x[1]-x[0]);
584 0 : alphY=alphX*Marg01;
585 : }
586 0 : if (abs(MargN1)>=MAXROOT) {
587 0 : betX=0.0;
588 0 : betY=sign(1.0,y[n]-y[n-1]);
589 : } else {
590 0 : betX=sign(sqrt(1.0/(1.0+MargN1*MargN1)),x[n]-x[n-1]);
591 0 : betY=betX*MargN1;
592 : }
593 : }
594 : } // switch MargCond
595 6 : if (MargCond==3) {
596 6 : Error=PeriodicSpline(n,T,x,bx,cx,dx);
597 6 : if (Error!=0) return(Error+4);
598 6 : Error=PeriodicSpline(n,T,y,by,cy,dy);
599 6 : if (Error!=0) return(Error+10);
600 : } else {
601 0 : Error=NaturalSpline(n,T,x,alphX,betX,MargCond,bx,cx,dx);
602 0 : if (Error!=0) return(Error+4);
603 0 : Error=NaturalSpline(n,T,y,alphY,betY,MargCond,by,cy,dy);
604 0 : if (Error!=0) return(Error+9);
605 : }
606 6 : return 0;
607 : }
608 :
609 6 : bool CalcSpline(Polygon& rPoly, bool Periodic, sal_uInt16& n,
610 : double*& ax, double*& ay, double*& bx, double*& by,
611 : double*& cx, double*& cy, double*& dx, double*& dy, double*& T)
612 : {
613 : sal_uInt8 Marg;
614 : double Marg01;
615 : double MargN1,MargN2;
616 : sal_uInt16 i;
617 6 : Point P0(-32768,-32768);
618 6 : Point Pt;
619 :
620 6 : n=rPoly.GetSize();
621 6 : ax=new double[rPoly.GetSize()+2];
622 6 : ay=new double[rPoly.GetSize()+2];
623 :
624 6 : n=0;
625 84 : for (i=0;i<rPoly.GetSize();i++) {
626 78 : Pt=rPoly.GetPoint(i);
627 78 : if (i==0 || Pt!=P0) {
628 78 : ax[n]=Pt.X();
629 78 : ay[n]=Pt.Y();
630 78 : n++;
631 78 : P0=Pt;
632 : }
633 : }
634 :
635 6 : if (Periodic) {
636 6 : Marg=3;
637 6 : ax[n]=ax[0];
638 6 : ay[n]=ay[0];
639 6 : n++;
640 : } else {
641 0 : Marg=2;
642 : }
643 :
644 6 : bx=new double[n+1];
645 6 : by=new double[n+1];
646 6 : cx=new double[n+1];
647 6 : cy=new double[n+1];
648 6 : dx=new double[n+1];
649 6 : dy=new double[n+1];
650 6 : T =new double[n+1];
651 :
652 6 : Marg01=0.0;
653 6 : MargN1=0.0;
654 6 : MargN2=0.0;
655 6 : if (n>0) n--; // correct n (number of partial polynoms)
656 :
657 6 : bool bRet = false;
658 6 : if ( ( Marg == 3 && n >= 3 ) || ( Marg == 2 && n >= 2 ) )
659 : {
660 6 : bRet = ParaSpline(n,ax,ay,Marg,Marg01,Marg01,MargN1,MargN2,false,T,bx,cx,dx,by,cy,dy) == 0;
661 : }
662 6 : if ( !bRet )
663 : {
664 0 : delete[] ax;
665 0 : delete[] ay;
666 0 : delete[] bx;
667 0 : delete[] by;
668 0 : delete[] cx;
669 0 : delete[] cy;
670 0 : delete[] dx;
671 0 : delete[] dy;
672 0 : delete[] T;
673 0 : n=0;
674 : }
675 6 : return bRet;
676 : }
677 :
678 6 : bool Spline2Poly(Polygon& rSpln, bool Periodic, Polygon& rPoly)
679 : {
680 6 : short MinKoord=-32000; // to prevent
681 6 : short MaxKoord=32000; // overflows
682 :
683 : double* ax; // coefficients of the polynoms
684 : double* ay;
685 : double* bx;
686 : double* by;
687 : double* cx;
688 : double* cy;
689 : double* dx;
690 : double* dy;
691 : double* tv;
692 :
693 : double Step; // stepsize for t
694 : double dt1,dt2,dt3; // delta t, y, ^3
695 : double t;
696 : bool bEnde; // partial polynom ended?
697 : sal_uInt16 n; // number of partial polynoms to draw
698 : sal_uInt16 i; // actual partial polynom
699 : bool bOk; // all still ok?
700 6 : sal_uInt16 PolyMax=16380; // max number of polygon points
701 : long x,y;
702 :
703 6 : bOk=CalcSpline(rSpln,Periodic,n,ax,ay,bx,by,cx,cy,dx,dy,tv);
704 6 : if (bOk) {
705 6 : Step =10;
706 :
707 6 : rPoly.SetSize(1);
708 6 : rPoly.SetPoint(Point(short(ax[0]),short(ay[0])),0); // first point
709 6 : i=0;
710 90 : while (i<n) { // draw n partial polynoms
711 78 : t=tv[i]+Step;
712 78 : bEnde=false;
713 21428 : while (!bEnde) { // extrapolate one partial polynom
714 21272 : bEnde=t>=tv[i+1];
715 21272 : if (bEnde) t=tv[i+1];
716 21272 : dt1=t-tv[i]; dt2=dt1*dt1; dt3=dt2*dt1;
717 21272 : x=long(ax[i]+bx[i]*dt1+cx[i]*dt2+dx[i]*dt3);
718 21272 : y=long(ay[i]+by[i]*dt1+cy[i]*dt2+dy[i]*dt3);
719 21272 : if (x<MinKoord) x=MinKoord; if (x>MaxKoord) x=MaxKoord;
720 21272 : if (y<MinKoord) y=MinKoord; if (y>MaxKoord) y=MaxKoord;
721 21272 : if (rPoly.GetSize()<PolyMax) {
722 21272 : rPoly.SetSize(rPoly.GetSize()+1);
723 21272 : rPoly.SetPoint(Point(short(x),short(y)),rPoly.GetSize()-1);
724 : } else {
725 0 : bOk=false; // error: polygon becomes to large
726 : }
727 21272 : t=t+Step;
728 : } // end of partial polynom
729 78 : i++; // next partial polynom
730 : }
731 6 : delete[] ax;
732 6 : delete[] ay;
733 6 : delete[] bx;
734 6 : delete[] by;
735 6 : delete[] cx;
736 6 : delete[] cy;
737 6 : delete[] dx;
738 6 : delete[] dy;
739 6 : delete[] tv;
740 6 : return bOk;
741 : } // end of if (bOk)
742 0 : rPoly.SetSize(0);
743 0 : return false;
744 : }
745 :
746 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|