Annotation of OpenXM_contrib/gnuplot/interpol.c, Revision 1.1.1.3
1.1 maekawa 1: #ifndef lint
1.1.1.3 ! ohara 2: static char *RCSid = "$Id: interpol.c,v 1.6.2.2 2002/02/01 18:41:05 broeker Exp $";
1.1 maekawa 3: #endif
4:
5: /* GNUPLOT - interpol.c */
6:
7: /*[
8: * Copyright 1986 - 1993, 1998 Thomas Williams, Colin Kelley
9: *
10: * Permission to use, copy, and distribute this software and its
11: * documentation for any purpose with or without fee is hereby granted,
12: * provided that the above copyright notice appear in all copies and
13: * that both that copyright notice and this permission notice appear
14: * in supporting documentation.
15: *
16: * Permission to modify the software is granted, but not the right to
17: * distribute the complete modified source code. Modifications are to
18: * be distributed as patches to the released version. Permission to
19: * distribute binaries produced by compiling modified sources is granted,
20: * provided you
21: * 1. distribute the corresponding source modifications from the
22: * released version in the form of a patch file along with the binaries,
23: * 2. add special version identification to distinguish your version
24: * in addition to the base release version number,
25: * 3. provide your name and address as the primary contact for the
26: * support of your modified version, and
27: * 4. retain our contact information in regard to use of the base
28: * software.
29: * Permission to distribute the released version of the source code along
30: * with corresponding source modifications in the form of a patch file is
31: * granted with same provisions 2 through 4 for binary distributions.
32: *
33: * This software is provided "as is" without express or implied warranty
34: * to the extent permitted by applicable law.
35: ]*/
36:
37:
38: /*
39: * C-Source file identification Header
40: *
41: * This file belongs to a project which is:
42: *
43: * done 1993 by MGR-Software, Asgard (Lars Hanke)
44: * written by Lars Hanke
45: *
46: * Contact me via:
47: *
48: * InterNet: mgr@asgard.bo.open.de
49: * FIDO: Lars Hanke @ 2:243/4802.22 (as long as they keep addresses)
50: *
51: **************************************************************************
52: *
53: * Project: gnuplot
54: * Module:
55: * File: interpol.c
56: *
57: * Revisor: Lars Hanke
58: * Revised: 26/09/93
59: * Revision: 1.0
60: *
61: **************************************************************************
62: *
63: * LEGAL
64: * This module is part of gnuplot and distributed under whatever terms
65: * gnuplot is or will be published, unless exclusive rights are claimed.
66: *
67: * DESCRIPTION
68: * Supplies 2-D data interpolation and approximation routines
69: *
70: * IMPORTS
71: * plot.h
72: * - cp_extend()
73: * - structs: curve_points, coordval, coordinate
74: *
75: * setshow.h
76: * - samples, is_log_x, base_log_x, xmin, xmax, autoscale_lx
77: * - plottypes
78: *
79: * proto.h
80: * - solve_tri_diag()
81: * - typedef tri_diag
82: *
83: * EXPORTS
84: * gen_interp()
85: * sort_points()
86: * cp_implode()
87: *
88: * BUGS and TODO
89: * I would really have liked to use Gershon Elbers contouring code for
90: * all the stuff done here, but I failed. So I used my own code.
91: * If somebody is able to consolidate Gershon's code for this purpose
92: * a lot of gnuplot users would be very happy - due to memory problems.
93: *
94: **************************************************************************
95: *
96: * HISTORY
97: * Changes:
98: * Nov 24, 1995 Markus Schuh (M.Schuh@meteo.uni-koeln.de):
99: * changed the algorithm for csplines
100: * added algorithm for approximation csplines
101: * copied point storage and range fix from plot2d.c
102: *
103: * Dec 12, 1995 David Denholm
104: * oops - at the time this is called, stored co-ords are
105: * internal (ie maybe log of data) but min/max are in
106: * user co-ordinates.
107: * Work with min and max of internal co-ords, and
108: * check at the end whether external min and max need to
109: * be increased. (since samples is typically 100 ; we
110: * dont want to take more logs than necessary)
111: * Also, need to take into account which axes are active
112: *
113: * Jun 30, 1996 Jens Emmerich
114: * implemented handling of UNDEFINED points
115: */
116:
117: #include "plot.h"
118: #include "setshow.h"
119:
120:
121: /* in order to support multiple axes, and to simplify ranging in
122: * parametric plots, we use arrays to store some things. For 2d plots,
123: * elements are z=0,y1=1,x1=2,z2=4,y2=5,x2=6 these are given symbolic
124: * names in plot.h
125: */
126:
127: extern double min_array[AXIS_ARRAY_SIZE];
128: extern double max_array[AXIS_ARRAY_SIZE];
129: extern int auto_array[AXIS_ARRAY_SIZE];
130: extern TBOOLEAN log_array[AXIS_ARRAY_SIZE];
131: extern double base_array[AXIS_ARRAY_SIZE];
132: extern double log_base_array[AXIS_ARRAY_SIZE];
133:
134:
135: #define Inc_c_token if (++c_token >= num_tokens) \
136: int_error ("Syntax error", c_token);
137:
138:
139: /*
140: * IMHO, code is getting too cluttered with repeated chunks of
141: * code. Some macros to simplify, I hope.
142: *
143: * do { } while(0) is comp.lang.c recommendation for complex macros
144: * also means that break can be specified as an action, and it will
145: *
146: */
147:
148:
149: /* store VALUE or log(VALUE) in STORE, set TYPE as appropriate
150: * Do OUT_ACTION or UNDEF_ACTION as appropriate
151: * adjust range provided type is INRANGE (ie dont adjust y if x is outrange
152: * VALUE must not be same as STORE
153: */
154:
155: #define STORE_AND_FIXUP_RANGE(STORE, VALUE, TYPE, MIN, MAX, AUTO, OUT_ACTION, UNDEF_ACTION)\
156: do { STORE=VALUE; \
157: if (TYPE != INRANGE) break; /* dont set y range if x is outrange, for example */ \
158: if ( VALUE<MIN ) { \
159: if (AUTO & 1) MIN=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; break; } \
160: } \
161: if ( VALUE>MAX ) {\
162: if (AUTO & 2) MAX=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; } \
163: } \
164: } while(0)
165:
166: #define UPDATE_RANGE(TEST,OLD,NEW,AXIS) \
167: do { if (TEST) { \
168: if (log_array[AXIS]) OLD = pow(base_array[AXIS], NEW); else OLD = NEW; \
169: } \
170: } while(0)
171:
172: /* use this instead empty macro arguments to work around NeXT cpp bug */
173: /* if this fails on any system, we might use ((void)0) */
174:
175: #define NOOP /* */
176:
177: #define spline_coeff_size 4
178: typedef double spline_coeff[spline_coeff_size];
179: typedef double five_diag[5];
180:
181: static int next_curve __PROTO((struct curve_points * plot, int *curve_start));
182: static int num_curves __PROTO((struct curve_points * plot));
183: static double *cp_binomial __PROTO((int points));
184: GP_INLINE static double s_pow __PROTO((double base, unsigned int exponent));
185: static void eval_bezier __PROTO((struct curve_points * cp, int first_point, int num_points, double sr, coordval * px, coordval * py, double *c));
186: static void do_bezier __PROTO((struct curve_points * cp, double *bc, int first_point, int num_points, struct coordinate * dest));
187: static int solve_five_diag __PROTO((five_diag m[], double r[], double x[], int n));
188: static spline_coeff *cp_approx_spline __PROTO((struct curve_points * plot, int first_point, int num_points));
189: static spline_coeff *cp_tridiag __PROTO((struct curve_points * plot, int first_point, int num_points));
190: static void do_cubic __PROTO((struct curve_points * plot, spline_coeff * sc, int first_point, int num_points, struct coordinate * dest));
1.1.1.3 ! ohara 191: int compare_points __PROTO((SORTFUNC_ARGS p1, SORTFUNC_ARGS p2));
1.1 maekawa 192:
193:
194: /*
195: * position curve_start to index the next non-UNDEFINDED point,
196: * start search at initial curve_start,
197: * return number of non-UNDEFINDED points from there on,
198: * if no more valid points are found, curve_start is set
199: * to plot->p_count and 0 is returned
200: */
201:
202: static int next_curve(plot, curve_start)
203: struct curve_points *plot;
204: int *curve_start;
205: {
206: int curve_length;
207:
208: /* Skip undefined points */
209: while (*curve_start < plot->p_count
210: && plot->points[*curve_start].type == UNDEFINED) {
211: (*curve_start)++;
212: };
213: curve_length = 0;
214: /* curve_length is first used as an offset, then the correkt # points */
215: while ((*curve_start) + curve_length < plot->p_count
216: && plot->points[(*curve_start) + curve_length].type != UNDEFINED) {
217: curve_length++;
218: };
219: return (curve_length);
220: }
221:
222:
223: /*
224: * determine the number of curves in plot->points, separated by
225: * UNDEFINED points
226: */
227:
228: static int num_curves(plot)
229: struct curve_points *plot;
230: {
231: int curves;
232: int first_point;
233: int num_points;
234:
235: first_point = 0;
236: curves = 0;
237: while ((num_points = next_curve(plot, &first_point)) > 0) {
238: curves++;
239: first_point += num_points;
240: }
241: return (curves);
242: }
243:
244:
245:
246: /*
247: * build up a cntr_struct list from curve_points
248: * this funtion is only used for the alternate entry point to
249: * Gershon's code and thus commented out
250: ***deleted***
251: */
252:
253:
1.1.1.2 maekawa 254: /* HBB 990205: rewrote the 'bezier' interpolation routine,
255: * to prevent numerical overflow and other undesirable things happening
256: * for large data files (num_data about 1000 or so), where binomial
257: * coefficients would explode, and powers of 'sr' (0 < sr < 1) become
258: * extremely small. Method used: compute logarithms of these
259: * extremely large and small numbers, and only go back to the
260: * real numbers once they've cancelled out each other, leaving
261: * a reasonable-sized one. */
262:
1.1 maekawa 263: /*
264: * cp_binomial() computes the binomial coefficients needed for BEZIER stuff
265: * and stores them into an array which is hooked to sdat.
266: * (MGR 1992)
267: */
268: static double *cp_binomial(points)
269: int points;
270: {
271: register double *coeff;
272: register int n, k;
273: int e;
274:
275: e = points; /* well we're going from k=0 to k=p_count-1 */
276: coeff = (double *) gp_alloc(e * sizeof(double), "bezier coefficients");
277:
278: n = points - 1;
279: e = n / 2;
1.1.1.2 maekawa 280: /* HBB 990205: calculate these in 'logarithmic space',
281: * as they become _very_ large, with growing n (4^n) */
282: coeff[0] = 0.0;
1.1 maekawa 283:
284: for (k = 0; k < e; k++) {
1.1.1.2 maekawa 285: coeff[k + 1] = coeff[k] + log(((double) (n-k)) / ((double) (k+1)));
1.1 maekawa 286: }
287:
288: for (k = n; k >= e; k--)
289: coeff[k] = coeff[n - k];
290:
291: return (coeff);
292: }
293:
294:
295: /* This is a subfunction of do_bezier() for BEZIER style computations.
296: * It is passed the stepration (STEP/MAXSTEPS) and the addresses of
297: * the double values holding the next x and y coordinates.
298: * (MGR 1992)
299: */
300:
301: /*
302: * well, this routine runs faster with the 68040 striptease FPU
303: * and it handles zeroes correctly - there had been some trouble with TC
304: */
305:
306: GP_INLINE static double s_pow(base, exponent)
307: double base;
308: unsigned int exponent;
309: {
310: double y;
311:
312: if (exponent == 0)
313: return (1.0);
314: if (base == 0.0)
315: return (0.0);
316:
317: /* consider i in binary = abcd
318: * x^i = x^(8a+4b+2c+d) = x^8a * x^4b * x^2b * x^d
319: * = a?x^2^2^2:1 * b?x^2^2:1 + ...
320: * so for each bit in exponent, square x, multiplying it into accumulator
321: *
322: */
323:
324: y = 1;
325: while (exponent) {
326: if (exponent & 1)
327: y *= base;
328: base *= base;
329: /* if exponent was signed, this could be trouble ! */
330: exponent >>= 1;
331: }
332:
333: return (y);
334: }
335:
1.1.1.2 maekawa 336:
1.1 maekawa 337: static void eval_bezier(cp, first_point, num_points, sr, px, py, c)
338: struct curve_points *cp;
339: int first_point; /* where to start in plot->points (to find x-range) */
340: int num_points; /* to determine end in plot->points */
341: double sr;
342: coordval *px;
343: coordval *py;
344: double *c;
345: {
346: unsigned int n = num_points - 1;
347: /* HBB 980308: added 'GPHUGE' tag for DOS */
348: struct coordinate GPHUGE *this_points;
349:
350: this_points = (cp->points) + first_point;
351:
352: if (sr == 0.0) {
353: *px = this_points[0].x;
354: *py = this_points[0].y;
355: } else if (sr == 1.0) {
356: *px = this_points[n].x;
357: *py = this_points[n].y;
358: } else {
1.1.1.2 maekawa 359: /* HBB 990205: do calculation in 'logarithmic space',
360: * to avoid over/underflow errors, which would exactly cancel
361: * out each other, anyway, in an exact calculation
362: */
1.1 maekawa 363: unsigned int i;
364: double lx = 0.0, ly = 0.0;
1.1.1.2 maekawa 365: double log_dsr_to_the_n = n * log(1 -sr);
366: double log_sr_over_dsr = log(sr) - log(1 - sr);
1.1 maekawa 367:
368: for (i = 0; i <= n; i++) {
1.1.1.2 maekawa 369: double u = exp(c[i] + log_dsr_to_the_n + i * log_sr_over_dsr);
370:
1.1 maekawa 371: lx += this_points[i].x * u;
372: ly += this_points[i].y * u;
373: }
374:
375: *px = lx;
376: *py = ly;
377: }
378: }
379:
380: /*
381: * generate a new set of coordinates representing the bezier curve and
382: * set it to the plot
383: */
384:
385: static void do_bezier(cp, bc, first_point, num_points, dest)
386: struct curve_points *cp;
387: double *bc;
388: int first_point; /* where to start in plot->points */
389: int num_points; /* to determine end in plot->points */
390: struct coordinate *dest; /* where to put the interpolated data */
391: {
392: int i;
393: coordval x, y;
394: int xaxis = cp->x_axis;
395: int yaxis = cp->y_axis;
396:
397: /* min and max in internal (eg logged) co-ordinates. We update
398: * these, then update the external extrema in user co-ordinates
399: * at the end.
400: */
401:
402: double ixmin, ixmax, iymin, iymax;
403: double sxmin, sxmax, symin, symax; /* starting values of above */
404:
405: if (log_array[xaxis]) {
406: ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
407: ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
408: } else {
409: ixmin = sxmin = min_array[xaxis];
410: ixmax = sxmax = max_array[xaxis];
411: }
412:
413: if (log_array[yaxis]) {
414: iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
415: iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
416: } else {
417: iymin = symin = min_array[yaxis];
418: iymax = symax = max_array[yaxis];
419: }
420:
421: for (i = 0; i < samples; i++) {
422: eval_bezier(cp, first_point, num_points, (double) i / (double) (samples - 1), &x, &y, bc);
423:
424: /* now we have to store the points and adjust the ranges */
425:
426: dest[i].type = INRANGE;
427: STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue);
428: STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP);
429:
430: dest[i].xlow = dest[i].xhigh = dest[i].x;
431: dest[i].ylow = dest[i].yhigh = dest[i].y;
432:
433: dest[i].z = -1;
434: }
435:
436: UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
437: UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
438: UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
439: UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
440: }
441:
442: /*
443: * call contouring routines -- main entry
444: */
445:
446: /*
447: * it should be like this, but it doesn't run. If you find out why,
448: * contact me: mgr@asgard.bo.open.de or Lars Hanke 2:243/4802.22@fidonet
449: *
450: * Well, all this had originally been inside contour.c, so maybe links
451: * to functions and of contour.c are broken.
452: * ***deleted***
453: * end of unused entry point to Gershon's code
454: *
455: */
456:
457: /*
458: * Solve five diagonal linear system equation. The five diagonal matrix is
459: * defined via matrix M, right side is r, and solution X i.e. M * X = R.
460: * Size of system given in n. Return TRUE if solution exist.
461: * G. Engeln-Muellges/ F.Reutter:
462: * "Formelsammlung zur Numerischen Mathematik mit Standard-FORTRAN-Programmen"
463: * ISBN 3-411-01677-9
464: *
465: * / m02 m03 m04 0 0 0 0 . . . \ / x0 \ / r0 \
466: * I m11 m12 m13 m14 0 0 0 . . . I I x1 I I r1 I
467: * I m20 m21 m22 m23 m24 0 0 . . . I * I x2 I = I r2 I
468: * I 0 m30 m31 m32 m33 m34 0 . . . I I x3 I I r3 I
469: * . . . . . . . . . . . .
470: * \ m(n-3)0 m(n-2)1 m(n-1)2 / \x(n-1)/ \r(n-1)/
471: *
472: */
473: static int solve_five_diag(m, r, x, n)
474: five_diag m[];
475: double r[], x[];
476: int n;
477: {
478: int i;
479: five_diag *hv;
480:
481: hv = (five_diag *) gp_alloc((n + 1) * sizeof(five_diag), "five_diag help vars");
482:
483: hv[0][0] = m[0][2];
484: if (hv[0][0] == 0) {
485: free(hv);
486: return FALSE;
487: }
488: hv[0][1] = m[0][3] / hv[0][0];
489: hv[0][2] = m[0][4] / hv[0][0];
490:
491: hv[1][3] = m[1][1];
492: hv[1][0] = m[1][2] - hv[1][3] * hv[0][1];
493: if (hv[1][0] == 0) {
494: free(hv);
495: return FALSE;
496: }
497: hv[1][1] = (m[1][3] - hv[1][3] * hv[0][2]) / hv[1][0];
498: hv[1][2] = m[1][4] / hv[1][0];
499:
500: for (i = 2; i <= n - 1; i++) {
501: hv[i][3] = m[i][1] - m[i][0] * hv[i - 2][1];
502: hv[i][0] = m[i][2] - m[i][0] * hv[i - 2][2] - hv[i][3] * hv[i - 1][1];
503: if (hv[i][0] == 0) {
504: free(hv);
505: return FALSE;
506: }
507: hv[i][1] = (m[i][3] - hv[i][3] * hv[i - 1][2]) / hv[i][0];
508: hv[i][2] = m[i][4] / hv[i][0];
509: }
510:
511: hv[0][4] = 0;
512: hv[1][4] = r[0] / hv[0][0];
513: for (i = 1; i <= n - 1; i++) {
514: hv[i + 1][4] = (r[i] - m[i][0] * hv[i - 1][4] - hv[i][3] * hv[i][4]) / hv[i][0];
515: }
516:
517: x[n - 1] = hv[n][4];
518: x[n - 2] = hv[n - 1][4] - hv[n - 2][1] * x[n - 1];
519: for (i = n - 3; i >= 0; i--)
520: x[i] = hv[i + 1][4] - hv[i][1] * x[i + 1] - hv[i][2] * x[i + 2];
521:
522: free(hv);
523: return TRUE;
524: }
525:
526:
527: /*
528: * Calculation of approximation cubic splines
529: * Input: x[i], y[i], weights z[i]
530: *
531: * Returns matrix of spline coefficients
532: */
533: static spline_coeff *cp_approx_spline(plot, first_point, num_points)
534: struct curve_points *plot;
535: int first_point; /* where to start in plot->points */
536: int num_points; /* to determine end in plot->points */
537: {
538: spline_coeff *sc;
539: five_diag *m;
540: int xaxis = plot->x_axis;
541: int yaxis = plot->y_axis;
542: double *r, *x, *h, *xp, *yp;
543: /* HBB 980308: added 'GPHUGE' tag */
544: struct coordinate GPHUGE *this_points;
545: int i;
546:
547: sc = (spline_coeff *) gp_alloc((num_points) * sizeof(spline_coeff),
548: "spline matrix");
549:
550: if (num_points < 4)
551: int_error("Can't calculate approximation splines, need at least 4 points", NO_CARET);
552:
553: this_points = (plot->points) + first_point;
554:
555: for (i = 0; i <= num_points - 1; i++)
556: if (this_points[i].z <= 0)
557: int_error("Can't calculate approximation splines, all weights have to be > 0", NO_CARET);
558:
559: m = (five_diag *) gp_alloc((num_points - 2) * sizeof(five_diag), "spline help matrix");
560:
561: r = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline right side");
562: x = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
563: h = (double *) gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
564:
565: xp = (double *) gp_alloc((num_points) * sizeof(double), "x pos");
566: yp = (double *) gp_alloc((num_points) * sizeof(double), "y pos");
567:
568: /* KB 981107: With logarithmic axis first convert back to linear scale */
569:
570: if (log_array[xaxis]) {
571: for (i = 0; i <= num_points - 1; i++)
572: xp[i] = exp(this_points[i].x * log_base_array[xaxis]);
573: } else {
574: for (i = 0; i <= num_points - 1; i++)
575: xp[i] = this_points[i].x;
576: }
577: if (log_array[yaxis]) {
578: for (i = 0; i <= num_points - 1; i++)
579: yp[i] = exp(this_points[i].y * log_base_array[yaxis]);
580: } else {
581: for (i = 0; i <= num_points - 1; i++)
582: yp[i] = this_points[i].y;
583: }
584:
585: for (i = 0; i <= num_points - 2; i++)
586: h[i] = xp[i + 1] - xp[i];
587:
588: /* set up the matrix and the vector */
589:
590: for (i = 0; i <= num_points - 3; i++) {
591: r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
592: - (yp[i + 1] - yp[i]) / h[i]);
593:
594: if (i < 2)
595: m[i][0] = 0;
596: else
597: m[i][0] = 6 / this_points[i].z / h[i - 1] / h[i];
598:
599: if (i < 1)
600: m[i][1] = 0;
601: else
602: m[i][1] = h[i] - 6 / this_points[i].z / h[i] * (1 / h[i - 1] + 1 / h[i])
603: - 6 / this_points[i + 1].z / h[i] * (1 / h[i] + 1 / h[i + 1]);
604:
605: m[i][2] = 2 * (h[i] + h[i + 1])
606: + 6 / this_points[i].z / h[i] / h[i]
607: + 6 / this_points[i + 1].z * (1 / h[i] + 1 / h[i + 1]) * (1 / h[i] + 1 / h[i + 1])
608: + 6 / this_points[i + 2].z / h[i + 1] / h[i + 1];
609:
610: if (i > num_points - 4)
611: m[i][3] = 0;
612: else
613: m[i][3] = h[i + 1] - 6 / this_points[i + 1].z / h[i + 1] * (1 / h[i] + 1 / h[i + 1])
614: - 6 / this_points[i + 2].z / h[i + 1] * (1 / h[i + 1] + 1 / h[i + 2]);
615:
616: if (i > num_points - 5)
617: m[i][4] = 0;
618: else
619: m[i][4] = 6 / this_points[i + 2].z / h[i + 1] / h[i + 2];
620: }
621:
622: /* solve the matrix */
623: if (!solve_five_diag(m, r, x, num_points - 2)) {
624: free(h);
625: free(x);
626: free(r);
627: free(m);
628: free(xp);
629: free(yp);
630: int_error("Can't calculate approximation splines", NO_CARET);
631: }
632: sc[0][2] = 0;
633: for (i = 1; i <= num_points - 2; i++)
634: sc[i][2] = x[i - 1];
635: sc[num_points - 1][2] = 0;
636:
637: sc[0][0] = yp[0] + 2 / this_points[0].z / h[0] * (sc[0][2] - sc[1][2]);
638: for (i = 1; i <= num_points - 2; i++)
639: sc[i][0] = yp[i] - 2 / this_points[i].z *
640: (sc[i - 1][2] / h[i - 1]
641: - sc[i][2] * (1 / h[i - 1] + 1 / h[i])
642: + sc[i + 1][2] / h[i]);
643: sc[num_points - 1][0] = yp[num_points - 1]
644: - 2 / this_points[num_points - 1].z / h[num_points - 2]
645: * (sc[num_points - 2][2] - sc[num_points - 1][2]);
646:
647: for (i = 0; i <= num_points - 2; i++) {
648: sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
649: - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
650: sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
651: }
652:
653: free(h);
654: free(x);
655: free(r);
656: free(m);
657: free(xp);
658: free(yp);
659:
660: return (sc);
661: }
662:
663:
664: /*
665: * Calculation of cubic splines
666: *
667: * This can be treated as a special case of approximation cubic splines, with
668: * all weights -> infinity.
669: *
670: * Returns matrix of spline coefficients
671: */
672: static spline_coeff *cp_tridiag(plot, first_point, num_points)
673: struct curve_points *plot;
674: int first_point, num_points;
675: {
676: spline_coeff *sc;
677: tri_diag *m;
678: int xaxis = plot->x_axis;
679: int yaxis = plot->y_axis;
680: double *r, *x, *h, *xp, *yp;
681: /* HBB 980308: added 'GPHUGE' tag */
682: struct coordinate GPHUGE *this_points;
683: int i;
684:
685: if (num_points < 3)
686: int_error("Can't calculate splines, need at least 3 points", NO_CARET);
687:
688: this_points = (plot->points) + first_point;
689:
690: sc = (spline_coeff *) gp_alloc((num_points) * sizeof(spline_coeff), "spline matrix");
691: m = (tri_diag *) gp_alloc((num_points - 2) * sizeof(tri_diag), "spline help matrix");
692:
693: r = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline right side");
694: x = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
695: h = (double *) gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
696:
697: xp = (double *) gp_alloc((num_points) * sizeof(double), "x pos");
698: yp = (double *) gp_alloc((num_points) * sizeof(double), "y pos");
699:
700: /* KB 981107: With logarithmic axis first convert back to linear scale */
701:
702: if (log_array[xaxis]) {
703: for (i = 0; i <= num_points - 1; i++)
704: xp[i] = exp(this_points[i].x * log_base_array[xaxis]);
705: } else {
706: for (i = 0; i <= num_points - 1; i++)
707: xp[i] = this_points[i].x;
708: }
709: if (log_array[yaxis]) {
710: for (i = 0; i <= num_points - 1; i++)
711: yp[i] = exp(this_points[i].y * log_base_array[yaxis]);
712: } else {
713: for (i = 0; i <= num_points - 1; i++)
714: yp[i] = this_points[i].y;
715: }
716:
717: for (i = 0; i <= num_points - 2; i++)
718: h[i] = xp[i + 1] - xp[i];
719:
720: /* set up the matrix and the vector */
721:
722: for (i = 0; i <= num_points - 3; i++) {
723: r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
724: - (yp[i + 1] - yp[i]) / h[i]);
725:
726: if (i < 1)
727: m[i][0] = 0;
728: else
729: m[i][0] = h[i];
730:
731: m[i][1] = 2 * (h[i] + h[i + 1]);
732:
733: if (i > num_points - 4)
734: m[i][2] = 0;
735: else
736: m[i][2] = h[i + 1];
737: }
738:
739: /* solve the matrix */
740: if (!solve_tri_diag(m, r, x, num_points - 2)) {
741: free(h);
742: free(x);
743: free(r);
744: free(m);
745: free(xp);
746: free(yp);
747: int_error("Can't calculate cubic splines", NO_CARET);
748: }
749: sc[0][2] = 0;
750: for (i = 1; i <= num_points - 2; i++)
751: sc[i][2] = x[i - 1];
752: sc[num_points - 1][2] = 0;
753:
754: for (i = 0; i <= num_points - 1; i++)
755: sc[i][0] = yp[i];
756:
757: for (i = 0; i <= num_points - 2; i++) {
758: sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
759: - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
760: sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
761: }
762:
763: free(h);
764: free(x);
765: free(r);
766: free(m);
767: free(xp);
768: free(yp);
769:
770: return (sc);
771: }
772:
773: static void do_cubic(plot, sc, first_point, num_points, dest)
774: struct curve_points *plot; /* still containes old plot->points */
775: spline_coeff *sc; /* generated by cp_tridiag */
776: int first_point; /* where to start in plot->points */
777: int num_points; /* to determine end in plot->points */
778: struct coordinate *dest; /* where to put the interpolated data */
779: {
780: double xdiff, temp, x, y;
781: int i, l;
782: int xaxis = plot->x_axis;
783: int yaxis = plot->y_axis;
784: /* HBB 980308: added 'GPHUGE' tag */
785: struct coordinate GPHUGE *this_points;
786:
787: /* min and max in internal (eg logged) co-ordinates. We update
788: * these, then update the external extrema in user co-ordinates
789: * at the end.
790: */
791:
792: double ixmin, ixmax, iymin, iymax;
793: double sxmin, sxmax, symin, symax; /* starting values of above */
794:
795: if (log_array[xaxis]) {
796: ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
797: ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
798: } else {
799: ixmin = sxmin = min_array[xaxis];
800: ixmax = sxmax = max_array[xaxis];
801: }
802:
803: if (log_array[yaxis]) {
804: iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
805: iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
806: } else {
807: iymin = symin = min_array[yaxis];
808: iymax = symax = max_array[yaxis];
809: }
810:
811:
812: this_points = (plot->points) + first_point;
813:
814: l = 0;
815:
816: xdiff = (this_points[num_points - 1].x - this_points[0].x) / (samples - 1);
817:
818: for (i = 0; i < samples; i++) {
819: x = this_points[0].x + i * xdiff;
820:
821: while ((x >= this_points[l + 1].x) && (l < num_points - 2))
822: l++;
823:
824: /* KB 981107: With logarithmic x axis the values were converted back to linear */
825: /* scale before calculating the coefficients. Use exponential for log x values. */
826:
827: if (log_array[xaxis]) {
828: temp = exp(x * log_base_array[xaxis]) - exp(this_points[l].x * log_base_array[xaxis]);
829: y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
830: } else {
831: temp = x - this_points[l].x;
832: y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
833: }
834: /* With logarithmic y axis, we need to convert from linear to log scale now. */
835: if (log_array[yaxis]) {
836: if (y > 0.)
837: y = log(y) / log_base_array[yaxis];
838: else
839: y = symin - (symax - symin);
840: }
841: dest[i].type = INRANGE;
842: STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue);
843: STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP);
844:
845: dest[i].xlow = dest[i].xhigh = dest[i].x;
846: dest[i].ylow = dest[i].yhigh = dest[i].y;
847:
848: dest[i].z = -1;
849:
850: }
851:
852: UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
853: UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
854: UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
855: UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
856:
857: }
858:
859: /*
860: * This is the main entry point used. As stated in the header, it is fine,
861: * but I'm not too happy with it.
862: */
863:
864: void gen_interp(plot)
865: struct curve_points *plot;
866: {
867:
868: spline_coeff *sc;
869: double *bc;
870: struct coordinate *new_points;
871: int i, curves;
872: int first_point, num_points;
873:
874: curves = num_curves(plot);
875: new_points = (struct coordinate *) gp_alloc((samples + 1) * curves * sizeof(struct coordinate), "interpolation table");
876:
877: first_point = 0;
878: for (i = 0; i < curves; i++) {
879: num_points = next_curve(plot, &first_point);
880: switch (plot->plot_smooth) {
881: case CSPLINES:
882: sc = cp_tridiag(plot, first_point, num_points);
883: do_cubic(plot, sc, first_point, num_points,
884: new_points + i * (samples + 1));
885: free(sc);
886: break;
887: case ACSPLINES:
888: sc = cp_approx_spline(plot, first_point, num_points);
889: do_cubic(plot, sc, first_point, num_points,
890: new_points + i * (samples + 1));
891: free(sc);
892: break;
893:
894: case BEZIER:
895: case SBEZIER:
896: bc = cp_binomial(num_points);
897: do_bezier(plot, bc, first_point, num_points,
898: new_points + i * (samples + 1));
899: free((char *) bc);
900: break;
901: default: /* keep gcc -Wall quiet */
902: ;
903: }
904: new_points[(i + 1) * (samples + 1) - 1].type = UNDEFINED;
905: first_point += num_points;
906: }
907:
908: free(plot->points);
909: plot->points = new_points;
910: plot->p_max = curves * (samples + 1);
911: plot->p_count = plot->p_max - 1;
912:
913: return;
914: }
915:
916: /*
917: * sort_points
918: *
919: * sort data succession for further evaluation by plot_splines, etc.
920: * This routine is mainly introduced for compilers *NOT* supporting the
921: * UNIX qsort() routine. You can then easily replace it by more convenient
922: * stuff for your compiler.
923: * (MGR 1992)
924: */
925:
1.1.1.3 ! ohara 926: int compare_points(argp1, argp2)
! 927: SORTFUNC_ARGS argp1;
! 928: SORTFUNC_ARGS argp2;
1.1 maekawa 929: {
1.1.1.3 ! ohara 930: const struct coordinate *p1 = argp1;
! 931: const struct coordinate *p2 = argp2;
! 932:
1.1 maekawa 933: if (p1->x > p2->x)
934: return (1);
935: if (p1->x < p2->x)
936: return (-1);
937: return (0);
938: }
939:
940: void sort_points(plot)
1.1.1.3 ! ohara 941: struct curve_points *plot;
1.1 maekawa 942: {
943: int first_point, num_points;
944:
945: first_point = 0;
946: while ((num_points = next_curve(plot, &first_point)) > 0) {
947: /* Sort this set of points, does qsort handle 1 point correctly? */
948: qsort((char *) (plot->points + first_point), num_points,
1.1.1.3 ! ohara 949: sizeof(struct coordinate), compare_points);
1.1 maekawa 950: first_point += num_points;
951: }
952: return;
953: }
954:
955:
956: /*
957: * cp_implode() if averaging is selected this function computes the new
958: * entries and shortens the whole thing to the necessary
959: * size
960: * MGR Addendum
961: */
962:
963: void cp_implode(cp)
964: struct curve_points *cp;
965: {
966: int first_point, num_points;
967: int i, j, k;
968: double x, y, sux, slx, suy, sly;
969: enum coord_type dot;
970:
971:
972: j = 0;
973: first_point = 0;
974: while ((num_points = next_curve(cp, &first_point)) > 0) {
975: k = 0;
976: for (i = first_point; i < first_point + num_points; i++) {
977: if (!k) {
978: x = cp->points[i].x;
979: y = cp->points[i].y;
980: sux = cp->points[i].xhigh;
981: slx = cp->points[i].xlow;
982: suy = cp->points[i].yhigh;
983: sly = cp->points[i].ylow;
984: dot = INRANGE;
985: if (cp->points[i].type != INRANGE)
986: dot = UNDEFINED; /* This means somthing other than usual *//* just signal to check if INRANGE */
987: k = 1;
988: } else if (cp->points[i].x == x) {
989: y += cp->points[i].y;
990: sux += cp->points[i].xhigh;
991: slx += cp->points[i].xlow;
992: suy += cp->points[i].yhigh;
993: sly += cp->points[i].ylow;
994: if (cp->points[i].type != INRANGE)
995: dot = UNDEFINED;
996: k++;
997: } else {
998: cp->points[j].x = x;
999: cp->points[j].y = y / (double) k;
1000: cp->points[j].xhigh = sux / (double) k;
1001: cp->points[j].xlow = slx / (double) k;
1002: cp->points[j].yhigh = suy / (double) k;
1003: cp->points[j].ylow = sly / (double) k;
1004: cp->points[j].type = INRANGE;
1005: if (dot != INRANGE) {
1006: if ((cp->points[j].x > xmax) || (cp->points[j].x < xmin))
1007: cp->points[j].type = OUTRANGE;
1008: else if (!autoscale_ly) {
1009: if ((cp->points[j].y > ymax) || (cp->points[j].y < ymin))
1010: cp->points[j].type = OUTRANGE;
1011: } else
1012: cp->points[j].type = INRANGE;
1013: }
1014: j++; /* next valid entry */
1015: k = 0; /* to read */
1016: i--; /* from this (-> last after for(;;)) entry */
1017: }
1018: }
1019: if (k) {
1020: cp->points[j].x = x;
1021: cp->points[j].y = y / (double) k;
1022: cp->points[j].xhigh = sux / (double) k;
1023: cp->points[j].xlow = slx / (double) k;
1024: cp->points[j].yhigh = suy / (double) k;
1025: cp->points[j].ylow = sly / (double) k;
1026: cp->points[j].type = INRANGE;
1027: if (dot != INRANGE) {
1028: if ((cp->points[j].x > xmax) || (cp->points[j].x < xmin))
1029: cp->points[j].type = OUTRANGE;
1030: else if (!autoscale_ly) {
1031: if ((cp->points[j].y > ymax) || (cp->points[j].y < ymin))
1032: cp->points[j].type = OUTRANGE;
1033: } else
1034: cp->points[j].type = INRANGE;
1035: }
1036: j++; /* next valid entry */
1037: }
1038: /* insert invalid point to separate curves */
1039: if (j < cp->p_count) {
1040: cp->points[j].type = UNDEFINED;
1041: j++;
1042: }
1043: first_point += num_points;
1044: } /* end while */
1045: cp->p_count = j;
1046: cp_extend(cp, j);
1047: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to interpol.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gnuplot