Annotation of OpenXM_contrib/gnuplot/interpol.c, Revision 1.1.1.1
1.1 maekawa 1: #ifndef lint
2: static char *RCSid = "$Id: interpol.c,v 1.29 1998/04/14 00:15:45 drd Exp $";
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));
191: static int compare_points __PROTO((struct coordinate * p1, struct coordinate * p2));
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:
254: /*
255: * cp_binomial() computes the binomial coefficients needed for BEZIER stuff
256: * and stores them into an array which is hooked to sdat.
257: * (MGR 1992)
258: */
259: static double *cp_binomial(points)
260: int points;
261: {
262: register double *coeff;
263: register int n, k;
264: int e;
265:
266: e = points; /* well we're going from k=0 to k=p_count-1 */
267: coeff = (double *) gp_alloc(e * sizeof(double), "bezier coefficients");
268:
269: n = points - 1;
270: e = n / 2;
271: coeff[0] = 1.0;
272:
273: for (k = 0; k < e; k++) {
274: coeff[k + 1] = coeff[k] * ((double) (n - k)) / ((double) (k + 1));
275: }
276:
277: for (k = n; k >= e; k--)
278: coeff[k] = coeff[n - k];
279:
280: return (coeff);
281: }
282:
283:
284: /* This is a subfunction of do_bezier() for BEZIER style computations.
285: * It is passed the stepration (STEP/MAXSTEPS) and the addresses of
286: * the double values holding the next x and y coordinates.
287: * (MGR 1992)
288: */
289:
290: /*
291: * well, this routine runs faster with the 68040 striptease FPU
292: * and it handles zeroes correctly - there had been some trouble with TC
293: */
294:
295: GP_INLINE static double s_pow(base, exponent)
296: double base;
297: unsigned int exponent;
298: {
299: double y;
300:
301: if (exponent == 0)
302: return (1.0);
303: if (base == 0.0)
304: return (0.0);
305:
306: /* consider i in binary = abcd
307: * x^i = x^(8a+4b+2c+d) = x^8a * x^4b * x^2b * x^d
308: * = a?x^2^2^2:1 * b?x^2^2:1 + ...
309: * so for each bit in exponent, square x, multiplying it into accumulator
310: *
311: */
312:
313: y = 1;
314: while (exponent) {
315: if (exponent & 1)
316: y *= base;
317: base *= base;
318: /* if exponent was signed, this could be trouble ! */
319: exponent >>= 1;
320: }
321:
322: return (y);
323: }
324:
325: static void eval_bezier(cp, first_point, num_points, sr, px, py, c)
326: struct curve_points *cp;
327: int first_point; /* where to start in plot->points (to find x-range) */
328: int num_points; /* to determine end in plot->points */
329: double sr;
330: coordval *px;
331: coordval *py;
332: double *c;
333: {
334: unsigned int n = num_points - 1;
335: /* HBB 980308: added 'GPHUGE' tag for DOS */
336: struct coordinate GPHUGE *this_points;
337:
338: this_points = (cp->points) + first_point;
339:
340: if (sr == 0.0) {
341: *px = this_points[0].x;
342: *py = this_points[0].y;
343: } else if (sr == 1.0) {
344: *px = this_points[n].x;
345: *py = this_points[n].y;
346: } else {
347: unsigned int i;
348: double lx = 0.0, ly = 0.0;
349: double sr_to_the_i = 1.0;
350: double dsr_to_the_di = s_pow(1 - sr, n);
351: double reciprocal_dsr = 1.0 / (1 - sr);
352:
353: for (i = 0; i <= n; i++) {
354: double u = c[i] * sr_to_the_i * s_pow(1 - sr, n - i);
355: lx += this_points[i].x * u;
356: ly += this_points[i].y * u;
357: sr_to_the_i *= sr;
358: dsr_to_the_di *= reciprocal_dsr;
359: }
360:
361: *px = lx;
362: *py = ly;
363: }
364: }
365:
366: /*
367: * generate a new set of coordinates representing the bezier curve and
368: * set it to the plot
369: */
370:
371: static void do_bezier(cp, bc, first_point, num_points, dest)
372: struct curve_points *cp;
373: double *bc;
374: int first_point; /* where to start in plot->points */
375: int num_points; /* to determine end in plot->points */
376: struct coordinate *dest; /* where to put the interpolated data */
377: {
378: int i;
379: coordval x, y;
380: int xaxis = cp->x_axis;
381: int yaxis = cp->y_axis;
382:
383: /* min and max in internal (eg logged) co-ordinates. We update
384: * these, then update the external extrema in user co-ordinates
385: * at the end.
386: */
387:
388: double ixmin, ixmax, iymin, iymax;
389: double sxmin, sxmax, symin, symax; /* starting values of above */
390:
391: if (log_array[xaxis]) {
392: ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
393: ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
394: } else {
395: ixmin = sxmin = min_array[xaxis];
396: ixmax = sxmax = max_array[xaxis];
397: }
398:
399: if (log_array[yaxis]) {
400: iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
401: iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
402: } else {
403: iymin = symin = min_array[yaxis];
404: iymax = symax = max_array[yaxis];
405: }
406:
407: for (i = 0; i < samples; i++) {
408: eval_bezier(cp, first_point, num_points, (double) i / (double) (samples - 1), &x, &y, bc);
409:
410: /* now we have to store the points and adjust the ranges */
411:
412: dest[i].type = INRANGE;
413: STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue);
414: STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP);
415:
416: dest[i].xlow = dest[i].xhigh = dest[i].x;
417: dest[i].ylow = dest[i].yhigh = dest[i].y;
418:
419: dest[i].z = -1;
420: }
421:
422: UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
423: UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
424: UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
425: UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
426: }
427:
428: /*
429: * call contouring routines -- main entry
430: */
431:
432: /*
433: * it should be like this, but it doesn't run. If you find out why,
434: * contact me: mgr@asgard.bo.open.de or Lars Hanke 2:243/4802.22@fidonet
435: *
436: * Well, all this had originally been inside contour.c, so maybe links
437: * to functions and of contour.c are broken.
438: * ***deleted***
439: * end of unused entry point to Gershon's code
440: *
441: */
442:
443: /*
444: * Solve five diagonal linear system equation. The five diagonal matrix is
445: * defined via matrix M, right side is r, and solution X i.e. M * X = R.
446: * Size of system given in n. Return TRUE if solution exist.
447: * G. Engeln-Muellges/ F.Reutter:
448: * "Formelsammlung zur Numerischen Mathematik mit Standard-FORTRAN-Programmen"
449: * ISBN 3-411-01677-9
450: *
451: * / m02 m03 m04 0 0 0 0 . . . \ / x0 \ / r0 \
452: * I m11 m12 m13 m14 0 0 0 . . . I I x1 I I r1 I
453: * I m20 m21 m22 m23 m24 0 0 . . . I * I x2 I = I r2 I
454: * I 0 m30 m31 m32 m33 m34 0 . . . I I x3 I I r3 I
455: * . . . . . . . . . . . .
456: * \ m(n-3)0 m(n-2)1 m(n-1)2 / \x(n-1)/ \r(n-1)/
457: *
458: */
459: static int solve_five_diag(m, r, x, n)
460: five_diag m[];
461: double r[], x[];
462: int n;
463: {
464: int i;
465: five_diag *hv;
466:
467: hv = (five_diag *) gp_alloc((n + 1) * sizeof(five_diag), "five_diag help vars");
468:
469: hv[0][0] = m[0][2];
470: if (hv[0][0] == 0) {
471: free(hv);
472: return FALSE;
473: }
474: hv[0][1] = m[0][3] / hv[0][0];
475: hv[0][2] = m[0][4] / hv[0][0];
476:
477: hv[1][3] = m[1][1];
478: hv[1][0] = m[1][2] - hv[1][3] * hv[0][1];
479: if (hv[1][0] == 0) {
480: free(hv);
481: return FALSE;
482: }
483: hv[1][1] = (m[1][3] - hv[1][3] * hv[0][2]) / hv[1][0];
484: hv[1][2] = m[1][4] / hv[1][0];
485:
486: for (i = 2; i <= n - 1; i++) {
487: hv[i][3] = m[i][1] - m[i][0] * hv[i - 2][1];
488: hv[i][0] = m[i][2] - m[i][0] * hv[i - 2][2] - hv[i][3] * hv[i - 1][1];
489: if (hv[i][0] == 0) {
490: free(hv);
491: return FALSE;
492: }
493: hv[i][1] = (m[i][3] - hv[i][3] * hv[i - 1][2]) / hv[i][0];
494: hv[i][2] = m[i][4] / hv[i][0];
495: }
496:
497: hv[0][4] = 0;
498: hv[1][4] = r[0] / hv[0][0];
499: for (i = 1; i <= n - 1; i++) {
500: hv[i + 1][4] = (r[i] - m[i][0] * hv[i - 1][4] - hv[i][3] * hv[i][4]) / hv[i][0];
501: }
502:
503: x[n - 1] = hv[n][4];
504: x[n - 2] = hv[n - 1][4] - hv[n - 2][1] * x[n - 1];
505: for (i = n - 3; i >= 0; i--)
506: x[i] = hv[i + 1][4] - hv[i][1] * x[i + 1] - hv[i][2] * x[i + 2];
507:
508: free(hv);
509: return TRUE;
510: }
511:
512:
513: /*
514: * Calculation of approximation cubic splines
515: * Input: x[i], y[i], weights z[i]
516: *
517: * Returns matrix of spline coefficients
518: */
519: static spline_coeff *cp_approx_spline(plot, first_point, num_points)
520: struct curve_points *plot;
521: int first_point; /* where to start in plot->points */
522: int num_points; /* to determine end in plot->points */
523: {
524: spline_coeff *sc;
525: five_diag *m;
526: int xaxis = plot->x_axis;
527: int yaxis = plot->y_axis;
528: double *r, *x, *h, *xp, *yp;
529: /* HBB 980308: added 'GPHUGE' tag */
530: struct coordinate GPHUGE *this_points;
531: int i;
532:
533: sc = (spline_coeff *) gp_alloc((num_points) * sizeof(spline_coeff),
534: "spline matrix");
535:
536: if (num_points < 4)
537: int_error("Can't calculate approximation splines, need at least 4 points", NO_CARET);
538:
539: this_points = (plot->points) + first_point;
540:
541: for (i = 0; i <= num_points - 1; i++)
542: if (this_points[i].z <= 0)
543: int_error("Can't calculate approximation splines, all weights have to be > 0", NO_CARET);
544:
545: m = (five_diag *) gp_alloc((num_points - 2) * sizeof(five_diag), "spline help matrix");
546:
547: r = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline right side");
548: x = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
549: h = (double *) gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
550:
551: xp = (double *) gp_alloc((num_points) * sizeof(double), "x pos");
552: yp = (double *) gp_alloc((num_points) * sizeof(double), "y pos");
553:
554: /* KB 981107: With logarithmic axis first convert back to linear scale */
555:
556: if (log_array[xaxis]) {
557: for (i = 0; i <= num_points - 1; i++)
558: xp[i] = exp(this_points[i].x * log_base_array[xaxis]);
559: } else {
560: for (i = 0; i <= num_points - 1; i++)
561: xp[i] = this_points[i].x;
562: }
563: if (log_array[yaxis]) {
564: for (i = 0; i <= num_points - 1; i++)
565: yp[i] = exp(this_points[i].y * log_base_array[yaxis]);
566: } else {
567: for (i = 0; i <= num_points - 1; i++)
568: yp[i] = this_points[i].y;
569: }
570:
571: for (i = 0; i <= num_points - 2; i++)
572: h[i] = xp[i + 1] - xp[i];
573:
574: /* set up the matrix and the vector */
575:
576: for (i = 0; i <= num_points - 3; i++) {
577: r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
578: - (yp[i + 1] - yp[i]) / h[i]);
579:
580: if (i < 2)
581: m[i][0] = 0;
582: else
583: m[i][0] = 6 / this_points[i].z / h[i - 1] / h[i];
584:
585: if (i < 1)
586: m[i][1] = 0;
587: else
588: m[i][1] = h[i] - 6 / this_points[i].z / h[i] * (1 / h[i - 1] + 1 / h[i])
589: - 6 / this_points[i + 1].z / h[i] * (1 / h[i] + 1 / h[i + 1]);
590:
591: m[i][2] = 2 * (h[i] + h[i + 1])
592: + 6 / this_points[i].z / h[i] / h[i]
593: + 6 / this_points[i + 1].z * (1 / h[i] + 1 / h[i + 1]) * (1 / h[i] + 1 / h[i + 1])
594: + 6 / this_points[i + 2].z / h[i + 1] / h[i + 1];
595:
596: if (i > num_points - 4)
597: m[i][3] = 0;
598: else
599: m[i][3] = h[i + 1] - 6 / this_points[i + 1].z / h[i + 1] * (1 / h[i] + 1 / h[i + 1])
600: - 6 / this_points[i + 2].z / h[i + 1] * (1 / h[i + 1] + 1 / h[i + 2]);
601:
602: if (i > num_points - 5)
603: m[i][4] = 0;
604: else
605: m[i][4] = 6 / this_points[i + 2].z / h[i + 1] / h[i + 2];
606: }
607:
608: /* solve the matrix */
609: if (!solve_five_diag(m, r, x, num_points - 2)) {
610: free(h);
611: free(x);
612: free(r);
613: free(m);
614: free(xp);
615: free(yp);
616: int_error("Can't calculate approximation splines", NO_CARET);
617: }
618: sc[0][2] = 0;
619: for (i = 1; i <= num_points - 2; i++)
620: sc[i][2] = x[i - 1];
621: sc[num_points - 1][2] = 0;
622:
623: sc[0][0] = yp[0] + 2 / this_points[0].z / h[0] * (sc[0][2] - sc[1][2]);
624: for (i = 1; i <= num_points - 2; i++)
625: sc[i][0] = yp[i] - 2 / this_points[i].z *
626: (sc[i - 1][2] / h[i - 1]
627: - sc[i][2] * (1 / h[i - 1] + 1 / h[i])
628: + sc[i + 1][2] / h[i]);
629: sc[num_points - 1][0] = yp[num_points - 1]
630: - 2 / this_points[num_points - 1].z / h[num_points - 2]
631: * (sc[num_points - 2][2] - sc[num_points - 1][2]);
632:
633: for (i = 0; i <= num_points - 2; i++) {
634: sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
635: - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
636: sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
637: }
638:
639: free(h);
640: free(x);
641: free(r);
642: free(m);
643: free(xp);
644: free(yp);
645:
646: return (sc);
647: }
648:
649:
650: /*
651: * Calculation of cubic splines
652: *
653: * This can be treated as a special case of approximation cubic splines, with
654: * all weights -> infinity.
655: *
656: * Returns matrix of spline coefficients
657: */
658: static spline_coeff *cp_tridiag(plot, first_point, num_points)
659: struct curve_points *plot;
660: int first_point, num_points;
661: {
662: spline_coeff *sc;
663: tri_diag *m;
664: int xaxis = plot->x_axis;
665: int yaxis = plot->y_axis;
666: double *r, *x, *h, *xp, *yp;
667: /* HBB 980308: added 'GPHUGE' tag */
668: struct coordinate GPHUGE *this_points;
669: int i;
670:
671: if (num_points < 3)
672: int_error("Can't calculate splines, need at least 3 points", NO_CARET);
673:
674: this_points = (plot->points) + first_point;
675:
676: sc = (spline_coeff *) gp_alloc((num_points) * sizeof(spline_coeff), "spline matrix");
677: m = (tri_diag *) gp_alloc((num_points - 2) * sizeof(tri_diag), "spline help matrix");
678:
679: r = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline right side");
680: x = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
681: h = (double *) gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
682:
683: xp = (double *) gp_alloc((num_points) * sizeof(double), "x pos");
684: yp = (double *) gp_alloc((num_points) * sizeof(double), "y pos");
685:
686: /* KB 981107: With logarithmic axis first convert back to linear scale */
687:
688: if (log_array[xaxis]) {
689: for (i = 0; i <= num_points - 1; i++)
690: xp[i] = exp(this_points[i].x * log_base_array[xaxis]);
691: } else {
692: for (i = 0; i <= num_points - 1; i++)
693: xp[i] = this_points[i].x;
694: }
695: if (log_array[yaxis]) {
696: for (i = 0; i <= num_points - 1; i++)
697: yp[i] = exp(this_points[i].y * log_base_array[yaxis]);
698: } else {
699: for (i = 0; i <= num_points - 1; i++)
700: yp[i] = this_points[i].y;
701: }
702:
703: for (i = 0; i <= num_points - 2; i++)
704: h[i] = xp[i + 1] - xp[i];
705:
706: /* set up the matrix and the vector */
707:
708: for (i = 0; i <= num_points - 3; i++) {
709: r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
710: - (yp[i + 1] - yp[i]) / h[i]);
711:
712: if (i < 1)
713: m[i][0] = 0;
714: else
715: m[i][0] = h[i];
716:
717: m[i][1] = 2 * (h[i] + h[i + 1]);
718:
719: if (i > num_points - 4)
720: m[i][2] = 0;
721: else
722: m[i][2] = h[i + 1];
723: }
724:
725: /* solve the matrix */
726: if (!solve_tri_diag(m, r, x, num_points - 2)) {
727: free(h);
728: free(x);
729: free(r);
730: free(m);
731: free(xp);
732: free(yp);
733: int_error("Can't calculate cubic splines", NO_CARET);
734: }
735: sc[0][2] = 0;
736: for (i = 1; i <= num_points - 2; i++)
737: sc[i][2] = x[i - 1];
738: sc[num_points - 1][2] = 0;
739:
740: for (i = 0; i <= num_points - 1; i++)
741: sc[i][0] = yp[i];
742:
743: for (i = 0; i <= num_points - 2; i++) {
744: sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
745: - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
746: sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
747: }
748:
749: free(h);
750: free(x);
751: free(r);
752: free(m);
753: free(xp);
754: free(yp);
755:
756: return (sc);
757: }
758:
759: static void do_cubic(plot, sc, first_point, num_points, dest)
760: struct curve_points *plot; /* still containes old plot->points */
761: spline_coeff *sc; /* generated by cp_tridiag */
762: int first_point; /* where to start in plot->points */
763: int num_points; /* to determine end in plot->points */
764: struct coordinate *dest; /* where to put the interpolated data */
765: {
766: double xdiff, temp, x, y;
767: int i, l;
768: int xaxis = plot->x_axis;
769: int yaxis = plot->y_axis;
770: /* HBB 980308: added 'GPHUGE' tag */
771: struct coordinate GPHUGE *this_points;
772:
773: /* min and max in internal (eg logged) co-ordinates. We update
774: * these, then update the external extrema in user co-ordinates
775: * at the end.
776: */
777:
778: double ixmin, ixmax, iymin, iymax;
779: double sxmin, sxmax, symin, symax; /* starting values of above */
780:
781: if (log_array[xaxis]) {
782: ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
783: ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
784: } else {
785: ixmin = sxmin = min_array[xaxis];
786: ixmax = sxmax = max_array[xaxis];
787: }
788:
789: if (log_array[yaxis]) {
790: iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
791: iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
792: } else {
793: iymin = symin = min_array[yaxis];
794: iymax = symax = max_array[yaxis];
795: }
796:
797:
798: this_points = (plot->points) + first_point;
799:
800: l = 0;
801:
802: xdiff = (this_points[num_points - 1].x - this_points[0].x) / (samples - 1);
803:
804: for (i = 0; i < samples; i++) {
805: x = this_points[0].x + i * xdiff;
806:
807: while ((x >= this_points[l + 1].x) && (l < num_points - 2))
808: l++;
809:
810: /* KB 981107: With logarithmic x axis the values were converted back to linear */
811: /* scale before calculating the coefficients. Use exponential for log x values. */
812:
813: if (log_array[xaxis]) {
814: temp = exp(x * log_base_array[xaxis]) - exp(this_points[l].x * log_base_array[xaxis]);
815: y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
816: } else {
817: temp = x - this_points[l].x;
818: y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
819: }
820: /* With logarithmic y axis, we need to convert from linear to log scale now. */
821: if (log_array[yaxis]) {
822: if (y > 0.)
823: y = log(y) / log_base_array[yaxis];
824: else
825: y = symin - (symax - symin);
826: }
827: dest[i].type = INRANGE;
828: STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue);
829: STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP);
830:
831: dest[i].xlow = dest[i].xhigh = dest[i].x;
832: dest[i].ylow = dest[i].yhigh = dest[i].y;
833:
834: dest[i].z = -1;
835:
836: }
837:
838: UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
839: UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
840: UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
841: UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
842:
843: }
844:
845: /*
846: * This is the main entry point used. As stated in the header, it is fine,
847: * but I'm not too happy with it.
848: */
849:
850: void gen_interp(plot)
851: struct curve_points *plot;
852: {
853:
854: spline_coeff *sc;
855: double *bc;
856: struct coordinate *new_points;
857: int i, curves;
858: int first_point, num_points;
859:
860: curves = num_curves(plot);
861: new_points = (struct coordinate *) gp_alloc((samples + 1) * curves * sizeof(struct coordinate), "interpolation table");
862:
863: first_point = 0;
864: for (i = 0; i < curves; i++) {
865: num_points = next_curve(plot, &first_point);
866: switch (plot->plot_smooth) {
867: case CSPLINES:
868: sc = cp_tridiag(plot, first_point, num_points);
869: do_cubic(plot, sc, first_point, num_points,
870: new_points + i * (samples + 1));
871: free(sc);
872: break;
873: case ACSPLINES:
874: sc = cp_approx_spline(plot, first_point, num_points);
875: do_cubic(plot, sc, first_point, num_points,
876: new_points + i * (samples + 1));
877: free(sc);
878: break;
879:
880: case BEZIER:
881: case SBEZIER:
882: bc = cp_binomial(num_points);
883: do_bezier(plot, bc, first_point, num_points,
884: new_points + i * (samples + 1));
885: free((char *) bc);
886: break;
887: default: /* keep gcc -Wall quiet */
888: ;
889: }
890: new_points[(i + 1) * (samples + 1) - 1].type = UNDEFINED;
891: first_point += num_points;
892: }
893:
894: free(plot->points);
895: plot->points = new_points;
896: plot->p_max = curves * (samples + 1);
897: plot->p_count = plot->p_max - 1;
898:
899: return;
900: }
901:
902: /*
903: * sort_points
904: *
905: * sort data succession for further evaluation by plot_splines, etc.
906: * This routine is mainly introduced for compilers *NOT* supporting the
907: * UNIX qsort() routine. You can then easily replace it by more convenient
908: * stuff for your compiler.
909: * (MGR 1992)
910: */
911:
912: static int compare_points(p1, p2)
913: struct coordinate *p1;
914: struct coordinate *p2;
915: {
916: if (p1->x > p2->x)
917: return (1);
918: if (p1->x < p2->x)
919: return (-1);
920: return (0);
921: }
922:
923: void sort_points(plot)
924: struct curve_points *plot;
925: {
926: int first_point, num_points;
927:
928: first_point = 0;
929: while ((num_points = next_curve(plot, &first_point)) > 0) {
930: /* Sort this set of points, does qsort handle 1 point correctly? */
931: qsort((char *) (plot->points + first_point), num_points,
932: sizeof(struct coordinate), (sortfunc) compare_points);
933: first_point += num_points;
934: }
935: return;
936: }
937:
938:
939: /*
940: * cp_implode() if averaging is selected this function computes the new
941: * entries and shortens the whole thing to the necessary
942: * size
943: * MGR Addendum
944: */
945:
946: void cp_implode(cp)
947: struct curve_points *cp;
948: {
949: int first_point, num_points;
950: int i, j, k;
951: double x, y, sux, slx, suy, sly;
952: enum coord_type dot;
953:
954:
955: j = 0;
956: first_point = 0;
957: while ((num_points = next_curve(cp, &first_point)) > 0) {
958: k = 0;
959: for (i = first_point; i < first_point + num_points; i++) {
960: if (!k) {
961: x = cp->points[i].x;
962: y = cp->points[i].y;
963: sux = cp->points[i].xhigh;
964: slx = cp->points[i].xlow;
965: suy = cp->points[i].yhigh;
966: sly = cp->points[i].ylow;
967: dot = INRANGE;
968: if (cp->points[i].type != INRANGE)
969: dot = UNDEFINED; /* This means somthing other than usual *//* just signal to check if INRANGE */
970: k = 1;
971: } else if (cp->points[i].x == x) {
972: y += cp->points[i].y;
973: sux += cp->points[i].xhigh;
974: slx += cp->points[i].xlow;
975: suy += cp->points[i].yhigh;
976: sly += cp->points[i].ylow;
977: if (cp->points[i].type != INRANGE)
978: dot = UNDEFINED;
979: k++;
980: } else {
981: cp->points[j].x = x;
982: cp->points[j].y = y / (double) k;
983: cp->points[j].xhigh = sux / (double) k;
984: cp->points[j].xlow = slx / (double) k;
985: cp->points[j].yhigh = suy / (double) k;
986: cp->points[j].ylow = sly / (double) k;
987: cp->points[j].type = INRANGE;
988: if (dot != INRANGE) {
989: if ((cp->points[j].x > xmax) || (cp->points[j].x < xmin))
990: cp->points[j].type = OUTRANGE;
991: else if (!autoscale_ly) {
992: if ((cp->points[j].y > ymax) || (cp->points[j].y < ymin))
993: cp->points[j].type = OUTRANGE;
994: } else
995: cp->points[j].type = INRANGE;
996: }
997: j++; /* next valid entry */
998: k = 0; /* to read */
999: i--; /* from this (-> last after for(;;)) entry */
1000: }
1001: }
1002: if (k) {
1003: cp->points[j].x = x;
1004: cp->points[j].y = y / (double) k;
1005: cp->points[j].xhigh = sux / (double) k;
1006: cp->points[j].xlow = slx / (double) k;
1007: cp->points[j].yhigh = suy / (double) k;
1008: cp->points[j].ylow = sly / (double) k;
1009: cp->points[j].type = INRANGE;
1010: if (dot != INRANGE) {
1011: if ((cp->points[j].x > xmax) || (cp->points[j].x < xmin))
1012: cp->points[j].type = OUTRANGE;
1013: else if (!autoscale_ly) {
1014: if ((cp->points[j].y > ymax) || (cp->points[j].y < ymin))
1015: cp->points[j].type = OUTRANGE;
1016: } else
1017: cp->points[j].type = INRANGE;
1018: }
1019: j++; /* next valid entry */
1020: }
1021: /* insert invalid point to separate curves */
1022: if (j < cp->p_count) {
1023: cp->points[j].type = UNDEFINED;
1024: j++;
1025: }
1026: first_point += num_points;
1027: } /* end while */
1028: cp->p_count = j;
1029: cp_extend(cp, j);
1030: }
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