2 static char *RCSid() { return RCSid("$Id: interpol.c,v 1.31 2004/07/25 12:25:01 broeker Exp $"); }
5 /* GNUPLOT - interpol.c */
8 * Copyright 1986 - 1993, 1998, 2004 Thomas Williams, Colin Kelley
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.
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,
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
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.
33 * This software is provided "as is" without express or implied warranty
34 * to the extent permitted by applicable law.
39 * C-Source file identification Header
41 * This file belongs to a project which is:
43 * done 1993 by MGR-Software, Asgard (Lars Hanke)
44 * written by Lars Hanke
48 * InterNet: mgr@asgard.bo.open.de
49 * FIDO: Lars Hanke @ 2:243/4802.22 (as long as they keep addresses)
51 **************************************************************************
61 **************************************************************************
64 * This module is part of gnuplot and distributed under whatever terms
65 * gnuplot is or will be published, unless exclusive rights are claimed.
68 * Supplies 2-D data interpolation and approximation routines
73 * - structs: curve_points, coordval, coordinate
76 * - samples, axis array[] variables
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.
94 **************************************************************************
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
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
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_1 is typically 100 ; we
110 * dont want to take more logs than necessary)
111 * Also, need to take into account which axes are active
113 * Jun 30, 1996 Jens Emmerich
114 * implemented handling of UNDEFINED points
117 #include "interpol.h"
122 #include "graphics.h"
125 /* #include "setshow.h" */
129 /* in order to support multiple axes, and to simplify ranging in
130 * parametric plots, we use arrays to store some things. For 2d plots,
131 * elements are z=0,y1=1,x1=2,z2=4,y2=5,x2=6 these are given symbolic
137 * IMHO, code is getting too cluttered with repeated chunks of
138 * code. Some macros to simplify, I hope.
140 * do { } while(0) is comp.lang.c recommendation for complex macros
141 * also means that break can be specified as an action, and it will
146 /* store VALUE or log(VALUE) in STORE, set TYPE as appropriate Do
147 * OUT_ACTION or UNDEF_ACTION as appropriate. Adjust range provided
148 * type is INRANGE (ie dont adjust y if x is outrange). VALUE must not
149 * be same as STORE */
150 /* FIXME 20010610: UNDEF_ACTION is completely unused ??? Furthermore,
151 * this is so similar to STORE_WITH_LOG_AND_UPDATE_RANGE() from axis.h
152 * that the two should probably be merged. */
153 #define STORE_AND_FIXUP_RANGE(store, value, type, min, max, auto, \
154 out_action, undef_action) \
157 if (type != INRANGE) \
158 break; /* don't set y range if x is outrange, for example */ \
159 if ((value) < (min)) { \
160 if ((auto) & AUTOSCALE_MIN) \
168 if ((value) > (max)) { \
169 if ((auto) & AUTOSCALE_MAX) \
178 #define UPDATE_RANGE(TEST,OLD,NEW,AXIS) \
181 (OLD) = AXIS_DE_LOG_VALUE(AXIS,NEW); \
184 #define spline_coeff_size 4
185 typedef double spline_coeff[spline_coeff_size];
186 typedef double five_diag[5];
188 static int next_curve __PROTO((struct curve_points * plot, int *curve_start));
189 static int num_curves __PROTO((struct curve_points * plot));
190 static double *cp_binomial __PROTO((int points));
191 static void eval_bezier __PROTO((struct curve_points * cp, int first_point,
192 int num_points, double sr, coordval * px,
193 coordval *py, double *c));
194 static void do_bezier __PROTO((struct curve_points * cp, double *bc, int first_point, int num_points, struct coordinate * dest));
195 static int solve_five_diag __PROTO((five_diag m[], double r[], double x[], int n));
196 static spline_coeff *cp_approx_spline __PROTO((struct curve_points * plot, int first_point, int num_points));
197 static spline_coeff *cp_tridiag __PROTO((struct curve_points * plot, int first_point, int num_points));
198 static void do_cubic __PROTO((struct curve_points * plot, spline_coeff * sc, int first_point, int num_points, struct coordinate * dest));
199 static void do_freq __PROTO((struct curve_points *plot, int first_point, int num_points));
200 int compare_points __PROTO((SORTFUNC_ARGS p1, SORTFUNC_ARGS p2));
204 * position curve_start to index the next non-UNDEFINDED point,
205 * start search at initial curve_start,
206 * return number of non-UNDEFINDED points from there on,
207 * if no more valid points are found, curve_start is set
208 * to plot->p_count and 0 is returned
212 next_curve(struct curve_points *plot, int *curve_start)
216 /* Skip undefined points */
217 while (*curve_start < plot->p_count
218 && plot->points[*curve_start].type == UNDEFINED) {
222 /* curve_length is first used as an offset, then the correct # points */
223 while ((*curve_start) + curve_length < plot->p_count
224 && plot->points[(*curve_start) + curve_length].type != UNDEFINED) {
227 return (curve_length);
232 * determine the number of curves in plot->points, separated by
237 num_curves(struct curve_points *plot)
245 while ((num_points = next_curve(plot, &first_point)) > 0) {
247 first_point += num_points;
255 * build up a cntr_struct list from curve_points
256 * this funtion is only used for the alternate entry point to
257 * Gershon's code and thus commented out
262 /* HBB 990205: rewrote the 'bezier' interpolation routine,
263 * to prevent numerical overflow and other undesirable things happening
264 * for large data files (num_data about 1000 or so), where binomial
265 * coefficients would explode, and powers of 'sr' (0 < sr < 1) become
266 * extremely small. Method used: compute logarithms of these
267 * extremely large and small numbers, and only go back to the
268 * real numbers once they've cancelled out each other, leaving
269 * a reasonable-sized one. */
272 * cp_binomial() computes the binomial coefficients needed for BEZIER stuff
273 * and stores them into an array which is hooked to sdat.
277 cp_binomial(int points)
283 e = points; /* well we're going from k=0 to k=p_count-1 */
284 coeff = gp_alloc(e * sizeof(double), "bezier coefficients");
288 /* HBB 990205: calculate these in 'logarithmic space',
289 * as they become _very_ large, with growing n (4^n) */
292 for (k = 0; k < e; k++) {
293 coeff[k + 1] = coeff[k] + log(((double) (n - k)) / ((double) (k + 1)));
296 for (k = n; k >= e; k--)
297 coeff[k] = coeff[n - k];
303 /* This is a subfunction of do_bezier() for BEZIER style computations.
304 * It is passed the stepration (STEP/MAXSTEPS) and the addresses of
305 * the double values holding the next x and y coordinates.
311 struct curve_points *cp,
312 int first_point, /* where to start in plot->points (to find x-range) */
313 int num_points, /* to determine end in plot->points */
314 double sr, /* position inside curve, range [0:1] */
315 coordval *px, /* OUTPUT: x and y */
317 double *c) /* Bezier coefficient array */
319 unsigned int n = num_points - 1;
320 struct coordinate GPHUGE *this_points;
322 this_points = (cp->points) + first_point;
325 *px = this_points[0].x;
326 *py = this_points[0].y;
327 } else if (sr == 1.0) {
328 *px = this_points[n].x;
329 *py = this_points[n].y;
331 /* HBB 990205: do calculation in 'logarithmic space',
332 * to avoid over/underflow errors, which would exactly cancel
333 * out each other, anyway, in an exact calculation
336 double lx = 0.0, ly = 0.0;
337 double log_dsr_to_the_n = n * log(1 - sr);
338 double log_sr_over_dsr = log(sr) - log(1 - sr);
340 for (i = 0; i <= n; i++) {
341 double u = exp(c[i] + log_dsr_to_the_n + i * log_sr_over_dsr);
343 lx += this_points[i].x * u;
344 ly += this_points[i].y * u;
353 * generate a new set of coordinates representing the bezier curve and
359 struct curve_points *cp,
360 double *bc, /* Bezier coefficient array */
361 int first_point, /* where to start in plot->points */
362 int num_points, /* to determine end in plot->points */
363 struct coordinate *dest) /* where to put the interpolated data */
368 /* min and max in internal (eg logged) co-ordinates. We update
369 * these, then update the external extrema in user co-ordinates
373 double ixmin, ixmax, iymin, iymax;
374 double sxmin, sxmax, symin, symax; /* starting values of above */
379 ixmin = sxmin = AXIS_LOG_VALUE(x_axis, X_AXIS.min);
380 ixmax = sxmax = AXIS_LOG_VALUE(x_axis, X_AXIS.max);
381 iymin = symin = AXIS_LOG_VALUE(y_axis, Y_AXIS.min);
382 iymax = symax = AXIS_LOG_VALUE(y_axis, Y_AXIS.max);
384 for (i = 0; i < samples_1; i++) {
385 eval_bezier(cp, first_point, num_points,
386 (double) i / (double) (samples_1 - 1),
389 /* now we have to store the points and adjust the ranges */
391 dest[i].type = INRANGE;
392 STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, X_AXIS.autoscale, NOOP, continue);
393 STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, Y_AXIS.autoscale, NOOP, NOOP);
395 dest[i].xlow = dest[i].xhigh = dest[i].x;
396 dest[i].ylow = dest[i].yhigh = dest[i].y;
401 UPDATE_RANGE(ixmax > sxmax, X_AXIS.max, ixmax, x_axis);
402 UPDATE_RANGE(ixmin < sxmin, X_AXIS.min, ixmin, x_axis);
403 UPDATE_RANGE(iymax > symax, Y_AXIS.max, iymax, y_axis);
404 UPDATE_RANGE(iymin < symin, Y_AXIS.min, iymin, y_axis);
408 * call contouring routines -- main entry
412 * it should be like this, but it doesn't run. If you find out why,
413 * contact me: mgr@asgard.bo.open.de or Lars Hanke 2:243/4802.22@fidonet
415 * Well, all this had originally been inside contour.c, so maybe links
416 * to functions and of contour.c are broken.
418 * end of unused entry point to Gershon's code
423 * Solve five diagonal linear system equation. The five diagonal matrix is
424 * defined via matrix M, right side is r, and solution X i.e. M * X = R.
425 * Size of system given in n. Return TRUE if solution exist.
426 * G. Engeln-Muellges/ F.Reutter:
427 * "Formelsammlung zur Numerischen Mathematik mit Standard-FORTRAN-Programmen"
430 * / m02 m03 m04 0 0 0 0 . . . \ / x0 \ / r0 \
431 * I m11 m12 m13 m14 0 0 0 . . . I I x1 I I r1 I
432 * I m20 m21 m22 m23 m24 0 0 . . . I * I x2 I = I r2 I
433 * I 0 m30 m31 m32 m33 m34 0 . . . I I x3 I I r3 I
434 * . . . . . . . . . . . .
435 * \ m(n-3)0 m(n-2)1 m(n-1)2 / \x(n-1)/ \r(n-1)/
439 solve_five_diag(five_diag m[], double r[], double x[], int n)
444 hv = gp_alloc((n + 1) * sizeof(five_diag), "five_diag help vars");
451 hv[0][1] = m[0][3] / hv[0][0];
452 hv[0][2] = m[0][4] / hv[0][0];
455 hv[1][0] = m[1][2] - hv[1][3] * hv[0][1];
460 hv[1][1] = (m[1][3] - hv[1][3] * hv[0][2]) / hv[1][0];
461 hv[1][2] = m[1][4] / hv[1][0];
463 for (i = 2; i < n; i++) {
464 hv[i][3] = m[i][1] - m[i][0] * hv[i - 2][1];
465 hv[i][0] = m[i][2] - m[i][0] * hv[i - 2][2] - hv[i][3] * hv[i - 1][1];
470 hv[i][1] = (m[i][3] - hv[i][3] * hv[i - 1][2]) / hv[i][0];
471 hv[i][2] = m[i][4] / hv[i][0];
475 hv[1][4] = r[0] / hv[0][0];
476 for (i = 1; i < n; i++) {
477 hv[i + 1][4] = (r[i] - m[i][0] * hv[i - 1][4] - hv[i][3] * hv[i][4]) / hv[i][0];
481 x[n - 2] = hv[n - 1][4] - hv[n - 2][1] * x[n - 1];
482 for (i = n - 3; i >= 0; i--)
483 x[i] = hv[i + 1][4] - hv[i][1] * x[i + 1] - hv[i][2] * x[i + 2];
491 * Calculation of approximation cubic splines
492 * Input: x[i], y[i], weights z[i]
494 * Returns matrix of spline coefficients
496 static spline_coeff *
498 struct curve_points *plot,
499 int first_point, /* where to start in plot->points */
500 int num_points) /* to determine end in plot->points */
504 double *r, *x, *h, *xp, *yp;
505 struct coordinate GPHUGE *this_points;
508 x_axis = plot->x_axis;
509 y_axis = plot->y_axis;
511 sc = gp_alloc((num_points) * sizeof(spline_coeff),
515 int_error(plot->token, "Can't calculate approximation splines, need at least 4 points");
517 this_points = (plot->points) + first_point;
519 for (i = 0; i < num_points; i++)
520 if (this_points[i].z <= 0)
521 int_error(plot->token, "Can't calculate approximation splines, all weights have to be > 0");
523 m = gp_alloc((num_points - 2) * sizeof(five_diag), "spline help matrix");
525 r = gp_alloc((num_points - 2) * sizeof(double), "spline right side");
526 x = gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
527 h = gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
529 xp = gp_alloc((num_points) * sizeof(double), "x pos");
530 yp = gp_alloc((num_points) * sizeof(double), "y pos");
532 /* KB 981107: With logarithmic axis first convert back to linear scale */
534 xp[0] = AXIS_DE_LOG_VALUE(x_axis, this_points[0].x);
535 yp[0] = AXIS_DE_LOG_VALUE(y_axis, this_points[0].y);
536 for (i = 1; i < num_points; i++) {
537 xp[i] = AXIS_DE_LOG_VALUE(x_axis, this_points[i].x);
538 yp[i] = AXIS_DE_LOG_VALUE(y_axis, this_points[i].y);
539 h[i - 1] = xp[i] - xp[i - 1];
542 /* set up the matrix and the vector */
544 for (i = 0; i <= num_points - 3; i++) {
545 r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
546 - (yp[i + 1] - yp[i]) / h[i]);
551 m[i][0] = 6 / this_points[i].z / h[i - 1] / h[i];
556 m[i][1] = h[i] - 6 / this_points[i].z / h[i] * (1 / h[i - 1] + 1 / h[i])
557 - 6 / this_points[i + 1].z / h[i] * (1 / h[i] + 1 / h[i + 1]);
559 m[i][2] = 2 * (h[i] + h[i + 1])
560 + 6 / this_points[i].z / h[i] / h[i]
561 + 6 / this_points[i + 1].z * (1 / h[i] + 1 / h[i + 1]) * (1 / h[i] + 1 / h[i + 1])
562 + 6 / this_points[i + 2].z / h[i + 1] / h[i + 1];
564 if (i > num_points - 4)
567 m[i][3] = h[i + 1] - 6 / this_points[i + 1].z / h[i + 1] * (1 / h[i] + 1 / h[i + 1])
568 - 6 / this_points[i + 2].z / h[i + 1] * (1 / h[i + 1] + 1 / h[i + 2]);
570 if (i > num_points - 5)
573 m[i][4] = 6 / this_points[i + 2].z / h[i + 1] / h[i + 2];
576 /* solve the matrix */
577 if (!solve_five_diag(m, r, x, num_points - 2)) {
584 int_error(plot->token, "Can't calculate approximation splines");
587 for (i = 1; i <= num_points - 2; i++)
589 sc[num_points - 1][2] = 0;
591 sc[0][0] = yp[0] + 2 / this_points[0].z / h[0] * (sc[0][2] - sc[1][2]);
592 for (i = 1; i <= num_points - 2; i++)
593 sc[i][0] = yp[i] - 2 / this_points[i].z *
594 (sc[i - 1][2] / h[i - 1]
595 - sc[i][2] * (1 / h[i - 1] + 1 / h[i])
596 + sc[i + 1][2] / h[i]);
597 sc[num_points - 1][0] = yp[num_points - 1]
598 - 2 / this_points[num_points - 1].z / h[num_points - 2]
599 * (sc[num_points - 2][2] - sc[num_points - 1][2]);
601 for (i = 0; i <= num_points - 2; i++) {
602 sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
603 - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
604 sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
619 * Calculation of cubic splines
621 * This can be treated as a special case of approximation cubic splines, with
622 * all weights -> infinity.
624 * Returns matrix of spline coefficients
626 static spline_coeff *
627 cp_tridiag(struct curve_points *plot, int first_point, int num_points)
631 double *r, *x, *h, *xp, *yp;
632 struct coordinate GPHUGE *this_points;
635 x_axis = plot->x_axis;
636 y_axis = plot->y_axis;
638 int_error(plot->token, "Can't calculate splines, need at least 3 points");
640 this_points = (plot->points) + first_point;
642 sc = gp_alloc((num_points) * sizeof(spline_coeff), "spline matrix");
643 m = gp_alloc((num_points - 2) * sizeof(tri_diag), "spline help matrix");
645 r = gp_alloc((num_points - 2) * sizeof(double), "spline right side");
646 x = gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
647 h = gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
649 xp = gp_alloc((num_points) * sizeof(double), "x pos");
650 yp = gp_alloc((num_points) * sizeof(double), "y pos");
652 /* KB 981107: With logarithmic axis first convert back to linear scale */
654 xp[0] = AXIS_DE_LOG_VALUE(x_axis,this_points[0].x);
655 yp[0] = AXIS_DE_LOG_VALUE(y_axis,this_points[0].y);
656 for (i = 1; i < num_points; i++) {
657 xp[i] = AXIS_DE_LOG_VALUE(x_axis,this_points[i].x);
658 yp[i] = AXIS_DE_LOG_VALUE(y_axis,this_points[i].y);
659 h[i - 1] = xp[i] - xp[i - 1];
662 /* set up the matrix and the vector */
664 for (i = 0; i <= num_points - 3; i++) {
665 r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
666 - (yp[i + 1] - yp[i]) / h[i]);
673 m[i][1] = 2 * (h[i] + h[i + 1]);
675 if (i > num_points - 4)
681 /* solve the matrix */
682 if (!solve_tri_diag(m, r, x, num_points - 2)) {
689 int_error(plot->token, "Can't calculate cubic splines");
692 for (i = 1; i <= num_points - 2; i++)
694 sc[num_points - 1][2] = 0;
696 for (i = 0; i <= num_points - 1; i++)
699 for (i = 0; i <= num_points - 2; i++) {
700 sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
701 - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
702 sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
717 struct curve_points *plot, /* still containes old plot->points */
718 spline_coeff *sc, /* generated by cp_tridiag */
719 int first_point, /* where to start in plot->points */
720 int num_points, /* to determine end in plot->points */
721 struct coordinate *dest) /* where to put the interpolated data */
723 double xdiff, temp, x, y;
724 double xstart, xend; /* Endpoints of the sampled x range */
726 struct coordinate GPHUGE *this_points;
728 /* min and max in internal (eg logged) co-ordinates. We update
729 * these, then update the external extrema in user co-ordinates
732 double ixmin, ixmax, iymin, iymax;
733 double sxmin, sxmax, symin, symax; /* starting values of above */
735 x_axis = plot->x_axis;
736 y_axis = plot->y_axis;
738 ixmin = sxmin = AXIS_LOG_VALUE(x_axis, X_AXIS.min);
739 ixmax = sxmax = AXIS_LOG_VALUE(x_axis, X_AXIS.max);
740 iymin = symin = AXIS_LOG_VALUE(y_axis, Y_AXIS.min);
741 iymax = symax = AXIS_LOG_VALUE(y_axis, Y_AXIS.max);
743 this_points = (plot->points) + first_point;
747 /* HBB 20010727: Sample only across the actual x range, not the
748 * full range of input data */
749 #if SAMPLE_CSPLINES_TO_FULL_RANGE
750 xstart = this_points[0].x;
751 xend = this_points[num_points - 1].x;
753 xstart = GPMAX(this_points[0].x, sxmin);
754 xend = GPMIN(this_points[num_points - 1].x, sxmax);
757 int_error(plot->token,
758 "Cannot smooth: no data within fixed xrange!");
760 xdiff = (xend - xstart) / (samples_1 - 1);
762 for (i = 0; i < samples_1; i++) {
763 x = xstart + i * xdiff;
765 /* Move forward to the spline interval this point is in */
766 while ((x >= this_points[l + 1].x) && (l < num_points - 2))
769 /* KB 981107: With logarithmic x axis the values were
770 * converted back to linear scale before calculating the
771 * coefficients. Use exponential for log x values. */
772 temp = AXIS_DE_LOG_VALUE(x_axis, x)
773 - AXIS_DE_LOG_VALUE(x_axis, this_points[l].x);
775 /* Evaluate cubic spline polynomial */
776 y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
778 /* With logarithmic y axis, we need to convert from linear to
782 y = AXIS_DO_LOG(y_axis, y);
784 y = symin - (symax - symin);
787 dest[i].type = INRANGE;
788 STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, X_AXIS.autoscale, NOOP, continue);
789 STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, Y_AXIS.autoscale, NOOP, NOOP);
791 dest[i].xlow = dest[i].xhigh = dest[i].x;
792 dest[i].ylow = dest[i].yhigh = dest[i].y;
798 UPDATE_RANGE(ixmax > sxmax, X_AXIS.max, ixmax, x_axis);
799 UPDATE_RANGE(ixmin < sxmin, X_AXIS.min, ixmin, x_axis);
800 UPDATE_RANGE(iymax > symax, Y_AXIS.max, iymax, y_axis);
801 UPDATE_RANGE(iymin < symin, Y_AXIS.min, iymin, y_axis);
807 * do_freq() is like the other smoothers only in that it
808 * needs to adjust the plot ranges. We don't have to copy
809 * approximated curves or anything like that.
814 struct curve_points *plot, /* still contains old plot->points */
815 int first_point, /* where to start in plot->points */
816 int num_points) /* to determine end in plot->points */
820 int x_axis = plot->x_axis;
821 int y_axis = plot->y_axis;
822 struct coordinate GPHUGE *this;
824 /* min and max in internal (eg logged) co-ordinates. We update
825 * these, then update the external extrema in user co-ordinates
829 double ixmin, ixmax, iymin, iymax;
830 double sxmin, sxmax, symin, symax; /* starting values of above */
832 ixmin = sxmin = AXIS_LOG_VALUE(x_axis, X_AXIS.min);
833 ixmax = sxmax = AXIS_LOG_VALUE(x_axis, X_AXIS.max);
834 iymin = symin = AXIS_LOG_VALUE(y_axis, Y_AXIS.min);
835 iymax = symax = AXIS_LOG_VALUE(y_axis, Y_AXIS.max);
837 this = (plot->points) + first_point;
839 for (i=0; i<num_points; i++) {
844 this[i].type = INRANGE;
846 STORE_AND_FIXUP_RANGE(this[i].x, x, this[i].type, ixmin, ixmax, X_AXIS.autoscale, NOOP, continue);
847 STORE_AND_FIXUP_RANGE(this[i].y, y, this[i].type, iymin, iymax, Y_AXIS.autoscale, NOOP, NOOP);
849 this[i].xlow = this[i].xhigh = this[i].x;
850 this[i].ylow = this[i].yhigh = this[i].y;
854 UPDATE_RANGE(ixmax > sxmax, X_AXIS.max, ixmax, x_axis);
855 UPDATE_RANGE(ixmin < sxmin, X_AXIS.min, ixmin, x_axis);
856 UPDATE_RANGE(iymax > symax, Y_AXIS.max, iymax, y_axis);
857 UPDATE_RANGE(iymin < symin, Y_AXIS.min, iymin, y_axis);
862 * Frequency plots have don't need new points allocated; we just need
863 * to adjust the plot ranges. Wedging this into gen_interp() would
864 * make that code even harder to read.
868 gen_interp_frequency(struct curve_points *plot)
872 int first_point, num_points;
874 curves = num_curves(plot);
877 for (i = 0; i < curves; i++) {
878 num_points = next_curve(plot, &first_point);
879 do_freq(plot, first_point, num_points);
880 first_point += num_points + 1;
886 * This is the main entry point used for everything except frequencies.
887 * As stated in the header, it is fine, but I'm not too happy with it.
891 gen_interp(struct curve_points *plot)
896 struct coordinate *new_points;
898 int first_point, num_points;
900 curves = num_curves(plot);
901 new_points = gp_alloc((samples_1 + 1) * curves * sizeof(struct coordinate),
902 "interpolation table");
905 for (i = 0; i < curves; i++) {
906 num_points = next_curve(plot, &first_point);
907 switch (plot->plot_smooth) {
908 case SMOOTH_CSPLINES:
909 sc = cp_tridiag(plot, first_point, num_points);
910 do_cubic(plot, sc, first_point, num_points,
911 new_points + i * (samples_1 + 1));
914 case SMOOTH_ACSPLINES:
915 sc = cp_approx_spline(plot, first_point, num_points);
916 do_cubic(plot, sc, first_point, num_points,
917 new_points + i * (samples_1 + 1));
923 bc = cp_binomial(num_points);
924 do_bezier(plot, bc, first_point, num_points,
925 new_points + i * (samples_1 + 1));
928 default: /* keep gcc -Wall quiet */
931 new_points[(i + 1) * (samples_1 + 1) - 1].type = UNDEFINED;
932 first_point += num_points;
936 plot->points = new_points;
937 plot->p_max = curves * (samples_1 + 1);
938 plot->p_count = plot->p_max - 1;
946 * sort data succession for further evaluation by plot_splines, etc.
947 * This routine is mainly introduced for compilers *NOT* supporting the
948 * UNIX qsort() routine. You can then easily replace it by more convenient
949 * stuff for your compiler.
953 /* HBB 20010720: To avoid undefined behaviour that would be caused by
954 * casting functions pointers around, changed arguments to what
955 * qsort() *really* wants */
956 /* HBB 20010720: removed 'static' to avoid HP-sUX gcc bug */
958 compare_points(SORTFUNC_ARGS arg1, SORTFUNC_ARGS arg2)
960 struct coordinate const *p1 = arg1;
961 struct coordinate const *p2 = arg2;
971 sort_points(struct curve_points *plot)
973 int first_point, num_points;
976 while ((num_points = next_curve(plot, &first_point)) > 0) {
977 /* Sort this set of points, does qsort handle 1 point correctly? */
978 /* HBB 20010720: removed casts -- they don't help a thing, but
979 * may hide problems */
980 qsort(plot->points + first_point, num_points,
981 sizeof(struct coordinate), compare_points);
982 first_point += num_points;
989 * cp_implode() if averaging is selected this function computes the new
990 * entries and shortens the whole thing to the necessary
996 cp_implode(struct curve_points *cp)
998 int first_point, num_points;
1000 double x = 0., y = 0., sux = 0., slx = 0., suy = 0., sly = 0.;
1001 /* int x_axis = cp->x_axis; */
1002 /* int y_axis = cp->y_axis; */
1003 TBOOLEAN all_inrange = FALSE;
1005 x_axis = cp->x_axis;
1006 y_axis = cp->y_axis;
1009 while ((num_points = next_curve(cp, &first_point)) > 0) {
1011 for (i = first_point; i < first_point + num_points; i++) {
1012 /* HBB 20020801: don't try to use undefined datapoints */
1013 if (cp->points[i].type == UNDEFINED)
1016 x = cp->points[i].x;
1017 y = cp->points[i].y;
1018 sux = cp->points[i].xhigh;
1019 slx = cp->points[i].xlow;
1020 suy = cp->points[i].yhigh;
1021 sly = cp->points[i].ylow;
1022 all_inrange = (cp->points[i].type == INRANGE);
1024 } else if (cp->points[i].x == x) {
1025 y += cp->points[i].y;
1026 sux += cp->points[i].xhigh;
1027 slx += cp->points[i].xlow;
1028 suy += cp->points[i].yhigh;
1029 sly += cp->points[i].ylow;
1030 if (cp->points[i].type != INRANGE)
1031 all_inrange = FALSE;
1034 cp->points[j].x = x;
1035 if ( cp->plot_smooth == SMOOTH_FREQUENCY )
1037 cp->points[j].y = y /= (double) k;
1038 cp->points[j].xhigh = sux / (double) k;
1039 cp->points[j].xlow = slx / (double) k;
1040 cp->points[j].yhigh = suy / (double) k;
1041 cp->points[j].ylow = sly / (double) k;
1042 /* HBB 20000405: I wanted to use STORE_AND_FIXUP_RANGE
1043 * here, but won't: it assumes we want to modify the
1044 * range, and that the range is given in 'input'
1045 * coordinates. For logarithmic axes, the overhead
1046 * would be larger than the possible gain, so write it
1047 * out explicitly, instead:
1049 cp->points[j].type = INRANGE;
1050 if (! all_inrange) {
1052 if (x <= -VERYLARGE) {
1053 cp->points[j].type = OUTRANGE;
1056 x = AXIS_UNDO_LOG(x_axis, x);
1058 if (((x < X_AXIS.min) && !(X_AXIS.autoscale & AUTOSCALE_MIN))
1059 || ((x > X_AXIS.max) && !(X_AXIS.autoscale & AUTOSCALE_MAX))) {
1060 cp->points[j].type = OUTRANGE;
1064 if (y <= -VERYLARGE) {
1065 cp->points[j].type = OUTRANGE;
1068 y = AXIS_UNDO_LOG(y_axis, y);
1070 if (((y < Y_AXIS.min) && !(Y_AXIS.autoscale & AUTOSCALE_MIN))
1071 || ((y > Y_AXIS.max) && !(Y_AXIS.autoscale & AUTOSCALE_MAX)))
1072 cp->points[j].type = OUTRANGE;
1075 } /* if(! all inrange) */
1077 j++; /* next valid entry */
1078 k = 0; /* to read */
1079 i--; /* from this (-> last after for(;;)) entry */
1080 } /* else (same x position) */
1081 } /* for(points in curve) */
1084 cp->points[j].x = x;
1085 if ( cp->plot_smooth == SMOOTH_FREQUENCY )
1087 cp->points[j].y = y /= (double) k;
1088 cp->points[j].xhigh = sux / (double) k;
1089 cp->points[j].xlow = slx / (double) k;
1090 cp->points[j].yhigh = suy / (double) k;
1091 cp->points[j].ylow = sly / (double) k;
1092 cp->points[j].type = INRANGE;
1093 if (! all_inrange) {
1095 if (x <= -VERYLARGE) {
1096 cp->points[j].type = OUTRANGE;
1099 x = AXIS_UNDO_LOG(x_axis, x);
1101 if (((x < X_AXIS.min) && !(X_AXIS.autoscale & AUTOSCALE_MIN))
1102 || ((x > X_AXIS.max) && !(X_AXIS.autoscale & AUTOSCALE_MAX))) {
1103 cp->points[j].type = OUTRANGE;
1107 if (y <= -VERYLARGE) {
1108 cp->points[j].type = OUTRANGE;
1111 y = AXIS_UNDO_LOG(y_axis, y);
1113 if (((y < Y_AXIS.min) && !(Y_AXIS.autoscale & AUTOSCALE_MIN))
1114 || ((y > Y_AXIS.max) && !(Y_AXIS.autoscale & AUTOSCALE_MAX)))
1115 cp->points[j].type = OUTRANGE;
1119 j++; /* next valid entry */
1121 /* insert invalid point to separate curves */
1122 if (j < cp->p_count) {
1123 cp->points[j].type = UNDEFINED;
1126 first_point += num_points;