1 /*****************************************************************************
2 ** Copyright (C) 1998-2001 Ljubomir Milanovic & Horst Wagner
3 ** Copyright (C) 1999-2006 Tijs Michels
4 ** This file is part of the g2 library
5 **
6 ** This library is free software; you can redistribute it and/or
7 ** modify it under the terms of the GNU Lesser General Public
8 ** License as published by the Free Software Foundation; either
9 ** version 2.1 of the License, or (at your option) any later version.
10 **
11 ** This library is distributed in the hope that it will be useful,
12 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
13 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 ** Lesser General Public License for more details.
15 **
16 ** You should have received a copy of the GNU Lesser General Public
17 ** License along with this library; if not, write to the Free Software
18 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 ******************************************************************************/
20 /*
21 * g2_splines.c
22 * Tijs Michels
23 * tijs@users.sourceforge.net
24 * 06/16/99 : initial release
25 * 19/02/06 : eliminated duplicates by using pointers to functions
26 */
27
28 #include <math.h>
29 #include <stdio.h>
30 #include "g2.h"
31 #include "g2_util.h"
32
33 /**
34 * \ingroup graphic
35 * \defgroup splines splines
36 */
37
38 typedef void calc_f( int, const double *, int, double *);
39 typedef void calc_d( int, const double *, double, double *);
40
41 static calc_f g2_c_spline;
42 static calc_f g2_c_b_spline;
43 static calc_d g2_c_raspln;
44 static calc_f g2_c_para_3;
45 static calc_f g2_c_para_5;
46
g2_p_spline(int id,int n,const double * points,int o,calc_f * f)47 static void g2_p_spline(int id, int n, const double *points, int o, calc_f *f)
48 {
49 const int m = (n-1)*o+1;
50 double * const sxy = (double *) g2_malloc(m*2*sizeof(double));
51
52 (*f)(n, points, m, sxy);
53 g2_poly_line(id, m, sxy);
54
55 g2_free(sxy);
56 }
57
g2_p_filled_spline(int id,int n,const double * points,int o,calc_f * f)58 static void g2_p_filled_spline(int id, int n, const double *points, int o, calc_f *f)
59 {
60 const int m = (n-1)*o+2;
61 double * const sxy = (double *) g2_malloc(m*2*sizeof(double));
62
63 (*f)(n, points, m, sxy);
64 sxy[m+m-2] = points[n+n-2];
65 sxy[m+m-1] = points[1];
66 g2_filled_polygon(id, m, sxy);
67 g2_free(sxy);
68 }
69
g2_split(int n,const double * points,double * x,double * y)70 static void g2_split(int n, const double *points, double *x, double *y)
71 {
72 int i;
73 for (i = 0; i < n; i++) {
74 x[i] = points[i+i];
75 y[i] = points[i+i+1];
76 }
77 }
78
79 #define eps 1.e-12
80
g2_c_spline(int n,const double * points,int m,double * sxy)81 void g2_c_spline(int n, const double *points, int m, double *sxy)
82
83 /*
84 * FUNCTIONAL DESCRIPTION:
85 *
86 * Compute a curve of m points (sx[j],sy[j])
87 * -- j being a positive integer < m --
88 * passing through the n data points (x[i],y[i])
89 * -- i being a positive integer < n --
90 * supplied by the user.
91 * The procedure to determine sy[j] involves
92 * Young's method of successive over-relaxation.
93 *
94 * FORMAL ARGUMENTS:
95 *
96 * n number of data points
97 * points data points (x[i],y[i])
98 * m number of interpolated points; m = (n-1)*o+1
99 * for o curve points for every data point
100 * sxy interpolated points (sx[j],sy[j])
101 *
102 * IMPLICIT INPUTS: NONE
103 * IMPLICIT OUTPUTS: NONE
104 * SIDE EFFECTS: NONE
105 *
106 * REFERENCES:
107 *
108 * 1. Ralston and Wilf, Mathematical Methods for Digital Computers,
109 * Vol. II, John Wiley and Sons, New York 1967, pp. 156-158.
110 * 2. Greville, T.N.E., Ed., Proceedings of An Advanced Seminar
111 * Conducted by the Mathematics Research Center, U.S. Army,
112 * University of Wisconsin, Madison. October 7-9, 1968. Theory
113 * and Applications of Spline Functions, Academic Press,
114 * New York / London 1969, pp. 156-167.
115 *
116 * AUTHORS:
117 *
118 * Josef Heinen 04/06/88 <J.Heinen@KFA-Juelich.de>
119 * Tijs Michels 06/16/99 <t.michels@vimec.nl>
120 */
121
122 {
123 int i, j;
124 double *x, *y, *g, *h;
125 double k, u, delta_g;
126
127 if (n < 3) {
128 fputs("\nERROR calling function \"g2_c_spline\":\n"
129 "number of data points input should be at least three\n", stderr);
130 return;
131 }
132 if ((m-1)%(n-1)) {
133 fputs("\nWARNING from function \"g2_c_spline\":\n"
134 "number of curve points output for every data point input "
135 "is not an integer\n", stderr);
136 }
137
138 x = (double *) g2_malloc(n*4*sizeof(double));
139 y = x + n;
140 g = y + n;
141 h = g + n; /* for the constant copy of g */
142 g2_split(n, points, x, y);
143
144 n--; /* last value index */
145 k = x[0]; /* look up once */
146 u = (x[n] - k) / (m - 1); /* calculate step outside loop */
147 for (j = 0; j < m; j++) sxy[j+j] = j * u + k; /* x-coordinates */
148
149 for (i = 1; i < n; i++) {
150 g[i] = 2. * ((y[i+1] - y[i]) / (x[i+1] - x[i]) -
151 (y[i] - y[i-1]) / (x[i] - x[i-1]))
152 / (x[i+1] - x[i-1]); /* whereas g[i] will later be changed repeatedly */
153 h[i] = 1.5 * g[i]; /* copy h[i] of g[i] will remain constant */
154 }
155
156 k = 0.;
157
158 do {
159 for (u = 0., i = 1; i < n; i++) {
160 delta_g = .5 * (x[i] - x[i-1]) / (x[i+1] - x[i-1]);
161 delta_g = (h[i] -
162 g[i] -
163 g[i-1] * delta_g - /* 8. - 4 * sqrt(3.) */
164 g[i+1] * (.5 - delta_g)) * 1.0717967697244907832;
165 g[i] += delta_g;
166
167 if (fabs(delta_g) > u) u = fabs(delta_g);
168 } /* On loop termination u holds the largest delta_g. */
169
170 if (k == 0.) k = u * eps;
171 /* Only executed once, at the end of pass one. So k preserves
172 * the largest delta_g of pass one, multiplied by eps.
173 */
174 } while (u > k);
175
176 m += m, i = 1, j = 0;
177 do {
178 u = sxy[j++]; /* x-coordinate */
179
180 while (x[i] < u) i++;
181
182 if (--i > n) i = n;
183
184 k = (u - x[i]) / (x[i+1] - x[i]); /* calculate outside loop */
185 sxy[j++] = y[i] +
186 (y[i+1] - y[i]) * k +
187 (u - x[i]) * (u - x[i+1]) *
188 ((2. - k) * g[i] +
189 (1. + k) * g[i+1]) / 6.; /* y-coordinate */
190 } while (j < m);
191 g2_free(x);
192 }
193
194 /**
195 *
196 * Using Young's method of successive over-relaxation,
197 * plot a spline curve with \a o interpolated points per data point.
198 * So the larger \a o, the more fluent the curve.
199 *
200 * \param dev device id
201 * \param n number of data points (not the size of buffer \a points)
202 * \param points buffer of \a n data points x1, y1, ... x\a n, y\a n
203 * \param o number of interpolated points per data point
204 *
205 * \ingroup splines
206 */
g2_spline(int dev,int n,double * points,int o)207 void g2_spline(int dev, int n, double *points, int o)
208 {
209 g2_p_spline(dev, n, points, o, g2_c_spline);
210 }
211
212 /**
213 *
214 * Using Young's method of successive over-relaxation,
215 * plot a filled spline curve with \a o interpolated points per data point.
216 * So the larger \a o, the more fluent the curve.
217 *
218 * \param dev device id
219 * \param n number of data points (not the size of buffer \a points)
220 * \param points buffer of \a n data points x1, y1, ... x\a n, y\a n
221 * \param o number of interpolated points per data point
222 *
223 * \ingroup splines
224 */
g2_filled_spline(int dev,int n,double * points,int o)225 void g2_filled_spline(int dev, int n, double *points, int o)
226 {
227 g2_p_filled_spline(dev, n, points, o, g2_c_spline);
228 }
229
g2_c_b_spline(int n,const double * points,int m,double * sxy)230 void g2_c_b_spline(int n, const double *points, int m, double *sxy)
231
232 /*
233 * g2_c_b_spline takes n input points. It uses parameter t
234 * to compute sx(t) and sy(t) respectively
235 */
236
237 {
238 int i, j;
239 double *x, *y;
240 double t, bl1, bl2, bl3, bl4;
241 double interval, xi_3, yi_3, xi, yi;
242
243 if (n < 3) {
244 fputs("\nERROR calling function \"g2_c_b_spline\":\n"
245 "number of data points input should be at least three\n", stderr);
246 return;
247 }
248 x = (double *) g2_malloc(n*2*sizeof(double));
249 y = x + n;
250 g2_split(n, points, x, y);
251
252 m--; /* last value index */
253 n--; /* last value index */
254 interval = (double)n / m;
255
256 for (m += m, i = 2, j = 0; i <= n+1; i++) {
257 if (i == 2) {
258 xi_3 = 2 * x[0] - x[1];
259 yi_3 = 2 * y[0] - y[1];
260 } else {
261 xi_3 = x[i-3];
262 yi_3 = y[i-3];
263 }
264 if (i == n+1) {
265 xi = 2 * x[n] - x[n-1];
266 yi = 2 * y[n] - y[n-1];
267 } else {
268 xi = x[i];
269 yi = y[i];
270 }
271
272 t = fmod(j * interval, 1.);
273
274 while (t < 1. && j < m) {
275 bl1 = (1. - t);
276 bl2 = t * t; /* t^2 */
277 bl4 = t * bl2; /* t^3 */
278 bl3 = bl4 - bl2;
279
280 bl1 = bl1 * bl1 * bl1;
281 bl2 = 3. * (bl3 - bl2) + 4.;
282 bl3 = 3. * ( t - bl3) + 1.;
283
284 sxy[j++] = (bl1 * xi_3 + bl2 * x[i-2] + bl3 * x[i-1] + bl4 * xi) / 6.; /* x-coordinate */
285 sxy[j++] = (bl1 * yi_3 + bl2 * y[i-2] + bl3 * y[i-1] + bl4 * yi) / 6.; /* y-coordinate */
286
287 t += interval;
288 }
289 }
290 sxy[m] = x[n];
291 sxy[m+1] = y[n];
292 g2_free(x);
293 }
294
295 /**
296 *
297 * Plot a b-spline curve with \a o interpolated points per data point.
298 * So the larger \a o, the more fluent the curve.
299 * For most averaging purposes, this is the right spline.
300 *
301 * \param dev device id
302 * \param n number of data points (not the size of buffer \a points)
303 * \param points buffer of \a n data points x1, y1, ... x\a n, y\a n
304 * \param o number of interpolated points per data point
305 *
306 * \ingroup splines
307 */
g2_b_spline(int dev,int n,double * points,int o)308 void g2_b_spline(int dev, int n, double *points, int o)
309 {
310 g2_p_spline(dev, n, points, o, g2_c_b_spline);
311 }
312
313 /**
314 *
315 * Plot a filled b-spline curve with \a o interpolated points per data point.
316 * So the larger \a o, the more fluent the curve.
317 * For most averaging purposes, this is the right spline.
318 *
319 * \param dev device id
320 * \param n number of data points (not the size of buffer \a points)
321 * \param points buffer of \a n data points x1, y1, ... x\a n, y\a n
322 * \param o number of interpolated points per data point
323 *
324 * \ingroup splines
325 */
g2_filled_b_spline(int dev,int n,double * points,int o)326 void g2_filled_b_spline(int dev, int n, double *points, int o)
327 {
328 g2_p_filled_spline(dev, n, points, o, g2_c_b_spline);
329 }
330
331 /*
332 * FUNCTION g2_c_raspln
333 *
334 * FUNCTIONAL DESCRIPTION:
335 *
336 * This function draws a piecewise cubic polynomial through
337 * the specified data points. The (n-1) cubic polynomials are
338 * basically parametric cubic Hermite polynomials through the
339 * n specified data points with tangent values at the data
340 * points determined by a weighted average of the slopes of
341 * the secant lines. A tension parameter "tn" is provided to
342 * adjust the length of the tangent vector at the data points.
343 * This allows the "roundness" of the curve to be adjusted.
344 * For further information and references on this technique see:
345 *
346 * D. Kochanek and R. Bartels, Interpolating Splines With Local
347 * Tension, Continuity and Bias Control, Computer Graphics,
348 * 18(1984)3.
349 *
350 * AUTHORS:
351 *
352 * Dennis Mikkelson distributed in GPLOT Jan 7, 1988 F77
353 * Tijs Michels t.michels@vimec.nl Jun 7, 1999 C
354 *
355 * FORMAL ARGUMENTS:
356 *
357 * n number of data points, n > 2
358 * points double array holding the x and y-coords of the data points
359 * tn double parameter in [0.0, 2.0]. When tn = 0.0,
360 * the curve through the data points is very rounded.
361 * As tn increases the curve is gradually pulled tighter.
362 * When tn = 2.0, the curve is essentially a polyline
363 * through the given data points.
364 * sxy double array holding the coords of the spline curve
365 *
366 * IMPLICIT INPUTS: NONE
367 * IMPLICIT OUTPUTS: NONE
368 * SIDE EFFECTS: NONE
369 */
370
371 #define nb 40
372 /*
373 * Number of straight connecting lines of which each polynomial consists.
374 * So between one data point and the next, (nb-1) points are placed.
375 */
376
g2_c_raspln(int n,const double * points,double tn,double * sxy)377 void g2_c_raspln(int n, const double *points, double tn, double *sxy)
378 {
379 int i, j;
380 double bias, tnFactor, tangentL1, tangentL2;
381 double D1x, D1y, D2x, D2y, t1x, t1y, t2x, t2y;
382 double h1[nb+1]; /* Values of the Hermite basis functions */
383 double h2[nb+1]; /* at nb+1 evenly spaced points in [0,1] */
384 double h3[nb+1];
385 double h4[nb+1];
386
387 double * const x = (double *) g2_malloc(n*2*sizeof(double));
388 double * const y = x + n;
389 g2_split(n, points, x, y);
390
391 /*
392 * First, store the values of the Hermite basis functions in a table h[ ]
393 * so no time is wasted recalculating them
394 */
395 for (i = 0; i < nb+1; i++) {
396 double t, tt, ttt;
397 t = (double) i / nb;
398 tt = t * t;
399 ttt = t * tt;
400 h1[i] = 2. * ttt - 3. * tt + 1.;
401 h2[i] = -2. * ttt + 3. * tt;
402 h3[i] = ttt - 2. * tt + t;
403 h4[i] = ttt - tt;
404 }
405
406 /*
407 * Set local tnFactor based on input parameter tn
408 */
409 if (tn <= 0.) {
410 tnFactor = 2.;
411 fputs("g2_c_raspln: Using Tension Factor 0.0: very rounded", stderr);
412 }
413 else if (tn >= 2.) {
414 tnFactor = 0.;
415 fputs("g2_c_raspln: Using Tension Factor 2.0: not rounded at all", stderr);
416 }
417 else tnFactor = 2. - tn;
418
419 D1x = D1y = 0.; /* first point has no preceding point */
420 for (j = 0; j < n - 2; j++) {
421 t1x = x[j+1] - x[j];
422 t1y = y[j+1] - y[j];
423 t2x = x[j+2] - x[j+1];
424 t2y = y[j+2] - y[j+1];
425 tangentL1 = t1x * t1x + t1y * t1y;
426 tangentL2 = t2x * t2x + t2y * t2y;
427 if (tangentL1 + tangentL2 == 0) bias = .5;
428 else bias = tangentL2 / (tangentL1 + tangentL2);
429 D2x = tnFactor * (bias * t1x + (1 - bias) * t2x);
430 D2y = tnFactor * (bias * t1y + (1 - bias) * t2y);
431 for (i = 0; i < nb; i++) {
432 sxy[2 * nb * j + i + i] =
433 h1[i] * x[j] + h2[i] * x[j+1] + h3[i] * D1x + h4[i] * D2x;
434 sxy[2 * nb * j + i + i + 1] =
435 h1[i] * y[j] + h2[i] * y[j+1] + h3[i] * D1y + h4[i] * D2y;
436 }
437 D1x = D2x; /* store as preceding point in */
438 D1y = D2y; /* the next pass */
439 }
440
441 /*
442 * Do the last subinterval as a special case since no point follows the
443 * last point
444 */
445 for (i = 0; i < nb+1; i++) {
446 sxy[2 * nb * (n-2) + i + i] =
447 h1[i] * x[n-2] + h2[i] * x[n-1] + h3[i] * D1x;
448 sxy[2 * nb * (n-2) + i + i + 1] =
449 h1[i] * y[n-2] + h2[i] * y[n-1] + h3[i] * D1y;
450 }
451 g2_free(x);
452 }
453
454 /**
455 *
456 * Plot a piecewise cubic polynomial with adjustable roundness
457 * through the given data points.
458 * Each Hermite polynomial between two data points is made up of 40 lines.
459 * Tension factor \a tn must be between 0.0 (very rounded)
460 * and 2.0 (not rounded at all, i.e. essentially a \ref g2_poly_line "polyline").
461 *
462 * \param dev device id
463 * \param n number of data points (not the size of buffer \a points)
464 * \param points buffer of \a n data points x1, y1, ... x\a n, y\a n
465 * \param tn tension factor in the range [0.0, 2.0]
466 *
467 * \ingroup splines
468 */
g2_raspln(int dev,int n,double * points,double tn)469 void g2_raspln(int dev, int n, double *points, double tn)
470 {
471 const int m = (n-1)*nb+1;
472 double * const sxy = (double *) g2_malloc(m*2*sizeof(double)); /* coords of the entire spline curve */
473
474 g2_c_raspln(n, points, tn, sxy);
475 g2_poly_line(dev, m, sxy);
476
477 g2_free(sxy);
478 }
479
480 /**
481 *
482 * Plot a filled piecewise cubic polynomial with adjustable roundness
483 * through the given data points.
484 * Each Hermite polynomial between two data points is made up of 40 lines.
485 * Tension factor \a tn must be between 0.0 (very rounded)
486 * and 2.0 (not rounded at all, i.e. essentially a \ref g2_poly_line "polyline").
487 *
488 * \param dev device id
489 * \param n number of data points (not the size of buffer \a points)
490 * \param points buffer of \a n data points x1, y1, ... x\a n, y\a n
491 * \param tn tension factor in the range [0.0, 2.0]
492 *
493 * \ingroup splines
494 */
g2_filled_raspln(int dev,int n,double * points,double tn)495 void g2_filled_raspln(int dev, int n, double *points, double tn)
496 {
497 const int m = (n-1)*nb+2;
498 double * const sxy = (double *) g2_malloc(m*2*sizeof(double)); /* coords of the entire spline curve */
499
500 g2_c_raspln(n, points, tn, sxy);
501 sxy[m+m-2] = points[n+n-2];
502 sxy[m+m-1] = points[1];
503 g2_filled_polygon(dev, m, sxy);
504
505 g2_free(sxy);
506 }
507
508 /* ---- And now for a rather different approach ---- */
509
510 /*
511 * FUNCTION g2_c_newton
512 *
513 * FUNCTIONAL DESCRIPTION:
514 *
515 * Use Newton's Divided Differences to calculate an interpolation
516 * polynomial through the specified data points.
517 * This function is called by
518 * g2_c_para_3 and
519 * g2_c_para_5.
520 *
521 * Dennis Mikkelson distributed in GPLOT Jan 5, 1988 F77
522 * Tijs Michels t.michels@vimec.nl Jun 16, 1999 C
523 *
524 * FORMAL ARGUMENTS:
525 *
526 * n number of entries in c1 and c2, 4 <= n <= MaxPts
527 * for para_3 (degree 3) n = 4
528 * for para_5 (degree 5) n = 6
529 * for para_i (degree i) n = (i + 1)
530 * c1 double array holding at most MaxPts values giving the
531 * first coords of the points to be interpolated
532 * c2 double array holding at most MaxPts values giving the
533 * second coords of the points to be interpolated
534 * o number of points at which the interpolation
535 * polynomial is to be evaluated
536 * xv double array holding o points at which to
537 * evaluate the interpolation polynomial
538 * yv double array holding upon return the values of the
539 * interpolation polynomial at the corresponding points in xv
540 *
541 * yv is the OUTPUT
542 *
543 * IMPLICIT INPUTS: NONE
544 * IMPLICIT OUTPUTS: NONE
545 * SIDE EFFECTS: NONE
546 */
547
548 #define MaxPts 21
549 #define xstr(s) __str(s)
550 #define __str(s) #s
551
552 /*
553 * Maximum number of data points allowed
554 * 21 would correspond to a polynomial of degree 20
555 */
556
g2_c_newton(int n,const double * c1,const double * c2,int o,const double * xv,double * yv)557 static void g2_c_newton(int n, const double *c1, const double *c2,
558 int o, const double *xv, double *yv)
559 {
560 int i, j;
561 double p, s;
562 double ddt[MaxPts][MaxPts]; /* Divided Difference Table */
563
564 if (n < 4) {
565 fputs("g2_c_newton: Error! Less than 4 points passed "
566 "to function g2_c_newton\n", stderr);
567 return;
568 }
569
570 if (n > MaxPts) {
571 fputs("g2_c_newton: Error! More than " xstr(MaxPts) " points passed "
572 "to function g2_c_newton\n", stderr);
573 return;
574 }
575
576 /* First, build the divided difference table */
577
578 for (i = 0; i < n; i++) ddt[i][0] = c2[i];
579 for (j = 1; j < n; j++) {
580 for (i = 0; i < n - j; i++)
581 ddt[i][j] = (ddt[i+1][j-1] - ddt[i][j-1]) / (c1[i+j] - c1[i]);
582 }
583
584 /* Next, evaluate the polynomial at the specified points */
585
586 for (i = 0; i < o; i++) {
587 for (p = 1., s = ddt[0][0], j = 1; j < n; j++) {
588 p *= xv[i] - c1[j-1];
589 s += p * ddt[0][j];
590 }
591 yv[i] = s;
592 }
593 }
594
595 /*
596 * FUNCTION: g2_c_para_3
597 *
598 * FUNCTIONAL DESCRIPTION:
599 *
600 * This function draws a piecewise parametric interpolation
601 * polynomial of degree 3 through the specified data points.
602 * The effect is similar to that obtained using DISSPLA to
603 * draw a curve after a call to the DISSPLA routine PARA3.
604 * The curve is parameterized using an approximation to the
605 * curve's arc length. The basic interpolation is done
606 * using function g2_c_newton.
607 *
608 * Dennis Mikkelson distributed in GPLOT Jan 7, 1988 F77
609 * Tijs Michels t.michels@vimec.nl Jun 17, 1999 C
610 *
611 * FORMAL ARGUMENTS:
612 *
613 * n number of data points through which to draw the curve
614 * points double array containing the x and y-coords of the data points
615 *
616 * IMPLICIT INPUTS: NONE
617 * IMPLICIT OUTPUTS: NONE
618 * SIDE EFFECTS: NONE
619 */
620
621 /*
622 * #undef nb
623 * #define nb 40
624 * Number of straight connecting lines of which each polynomial consists.
625 * So between one data point and the next, (nb-1) points are placed.
626 */
627
g2_c_para_3(int n,const double * points,int m,double * sxy)628 void g2_c_para_3(int n, const double *points, int m, double *sxy)
629 {
630 #define dgr (3+1)
631 #define nb2 (nb*2)
632 int i, j;
633 double x1t, y1t;
634 double o, step;
635 double X[nb2]; /* x-coords of the current curve piece */
636 double Y[nb2]; /* y-coords of the current curve piece */
637 double t[dgr]; /* data point parameter values */
638 double Xpts[dgr]; /* x-coords data point subsection */
639 double Ypts[dgr]; /* y-coords data point subsection */
640 double s[nb2]; /* parameter values at which to interpolate */
641 m = m; /* dummy argument; stop compiler complaints */
642
643 /* Do first TWO subintervals first */
644
645 g2_split(dgr, points, Xpts, Ypts);
646
647 t[0] = 0.;
648 for (i = 1; i < dgr; i++) {
649 x1t = Xpts[i] - Xpts[i-1];
650 y1t = Ypts[i] - Ypts[i-1];
651 t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);
652 }
653
654 step = t[2] / nb2;
655 for (i = 0; i < nb2; i++) s[i] = i * step;
656
657 g2_c_newton(dgr, t, Xpts, nb2, s, X);
658 g2_c_newton(dgr, t, Ypts, nb2, s, Y);
659 for (i = 0; i < nb2; i++) {
660 sxy[i+i] = X[i];
661 sxy[i+i+1] = Y[i];
662 }
663
664 /* Next, do later central subintervals */
665
666 for (j = 1; j < n - dgr + 1; j++) {
667 g2_split(dgr, points + j + j, Xpts, Ypts);
668
669 for (i = 1; i < dgr; i++) {
670 x1t = Xpts[i] - Xpts[i-1];
671 y1t = Ypts[i] - Ypts[i-1];
672 t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);
673 }
674
675 o = t[1]; /* look up once */
676 step = (t[2] - o) / nb;
677 for (i = 0; i < nb; i++) s[i] = i * step + o;
678
679 g2_c_newton(dgr, t, Xpts, nb, s, X);
680 g2_c_newton(dgr, t, Ypts, nb, s, Y);
681
682 for (i = 0; i < nb; i++) {
683 sxy[(j + 1) * nb2 + i + i] = X[i];
684 sxy[(j + 1) * nb2 + i + i + 1] = Y[i];
685 }
686 }
687
688 /* Now do last subinterval */
689
690 o = t[2];
691 step = (t[3] - o) / nb;
692 for (i = 0; i < nb; i++) s[i] = i * step + o;
693
694 g2_c_newton(dgr, t, Xpts, nb, s, X);
695 g2_c_newton(dgr, t, Ypts, nb, s, Y);
696
697 for (i = 0; i < nb; i++) {
698 sxy[(n - dgr + 2) * nb2 + i + i] = X[i];
699 sxy[(n - dgr + 2) * nb2 + i + i + 1] = Y[i];
700 }
701 sxy[(n - 1) * nb2] = points[n+n-2];
702 sxy[(n - 1) * nb2 + 1] = points[n+n-1];
703 }
704
705 /**
706 *
707 * Using Newton's Divided Differences method, plot a piecewise
708 * parametric interpolation polynomial of degree 3
709 * through the given data points.
710 *
711 * \param dev device id
712 * \param n number of data points (not the size of buffer \a points)
713 * \param points buffer of \a n data points x1, y1, ... x\a n, y\a n
714 *
715 * \ingroup splines
716 */
g2_para_3(int dev,int n,double * points)717 void g2_para_3(int dev, int n, double *points)
718 {
719 g2_p_spline(dev, n, points, nb, g2_c_para_3);
720 }
721
722 /**
723 *
724 * Using Newton's Divided Differences method, plot a filled piecewise
725 * parametric interpolation polynomial of degree 3
726 * through the given data points.
727 *
728 * \param dev device id
729 * \param n number of data points (not the size of buffer \a points)
730 * \param points buffer of \a n data points x1, y1, ... x\a n, y\a n
731 *
732 * \ingroup splines
733 */
g2_filled_para_3(int dev,int n,double * points)734 void g2_filled_para_3(int dev, int n, double *points)
735 {
736 g2_p_filled_spline(dev, n, points, nb, g2_c_para_3);
737 }
738
739 /*
740 * FUNCTION: g2_c_para_5
741 *
742 * As g2_c_para_3, but now plot a polynomial of degree 5
743 */
744
745 /*
746 * Number of straight connecting lines of which each polynomial consists.
747 * So between one data point and the next, (nb-1) points are placed.
748 */
749
g2_c_para_5(int n,const double * points,int m,double * sxy)750 void g2_c_para_5(int n, const double *points, int m, double *sxy)
751 {
752 #undef dgr
753 #define dgr (5+1)
754 #define nb3 (nb*3)
755 int i, j;
756 double x1t, y1t;
757 double o, step;
758 double X[nb3]; /* x-coords of the current curve piece */
759 double Y[nb3]; /* y-coords of the current curve piece */
760 double t[dgr]; /* data point parameter values */
761 double Xpts[dgr]; /* x-coords data point subsection */
762 double Ypts[dgr]; /* y-coords data point subsection */
763 double s[nb3]; /* parameter values at which to interpolate */
764 m = m; /* dummy argument; stop compiler complaints */
765
766 /* Do first THREE subintervals first */
767
768 g2_split(dgr, points, Xpts, Ypts);
769
770 t[0] = 0.;
771 for (i = 1; i < dgr; i++) {
772 x1t = Xpts[i] - Xpts[i-1];
773 y1t = Ypts[i] - Ypts[i-1];
774 t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);
775 }
776
777 step = t[3] / nb3;
778 for (i = 0; i < nb3; i++) s[i] = i * step;
779
780 g2_c_newton(dgr, t, Xpts, nb3, s, X);
781 g2_c_newton(dgr, t, Ypts, nb3, s, Y);
782 for (i = 0; i < nb3; i++) {
783 sxy[i+i] = X[i];
784 sxy[i+i+1] = Y[i];
785 }
786
787 /* Next, do later central subintervals */
788
789 for (j = 1; j < n - dgr + 1; j++) {
790 g2_split(dgr, points + j + j, Xpts, Ypts);
791
792 for (i = 1; i < dgr; i++) {
793 x1t = Xpts[i] - Xpts[i-1];
794 y1t = Ypts[i] - Ypts[i-1];
795 t[i] = t[i-1] + sqrt(x1t * x1t + y1t * y1t);
796 }
797
798 o = t[2]; /* look up once */
799 step = (t[3] - o) / nb;
800 for (i = 0; i < nb; i++) s[i] = i * step + o;
801
802 g2_c_newton(dgr, t, Xpts, nb, s, X);
803 g2_c_newton(dgr, t, Ypts, nb, s, Y);
804
805 for (i = 0; i < nb; i++) {
806 sxy[(j + 2) * nb2 + i + i] = X[i];
807 sxy[(j + 2) * nb2 + i + i + 1] = Y[i];
808 }
809 }
810
811 /* Now do last TWO subinterval */
812
813 o = t[3];
814 step = (t[5] - o) / nb2;
815 for (i = 0; i < nb2; i++) s[i] = i * step + o;
816
817 g2_c_newton(dgr, t, Xpts, nb2, s, X);
818 g2_c_newton(dgr, t, Ypts, nb2, s, Y);
819
820 for (i = 0; i < nb2; i++) {
821 sxy[(n - dgr + 3) * nb2 + i + i] = X[i];
822 sxy[(n - dgr + 3) * nb2 + i + i + 1] = Y[i];
823 }
824 sxy[(n - 1) * nb2] = points[n+n-2];
825 sxy[(n - 1) * nb2 + 1] = points[n+n-1];
826 }
827
828 /**
829 *
830 * Using Newton's Divided Differences method, plot a piecewise
831 * parametric interpolation polynomial of degree 5
832 * through the given data points.
833 *
834 * \param dev device id
835 * \param n number of data points (not the size of buffer \a points)
836 * \param points buffer of \a n data points x1, y1, ... x\a n, y\a n
837 *
838 * \ingroup splines
839 */
g2_para_5(int dev,int n,double * points)840 void g2_para_5(int dev, int n, double *points)
841 {
842 g2_p_spline(dev, n, points, nb, g2_c_para_5);
843 }
844
845 /**
846 *
847 * Using Newton's Divided Differences method, plot a filled piecewise
848 * parametric interpolation polynomial of degree 5
849 * through the given data points.
850 *
851 * \param dev device id
852 * \param n number of data points (not the size of buffer \a points)
853 * \param points buffer of \a n data points x1, y1, ... x\a n, y\a n
854 *
855 * \ingroup splines
856 */
g2_filled_para_5(int dev,int n,double * points)857 void g2_filled_para_5(int dev, int n, double *points)
858 {
859 g2_p_filled_spline(dev, n, points, nb, g2_c_para_5);
860 }
861