source: trunk/kitgen/8.x/blt/generic/bltSpline.c@ 176

Last change on this file since 176 was 175, checked in by demin, 12 years ago

initial commit

File size: 38.8 KB
Line 
1#include "bltInt.h"
2
3typedef double TriDiagonalMatrix[3];
4typedef struct {
5 double b, c, d;
6} Cubic2D;
7
8typedef struct {
9 double b, c, d, e, f;
10} Quint2D;
11
12/*
13 * Quadratic spline parameters
14 */
15#define E1 param[0]
16#define E2 param[1]
17#define V1 param[2]
18#define V2 param[3]
19#define W1 param[4]
20#define W2 param[5]
21#define Z1 param[6]
22#define Z2 param[7]
23#define Y1 param[8]
24#define Y2 param[9]
25
26static Tcl_CmdProc SplineCmd;
27
28/*
29 * -----------------------------------------------------------------------
30 *
31 * Search --
32 *
33 * Conducts a binary search for a value. This routine is called
34 * only if key is between x(0) and x(len - 1).
35 *
36 * Results:
37 * Returns the index of the largest value in xtab for which
38 * x[i] < key.
39 *
40 * -----------------------------------------------------------------------
41 */
42static int
43Search(points, nPoints, key, foundPtr)
44 Point2D points[]; /* Contains the abscissas of the data
45 * points of interpolation. */
46 int nPoints; /* Dimension of x. */
47 double key; /* Value whose relative position in
48 * x is to be located. */
49 int *foundPtr; /* (out) Returns 1 if s is found in
50 * x and 0 otherwise. */
51{
52 int high, low, mid;
53
54 low = 0;
55 high = nPoints - 1;
56
57 while (high >= low) {
58 mid = (high + low) / 2;
59 if (key > points[mid].x) {
60 low = mid + 1;
61 } else if (key < points[mid].x) {
62 high = mid - 1;
63 } else {
64 *foundPtr = 1;
65 return mid;
66 }
67 }
68 *foundPtr = 0;
69 return low;
70}
71
72/*
73 *-----------------------------------------------------------------------
74 *
75 * QuadChoose --
76 *
77 * Determines the case needed for the computation of the parame-
78 * ters of the quadratic spline.
79 *
80 * Results:
81 * Returns a case number (1-4) which controls how the parameters
82 * of the quadratic spline are evaluated.
83 *
84 *-----------------------------------------------------------------------
85 */
86static int
87QuadChoose(p, q, m1, m2, epsilon)
88 Point2D *p; /* Coordinates of one of the points of
89 * interpolation */
90 Point2D *q; /* Coordinates of one of the points of
91 * interpolation */
92 double m1; /* Derivative condition at point P */
93 double m2; /* Derivative condition at point Q */
94 double epsilon; /* Error tolerance used to distinguish
95 * cases when m1 or m2 is relatively
96 * close to the slope or twice the
97 * slope of the line segment joining
98 * the points P and Q. If
99 * epsilon is not 0.0, then epsilon
100 * should be greater than or equal to
101 * machine epsilon. */
102{
103 double slope;
104
105 /* Calculate the slope of the line joining P and Q. */
106 slope = (q->y - p->y) / (q->x - p->x);
107
108 if (slope != 0.0) {
109 double relerr;
110 double mref, mref1, mref2, prod1, prod2;
111
112 prod1 = slope * m1;
113 prod2 = slope * m2;
114
115 /* Find the absolute values of the slopes slope, m1, and m2. */
116 mref = FABS(slope);
117 mref1 = FABS(m1);
118 mref2 = FABS(m2);
119
120 /*
121 * If the relative deviation of m1 or m2 from slope is less than
122 * epsilon, then choose case 2 or case 3.
123 */
124 relerr = epsilon * mref;
125 if ((FABS(slope - m1) > relerr) && (FABS(slope - m2) > relerr) &&
126 (prod1 >= 0.0) && (prod2 >= 0.0)) {
127 double prod;
128
129 prod = (mref - mref1) * (mref - mref2);
130 if (prod < 0.0) {
131 /*
132 * l1, the line through (x1,y1) with slope m1, and l2,
133 * the line through (x2,y2) with slope m2, intersect
134 * at a point whose abscissa is between x1 and x2.
135 * The abscissa becomes a knot of the spline.
136 */
137 return 1;
138 }
139 if (mref1 > (mref * 2.0)) {
140 if (mref2 <= ((2.0 - epsilon) * mref)) {
141 return 3;
142 }
143 } else if (mref2 <= (mref * 2.0)) {
144 /*
145 * Both l1 and l2 cross the line through
146 * (x1+x2)/2.0,y1 and (x1+x2)/2.0,y2, which is the
147 * midline of the rectangle formed by P and Q or both
148 * m1 and m2 have signs different than the sign of
149 * slope, or one of m1 and m2 has opposite sign from
150 * slope and l1 and l2 intersect to the left of x1 or
151 * to the right of x2. The point (x1+x2)/2. is a knot
152 * of the spline.
153 */
154 return 2;
155 } else if (mref1 <= ((2.0 - epsilon) * mref)) {
156 /*
157 * In cases 3 and 4, sign(m1)=sign(m2)=sign(slope).
158 * Either l1 or l2 crosses the midline, but not both.
159 * Choose case 4 if mref1 is greater than
160 * (2.-epsilon)*mref; otherwise, choose case 3.
161 */
162 return 3;
163 }
164 /*
165 * If neither l1 nor l2 crosses the midline, the spline
166 * requires two knots between x1 and x2.
167 */
168 return 4;
169 } else {
170 /*
171 * The sign of at least one of the slopes m1 or m2 does not
172 * agree with the sign of *slope*.
173 */
174 if ((prod1 < 0.0) && (prod2 < 0.0)) {
175 return 2;
176 } else if (prod1 < 0.0) {
177 if (mref2 > ((epsilon + 1.0) * mref)) {
178 return 1;
179 } else {
180 return 2;
181 }
182 } else if (mref1 > ((epsilon + 1.0) * mref)) {
183 return 1;
184 } else {
185 return 2;
186 }
187 }
188 } else if ((m1 * m2) >= 0.0) {
189 return 2;
190 } else {
191 return 1;
192 }
193}
194
195/*
196 * -----------------------------------------------------------------------
197 *
198 * QuadCases --
199 *
200 * Computes the knots and other parameters of the spline on the
201 * interval PQ.
202 *
203 *
204 * On input--
205 *
206 * P and Q are the coordinates of the points of interpolation.
207 *
208 * m1 is the slope at P.
209 *
210 * m2 is the slope at Q.
211 *
212 * ncase controls the number and location of the knots.
213 *
214 *
215 * On output--
216 *
217 * (v1,v2),(w1,w2),(z1,z2), and (e1,e2) are the coordinates of
218 * the knots and other parameters of the spline on P.
219 * (e1,e2) and Q are used only if ncase=4.
220 *
221 * -----------------------------------------------------------------------
222 */
223static void
224QuadCases(p, q, m1, m2, param, which)
225 Point2D *p, *q;
226 double m1, m2;
227 double param[];
228 int which;
229{
230 if ((which == 3) || (which == 4)) { /* Parameters used in both 3 and 4 */
231 double mbar1, mbar2, mbar3, c1, d1, h1, j1, k1;
232
233 c1 = p->x + (q->y - p->y) / m1;
234 d1 = q->x + (p->y - q->y) / m2;
235 h1 = c1 * 2.0 - p->x;
236 j1 = d1 * 2.0 - q->x;
237 mbar1 = (q->y - p->y) / (h1 - p->x);
238 mbar2 = (p->y - q->y) / (j1 - q->x);
239
240 if (which == 4) { /* Case 4. */
241 Y1 = (p->x + c1) / 2.0;
242 V1 = (p->x + Y1) / 2.0;
243 V2 = m1 * (V1 - p->x) + p->y;
244 Z1 = (d1 + q->x) / 2.0;
245 W1 = (q->x + Z1) / 2.0;
246 W2 = m2 * (W1 - q->x) + q->y;
247 mbar3 = (W2 - V2) / (W1 - V1);
248 Y2 = mbar3 * (Y1 - V1) + V2;
249 Z2 = mbar3 * (Z1 - V1) + V2;
250 E1 = (Y1 + Z1) / 2.0;
251 E2 = mbar3 * (E1 - V1) + V2;
252 } else { /* Case 3. */
253 k1 = (p->y - q->y + q->x * mbar2 - p->x * mbar1) / (mbar2 - mbar1);
254 if (FABS(m1) > FABS(m2)) {
255 Z1 = (k1 + p->x) / 2.0;
256 } else {
257 Z1 = (k1 + q->x) / 2.0;
258 }
259 V1 = (p->x + Z1) / 2.0;
260 V2 = p->y + m1 * (V1 - p->x);
261 W1 = (q->x + Z1) / 2.0;
262 W2 = q->y + m2 * (W1 - q->x);
263 Z2 = V2 + (W2 - V2) / (W1 - V1) * (Z1 - V1);
264 }
265 } else if (which == 2) { /* Case 2. */
266 Z1 = (p->x + q->x) / 2.0;
267 V1 = (p->x + Z1) / 2.0;
268 V2 = p->y + m1 * (V1 - p->x);
269 W1 = (Z1 + q->x) / 2.0;
270 W2 = q->y + m2 * (W1 - q->x);
271 Z2 = (V2 + W2) / 2.0;
272 } else { /* Case 1. */
273 double ztwo;
274
275 Z1 = (p->y - q->y + m2 * q->x - m1 * p->x) / (m2 - m1);
276 ztwo = p->y + m1 * (Z1 - p->x);
277 V1 = (p->x + Z1) / 2.0;
278 V2 = (p->y + ztwo) / 2.0;
279 W1 = (Z1 + q->x) / 2.0;
280 W2 = (ztwo + q->y) / 2.0;
281 Z2 = V2 + (W2 - V2) / (W1 - V1) * (Z1 - V1);
282 }
283}
284
285static int
286QuadSelect(p, q, m1, m2, epsilon, param)
287 Point2D *p, *q;
288 double m1, m2;
289 double epsilon;
290 double param[];
291{
292 int ncase;
293
294 ncase = QuadChoose(p, q, m1, m2, epsilon);
295 QuadCases(p, q, m1, m2, param, ncase);
296 return ncase;
297}
298
299/*
300 * -----------------------------------------------------------------------
301 *
302 * QuadGetImage --
303 *
304 * -----------------------------------------------------------------------
305 */
306INLINE static double
307QuadGetImage(p1, p2, p3, x1, x2, x3)
308 double p1, p2, p3;
309 double x1, x2, x3;
310{
311 double A, B, C;
312 double y;
313
314 A = x1 - x2;
315 B = x2 - x3;
316 C = x1 - x3;
317
318 y = (p1 * (A * A) + p2 * 2.0 * B * A + p3 * (B * B)) / (C * C);
319 return y;
320}
321
322/*
323 * -----------------------------------------------------------------------
324 *
325 * QuadSpline --
326 *
327 * Finds the image of a point in x.
328 *
329 * On input
330 *
331 * x Contains the value at which the spline is evaluated.
332 * leftX, leftY
333 * Coordinates of the left-hand data point used in the
334 * evaluation of x values.
335 * rightX, rightY
336 * Coordinates of the right-hand data point used in the
337 * evaluation of x values.
338 * Z1, Z2, Y1, Y2, E2, W2, V2
339 * Parameters of the spline.
340 * ncase Controls the evaluation of the spline by indicating
341 * whether one or two knots were placed in the interval
342 * (xtabs,xtabs1).
343 *
344 * Results:
345 * The image of the spline at x.
346 *
347 * -----------------------------------------------------------------------
348 */
349static void
350QuadSpline(intp, left, right, param, ncase)
351 Point2D *intp; /* Value at which spline is evaluated */
352 Point2D *left; /* Point to the left of the data point to
353 * be evaluated */
354 Point2D *right; /* Point to the right of the data point to
355 * be evaluated */
356 double param[]; /* Parameters of the spline */
357 int ncase; /* Controls the evaluation of the
358 * spline by indicating whether one or
359 * two knots were placed in the
360 * interval (leftX,rightX) */
361{
362 double y;
363
364 if (ncase == 4) {
365 /*
366 * Case 4: More than one knot was placed in the interval.
367 */
368
369 /*
370 * Determine the location of data point relative to the 1st knot.
371 */
372 if (Y1 > intp->x) {
373 y = QuadGetImage(left->y, V2, Y2, Y1, intp->x, left->x);
374 } else if (Y1 < intp->x) {
375 /*
376 * Determine the location of the data point relative to
377 * the 2nd knot.
378 */
379 if (Z1 > intp->x) {
380 y = QuadGetImage(Y2, E2, Z2, Z1, intp->x, Y1);
381 } else if (Z1 < intp->x) {
382 y = QuadGetImage(Z2, W2, right->y, right->x, intp->x, Z1);
383 } else {
384 y = Z2;
385 }
386 } else {
387 y = Y2;
388 }
389 } else {
390
391 /*
392 * Cases 1, 2, or 3:
393 *
394 * Determine the location of the data point relative to the
395 * knot.
396 */
397 if (Z1 < intp->x) {
398 y = QuadGetImage(Z2, W2, right->y, right->x, intp->x, Z1);
399 } else if (Z1 > intp->x) {
400 y = QuadGetImage(left->y, V2, Z2, Z1, intp->x, left->x);
401 } else {
402 y = Z2;
403 }
404 }
405 intp->y = y;
406}
407
408/*
409 * -----------------------------------------------------------------------
410 *
411 * QuadSlopes --
412 *
413 * Calculates the derivative at each of the data points. The
414 * slopes computed will insure that an osculatory quadratic
415 * spline will have one additional knot between two adjacent
416 * points of interpolation. Convexity and monotonicity are
417 * preserved wherever these conditions are compatible with the
418 * data.
419 *
420 * Results:
421 * The output array "m" is filled with the derivates at each
422 * data point.
423 *
424 * -----------------------------------------------------------------------
425 */
426static void
427QuadSlopes(points, m, nPoints)
428 Point2D points[];
429 double *m; /* (out) To be filled with the first
430 * derivative at each data point. */
431 int nPoints; /* Number of data points (dimension of
432 * x, y, and m). */
433{
434 double xbar, xmid, xhat, ydif1, ydif2;
435 double yxmid;
436 double m1, m2;
437 double m1s, m2s;
438 register int i, n, l;
439
440 m1s = m2s = m1 = m2 = 0;
441 for (l = 0, i = 1, n = 2; i < (nPoints - 1); l++, i++, n++) {
442 /*
443 * Calculate the slopes of the two lines joining three
444 * consecutive data points.
445 */
446 ydif1 = points[i].y - points[l].y;
447 ydif2 = points[n].y - points[i].y;
448 m1 = ydif1 / (points[i].x - points[l].x);
449 m2 = ydif2 / (points[n].x - points[i].x);
450 if (i == 1) {
451 m1s = m1, m2s = m2; /* Save slopes of starting point */
452 }
453 /*
454 * If one of the preceding slopes is zero or if they have opposite
455 * sign, assign the value zero to the derivative at the middle
456 * point.
457 */
458 if ((m1 == 0.0) || (m2 == 0.0) || ((m1 * m2) <= 0.0)) {
459 m[i] = 0.0;
460 } else if (FABS(m1) > FABS(m2)) {
461 /*
462 * Calculate the slope by extending the line with slope m1.
463 */
464 xbar = ydif2 / m1 + points[i].x;
465 xhat = (xbar + points[n].x) / 2.0;
466 m[i] = ydif2 / (xhat - points[i].x);
467 } else {
468 /*
469 * Calculate the slope by extending the line with slope m2.
470 */
471 xbar = -ydif1 / m2 + points[i].x;
472 xhat = (points[l].x + xbar) / 2.0;
473 m[i] = ydif1 / (points[i].x - xhat);
474 }
475 }
476
477 /* Calculate the slope at the last point, x(n). */
478 i = nPoints - 2;
479 n = nPoints - 1;
480 if ((m1 * m2) < 0.0) {
481 m[n] = m2 * 2.0;
482 } else {
483 xmid = (points[i].x + points[n].x) / 2.0;
484 yxmid = m[i] * (xmid - points[i].x) + points[i].y;
485 m[n] = (points[n].y - yxmid) / (points[n].x - xmid);
486 if ((m[n] * m2) < 0.0) {
487 m[n] = 0.0;
488 }
489 }
490
491 /* Calculate the slope at the first point, x(0). */
492 if ((m1s * m2s) < 0.0) {
493 m[0] = m1s * 2.0;
494 } else {
495 xmid = (points[0].x + points[1].x) / 2.0;
496 yxmid = m[1] * (xmid - points[1].x) + points[1].y;
497 m[0] = (yxmid - points[0].y) / (xmid - points[0].x);
498 if ((m[0] * m1s) < 0.0) {
499 m[0] = 0.0;
500 }
501 }
502
503}
504
505/*
506 * -----------------------------------------------------------------------
507 *
508 * QuadEval --
509 *
510 * QuadEval controls the evaluation of an osculatory quadratic
511 * spline. The user may provide his own slopes at the points of
512 * interpolation or use the subroutine 'QuadSlopes' to calculate
513 * slopes which are consistent with the shape of the data.
514 *
515 * ON INPUT--
516 * intpPts must be a nondecreasing vector of points at which the
517 * spline will be evaluated.
518 * origPts contains the abscissas of the data points to be
519 * interpolated. xtab must be increasing.
520 * y contains the ordinates of the data points to be
521 * interpolated.
522 * m contains the slope of the spline at each point of
523 * interpolation.
524 * nPoints number of data points (dimension of xtab and y).
525 * numEval is the number of points of evaluation (dimension of
526 * xval and yval).
527 * epsilon is a relative error tolerance used in subroutine
528 * 'QuadChoose' to distinguish the situation m(i) or
529 * m(i+1) is relatively close to the slope or twice
530 * the slope of the linear segment between xtab(i) and
531 * xtab(i+1). If this situation occurs, roundoff may
532 * cause a change in convexity or monotonicity of the
533 * resulting spline and a change in the case number
534 * provided by 'QuadChoose'. If epsilon is not equal to zero,
535 * then epsilon should be greater than or equal to machine
536 * epsilon.
537 * ON OUTPUT--
538 * yval contains the images of the points in xval.
539 * err is one of the following error codes:
540 * 0 - QuadEval ran normally.
541 * 1 - xval(i) is less than xtab(1) for at least one
542 * i or xval(i) is greater than xtab(num) for at
543 * least one i. QuadEval will extrapolate to provide
544 * function values for these abscissas.
545 * 2 - xval(i+1) < xval(i) for some i.
546 *
547 *
548 * QuadEval calls the following subroutines or functions:
549 * Search
550 * QuadCases
551 * QuadChoose
552 * QuadSpline
553 * -----------------------------------------------------------------------
554 */
555static int
556QuadEval(origPts, nOrigPts, intpPts, nIntpPts, m, epsilon)
557 Point2D origPts[];
558 int nOrigPts;
559 Point2D intpPts[];
560 int nIntpPts;
561 double *m; /* Slope of the spline at each point
562 * of interpolation. */
563 double epsilon; /* Relative error tolerance (see choose) */
564{
565 int error;
566 register int i, j;
567 double param[10];
568 int ncase;
569 int start, end;
570 int l, p;
571 register int n;
572 int found;
573
574 /* Initialize indices and set error result */
575 error = 0;
576 l = nOrigPts - 1;
577 p = l - 1;
578 ncase = 1;
579
580 /*
581 * Determine if abscissas of new vector are non-decreasing.
582 */
583 for (j = 1; j < nIntpPts; j++) {
584 if (intpPts[j].x < intpPts[j - 1].x) {
585 return 2;
586 }
587 }
588 /*
589 * Determine if any of the points in xval are LESS than the
590 * abscissa of the first data point.
591 */
592 for (start = 0; start < nIntpPts; start++) {
593 if (intpPts[start].x >= origPts[0].x) {
594 break;
595 }
596 }
597 /*
598 * Determine if any of the points in xval are GREATER than the
599 * abscissa of the l data point.
600 */
601 for (end = nIntpPts - 1; end >= 0; end--) {
602 if (intpPts[end].x <= origPts[l].x) {
603 break;
604 }
605 }
606
607 if (start > 0) {
608 error = 1; /* Set error value to indicate that
609 * extrapolation has occurred. */
610 /*
611 * Calculate the images of points of evaluation whose abscissas
612 * are less than the abscissa of the first data point.
613 */
614 ncase = QuadSelect(origPts, origPts + 1, m[0], m[1], epsilon, param);
615 for (j = 0; j < (start - 1); j++) {
616 QuadSpline(intpPts + j, origPts, origPts + 1, param, ncase);
617 }
618 if (nIntpPts == 1) {
619 return error;
620 }
621 }
622 if ((nIntpPts == 1) && (end != (nIntpPts - 1))) {
623 goto noExtrapolation;
624 }
625
626 /*
627 * Search locates the interval in which the first in-range
628 * point of evaluation lies.
629 */
630
631 i = Search(origPts, nOrigPts, intpPts[start].x, &found);
632
633 n = i + 1;
634 if (n >= nOrigPts) {
635 n = nOrigPts - 1;
636 i = nOrigPts - 2;
637 }
638 /*
639 * If the first in-range point of evaluation is equal to one
640 * of the data points, assign the appropriate value from y.
641 * Continue until a point of evaluation is found which is not
642 * equal to a data point.
643 */
644 if (found) {
645 do {
646 intpPts[start].y = origPts[i].y;
647 start++;
648 if (start >= nIntpPts) {
649 return error;
650 }
651 } while (intpPts[start - 1].x == intpPts[start].x);
652
653 for (;;) {
654 if (intpPts[start].x < origPts[n].x) {
655 break; /* Break out of for-loop */
656 }
657 if (intpPts[start].x == origPts[n].x) {
658 do {
659 intpPts[start].y = origPts[n].y;
660 start++;
661 if (start >= nIntpPts) {
662 return error;
663 }
664 } while (intpPts[start].x == intpPts[start - 1].x);
665 }
666 i++;
667 n++;
668 }
669 }
670 /*
671 * Calculate the images of all the points which lie within
672 * range of the data.
673 */
674 if ((i > 0) || (error != 1)) {
675 ncase = QuadSelect(origPts + i, origPts + n, m[i], m[n],
676 epsilon, param);
677 }
678 for (j = start; j <= end; j++) {
679 /*
680 * If xx(j) - x(n) is negative, do not recalculate
681 * the parameters for this section of the spline since
682 * they are already known.
683 */
684 if (intpPts[j].x == origPts[n].x) {
685 intpPts[j].y = origPts[n].y;
686 continue;
687 } else if (intpPts[j].x > origPts[n].x) {
688 double delta;
689
690 /* Determine that the routine is in the correct part of
691 the spline. */
692 do {
693 i++, n++;
694 delta = intpPts[j].x - origPts[n].x;
695 } while (delta > 0.0);
696
697 if (delta < 0.0) {
698 ncase = QuadSelect(origPts + i, origPts + n, m[i],
699 m[n], epsilon, param);
700 } else if (delta == 0.0) {
701 intpPts[j].y = origPts[n].y;
702 continue;
703 }
704 }
705 QuadSpline(intpPts + j, origPts + i, origPts + n, param, ncase);
706 }
707
708 if (end == (nIntpPts - 1)) {
709 return error;
710 }
711 if ((n == l) && (intpPts[end].x != origPts[l].x)) {
712 goto noExtrapolation;
713 }
714
715 error = 1; /* Set error value to indicate that
716 * extrapolation has occurred. */
717 ncase = QuadSelect(origPts + p, origPts + l, m[p], m[l], epsilon, param);
718
719 noExtrapolation:
720 /*
721 * Calculate the images of the points of evaluation whose
722 * abscissas are greater than the abscissa of the last data point.
723 */
724 for (j = (end + 1); j < nIntpPts; j++) {
725 QuadSpline(intpPts + j, origPts + p, origPts + l, param, ncase);
726 }
727 return error;
728}
729
730/*
731 * -----------------------------------------------------------------------
732 *
733 * Shape preserving quadratic splines
734 * by D.F.Mcallister & J.A.Roulier
735 * Coded by S.L.Dodd & M.Roulier
736 * N.C.State University
737 *
738 * -----------------------------------------------------------------------
739 */
740/*
741 * Driver routine for quadratic spline package
742 * On input--
743 * X,Y Contain n-long arrays of data (x is increasing)
744 * XM Contains m-long array of x values (increasing)
745 * eps Relative error tolerance
746 * n Number of input data points
747 * m Number of output data points
748 * On output--
749 * work Contains the value of the first derivative at each data point
750 * ym Contains the interpolated spline value at each data point
751 */
752int
753Blt_QuadraticSpline(origPts, nOrigPts, intpPts, nIntpPts)
754 Point2D origPts[];
755 int nOrigPts;
756 Point2D intpPts[];
757 int nIntpPts;
758{
759 double epsilon;
760 double *work;
761 int result;
762
763 work = Blt_Malloc(nOrigPts * sizeof(double));
764 assert(work);
765
766 epsilon = 0.0; /* TBA: adjust error via command-line option */
767 /* allocate space for vectors used in calculation */
768 QuadSlopes(origPts, work, nOrigPts);
769 result = QuadEval(origPts, nOrigPts, intpPts, nIntpPts, work, epsilon);
770 Blt_Free(work);
771 if (result > 1) {
772 return FALSE;
773 }
774 return TRUE;
775}
776
777/*
778 * ------------------------------------------------------------------------
779 *
780 * Reference:
781 * Numerical Analysis, R. Burden, J. Faires and A. Reynolds.
782 * Prindle, Weber & Schmidt 1981 pp 112
783 *
784 * Parameters:
785 * origPts - vector of points, assumed to be sorted along x.
786 * intpPts - vector of new points.
787 *
788 * ------------------------------------------------------------------------
789 */
790int
791Blt_NaturalSpline(origPts, nOrigPts, intpPts, nIntpPts)
792 Point2D origPts[];
793 int nOrigPts;
794 Point2D intpPts[];
795 int nIntpPts;
796{
797 Cubic2D *eq;
798 Point2D *iPtr, *endPtr;
799 TriDiagonalMatrix *A;
800 double *dx; /* vector of deltas in x */
801 double x, dy, alpha;
802 int isKnot;
803 register int i, j, n;
804
805 dx = Blt_Malloc(sizeof(double) * nOrigPts);
806 /* Calculate vector of differences */
807 for (i = 0, j = 1; j < nOrigPts; i++, j++) {
808 dx[i] = origPts[j].x - origPts[i].x;
809 if (dx[i] < 0.0) {
810 return 0;
811 }
812 }
813 n = nOrigPts - 1; /* Number of intervals. */
814 A = Blt_Malloc(sizeof(TriDiagonalMatrix) * nOrigPts);
815 if (A == NULL) {
816 Blt_Free(dx);
817 return 0;
818 }
819 /* Vectors to solve the tridiagonal matrix */
820 A[0][0] = A[n][0] = 1.0;
821 A[0][1] = A[n][1] = 0.0;
822 A[0][2] = A[n][2] = 0.0;
823
824 /* Calculate the intermediate results */
825 for (i = 0, j = 1; j < n; j++, i++) {
826 alpha = 3.0 * ((origPts[j + 1].y / dx[j]) - (origPts[j].y / dx[i]) -
827 (origPts[j].y / dx[j]) + (origPts[i].y / dx[i]));
828 A[j][0] = 2 * (dx[j] + dx[i]) - dx[i] * A[i][1];
829 A[j][1] = dx[j] / A[j][0];
830 A[j][2] = (alpha - dx[i] * A[i][2]) / A[j][0];
831 }
832 eq = Blt_Malloc(sizeof(Cubic2D) * nOrigPts);
833
834 if (eq == NULL) {
835 Blt_Free(A);
836 Blt_Free(dx);
837 return FALSE;
838 }
839 eq[0].c = eq[n].c = 0.0;
840 for (j = n, i = n - 1; i >= 0; i--, j--) {
841 eq[i].c = A[i][2] - A[i][1] * eq[j].c;
842 dy = origPts[i+1].y - origPts[i].y;
843 eq[i].b = (dy) / dx[i] - dx[i] * (eq[j].c + 2.0 * eq[i].c) / 3.0;
844 eq[i].d = (eq[j].c - eq[i].c) / (3.0 * dx[i]);
845 }
846 Blt_Free(A);
847 Blt_Free(dx);
848
849 endPtr = intpPts + nIntpPts;
850 /* Now calculate the new values */
851 for (iPtr = intpPts; iPtr < endPtr; iPtr++) {
852 iPtr->y = 0.0;
853 x = iPtr->x;
854
855 /* Is it outside the interval? */
856 if ((x < origPts[0].x) || (x > origPts[n].x)) {
857 continue;
858 }
859 /* Search for the interval containing x in the point array */
860 i = Search(origPts, nOrigPts, x, &isKnot);
861 if (isKnot) {
862 iPtr->y = origPts[i].y;
863 } else {
864 i--;
865 x -= origPts[i].x;
866 iPtr->y = origPts[i].y +
867 x * (eq[i].b + x * (eq[i].c + x * eq[i].d));
868 }
869 }
870 Blt_Free(eq);
871 return TRUE;
872}
873
874static Blt_OpSpec splineOps[] =
875{
876 {"natural", 1, (Blt_Op)Blt_NaturalSpline, 6, 6, "x y splx sply",},
877 {"quadratic", 1, (Blt_Op)Blt_QuadraticSpline, 6, 6, "x y splx sply",},
878};
879static int nSplineOps = sizeof(splineOps) / sizeof(Blt_OpSpec);
880
881/*ARGSUSED*/
882static int
883SplineCmd(clientData, interp, argc, argv)
884 ClientData clientData; /* Not used. */
885 Tcl_Interp *interp;
886 int argc;
887 char **argv;
888{
889 Blt_Op proc;
890 Blt_Vector *x, *y, *splX, *splY;
891 double *xArr, *yArr;
892 register int i;
893 Point2D *origPts, *intpPts;
894 int nOrigPts, nIntpPts;
895
896 proc = Blt_GetOp(interp, nSplineOps, splineOps, BLT_OP_ARG1, argc, argv,0);
897 if (proc == NULL) {
898 return TCL_ERROR;
899 }
900 if ((Blt_GetVector(interp, argv[2], &x) != TCL_OK) ||
901 (Blt_GetVector(interp, argv[3], &y) != TCL_OK) ||
902 (Blt_GetVector(interp, argv[4], &splX) != TCL_OK)) {
903 return TCL_ERROR;
904 }
905 nOrigPts = Blt_VecLength(x);
906 if (nOrigPts < 3) {
907 Tcl_AppendResult(interp, "length of vector \"", argv[2], "\" is < 3",
908 (char *)NULL);
909 return TCL_ERROR;
910 }
911 for (i = 1; i < nOrigPts; i++) {
912 if (Blt_VecData(x)[i] < Blt_VecData(x)[i - 1]) {
913 Tcl_AppendResult(interp, "x vector \"", argv[2],
914 "\" must be monotonically increasing", (char *)NULL);
915 return TCL_ERROR;
916 }
917 }
918 /* Check that all the data points aren't the same. */
919 if (Blt_VecData(x)[i - 1] <= Blt_VecData(x)[0]) {
920 Tcl_AppendResult(interp, "x vector \"", argv[2],
921 "\" must be monotonically increasing", (char *)NULL);
922 return TCL_ERROR;
923 }
924 if (nOrigPts != Blt_VecLength(y)) {
925 Tcl_AppendResult(interp, "vectors \"", argv[2], "\" and \"", argv[3],
926 " have different lengths", (char *)NULL);
927 return TCL_ERROR;
928 }
929 nIntpPts = Blt_VecLength(splX);
930 if (Blt_GetVector(interp, argv[5], &splY) != TCL_OK) {
931 /*
932 * If the named vector to hold the ordinates of the spline
933 * doesn't exist, create one the same size as the vector
934 * containing the abscissas.
935 */
936 if (Blt_CreateVector(interp, argv[5], nIntpPts, &splY) != TCL_OK) {
937 return TCL_ERROR;
938 }
939 } else if (nIntpPts != Blt_VecLength(splY)) {
940 /*
941 * The x and y vectors differ in size. Make the number of ordinates
942 * the same as the number of abscissas.
943 */
944 if (Blt_ResizeVector(splY, nIntpPts) != TCL_OK) {
945 return TCL_ERROR;
946 }
947 }
948 origPts = Blt_Malloc(sizeof(Point2D) * nOrigPts);
949 if (origPts == NULL) {
950 Tcl_AppendResult(interp, "can't allocate \"", Blt_Itoa(nOrigPts),
951 "\" points", (char *)NULL);
952 return TCL_ERROR;
953 }
954 intpPts = Blt_Malloc(sizeof(Point2D) * nIntpPts);
955 if (intpPts == NULL) {
956 Tcl_AppendResult(interp, "can't allocate \"", Blt_Itoa(nIntpPts),
957 "\" points", (char *)NULL);
958 Blt_Free(origPts);
959 return TCL_ERROR;
960 }
961 xArr = Blt_VecData(x);
962 yArr = Blt_VecData(y);
963 for (i = 0; i < nOrigPts; i++) {
964 origPts[i].x = xArr[i];
965 origPts[i].y = yArr[i];
966 }
967 xArr = Blt_VecData(splX);
968 yArr = Blt_VecData(splY);
969 for (i = 0; i < nIntpPts; i++) {
970 intpPts[i].x = xArr[i];
971 intpPts[i].y = yArr[i];
972 }
973 if (!(*proc) (origPts, nOrigPts, intpPts, nIntpPts)) {
974 Tcl_AppendResult(interp, "error generating spline for \"",
975 Blt_NameOfVector(splY), "\"", (char *)NULL);
976 Blt_Free(origPts);
977 Blt_Free(intpPts);
978 return TCL_ERROR;
979 }
980 yArr = Blt_VecData(splY);
981 for (i = 0; i < nIntpPts; i++) {
982 yArr[i] = intpPts[i].y;
983 }
984 Blt_Free(origPts);
985 Blt_Free(intpPts);
986
987 /* Finally update the vector. The size of the vector hasn't
988 * changed, just the data. Reset the vector using TCL_STATIC to
989 * indicate this. */
990 if (Blt_ResetVector(splY, Blt_VecData(splY), Blt_VecLength(splY),
991 Blt_VecSize(splY), TCL_STATIC) != TCL_OK) {
992 return TCL_ERROR;
993 }
994 return TCL_OK;
995}
996
997int
998Blt_SplineInit(interp)
999 Tcl_Interp *interp;
1000{
1001 static Blt_CmdSpec cmdSpec =
1002 {"spline", SplineCmd,};
1003
1004 if (Blt_InitCmd(interp, "blt", &cmdSpec) == NULL) {
1005 return TCL_ERROR;
1006 }
1007 return TCL_OK;
1008}
1009
1010
1011#define SQR(x) ((x)*(x))
1012
1013typedef struct {
1014 double t; /* Arc length of interval. */
1015 double x; /* 2nd derivative of X with respect to T */
1016 double y; /* 2nd derivative of Y with respect to T */
1017} CubicSpline;
1018
1019
1020/*
1021 * The following two procedures solve the special linear system which arise
1022 * in cubic spline interpolation. If x is assumed cyclic ( x[i]=x[n+i] ) the
1023 * equations can be written as (i=0,1,...,n-1):
1024 * m[i][0] * x[i-1] + m[i][1] * x[i] + m[i][2] * x[i+1] = b[i] .
1025 * In matrix notation one gets A * x = b, where the matrix A is tridiagonal
1026 * with additional elements in the upper right and lower left position:
1027 * A[i][0] = A_{i,i-1} for i=1,2,...,n-1 and m[0][0] = A_{0,n-1} ,
1028 * A[i][1] = A_{i, i } for i=0,1,...,n-1
1029 * A[i][2] = A_{i,i+1} for i=0,1,...,n-2 and m[n-1][2] = A_{n-1,0}.
1030 * A should be symmetric (A[i+1][0] == A[i][2]) and positive definite.
1031 * The size of the system is given in n (n>=1).
1032 *
1033 * In the first procedure the Cholesky decomposition A = C^T * D * C
1034 * (C is upper triangle with unit diagonal, D is diagonal) is calculated.
1035 * Return TRUE if decomposition exist.
1036 */
1037static int
1038SolveCubic1(A, n)
1039 TriDiagonalMatrix A[];
1040 int n;
1041{
1042 int i;
1043 double m_ij, m_n, m_nn, d;
1044
1045 if (n < 1) {
1046 return FALSE; /* Dimension should be at least 1 */
1047 }
1048 d = A[0][1]; /* D_{0,0} = A_{0,0} */
1049 if (d <= 0.0) {
1050 return FALSE; /* A (or D) should be positive definite */
1051 }
1052 m_n = A[0][0]; /* A_{0,n-1} */
1053 m_nn = A[n - 1][1]; /* A_{n-1,n-1} */
1054 for (i = 0; i < n - 2; i++) {
1055 m_ij = A[i][2]; /* A_{i,1} */
1056 A[i][2] = m_ij / d; /* C_{i,i+1} */
1057 A[i][0] = m_n / d; /* C_{i,n-1} */
1058 m_nn -= A[i][0] * m_n; /* to get C_{n-1,n-1} */
1059 m_n = -A[i][2] * m_n; /* to get C_{i+1,n-1} */
1060 d = A[i + 1][1] - A[i][2] * m_ij; /* D_{i+1,i+1} */
1061 if (d <= 0.0) {
1062 return FALSE; /* Elements of D should be positive */
1063 }
1064 A[i + 1][1] = d;
1065 }
1066 if (n >= 2) { /* Complete last column */
1067 m_n += A[n - 2][2]; /* add A_{n-2,n-1} */
1068 A[n - 2][0] = m_n / d; /* C_{n-2,n-1} */
1069 A[n - 1][1] = d = m_nn - A[n - 2][0] * m_n; /* D_{n-1,n-1} */
1070 if (d <= 0.0) {
1071 return FALSE;
1072 }
1073 }
1074 return TRUE;
1075}
1076
1077/*
1078 * The second procedure solves the linear system, with the Cholesky
1079 * decomposition calculated above (in m[][]) and the right side b given
1080 * in x[]. The solution x overwrites the right side in x[].
1081 */
1082static void
1083SolveCubic2(A, spline, nIntervals)
1084 TriDiagonalMatrix A[];
1085 CubicSpline spline[];
1086 int nIntervals;
1087{
1088 int i;
1089 double x, y;
1090 int n, m;
1091
1092 n = nIntervals - 2;
1093 m = nIntervals - 1;
1094
1095 /* Division by transpose of C : b = C^{-T} * b */
1096 x = spline[m].x;
1097 y = spline[m].y;
1098 for (i = 0; i < n; i++) {
1099 spline[i + 1].x -= A[i][2] * spline[i].x; /* C_{i,i+1} * x(i) */
1100 spline[i + 1].y -= A[i][2] * spline[i].y; /* C_{i,i+1} * x(i) */
1101 x -= A[i][0] * spline[i].x; /* C_{i,n-1} * x(i) */
1102 y -= A[i][0] * spline[i].y; /* C_{i,n-1} * x(i) */
1103 }
1104 if (n >= 0) {
1105 /* C_{n-2,n-1} * x_{n-1} */
1106 spline[m].x = x - A[n][0] * spline[n].x;
1107 spline[m].y = y - A[n][0] * spline[n].y;
1108 }
1109 /* Division by D: b = D^{-1} * b */
1110 for (i = 0; i < nIntervals; i++) {
1111 spline[i].x /= A[i][1];
1112 spline[i].y /= A[i][1];
1113 }
1114
1115 /* Division by C: b = C^{-1} * b */
1116 x = spline[m].x;
1117 y = spline[m].y;
1118 if (n >= 0) {
1119 /* C_{n-2,n-1} * x_{n-1} */
1120 spline[n].x -= A[n][0] * x;
1121 spline[n].y -= A[n][0] * y;
1122 }
1123 for (i = (n - 1); i >= 0; i--) {
1124 /* C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1} */
1125 spline[i].x -= A[i][2] * spline[i + 1].x + A[i][0] * x;
1126 spline[i].y -= A[i][2] * spline[i + 1].y + A[i][0] * y;
1127 }
1128}
1129
1130/*
1131 * Find second derivatives (x''(t_i),y''(t_i)) of cubic spline interpolation
1132 * through list of points (x_i,y_i). The parameter t is calculated as the
1133 * length of the linear stroke. The number of points must be at least 3.
1134 * Note: For CLOSED_CONTOURs the first and last point must be equal.
1135 */
1136static CubicSpline *
1137CubicSlopes(points, nPoints, isClosed, unitX, unitY)
1138 Point2D points[];
1139 int nPoints; /* Number of points (nPoints>=3) */
1140 int isClosed; /* CLOSED_CONTOUR or OPEN_CONTOUR */
1141 double unitX, unitY; /* Unit length in x and y (norm=1) */
1142{
1143 CubicSpline *spline;
1144 register CubicSpline *s1, *s2;
1145 int n, i;
1146 double norm, dx, dy;
1147 TriDiagonalMatrix *A; /* The tri-diagonal matrix is saved here. */
1148
1149 spline = Blt_Malloc(sizeof(CubicSpline) * nPoints);
1150 if (spline == NULL) {
1151 return NULL;
1152 }
1153 A = Blt_Malloc(sizeof(TriDiagonalMatrix) * nPoints);
1154 if (A == NULL) {
1155 Blt_Free(spline);
1156 return NULL;
1157 }
1158 /*
1159 * Calculate first differences in (dxdt2[i], y[i]) and interval lengths
1160 * in dist[i]:
1161 */
1162 s1 = spline;
1163 for (i = 0; i < nPoints - 1; i++) {
1164 s1->x = points[i+1].x - points[i].x;
1165 s1->y = points[i+1].y - points[i].y;
1166
1167 /*
1168 * The Norm of a linear stroke is calculated in "normal coordinates"
1169 * and used as interval length:
1170 */
1171 dx = s1->x / unitX;
1172 dy = s1->y / unitY;
1173 s1->t = sqrt(dx * dx + dy * dy);
1174
1175 s1->x /= s1->t; /* first difference, with unit norm: */
1176 s1->y /= s1->t; /* || (dxdt2[i], y[i]) || = 1 */
1177 s1++;
1178 }
1179
1180 /*
1181 * Setup linear System: Ax = b
1182 */
1183 n = nPoints - 2; /* Without first and last point */
1184 if (isClosed) {
1185 /* First and last points must be equal for CLOSED_CONTOURs */
1186 spline[nPoints - 1].t = spline[0].t;
1187 spline[nPoints - 1].x = spline[0].x;
1188 spline[nPoints - 1].y = spline[0].y;
1189 n++; /* Add last point (= first point) */
1190 }
1191 s1 = spline, s2 = s1 + 1;
1192 for (i = 0; i < n; i++) {
1193 /* Matrix A, mainly tridiagonal with cyclic second index
1194 ("j = j+n mod n")
1195 */
1196 A[i][0] = s1->t; /* Off-diagonal element A_{i,i-1} */
1197 A[i][1] = 2.0 * (s1->t + s2->t); /* A_{i,i} */
1198 A[i][2] = s2->t; /* Off-diagonal element A_{i,i+1} */
1199
1200 /* Right side b_x and b_y */
1201 s1->x = (s2->x - s1->x) * 6.0;
1202 s1->y = (s2->y - s1->y) * 6.0;
1203
1204 /*
1205 * If the linear stroke shows a cusp of more than 90 degree,
1206 * the right side is reduced to avoid oscillations in the
1207 * spline:
1208 */
1209 /*
1210 * The Norm of a linear stroke is calculated in "normal coordinates"
1211 * and used as interval length:
1212 */
1213 dx = s1->x / unitX;
1214 dy = s1->y / unitY;
1215 norm = sqrt(dx * dx + dy * dy) / 8.5;
1216 if (norm > 1.0) {
1217 /* The first derivative will not be continuous */
1218 s1->x /= norm;
1219 s1->y /= norm;
1220 }
1221 s1++, s2++;
1222 }
1223
1224 if (!isClosed) {
1225 /* Third derivative is set to zero at both ends */
1226 A[0][1] += A[0][0]; /* A_{0,0} */
1227 A[0][0] = 0.0; /* A_{0,n-1} */
1228 A[n-1][1] += A[n-1][2]; /* A_{n-1,n-1} */
1229 A[n-1][2] = 0.0; /* A_{n-1,0} */
1230 }
1231 /* Solve linear systems for dxdt2[] and y[] */
1232
1233 if (SolveCubic1(A, n)) { /* Cholesky decomposition */
1234 SolveCubic2(A, spline, n); /* A * dxdt2 = b_x */
1235 } else { /* Should not happen, but who knows ... */
1236 Blt_Free(A);
1237 Blt_Free(spline);
1238 return NULL;
1239 }
1240 /* Shift all second derivatives one place right and update the ends. */
1241 s2 = spline + n, s1 = s2 - 1;
1242 for (/* empty */; s2 > spline; s2--, s1--) {
1243 s2->x = s1->x;
1244 s2->y = s1->y;
1245 }
1246 if (isClosed) {
1247 spline[0].x = spline[n].x;
1248 spline[0].y = spline[n].y;
1249 } else {
1250 /* Third derivative is 0.0 for the first and last interval. */
1251 spline[0].x = spline[1].x;
1252 spline[0].y = spline[1].y;
1253 spline[n + 1].x = spline[n].x;
1254 spline[n + 1].y = spline[n].y;
1255 }
1256 Blt_Free( A);
1257 return spline;
1258}
1259
1260
1261/*
1262 * Calculate interpolated values of the spline function (defined via p_cntr
1263 * and the second derivatives dxdt2[] and dydt2[]). The number of tabulated
1264 * values is n. On an equidistant grid n_intpol values are calculated.
1265 */
1266static int
1267CubicEval(origPts, nOrigPts, intpPts, nIntpPts, spline)
1268 Point2D origPts[];
1269 int nOrigPts;
1270 Point2D intpPts[];
1271 int nIntpPts;
1272 CubicSpline spline[];
1273{
1274 double t, tSkip, tMax;
1275 Point2D p, q;
1276 double d, hx, dx0, dx01, hy, dy0, dy01;
1277 register int i, j, count;
1278
1279 /* Sum the lengths of all the segments (intervals). */
1280 tMax = 0.0;
1281 for (i = 0; i < nOrigPts - 1; i++) {
1282 tMax += spline[i].t;
1283 }
1284
1285 /* Need a better way of doing this... */
1286
1287 /* The distance between interpolated points */
1288 tSkip = (1. - 1e-7) * tMax / (nIntpPts - 1);
1289
1290 t = 0.0; /* Spline parameter value. */
1291 q = origPts[0];
1292 count = 0;
1293
1294 intpPts[count++] = q; /* First point. */
1295 t += tSkip;
1296
1297 for (i = 0, j = 1; j < nOrigPts; i++, j++) {
1298 d = spline[i].t; /* Interval length */
1299 p = q;
1300 q = origPts[i+1];
1301 hx = (q.x - p.x) / d;
1302 hy = (q.y - p.y) / d;
1303 dx0 = (spline[j].x + 2 * spline[i].x) / 6.0;
1304 dy0 = (spline[j].y + 2 * spline[i].y) / 6.0;
1305 dx01 = (spline[j].x - spline[i].x) / (6.0 * d);
1306 dy01 = (spline[j].y - spline[i].y) / (6.0 * d);
1307 while (t <= spline[i].t) { /* t in current interval ? */
1308 p.x += t * (hx + (t - d) * (dx0 + t * dx01));
1309 p.y += t * (hy + (t - d) * (dy0 + t * dy01));
1310 intpPts[count++] = p;
1311 t += tSkip;
1312 }
1313 /* Parameter t relative to start of next interval */
1314 t -= spline[i].t;
1315 }
1316 return count;
1317}
1318
1319/*
1320 * Generate a cubic spline curve through the points (x_i,y_i) which are
1321 * stored in the linked list p_cntr.
1322 * The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the
1323 * parameter t is the length of the linear stroke.
1324 */
1325int
1326Blt_NaturalParametricSpline(origPts, nOrigPts, extsPtr, isClosed,
1327 intpPts, nIntpPts)
1328 Point2D origPts[];
1329 int nOrigPts;
1330 Extents2D *extsPtr;
1331 int isClosed;
1332 Point2D *intpPts;
1333 int nIntpPts;
1334{
1335 double unitX, unitY; /* To define norm (x,y)-plane */
1336 CubicSpline *spline;
1337 int result;
1338
1339 if (nOrigPts < 3) {
1340 return 0;
1341 }
1342 if (isClosed) {
1343 origPts[nOrigPts].x = origPts[0].x;
1344 origPts[nOrigPts].y = origPts[0].y;
1345 nOrigPts++;
1346 }
1347 /* Width and height of the grid is used at unit length (2d-norm) */
1348 unitX = extsPtr->right - extsPtr->left;
1349 unitY = extsPtr->bottom - extsPtr->top;
1350
1351 if (unitX < FLT_EPSILON) {
1352 unitX = FLT_EPSILON;
1353 }
1354 if (unitY < FLT_EPSILON) {
1355 unitY = FLT_EPSILON;
1356 }
1357 /* Calculate parameters for cubic spline:
1358 * t = arc length of interval.
1359 * dxdt2 = second derivatives of x with respect to t,
1360 * dydt2 = second derivatives of y with respect to t,
1361 */
1362 spline = CubicSlopes(origPts, nOrigPts, isClosed, unitX, unitY);
1363 if (spline == NULL) {
1364 return 0;
1365 }
1366 result= CubicEval(origPts, nOrigPts, intpPts, nIntpPts, spline);
1367 Blt_Free(spline);
1368 return result;
1369}
1370
1371static void
1372CatromCoeffs(p, a, b, c, d)
1373 Point2D *p;
1374 Point2D *a, *b, *c, *d;
1375{
1376 a->x = -p[0].x + 3.0 * p[1].x - 3.0 * p[2].x + p[3].x;
1377 b->x = 2.0 * p[0].x - 5.0 * p[1].x + 4.0 * p[2].x - p[3].x;
1378 c->x = -p[0].x + p[2].x;
1379 d->x = 2.0 * p[1].x;
1380 a->y = -p[0].y + 3.0 * p[1].y - 3.0 * p[2].y + p[3].y;
1381 b->y = 2.0 * p[0].y - 5.0 * p[1].y + 4.0 * p[2].y - p[3].y;
1382 c->y = -p[0].y + p[2].y;
1383 d->y = 2.0 * p[1].y;
1384}
1385
1386/*
1387 *----------------------------------------------------------------------
1388 *
1389 * Blt_ParametricCatromSpline --
1390 *
1391 * Computes a spline based upon the data points, returning a new
1392 * (larger) coordinate array or points.
1393 *
1394 * Results:
1395 * None.
1396 *
1397 *----------------------------------------------------------------------
1398 */
1399int
1400Blt_CatromParametricSpline(points, nPoints, intpPts, nIntpPts)
1401 Point2D *points;
1402 int nPoints;
1403 Point2D *intpPts;
1404 int nIntpPts;
1405{
1406 register int i;
1407 Point2D *origPts;
1408 double t;
1409 int interval;
1410 Point2D a, b, c, d;
1411
1412 assert(nPoints > 0);
1413 /*
1414 * The spline is computed in screen coordinates instead of data
1415 * points so that we can select the abscissas of the interpolated
1416 * points from each pixel horizontally across the plotting area.
1417 */
1418 origPts = Blt_Malloc((nPoints + 4) * sizeof(Point2D));
1419 memcpy(origPts + 1, points, sizeof(Point2D) * nPoints);
1420
1421 origPts[0] = origPts[1];
1422 origPts[nPoints + 2] = origPts[nPoints + 1] = origPts[nPoints];
1423
1424 for (i = 0; i < nIntpPts; i++) {
1425 interval = (int)intpPts[i].x;
1426 t = intpPts[i].y;
1427 assert(interval < nPoints);
1428 CatromCoeffs(origPts + interval, &a, &b, &c, &d);
1429 intpPts[i].x = (d.x + t * (c.x + t * (b.x + t * a.x))) / 2.0;
1430 intpPts[i].y = (d.y + t * (c.y + t * (b.y + t * a.y))) / 2.0;
1431 }
1432 Blt_Free(origPts);
1433 return 1;
1434}
Note: See TracBrowser for help on using the repository browser.