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 6 : 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 6 : if ( n < 3 ) return(1);
332 :
333 6 : if (!rep) /* If rep = false, */
334 : { /* calculate decomposition */
335 6 : lower[0] = upper[n-1] = 0.0; /* of the matrix. */
336 :
337 6 : if ( fabs (diag[0]) < MACH_EPS ) return(2);
338 : /* Do not decompose if the */
339 6 : temp = 1.0 / diag[0]; /* value of a diagonal */
340 6 : upper[0] *= temp; /* element is smaller then */
341 6 : ricol[0] *= temp; /* MACH_EPS */
342 :
343 66 : for (i = 1; i < n-2; i++)
344 60 : { diag[i] -= lower[i] * upper[i-1];
345 60 : if ( fabs(diag[i]) < MACH_EPS ) return(2);
346 60 : temp = 1.0 / diag[i];
347 60 : upper[i] *= temp;
348 60 : ricol[i] = -lower[i] * ricol[i-1] * temp;
349 : }
350 :
351 6 : diag[n-2] -= lower[n-2] * upper[n-3];
352 6 : if ( fabs(diag[n-2]) < MACH_EPS ) return(2);
353 :
354 66 : for (i = 1; i < n-2; i++)
355 60 : lowrow[i] = -lowrow[i-1] * upper[i-1];
356 :
357 6 : lower[n-1] -= lowrow[n-3] * upper[n-3];
358 6 : upper[n-2] = ( upper[n-2] - lower[n-2] * ricol[n-3] ) / diag[n-2];
359 :
360 72 : for (temp = 0.0, i = 0; i < n-2; i++)
361 66 : temp -= lowrow[i] * ricol[i];
362 6 : diag[n-1] += temp - lower[n-1] * upper[n-2];
363 :
364 6 : if ( fabs(diag[n-1]) < MACH_EPS ) return(2);
365 : }
366 :
367 6 : b[0] /= diag[0]; /* forward elimination */
368 72 : for (i = 1; i < n-1; i++)
369 66 : b[i] = ( b[i] - b[i-1] * lower[i] ) / diag[i];
370 :
371 72 : for (temp = 0.0, i = 0; i < n-2; i++)
372 66 : temp -= lowrow[i] * b[i];
373 :
374 6 : b[n-1] = ( b[n-1] + temp - lower[n-1] * b[n-2] ) / diag[n-1];
375 :
376 6 : b[n-2] -= b[n-1] * upper[n-2]; /* backward elimination */
377 72 : for (j = n-3; j >= 0; j--) {
378 66 : i=j;
379 66 : b[i] -= upper[i] * b[i+1] + ricol[i] * b[n-1];
380 : }
381 6 : return(0);
382 : }
383 :
384 : /*------------------ END of CYCLIC TRIDIAGONAL ---------------------*/
385 :
386 : } // extern "C"
387 :
388 : /*************************************************************************
389 : |*
390 : |* NaturalSpline()
391 : |*
392 : |* Description calculates the coefficients of natural
393 : |* cubic splines with n intervals.
394 : |*
395 : *************************************************************************/
396 :
397 0 : sal_uInt16 NaturalSpline(sal_uInt16 n, double* x, double* y,
398 : double Marg0, double MargN,
399 : sal_uInt8 MargCond,
400 : double* b, double* c, double* d)
401 : {
402 : sal_uInt16 i;
403 0 : boost::scoped_array<double> a;
404 0 : boost::scoped_array<double> h;
405 : sal_uInt16 error;
406 :
407 0 : if (n<2) return 1;
408 0 : if ( (MargCond & ~3) ) return 2;
409 0 : a.reset(new double[n+1]);
410 0 : h.reset(new double[n+1]);
411 0 : for (i=0;i<n;i++) {
412 0 : h[i]=x[i+1]-x[i];
413 0 : if (h[i]<=0.0) return 1;
414 : }
415 0 : for (i=0;i<n-1;i++) {
416 0 : a[i]=3.0*((y[i+2]-y[i+1])/h[i+1]-(y[i+1]-y[i])/h[i]);
417 0 : b[i]=h[i];
418 0 : c[i]=h[i+1];
419 0 : d[i]=2.0*(h[i]+h[i+1]);
420 : }
421 0 : switch (MargCond) {
422 : case 0: {
423 0 : if (n==2) {
424 0 : a[0]=a[0]/3.0;
425 0 : d[0]=d[0]*0.5;
426 : } else {
427 0 : a[0] =a[0]*h[1]/(h[0]+h[1]);
428 0 : a[n-2]=a[n-2]*h[n-2]/(h[n-1]+h[n-2]);
429 0 : d[0] =d[0]-h[0];
430 0 : d[n-2]=d[n-2]-h[n-1];
431 0 : c[0] =c[0]-h[0];
432 0 : b[n-2]=b[n-2]-h[n-1];
433 : }
434 : }
435 : case 1: {
436 0 : a[0] =a[0]-1.5*((y[1]-y[0])/h[0]-Marg0);
437 0 : a[n-2]=a[n-2]-1.5*(MargN-(y[n]-y[n-1])/h[n-1]);
438 0 : d[0] =d[0]-h[0]*0.5;
439 0 : d[n-2]=d[n-2]-h[n-1]*0.5;
440 : }
441 : case 2: {
442 0 : a[0] =a[0]-h[0]*Marg0*0.5;
443 0 : a[n-2]=a[n-2]-h[n-1]*MargN*0.5;
444 : }
445 : case 3: {
446 0 : a[0] =a[0]+Marg0*h[0]*h[0]*0.5;
447 0 : a[n-2]=a[n-2]-MargN*h[n-1]*h[n-1]*0.5;
448 0 : d[0] =d[0]+h[0];
449 0 : d[n-2]=d[n-2]+h[n-1];
450 : }
451 : } // switch MargCond
452 0 : if (n==2) {
453 0 : c[1]=a[0]/d[0];
454 : } else {
455 0 : error=TriDiagGS(false,n-1,b,d,c,a.get());
456 0 : if (error!=0) return error+2;
457 0 : for (i=0;i<n-1;i++) c[i+1]=a[i];
458 : }
459 0 : switch (MargCond) {
460 : case 0: {
461 0 : if (n==2) {
462 0 : c[2]=c[1];
463 0 : c[0]=c[1];
464 : } else {
465 0 : c[0]=c[1]+h[0]*(c[1]-c[2])/h[1];
466 0 : c[n]=c[n-1]+h[n-1]*(c[n-1]-c[n-2])/h[n-2];
467 : }
468 : }
469 : case 1: {
470 0 : c[0]=1.5*((y[1]-y[0])/h[0]-Marg0);
471 0 : c[0]=(c[0]-c[1]*h[0]*0.5)/h[0];
472 0 : c[n]=1.5*((y[n]-y[n-1])/h[n-1]-MargN);
473 0 : c[n]=(c[n]-c[n-1]*h[n-1]*0.5)/h[n-1];
474 : }
475 : case 2: {
476 0 : c[0]=Marg0*0.5;
477 0 : c[n]=MargN*0.5;
478 : }
479 : case 3: {
480 0 : c[0]=c[1]-Marg0*h[0]*0.5;
481 0 : c[n]=c[n-1]+MargN*h[n-1]*0.5;
482 : }
483 : } // switch MargCond
484 0 : for (i=0;i<n;i++) {
485 0 : b[i]=(y[i+1]-y[i])/h[i]-h[i]*(c[i+1]+2.0*c[i])/3.0;
486 0 : d[i]=(c[i+1]-c[i])/(3.0*h[i]);
487 : }
488 0 : return 0;
489 : }
490 :
491 : /*************************************************************************
492 : |*
493 : |* PeriodicSpline()
494 : |*
495 : |* Description calculates the coefficients of periodical
496 : |* cubic splines with n intervals.
497 : |*
498 : *************************************************************************/
499 :
500 6 : sal_uInt16 PeriodicSpline(sal_uInt16 n, double* x, double* y,
501 : double* b, double* c, double* d)
502 : { // array dimensions should range from [0..n]!
503 : sal_uInt16 Error;
504 : sal_uInt16 i,im1,nm1; //integer
505 : double hr,hl;
506 6 : boost::scoped_array<double> a;
507 12 : boost::scoped_array<double> lowrow;
508 12 : boost::scoped_array<double> ricol;
509 :
510 6 : if (n<2) return 4;
511 6 : nm1=n-1;
512 6 : for (i=0;i<=nm1;i++) if (x[i+1]<=x[i]) return 2; // should be strictly monotonically decreasing!
513 6 : if (y[n]!=y[0]) return 3; // begin and end should be equal!
514 :
515 6 : a.reset(new double[n+1]);
516 6 : lowrow.reset(new double[n+1]);
517 6 : ricol.reset(new double[n+1]);
518 :
519 6 : if (n==2) {
520 0 : c[1]=3.0*((y[2]-y[1])/(x[2]-x[1]));
521 0 : c[1]=c[1]-3.0*((y[i]-y[0])/(x[1]-x[0]));
522 0 : c[1]=c[1]/(x[2]-x[0]);
523 0 : c[2]=-c[1];
524 : } else {
525 78 : for (i=1;i<=nm1;i++) {
526 72 : im1=i-1;
527 72 : hl=x[i]-x[im1];
528 72 : hr=x[i+1]-x[i];
529 72 : b[im1]=hl;
530 72 : d[im1]=2.0*(hl+hr);
531 72 : c[im1]=hr;
532 72 : a[im1]=3.0*((y[i+1]-y[i])/hr-(y[i]-y[im1])/hl);
533 : }
534 6 : hl=x[n]-x[nm1];
535 6 : hr=x[1]-x[0];
536 6 : b[nm1]=hl;
537 6 : d[nm1]=2.0*(hl+hr);
538 6 : lowrow[0]=hr;
539 6 : ricol[0]=hr;
540 6 : a[nm1]=3.0*((y[1]-y[0])/hr-(y[n]-y[nm1])/hl);
541 6 : Error=ZyklTriDiagGS(false,n,b,d,c,lowrow.get(),ricol.get(),a.get());
542 6 : if ( Error != 0 )
543 : {
544 0 : return(Error+4);
545 : }
546 6 : for (i=0;i<=nm1;i++) c[i+1]=a[i];
547 : }
548 6 : c[0]=c[n];
549 84 : for (i=0;i<=nm1;i++) {
550 78 : hl=x[i+1]-x[i];
551 78 : b[i]=(y[i+1]-y[i])/hl;
552 78 : b[i]=b[i]-hl*(c[i+1]+2.0*c[i])/3.0;
553 78 : d[i]=(c[i+1]-c[i])/hl/3.0;
554 : }
555 12 : return 0;
556 : }
557 :
558 : /*************************************************************************
559 : |*
560 : |* ParaSpline()
561 : |*
562 : |* Description calculate the coefficients of parametric
563 : |* natural of periodical cubic splines with n intervals
564 : |*
565 : *************************************************************************/
566 :
567 3 : sal_uInt16 ParaSpline(sal_uInt16 n, double* x, double* y, sal_uInt8 MargCond,
568 : double Marg01, double Marg02,
569 : double MargN1, double MargN2,
570 : bool CondT, double* T,
571 : double* bx, double* cx, double* dx,
572 : double* by, double* cy, double* dy)
573 : {
574 : sal_uInt16 Error;
575 : sal_uInt16 i;
576 : double deltX,deltY,delt,
577 3 : alphX = 0,alphY = 0,
578 3 : betX = 0,betY = 0;
579 :
580 3 : if (n<2) return 1;
581 3 : if ((MargCond & ~3) && (MargCond != 4)) return 2; // invalid boundary condition
582 3 : if (!CondT) {
583 3 : T[0]=0.0;
584 42 : for (i=0;i<n;i++) {
585 39 : deltX=x[i+1]-x[i]; deltY=y[i+1]-y[i];
586 39 : delt =deltX*deltX+deltY*deltY;
587 39 : if (delt<=0.0) return 3; // two identical adjacent points!
588 39 : T[i+1]=T[i]+sqrt(delt);
589 : }
590 : }
591 3 : switch (MargCond) {
592 0 : case 0: break;
593 : case 1: case 2: {
594 0 : alphX=Marg01; betX=MargN1;
595 0 : alphY=Marg02; betY=MargN2;
596 0 : } break;
597 : case 3: {
598 3 : if (x[n]!=x[0]) return 3;
599 3 : if (y[n]!=y[0]) return 4;
600 3 : } break;
601 : case 4: {
602 0 : if (abs(Marg01)>=MAXROOT) {
603 0 : alphX=0.0;
604 0 : alphY=sign(1.0,y[1]-y[0]);
605 : } else {
606 0 : alphX=sign(sqrt(1.0/(1.0+Marg01*Marg01)),x[1]-x[0]);
607 0 : alphY=alphX*Marg01;
608 : }
609 0 : if (abs(MargN1)>=MAXROOT) {
610 0 : betX=0.0;
611 0 : betY=sign(1.0,y[n]-y[n-1]);
612 : } else {
613 0 : betX=sign(sqrt(1.0/(1.0+MargN1*MargN1)),x[n]-x[n-1]);
614 0 : betY=betX*MargN1;
615 : }
616 : }
617 : } // switch MargCond
618 3 : if (MargCond==3) {
619 3 : Error=PeriodicSpline(n,T,x,bx,cx,dx);
620 3 : if (Error!=0) return(Error+4);
621 3 : Error=PeriodicSpline(n,T,y,by,cy,dy);
622 3 : if (Error!=0) return(Error+10);
623 : } else {
624 0 : Error=NaturalSpline(n,T,x,alphX,betX,MargCond,bx,cx,dx);
625 0 : if (Error!=0) return(Error+4);
626 0 : Error=NaturalSpline(n,T,y,alphY,betY,MargCond,by,cy,dy);
627 0 : if (Error!=0) return(Error+9);
628 : }
629 3 : return 0;
630 : }
631 :
632 : /*************************************************************************
633 : |*
634 : |* CalcSpline()
635 : |*
636 : |* Description Calculates the coefficients of parametrised
637 : |* natural or periodic cubic polynom-splines.
638 : |* The corner points of the polygon passed are used
639 : |* as support points. n returns the number of partial polynoms.
640 : |* This function returns TRUE if no error occured.
641 : |* Only in this case memory for the coefficient array
642 : |* has been allocated, which can be freed by the caller
643 : |* using a delete.
644 : *************************************************************************/
645 :
646 3 : bool CalcSpline(Polygon& rPoly, bool Periodic, sal_uInt16& n,
647 : double*& ax, double*& ay, double*& bx, double*& by,
648 : double*& cx, double*& cy, double*& dx, double*& dy, double*& T)
649 : {
650 : sal_uInt8 Marg;
651 : double Marg01;
652 : double MargN1,MargN2;
653 : sal_uInt16 i;
654 3 : Point P0(-32768,-32768);
655 3 : Point Pt;
656 :
657 3 : n=rPoly.GetSize();
658 3 : ax=new double[rPoly.GetSize()+2];
659 3 : ay=new double[rPoly.GetSize()+2];
660 :
661 3 : n=0;
662 42 : for (i=0;i<rPoly.GetSize();i++) {
663 39 : Pt=rPoly.GetPoint(i);
664 39 : if (i==0 || Pt!=P0) {
665 39 : ax[n]=Pt.X();
666 39 : ay[n]=Pt.Y();
667 39 : n++;
668 39 : P0=Pt;
669 : }
670 : }
671 :
672 3 : if (Periodic) {
673 3 : Marg=3;
674 3 : ax[n]=ax[0];
675 3 : ay[n]=ay[0];
676 3 : n++;
677 : } else {
678 0 : Marg=2;
679 : }
680 :
681 3 : bx=new double[n+1];
682 3 : by=new double[n+1];
683 3 : cx=new double[n+1];
684 3 : cy=new double[n+1];
685 3 : dx=new double[n+1];
686 3 : dy=new double[n+1];
687 3 : T =new double[n+1];
688 :
689 3 : Marg01=0.0;
690 3 : MargN1=0.0;
691 3 : MargN2=0.0;
692 3 : if (n>0) n--; // correct n (number of partial polynoms)
693 :
694 3 : bool bRet = false;
695 3 : if ( ( Marg == 3 && n >= 3 ) || ( Marg == 2 && n >= 2 ) )
696 : {
697 3 : bRet = ParaSpline(n,ax,ay,Marg,Marg01,Marg01,MargN1,MargN2,false,T,bx,cx,dx,by,cy,dy) == 0;
698 : }
699 3 : if ( !bRet )
700 : {
701 0 : delete[] ax;
702 0 : delete[] ay;
703 0 : delete[] bx;
704 0 : delete[] by;
705 0 : delete[] cx;
706 0 : delete[] cy;
707 0 : delete[] dx;
708 0 : delete[] dy;
709 0 : delete[] T;
710 0 : n=0;
711 : }
712 3 : return bRet;
713 : }
714 :
715 : /*************************************************************************
716 : |*
717 : |* Spline2Poly()
718 : |*
719 : |* Description converts a parametrised cubic spline (natural
720 : |* or periodic) to an approximate polygon.
721 : |* The function returns false, if an error occured
722 : |* during the calculation of the coefficients or
723 : |* the polygon became too large (>PolyMax=16380).
724 : |* In the first case the polygon has 0, in the
725 : |* second case PolyMax points.
726 : |* To prevent coordinate overflows we limit
727 : |* them to +/-32000.
728 : |*
729 : *************************************************************************/
730 3 : bool Spline2Poly(Polygon& rSpln, bool Periodic, Polygon& rPoly)
731 : {
732 3 : short MinKoord=-32000; // to prevent
733 3 : short MaxKoord=32000; // overflows
734 :
735 : double* ax; // coefficients of the polynoms
736 : double* ay;
737 : double* bx;
738 : double* by;
739 : double* cx;
740 : double* cy;
741 : double* dx;
742 : double* dy;
743 : double* tv;
744 :
745 : double Step; // stepsize for t
746 : double dt1,dt2,dt3; // delta t, y, ^3
747 : double t;
748 : bool bEnde; // partial polynom ended?
749 : sal_uInt16 n; // number of partial polynoms to draw
750 : sal_uInt16 i; // actual partial polynom
751 : bool bOk; // all still ok?
752 3 : sal_uInt16 PolyMax=16380; // max number of polygon points
753 : long x,y;
754 :
755 3 : bOk=CalcSpline(rSpln,Periodic,n,ax,ay,bx,by,cx,cy,dx,dy,tv);
756 3 : if (bOk) {
757 3 : Step =10;
758 :
759 3 : rPoly.SetSize(1);
760 3 : rPoly.SetPoint(Point(short(ax[0]),short(ay[0])),0); // first point
761 3 : i=0;
762 45 : while (i<n) { // draw n partial polynoms
763 39 : t=tv[i]+Step;
764 39 : bEnde=false;
765 10714 : while (!bEnde) { // extrapolate one partial polynom
766 10636 : bEnde=t>=tv[i+1];
767 10636 : if (bEnde) t=tv[i+1];
768 10636 : dt1=t-tv[i]; dt2=dt1*dt1; dt3=dt2*dt1;
769 10636 : x=long(ax[i]+bx[i]*dt1+cx[i]*dt2+dx[i]*dt3);
770 10636 : y=long(ay[i]+by[i]*dt1+cy[i]*dt2+dy[i]*dt3);
771 10636 : if (x<MinKoord) x=MinKoord; if (x>MaxKoord) x=MaxKoord;
772 10636 : if (y<MinKoord) y=MinKoord; if (y>MaxKoord) y=MaxKoord;
773 10636 : if (rPoly.GetSize()<PolyMax) {
774 10636 : rPoly.SetSize(rPoly.GetSize()+1);
775 10636 : rPoly.SetPoint(Point(short(x),short(y)),rPoly.GetSize()-1);
776 : } else {
777 0 : bOk=false; // error: polygon becomes to large
778 : }
779 10636 : t=t+Step;
780 : } // end of partial polynom
781 39 : i++; // next partial polynom
782 : }
783 3 : delete[] ax;
784 3 : delete[] ay;
785 3 : delete[] bx;
786 3 : delete[] by;
787 3 : delete[] cx;
788 3 : delete[] cy;
789 3 : delete[] dx;
790 3 : delete[] dy;
791 3 : delete[] tv;
792 3 : return bOk;
793 : } // end of if (bOk)
794 0 : rPoly.SetSize(0);
795 0 : return false;
796 : }
797 :
798 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
