2 static char *RCSid() { return RCSid("$Id: contour.c,v 1.27 2005/08/07 09:43:28 mikulik Exp $"); }
5 /* GNUPLOT - contour.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.
44 * Improvements to the numerical algorithms:
45 * Hans-Martin Keller, 1995,1997 (hkeller@gwdg.de)
53 /* #include "setshow.h" */
55 /* exported variables (to be handled by the 'set' and friends): */
57 char contour_format[32] = "%8.3g"; /* format for contour key entries */
58 t_contour_kind contour_kind = CONTOUR_KIND_LINEAR;
59 t_contour_levels_kind contour_levels_kind = LEVELS_AUTO;
60 int contour_levels = DEFAULT_CONTOUR_LEVELS;
61 int contour_order = DEFAULT_CONTOUR_ORDER;
62 int contour_pts = DEFAULT_NUM_APPROX_PTS;
64 /* storage for z levels to draw contours at */
65 dynarray dyn_contour_levels_list;
67 /* position of edge in mesh */
68 typedef enum en_edge_position {
75 /* FIXME HBB 2000052: yet another local copy of 'epsilon'. Why? */
76 #define EPSILON 1e-5 /* Used to decide if two float are equal. */
83 #define MAX_POINTS_PER_CNTR 100
85 #define SQR(x) ((x) * (x))
87 typedef struct edge_struct {
88 struct poly_struct *poly[2]; /* Each edge belongs to up to 2 polygons */
89 struct coordinate GPHUGE *vertex[2]; /* The two extreme points of this edge. */
90 struct edge_struct *next; /* To chain lists */
91 TBOOLEAN is_active; /* is edge is 'active' at certain Z level? */
92 t_edge_position position; /* position of edge in mesh */
95 typedef struct poly_struct {
96 struct edge_struct *edge[3]; /* As we do triangolation here... */
97 struct poly_struct *next; /* To chain lists. */
100 /* Contours are saved using this struct list. */
101 typedef struct cntr_struct {
102 double X, Y; /* The coordinates of this vertex. */
103 struct cntr_struct *next; /* To chain lists. */
106 static struct gnuplot_contours *contour_list = NULL;
107 static double crnt_cntr[MAX_POINTS_PER_CNTR * 2];
108 static int crnt_cntr_pt_index = 0;
109 static double contour_level = 0.0;
111 /* Linear, Cubic interp., Bspline: */
112 static t_contour_kind interp_kind = CONTOUR_KIND_LINEAR;
114 static double x_min, y_min, z_min; /* Minimum values of x, y, and z */
115 static double x_max, y_max, z_max; /* Maximum values of x, y, and z */
117 static void add_cntr_point __PROTO((double x, double y));
118 static void end_crnt_cntr __PROTO((void));
119 static void gen_contours __PROTO((edge_struct *p_edges, double z_level,
120 double xx_min, double xx_max,
121 double yy_min, double yy_max));
122 static int update_all_edges __PROTO((edge_struct *p_edges,
124 static cntr_struct *gen_one_contour __PROTO((edge_struct *p_edges,
126 TBOOLEAN *contr_isclosed,
128 static cntr_struct *trace_contour __PROTO((edge_struct *pe_start,
131 TBOOLEAN contr_isclosed));
132 static cntr_struct *update_cntr_pt __PROTO((edge_struct *p_edge,
134 static int fuzzy_equal __PROTO((cntr_struct *p_cntr1,
135 cntr_struct *p_cntr2));
138 static void gen_triangle __PROTO((int num_isolines,
139 struct iso_curve *iso_lines,
140 poly_struct **p_polys,
141 edge_struct **p_edges));
142 static void calc_min_max __PROTO((int num_isolines,
143 struct iso_curve *iso_lines,
144 double *xx_min, double *yy_min,
146 double *xx_max, double *yy_max,
148 static edge_struct *add_edge __PROTO((struct coordinate GPHUGE *point0,
149 struct coordinate GPHUGE *point1,
152 edge_struct **pe_tail));
153 static poly_struct *add_poly __PROTO((edge_struct *edge0,
156 poly_struct **p_poly,
157 poly_struct **pp_tail));
159 static void put_contour __PROTO((cntr_struct *p_cntr,
160 double xx_min, double xx_max,
161 double yy_min, double yy_max,
162 TBOOLEAN contr_isclosed));
163 static void put_contour_nothing __PROTO((cntr_struct *p_cntr));
164 static int chk_contour_kind __PROTO((cntr_struct *p_cntr,
165 TBOOLEAN contr_isclosed));
166 static void put_contour_cubic __PROTO((cntr_struct *p_cntr,
167 double xx_min, double xx_max,
168 double yy_min, double yy_max,
169 TBOOLEAN contr_isclosed));
170 static void put_contour_bspline __PROTO((cntr_struct *p_cntr,
171 TBOOLEAN contr_isclosed));
172 static void free_contour __PROTO((cntr_struct *p_cntr));
173 static int count_contour __PROTO((cntr_struct *p_cntr));
174 static int gen_cubic_spline __PROTO((int num_pts, cntr_struct *p_cntr,
175 double d2x[], double d2y[],
177 TBOOLEAN contr_isclosed,
178 double unit_x, double unit_y));
179 static void intp_cubic_spline __PROTO((int n, cntr_struct *p_cntr,
180 double d2x[], double d2y[],
181 double delta_t[], int n_intpol));
182 static int solve_cubic_1 __PROTO((tri_diag m[], int n));
183 static void solve_cubic_2 __PROTO((tri_diag m[], double x[], int n));
184 static void gen_bspline_approx __PROTO((cntr_struct *p_cntr,
185 int num_of_points, int order,
186 TBOOLEAN contr_isclosed));
187 static void eval_bspline __PROTO((double t, cntr_struct *p_cntr,
188 int num_of_points, int order, int j,
189 TBOOLEAN contr_isclosed, double *x,
191 static double fetch_knot __PROTO((TBOOLEAN contr_isclosed, int num_of_points,
195 * Entry routine to this whole set of contouring module.
197 struct gnuplot_contours *
198 contour(int num_isolines, struct iso_curve *iso_lines)
201 int num_of_z_levels; /* # Z contour levels. */
202 poly_struct *p_polys, *p_poly;
203 edge_struct *p_edges, *p_edge;
204 double z = 0, dz = 0;
205 struct gnuplot_contours *save_contour_list;
207 /* HBB FIXME 20050804: The number of contour_levels as set by 'set
208 * cnrparam lev inc a,b,c' is almost certainly wrong if z axis is
210 num_of_z_levels = contour_levels;
211 interp_kind = contour_kind;
216 * Calculate min/max values :
218 calc_min_max(num_isolines, iso_lines,
219 &x_min, &y_min, &z_min, &x_max, &y_max, &z_max);
222 * Generate list of edges (p_edges) and list of triangles (p_polys):
224 gen_triangle(num_isolines, iso_lines, &p_polys, &p_edges);
225 crnt_cntr_pt_index = 0;
227 if (contour_levels_kind == LEVELS_AUTO) {
228 dz = fabs(z_max - z_min);
230 return NULL; /* empty z range ? */
231 /* Find a tic step that will generate approximately the
232 * desired number of contour levels. The "* 2" is historical.
234 dz = quantize_normal_tics(dz, ((int) contour_levels + 1) * 2);
235 z = floor(z_min / dz) * dz;
236 num_of_z_levels = (int) floor((z_max - z) / dz);
238 for (i = 0; i < num_of_z_levels; i++) {
239 switch (contour_levels_kind) {
243 case LEVELS_INCREMENTAL:
244 z = AXIS_LOG_VALUE(FIRST_Z_AXIS, contour_levels_list[0]) +
245 i * AXIS_LOG_VALUE(FIRST_Z_AXIS, contour_levels_list[1]);
247 case LEVELS_DISCRETE:
248 z = AXIS_LOG_VALUE(FIRST_Z_AXIS, contour_levels_list[i]);
252 save_contour_list = contour_list;
253 gen_contours(p_edges, z, x_min, x_max, y_min, y_max);
254 if (contour_list != save_contour_list) {
255 contour_list->isNewLevel = 1;
256 sprintf(contour_list->label, contour_format, AXIS_DE_LOG_VALUE(FIRST_Z_AXIS,z));
261 /* Free all contouring related temporary data. */
263 p_poly = p_polys->next;
268 p_edge = p_edges->next;
277 * Adds another point to the currently build contour.
280 add_cntr_point(double x, double y)
284 if (crnt_cntr_pt_index >= MAX_POINTS_PER_CNTR - 1) {
285 index = crnt_cntr_pt_index - 1;
287 crnt_cntr[0] = crnt_cntr[index * 2];
288 crnt_cntr[1] = crnt_cntr[index * 2 + 1];
289 crnt_cntr_pt_index = 1; /* Keep the last point as first of this one. */
291 crnt_cntr[crnt_cntr_pt_index * 2] = x;
292 crnt_cntr[crnt_cntr_pt_index * 2 + 1] = y;
293 crnt_cntr_pt_index++;
297 * Done with current contour - create gnuplot data structure for it.
303 struct gnuplot_contours *cntr =
304 gp_alloc(sizeof(struct gnuplot_contours), "gnuplot_contour");
306 gp_alloc(sizeof(struct coordinate) * crnt_cntr_pt_index,
309 for (i = 0; i < crnt_cntr_pt_index; i++) {
310 cntr->coords[i].x = crnt_cntr[i * 2];
311 cntr->coords[i].y = crnt_cntr[i * 2 + 1];
312 cntr->coords[i].z = contour_level;
314 cntr->num_pts = crnt_cntr_pt_index;
316 cntr->next = contour_list;
318 contour_list->isNewLevel = 0;
320 crnt_cntr_pt_index = 0;
324 * Generates all contours by tracing the intersecting triangles.
328 edge_struct *p_edges,
330 double xx_min, double xx_max,
331 double yy_min, double yy_max)
333 int num_active; /* Number of edges marked ACTIVE. */
334 TBOOLEAN contr_isclosed; /* Is this contour a closed line? */
337 num_active = update_all_edges(p_edges, z_level); /* Do pass 1. */
339 contr_isclosed = FALSE; /* Start to look for contour on boundaries. */
341 while (num_active > 0) { /* Do Pass 2. */
342 /* Generate One contour (and update NumActive as needed): */
343 p_cntr = gen_one_contour(p_edges, z_level, &contr_isclosed, &num_active);
344 /* Emit it in requested format: */
345 put_contour(p_cntr, xx_min, xx_max, yy_min, yy_max, contr_isclosed);
350 * Does pass 1, or marks the edges which are active (crosses this z_level)
351 * Returns number of active edges (marked ACTIVE).
354 update_all_edges(edge_struct *p_edges, double z_level)
359 /* use the same test at both vertices to avoid roundoff errors */
360 if ((p_edges->vertex[0]->z >= z_level) !=
361 (p_edges->vertex[1]->z >= z_level)) {
362 p_edges->is_active = TRUE;
365 p_edges->is_active = FALSE;
366 p_edges = p_edges->next;
373 * Does pass 2, or find one complete contour out of the triangulation
376 * Returns a pointer to the contour (as linked list), contr_isclosed
377 * tells if the contour is a closed line or not, and num_active is
382 edge_struct *p_edges, /* list of edges input */
383 double z_level, /* Z level of contour input */
384 TBOOLEAN *contr_isclosed, /* open or closed contour, in/out */
385 int *num_active) /* number of active edges in/out */
387 edge_struct *pe_temp;
389 if (! *contr_isclosed) {
390 /* Look for something to start with on boundary: */
393 if (pe_temp->is_active && (pe_temp->position == BOUNDARY))
395 pe_temp = pe_temp->next;
398 *contr_isclosed = TRUE; /* No more contours on boundary. */
400 return trace_contour(pe_temp, z_level, num_active, *contr_isclosed);
403 if (*contr_isclosed) {
404 /* Look for something to start with inside: */
407 if (pe_temp->is_active && (pe_temp->position != BOUNDARY))
409 pe_temp = pe_temp->next;
413 fprintf(stderr, "gen_one_contour: no contour found\n");
416 *contr_isclosed = TRUE;
417 return trace_contour(pe_temp, z_level, num_active, *contr_isclosed);
420 return NULL; /* We should never be here, but lint... */
424 * Search the data base along a contour starts at the edge pe_start until
425 * a boundary edge is detected or until we close the loop back to pe_start.
426 * Returns a linked list of all the points on the contour
427 * Also decreases num_active by the number of points on contour.
431 edge_struct *pe_start, /* edge to start contour input */
432 double z_level, /* Z level of contour input */
433 int *num_active, /* number of active edges in/out */
434 TBOOLEAN contr_isclosed) /* open or closed contour line (input) */
436 cntr_struct *p_cntr, *pc_tail;
437 edge_struct *p_edge, *p_next_edge;
438 poly_struct *p_poly, *PLastpoly = NULL;
441 p_edge = pe_start; /* first edge to start contour */
443 /* Generate the header of the contour - the point on pe_start. */
444 if (! contr_isclosed) {
445 pe_start->is_active = FALSE;
448 if (p_edge->poly[0] || p_edge->poly[1]) { /* more than one point */
450 p_cntr = pc_tail = update_cntr_pt(pe_start, z_level); /* first point */
453 /* Find polygon to continue (Not where we came from - PLastpoly): */
454 if (p_edge->poly[0] == PLastpoly)
455 p_poly = p_edge->poly[1];
457 p_poly = p_edge->poly[0];
458 p_next_edge = NULL; /* In case of error, remains NULL. */
459 for (i = 0; i < 3; i++) /* Test the 3 edges of the polygon: */
460 if (p_poly->edge[i] != p_edge)
461 if (p_poly->edge[i]->is_active)
462 p_next_edge = p_poly->edge[i];
463 if (!p_next_edge) { /* Error exit */
464 pc_tail->next = NULL;
465 free_contour(p_cntr);
466 fprintf(stderr, "trace_contour: unexpected end of contour\n");
469 p_edge = p_next_edge;
471 p_edge->is_active = FALSE;
474 /* Do not allocate contour points on diagonal edges */
475 if (p_edge->position != DIAGONAL) {
477 pc_tail->next = update_cntr_pt(p_edge, z_level);
479 /* Remove nearby points */
480 if (fuzzy_equal(pc_tail, pc_tail->next)) {
484 pc_tail = pc_tail->next;
486 } while ((p_edge != pe_start) && (p_edge->position != BOUNDARY));
488 pc_tail->next = NULL;
490 /* For closed contour the first and last point should be equal */
491 if (pe_start == p_edge) {
492 (p_cntr->X) = (pc_tail->X);
493 (p_cntr->Y) = (pc_tail->Y);
495 } else { /* only one point, forget it */
503 * Allocates one contour location and update it to to correct position
504 * according to z_level and edge p_edge.
507 update_cntr_pt(edge_struct *p_edge, double z_level)
512 t = (z_level - p_edge->vertex[0]->z) /
513 (p_edge->vertex[1]->z - p_edge->vertex[0]->z);
515 /* test if t is out of interval [0:1] (should not happen but who knows ...) */
516 t = (t < 0.0 ? 0.0 : t);
517 t = (t > 1.0 ? 1.0 : t);
519 p_cntr = gp_alloc(sizeof(cntr_struct), "contour cntr_struct");
521 p_cntr->X = p_edge->vertex[1]->x * t +
522 p_edge->vertex[0]->x * (1 - t);
523 p_cntr->Y = p_edge->vertex[1]->y * t +
524 p_edge->vertex[0]->y * (1 - t);
528 /* Simple routine to decide if two contour points are equal by
529 * calculating the relative error (< EPSILON). */
530 /* HBB 20010121: don't use absolute value 'zero' to compare to data
533 fuzzy_equal(cntr_struct *p_cntr1, cntr_struct *p_cntr2)
535 double unit_x, unit_y;
536 unit_x = fabs(x_max - x_min); /* reference */
537 unit_y = fabs(y_max - y_min);
538 return ((fabs(p_cntr1->X - p_cntr2->X) < unit_x * EPSILON)
539 && (fabs(p_cntr1->Y - p_cntr2->Y) < unit_y * EPSILON));
543 * Generate the triangles.
544 * Returns the lists (edges & polys) via pointers to their heads.
548 int num_isolines, /* number of iso-lines input */
549 struct iso_curve *iso_lines, /* iso-lines input */
550 poly_struct **p_polys, /* list of polygons output */
551 edge_struct **p_edges) /* list of edges output */
553 int i, j, grid_x_max = iso_lines->p_count;
554 edge_struct *p_edge1, *p_edge2, *edge0, *edge1, *edge2, *pe_tail,
556 poly_struct *pp_tail, *lower_tri, *upper_tri;
557 /* HBB 980308: need to tag *each* of them as GPHUGE! */
558 struct coordinate GPHUGE *p_vrtx1, GPHUGE * p_vrtx2;
560 (*p_polys) = pp_tail = NULL; /* clear lists */
561 (*p_edges) = pe_tail = NULL;
563 p_vrtx1 = iso_lines->points; /* first row of vertices */
564 p_edge1 = pe_tail = NULL; /* clear list of edges */
566 /* Generate edges of first row */
567 for (j = 0; j < grid_x_max - 1; j++)
568 add_edge(p_vrtx1 + j, p_vrtx1 + j + 1, &p_edge1, &pe_tail);
570 (*p_edges) = p_edge1; /* update main list */
574 * Combines vertices to edges and edges to triangles:
575 * ==================================================
576 * The edges are stored in the edge list, referenced by p_edges
577 * (pe_tail points on last edge).
579 * Temporary pointers:
580 * 1. p_edge2: Top horizontal edge list: +-----------------------+ 2
581 * 2. p_tail : end of middle edge list: |\ |\ |\ |\ |\ |\ |
582 * | \| \| \| \| \| \|
583 * 3. p_edge1: Bottom horizontal edge list: +-----------------------+ 1
585 * pe_tail2 : end of list beginning at p_edge2
586 * pe_temp : position inside list beginning at p_edge1
587 * p_edges : head of the master edge list (part of our output)
588 * p_vrtx1 : start of lower row of input vertices
589 * p_vrtx2 : start of higher row of input vertices
591 * The routine generates two triangle Lower Upper 1
592 * upper one and lower one: | \ ----
593 * (Nums. are edges order in polys) 0| \1 0\ |2
594 * The polygons are stored in the polygon ---- \ |
595 * list (*p_polys) (pp_tail points on 2
599 * In addition, the edge lists are updated - | \ 0 |
600 * each edge has two pointers on the two | \ |
601 * (one active if boundary) polygons which 0|1 0\1 0|1
602 * uses it. These two pointer to polygons | \ |
603 * are named: poly[0], poly[1]. The diagram | 1 \ |
604 * on the right show how they are used for the -----------
605 * upper and lower polygons (INNER_MESH polygons only). 0
608 for (i = 1; i < num_isolines; i++) {
609 /* Read next column and gen. polys. */
610 iso_lines = iso_lines->next;
612 p_vrtx2 = iso_lines->points; /* next row of vertices */
613 p_edge2 = pe_tail2 = NULL; /* clear top horizontal list */
614 pe_temp = p_edge1; /* pointer in bottom list */
617 * Generate edges and triagles for next row:
620 /* generate first vertical edge */
621 edge2 = add_edge(p_vrtx1, p_vrtx2, p_edges, &pe_tail);
623 for (j = 0; j < grid_x_max - 1; j++) {
625 /* copy vertical edge for lower triangle */
628 if (pe_temp && pe_temp->vertex[0] == p_vrtx1 + j) {
629 /* test lower edge */
631 pe_temp = pe_temp->next;
633 edge2 = NULL; /* edge is undefined */
636 /* generate diagonal edge */
637 edge1 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j, p_edges, &pe_tail);
639 edge1->position = DIAGONAL;
641 /* generate lower triangle */
642 lower_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
644 /* copy diagonal edge for upper triangle */
647 /* generate upper edge */
648 edge1 = add_edge(p_vrtx2 + j, p_vrtx2 + j + 1, &p_edge2, &pe_tail2);
650 /* generate vertical edge */
651 edge2 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j + 1, p_edges, &pe_tail);
653 /* generate upper triangle */
654 upper_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
658 /* HBB 19991130 bugfix: if p_edge2 list is empty,
659 * don't change p_edges list! Crashes by access
660 * to NULL pointer pe_tail, the second time through,
662 if ((*p_edges)) { /* Chain new edges to main list. */
663 pe_tail->next = p_edge2;
666 (*p_edges) = p_edge2;
671 /* this row finished, move list heads up one row: */
676 /* Update the boundary flag, saved in each edge, and update indexes: */
678 pe_temp = (*p_edges);
681 if ((!(pe_temp->poly[0])) || (!(pe_temp->poly[1])))
682 (pe_temp->position) = BOUNDARY;
683 pe_temp = pe_temp->next;
688 * Calculate minimum and maximum values
692 int num_isolines, /* number of iso-lines input */
693 struct iso_curve *iso_lines, /* iso-lines input */
694 double *xx_min, double *yy_min, double *zz_min,
695 double *xx_max, double *yy_max, double *zz_max) /* min/max values in/out */
697 int i, j, grid_x_max;
698 struct coordinate GPHUGE *vertex;
700 grid_x_max = iso_lines->p_count; /* number of vertices per iso_line */
702 (*xx_min) = (*yy_min) = (*zz_min) = VERYLARGE; /* clear min/max values */
703 (*xx_max) = (*yy_max) = (*zz_max) = -VERYLARGE;
705 for (j = 0; j < num_isolines; j++) {
707 vertex = iso_lines->points;
709 for (i = 0; i < grid_x_max; i++) {
710 if (vertex[i].type != UNDEFINED) {
711 if (vertex[i].x > (*xx_max))
712 (*xx_max) = vertex[i].x;
713 if (vertex[i].y > (*yy_max))
714 (*yy_max) = vertex[i].y;
715 if (vertex[i].z > (*zz_max))
716 (*zz_max) = vertex[i].z;
717 if (vertex[i].x < (*xx_min))
718 (*xx_min) = vertex[i].x;
719 if (vertex[i].y < (*yy_min))
720 (*yy_min) = vertex[i].y;
721 if (vertex[i].z < (*zz_min))
722 (*zz_min) = vertex[i].z;
725 iso_lines = iso_lines->next;
727 /* HBB 20000426: this code didn't take into account that axes might
730 /* HBB 20001220: DON'T. The values are actually already stored
731 * logarithmized, as should be! */
732 axis_unlog_interval(FIRST_X_AXIS, xx_min, xx_max, 0);
733 axis_unlog_interval(FIRST_Y_AXIS, yy_min, yy_max, 0);
734 axis_unlog_interval(FIRST_Z_AXIS, zz_min, zz_max, 0);
738 * fprintf(stderr," x: %g, %g\n", (*xx_min), (*xx_max));
739 * fprintf(stderr," y: %g, %g\n", (*yy_min), (*yy_max));
740 * fprintf(stderr," z: %g, %g\n", (*zz_min), (*zz_max));
745 * Generate new edge and append it to list, but only if both vertices are
746 * defined. The list is referenced by p_edge and pe_tail (p_edge points on
747 * first edge and pe_tail on last one).
748 * Note, the list may be empty (pe_edge==pe_tail==NULL) on entry and exit.
752 struct coordinate GPHUGE *point0, /* 2 vertices input */
753 struct coordinate GPHUGE *point1,
754 edge_struct **p_edge, /* pointers to edge list in/out */
755 edge_struct **pe_tail)
757 edge_struct *pe_temp = NULL;
760 if (point0->type == INRANGE && point1->type == INRANGE)
762 if (point0->type != UNDEFINED && point1->type != UNDEFINED)
765 pe_temp = gp_alloc(sizeof(edge_struct), "contour edge");
767 pe_temp->poly[0] = NULL; /* clear links */
768 pe_temp->poly[1] = NULL;
769 pe_temp->vertex[0] = point0; /* First vertex of edge. */
770 pe_temp->vertex[1] = point1; /* Second vertex of edge. */
771 pe_temp->next = NULL;
772 pe_temp->position = INNER_MESH; /* default position in mesh */
775 (*pe_tail)->next = pe_temp; /* Stick new record as last one. */
777 (*p_edge) = pe_temp; /* start new list if empty */
779 (*pe_tail) = pe_temp; /* continue to last record. */
782 return pe_temp; /* returns NULL, if no edge allocated */
786 * Generate new triangle and append it to list, but only if all edges are defined.
787 * The list is referenced by p_poly and pp_tail (p_poly points on first ploygon
788 * and pp_tail on last one).
789 * Note, the list may be empty (pe_ploy==pp_tail==NULL) on entry and exit.
795 edge_struct *edge2, /* 3 edges input */
796 poly_struct **p_poly,
797 poly_struct **pp_tail) /* pointers to polygon list in/out */
799 poly_struct *pp_temp = NULL;
801 if (edge0 && edge1 && edge2) {
802 pp_temp = gp_alloc(sizeof(poly_struct), "contour polygon");
804 pp_temp->edge[0] = edge0; /* First edge of triangle */
805 pp_temp->edge[1] = edge1; /* Second one */
806 pp_temp->edge[2] = edge2; /* Third one */
807 pp_temp->next = NULL;
809 if (edge0->poly[0]) /* update edge0 */
810 edge0->poly[1] = pp_temp;
812 edge0->poly[0] = pp_temp;
814 if (edge1->poly[0]) /* update edge1 */
815 edge1->poly[1] = pp_temp;
817 edge1->poly[0] = pp_temp;
819 if (edge2->poly[0]) /* update edge2 */
820 edge2->poly[1] = pp_temp;
822 edge2->poly[0] = pp_temp;
824 if ((*pp_tail)) /* Stick new record as last one. */
825 (*pp_tail)->next = pp_temp;
827 (*p_poly) = pp_temp; /* start new list if empty */
829 (*pp_tail) = pp_temp; /* continue to last record. */
832 return pp_temp; /* returns NULL, if no edge allocated */
838 * Calls the (hopefully) desired interpolation/approximation routine.
842 cntr_struct *p_cntr, /* contour structure input */
843 double xx_min, double xx_max,
844 double yy_min, double yy_max, /* minimum/maximum values input */
845 TBOOLEAN contr_isclosed) /* contour line closed? (input) */
849 return; /* Nothing to do if it is empty contour. */
851 switch (interp_kind) {
852 case CONTOUR_KIND_LINEAR: /* No interpolation/approximation. */
853 put_contour_nothing(p_cntr);
855 case CONTOUR_KIND_CUBIC_SPL: /* Cubic spline interpolation. */
856 put_contour_cubic(p_cntr, xx_min, xx_max, yy_min, yy_max,
857 chk_contour_kind(p_cntr, contr_isclosed));
860 case CONTOUR_KIND_BSPLINE: /* Bspline approximation. */
861 put_contour_bspline(p_cntr,
862 chk_contour_kind(p_cntr, contr_isclosed));
865 free_contour(p_cntr);
869 * Simply puts contour coordinates in order with no interpolation or
873 put_contour_nothing(cntr_struct *p_cntr)
876 add_cntr_point(p_cntr->X, p_cntr->Y);
877 p_cntr = p_cntr->next;
883 * for some reason contours are never flagged as 'isclosed'
884 * if first point == last point, set flag accordingly
889 chk_contour_kind(cntr_struct *p_cntr, TBOOLEAN contr_isclosed)
891 cntr_struct *pc_tail = NULL;
892 TBOOLEAN current_contr_isclosed;
894 current_contr_isclosed = contr_isclosed;
896 if (! contr_isclosed) {
898 while (pc_tail->next)
899 pc_tail = pc_tail->next; /* Find last point. */
901 /* test if first and last point are equal */
902 if (fuzzy_equal(pc_tail, p_cntr))
903 current_contr_isclosed = TRUE;
905 return (current_contr_isclosed);
909 * Generate a cubic spline curve through the points (x_i,y_i) which are
910 * stored in the linked list p_cntr.
911 * The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the
912 * parameter t is the length of the linear stroke.
917 double xx_min, double xx_max,
918 double yy_min, double yy_max,
919 TBOOLEAN contr_isclosed)
921 int num_pts, num_intpol;
922 double unit_x, unit_y; /* To define norm (x,y)-plane */
923 double *delta_t; /* Interval length t_{i+1}-t_i */
924 double *d2x, *d2y; /* Second derivatives x''(t_i), y''(t_i) */
925 cntr_struct *pc_tail;
927 num_pts = count_contour(p_cntr); /* Number of points in contour. */
929 pc_tail = p_cntr; /* Find last point. */
930 while (pc_tail->next)
931 pc_tail = pc_tail->next;
933 if (contr_isclosed) {
934 /* Test if first and last point are equal (should be) */
935 if (!fuzzy_equal(pc_tail, p_cntr)) {
936 pc_tail->next = p_cntr; /* Close contour list - make it circular. */
940 delta_t = gp_alloc(num_pts * sizeof(double), "contour delta_t");
941 d2x = gp_alloc(num_pts * sizeof(double), "contour d2x");
942 d2y = gp_alloc(num_pts * sizeof(double), "contour d2y");
944 /* Width and height of the grid is used as a unit length (2d-norm) */
945 unit_x = xx_max - xx_min;
946 unit_y = yy_max - yy_min;
947 /* FIXME HBB 20010121: 'zero' should not be used as an absolute
948 * figure to compare to data */
949 unit_x = (unit_x > zero ? unit_x : zero); /* should not be zero */
950 unit_y = (unit_y > zero ? unit_y : zero);
954 * Calculate second derivatives d2x[], d2y[] and interval lengths delta_t[]:
956 if (!gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t,
957 contr_isclosed, unit_x, unit_y)) {
962 pc_tail->next = NULL; /* Un-circular list */
966 /* If following (num_pts > 1) is TRUE then exactly 2 points in contour. */
967 else if (num_pts > 1) {
968 /* set all second derivatives to zero, interval length to 1 */
974 } else { /* Only one point ( ?? ) - ignore it. */
979 pc_tail->next = NULL; /* Un-circular list */
983 /* Calculate "num_intpol" interpolated values */
984 num_intpol = 1 + (num_pts - 1) * contour_pts; /* global: contour_pts */
985 intp_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, num_intpol);
992 pc_tail->next = NULL; /* Un-circular list */
999 * Find Bspline approximation for this data set.
1000 * Uses global variable contour_pts to determine number of samples per
1001 * interval, where the knot vector intervals are assumed to be uniform, and
1002 * global variable contour_order for the order of Bspline to use.
1005 put_contour_bspline(cntr_struct *p_cntr, TBOOLEAN contr_isclosed)
1008 int order = contour_order - 1;
1010 num_pts = count_contour(p_cntr); /* Number of points in contour. */
1012 return; /* Can't do nothing if empty or one points! */
1013 /* Order must be less than number of points in curve - fix it if needed. */
1014 if (order > num_pts - 1)
1015 order = num_pts - 1;
1017 gen_bspline_approx(p_cntr, num_pts, order, contr_isclosed);
1022 * Free all elements in the contour list.
1025 free_contour(cntr_struct *p_cntr)
1027 cntr_struct *pc_temp;
1031 p_cntr = p_cntr->next;
1037 * Counts number of points in contour.
1040 count_contour(cntr_struct *p_cntr)
1046 p_cntr = p_cntr->next;
1052 * Find second derivatives (x''(t_i),y''(t_i)) of cubic spline interpolation
1053 * through list of points (x_i,y_i). The parameter t is calculated as the
1054 * length of the linear stroke. The number of points must be at least 3.
1055 * Note: For closed contours the first and last point must be equal.
1059 int num_pts, /* Number of points (num_pts>=3), input */
1060 cntr_struct *p_cntr, /* List of points (x(t_i),y(t_i)), input */
1061 double d2x[], double d2y[], /* Second derivatives (x''(t_i),y''(t_i)), output */
1062 double delta_t[], /* List of interval lengths t_{i+1}-t_{i}, output */
1063 TBOOLEAN contr_isclosed, /* Closed or open contour?, input */
1064 double unit_x, double unit_y) /* Unit length in x and y (norm=1), input */
1068 tri_diag *m; /* The tri-diagonal matrix is saved here. */
1069 cntr_struct *pc_temp;
1071 m = gp_alloc(num_pts * sizeof(tri_diag), "contour tridiag m");
1074 * Calculate first differences in (d2x[i], d2y[i]) and interval lengths
1078 for (i = 0; i < num_pts - 1; i++) {
1079 d2x[i] = pc_temp->next->X - pc_temp->X;
1080 d2y[i] = pc_temp->next->Y - pc_temp->Y;
1082 * The norm of a linear stroke is calculated in "normal coordinates"
1083 * and used as interval length:
1085 delta_t[i] = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y));
1087 d2x[i] /= delta_t[i]; /* first difference, with unit norm: */
1088 d2y[i] /= delta_t[i]; /* || (d2x[i], d2y[i]) || = 1 */
1090 pc_temp = pc_temp->next;
1094 * Setup linear system: m * x = b
1096 n = num_pts - 2; /* Without first and last point */
1097 if (contr_isclosed) {
1098 /* First and last points must be equal for closed contours */
1099 delta_t[num_pts - 1] = delta_t[0];
1100 d2x[num_pts - 1] = d2x[0];
1101 d2y[num_pts - 1] = d2y[0];
1102 n++; /* Add last point (= first point) */
1104 for (i = 0; i < n; i++) {
1105 /* Matrix M, mainly tridiagonal with cyclic second index ("j = j+n mod n") */
1106 m[i][0] = delta_t[i]; /* Off-diagonal element M_{i,i-1} */
1107 m[i][1] = 2. * (delta_t[i] + delta_t[i + 1]); /* M_{i,i} */
1108 m[i][2] = delta_t[i + 1]; /* Off-diagonal element M_{i,i+1} */
1110 /* Right side b_x and b_y */
1111 d2x[i] = (d2x[i + 1] - d2x[i]) * 6.;
1112 d2y[i] = (d2y[i + 1] - d2y[i]) * 6.;
1115 * If the linear stroke shows a cusps of more than 90 degree, the right
1116 * side is reduced to avoid oscillations in the spline:
1118 norm = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y)) / 8.5;
1123 /* The first derivative will not be continuous */
1127 if (!contr_isclosed) {
1128 /* Third derivative is set to zero at both ends */
1129 m[0][1] += m[0][0]; /* M_{0,0} */
1130 m[0][0] = 0.; /* M_{0,n-1} */
1131 m[n - 1][1] += m[n - 1][2]; /* M_{n-1,n-1} */
1132 m[n - 1][2] = 0.; /* M_{n-1,0} */
1134 /* Solve linear systems for d2x[] and d2y[] */
1137 if (solve_cubic_1(m, n)) { /* Calculate Cholesky decomposition */
1138 solve_cubic_2(m, d2x, n); /* solve M * d2x = b_x */
1139 solve_cubic_2(m, d2y, n); /* solve M * d2y = b_y */
1141 } else { /* Should not happen, but who knows ... */
1146 /* Shift all second derivatives one place right and abdate end points */
1147 for (i = n; i > 0; i--) {
1148 d2x[i] = d2x[i - 1];
1149 d2y[i] = d2y[i - 1];
1151 if (contr_isclosed) {
1155 d2x[0] = d2x[1]; /* Third derivative is zero in */
1156 d2y[0] = d2y[1]; /* first and last interval */
1157 d2x[n + 1] = d2x[n];
1158 d2y[n + 1] = d2y[n];
1166 * Calculate interpolated values of the spline function (defined via p_cntr
1167 * and the second derivatives d2x[] and d2y[]). The number of tabulated
1168 * values is n. On an equidistant grid n_intpol values are calculated.
1173 cntr_struct *p_cntr,
1174 double d2x[], double d2y[], double delta_t[],
1177 double t, t_skip, t_max;
1178 double x0, x1, x, y0, y1, y;
1179 double d, hx, dx0, dx01, hy, dy0, dy01;
1182 /* The length of the total interval */
1184 for (i = 0; i < n - 1; i++)
1185 t_max += delta_t[i];
1187 /* The distance between interpolated points */
1188 t_skip = (1. - 1e-7) * t_max / (n_intpol - 1);
1190 t = 0.; /* Parameter value */
1193 add_cntr_point(x1, y1); /* First point. */
1196 for (i = 0; i < n - 1; i++) {
1197 p_cntr = p_cntr->next;
1199 d = delta_t[i]; /* Interval length */
1206 dx0 = (d2x[i + 1] + 2 * d2x[i]) / 6.;
1207 dy0 = (d2y[i + 1] + 2 * d2y[i]) / 6.;
1208 dx01 = (d2x[i + 1] - d2x[i]) / (6. * d);
1209 dy01 = (d2y[i + 1] - d2y[i]) / (6. * d);
1210 while (t <= delta_t[i]) { /* t in current interval ? */
1211 x = x0 + t * (hx + (t - d) * (dx0 + t * dx01));
1212 y = y0 + t * (hy + (t - d) * (dy0 + t * dy01));
1213 add_cntr_point(x, y); /* next point. */
1216 t -= delta_t[i]; /* Parameter t relative to start of next interval */
1221 * The following two procedures solve the special linear system which arise
1222 * in cubic spline interpolation. If x is assumed cyclic ( x[i]=x[n+i] ) the
1223 * equations can be written as (i=0,1,...,n-1):
1224 * m[i][0] * x[i-1] + m[i][1] * x[i] + m[i][2] * x[i+1] = b[i] .
1225 * In matrix notation one gets M * x = b, where the matrix M is tridiagonal
1226 * with additional elements in the upper right and lower left position:
1227 * m[i][0] = M_{i,i-1} for i=1,2,...,n-1 and m[0][0] = M_{0,n-1} ,
1228 * m[i][1] = M_{i, i } for i=0,1,...,n-1
1229 * m[i][2] = M_{i,i+1} for i=0,1,...,n-2 and m[n-1][2] = M_{n-1,0}.
1230 * M should be symmetric (m[i+1][0]=m[i][2]) and positiv definite.
1231 * The size of the system is given in n (n>=1).
1233 * In the first procedure the Cholesky decomposition M = C^T * D * C
1234 * (C is upper triangle with unit diagonal, D is diagonal) is calculated.
1235 * Return TRUE if decomposition exist.
1238 solve_cubic_1(tri_diag m[], int n)
1241 double m_ij, m_n, m_nn, d;
1244 return FALSE; /* Dimension should be at least 1 */
1246 d = m[0][1]; /* D_{0,0} = M_{0,0} */
1248 return FALSE; /* M (or D) should be positiv definite */
1249 m_n = m[0][0]; /* M_{0,n-1} */
1250 m_nn = m[n - 1][1]; /* M_{n-1,n-1} */
1251 for (i = 0; i < n - 2; i++) {
1252 m_ij = m[i][2]; /* M_{i,1} */
1253 m[i][2] = m_ij / d; /* C_{i,i+1} */
1254 m[i][0] = m_n / d; /* C_{i,n-1} */
1255 m_nn -= m[i][0] * m_n; /* to get C_{n-1,n-1} */
1256 m_n = -m[i][2] * m_n; /* to get C_{i+1,n-1} */
1257 d = m[i + 1][1] - m[i][2] * m_ij; /* D_{i+1,i+1} */
1259 return FALSE; /* Elements of D should be positiv */
1262 if (n >= 2) { /* Complete last column */
1263 m_n += m[n - 2][2]; /* add M_{n-2,n-1} */
1264 m[n - 2][0] = m_n / d; /* C_{n-2,n-1} */
1265 m[n - 1][1] = d = m_nn - m[n - 2][0] * m_n; /* D_{n-1,n-1} */
1273 * The second procedure solves the linear system, with the Choleky
1274 * decomposition calculated above (in m[][]) and the right side b given
1275 * in x[]. The solution x overwrites the right side in x[].
1278 solve_cubic_2(tri_diag m[], double x[], int n)
1283 /* Division by transpose of C : b = C^{-T} * b */
1285 for (i = 0; i < n - 2; i++) {
1286 x[i + 1] -= m[i][2] * x[i]; /* C_{i,i+1} * x_{i} */
1287 x_n -= m[i][0] * x[i]; /* C_{i,n-1} * x_{i} */
1290 x[n - 1] = x_n - m[n - 2][0] * x[n - 2]; /* C_{n-2,n-1} * x_{n-1} */
1292 /* Division by D: b = D^{-1} * b */
1293 for (i = 0; i < n; i++)
1296 /* Division by C: b = C^{-1} * b */
1299 x[n - 2] -= m[n - 2][0] * x_n; /* C_{n-2,n-1} * x_{n-1} */
1300 for (i = n - 3; i >= 0; i--) {
1301 /* C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1} */
1302 x[i] -= m[i][2] * x[i + 1] + m[i][0] * x_n;
1308 * Solve tri diagonal linear system equation. The tri diagonal matrix is
1309 * defined via matrix M, right side is r, and solution X i.e. M * X = R.
1310 * Size of system given in n. Return TRUE if solution exist.
1312 /* not used any more in "contour.c", but in "spline.c" (21. Dec. 1995) ! */
1315 solve_tri_diag(tri_diag m[], double r[], double x[], int n)
1320 for (i = 1; i < n; i++) { /* Eliminate element m[i][i-1] (lower diagonal). */
1321 if (m[i - 1][1] == 0)
1323 t = m[i][0] / m[i - 1][1]; /* Find ratio between the two lines. */
1324 /* m[i][0] = m[i][0] - m[i-1][1] * t; */
1325 /* m[i][0] is not used any more (and set to 0 in the above line) */
1326 m[i][1] = m[i][1] - m[i - 1][2] * t;
1327 r[i] = r[i] - r[i - 1] * t;
1329 /* Now do back subtitution - update the solution vector X: */
1330 if (m[n - 1][1] == 0)
1332 x[n - 1] = r[n - 1] / m[n - 1][1]; /* Find last element. */
1333 for (i = n - 2; i >= 0; i--) {
1336 x[i] = (r[i] - x[i + 1] * m[i][2]) / m[i][1];
1342 * Generate a Bspline curve defined by all the points given in linked list p:
1343 * Algorithm: using deBoor algorithm
1344 * Note: if Curvekind is open contour than Open end knot vector is assumed,
1345 * else (closed contour) Float end knot vector is assumed.
1346 * It is assumed that num_of_points is at least 2, and order of Bspline is less
1347 * than num_of_points!
1351 cntr_struct *p_cntr,
1354 TBOOLEAN contr_isclosed)
1356 int knot_index = 0, pts_count = 1;
1357 double dt, t, next_t, t_min, t_max, x, y;
1358 cntr_struct *pc_temp = p_cntr, *pc_tail = NULL;
1360 /* If the contour is Closed one we must update few things:
1361 * 1. Make the list temporary circular, so we can close the contour.
1362 * 2. Update num_of_points - increase it by "order-1" so contour will be
1363 * closed. This will evaluate order more sections to close it!
1365 if (contr_isclosed) {
1367 while (pc_tail->next)
1368 pc_tail = pc_tail->next; /* Find last point. */
1370 /* test if first and last point are equal */
1371 if (fuzzy_equal(pc_tail, p_cntr)) {
1372 /* Close contour list - make it circular. */
1373 pc_tail->next = p_cntr->next;
1374 num_of_points += order - 1;
1376 pc_tail->next = p_cntr;
1377 num_of_points += order;
1380 /* Find first (t_min) and last (t_max) t value to eval: */
1381 t = t_min = fetch_knot(contr_isclosed, num_of_points, order, order);
1382 t_max = fetch_knot(contr_isclosed, num_of_points, order, num_of_points);
1383 next_t = t_min + 1.0;
1385 dt = 1.0 / contour_pts; /* Number of points per one section. */
1390 pc_temp = pc_temp->next; /* Next order ctrl. pt. to blend. */
1394 eval_bspline(t, pc_temp, num_of_points, order, knot_index,
1395 contr_isclosed, &x, &y); /* Next pt. */
1396 add_cntr_point(x, y);
1398 /* As we might have some real number round off problems we do */
1399 /* the last point outside the loop */
1400 if (pts_count == contour_pts * (num_of_points - order) + 1)
1405 /* Now do the last point */
1406 eval_bspline(t_max - EPSILON, pc_temp, num_of_points, order, knot_index,
1407 contr_isclosed, &x, &y);
1408 add_cntr_point(x, y); /* Complete the contour. */
1410 if (contr_isclosed) /* Update list - un-circular it. */
1411 pc_tail->next = NULL;
1415 * The routine to evaluate the B-spline value at point t using knot vector
1416 * from function fetch_knot(), and the control points p_cntr.
1417 * Returns (x, y) of approximated B-spline. Note that p_cntr points on the
1418 * first control point to blend with. The B-spline is of order order.
1423 cntr_struct *p_cntr,
1424 int num_of_points, int order, int j,
1425 TBOOLEAN contr_isclosed,
1426 double *x, double *y)
1429 double ti, tikp, *dx, *dy; /* Copy p_cntr into it to make it faster. */
1431 dx = gp_alloc((order + j) * sizeof(double), "contour b_spline");
1432 dy = gp_alloc((order + j) * sizeof(double), "contour b_spline");
1434 /* Set the dx/dy - [0] iteration step, control points (p==0 iterat.): */
1435 for (i = j - order; i <= j; i++) {
1438 p_cntr = p_cntr->next;
1441 for (p = 1; p <= order; p++) { /* Iteration (b-spline level) counter. */
1442 for (i = j; i >= j - order + p; i--) { /* Control points indexing. */
1443 ti = fetch_knot(contr_isclosed, num_of_points, order, i);
1444 tikp = fetch_knot(contr_isclosed, num_of_points, order, i + order + 1 - p);
1445 if (ti == tikp) { /* Should not be a problems but how knows... */
1447 dx[i] = dx[i] * (t - ti) / (tikp - ti) + /* Calculate x. */
1448 dx[i - 1] * (tikp - t) / (tikp - ti);
1449 dy[i] = dy[i] * (t - ti) / (tikp - ti) + /* Calculate y. */
1450 dy[i - 1] * (tikp - t) / (tikp - ti);
1461 * Routine to get the i knot from uniform knot vector. The knot vector
1462 * might be float (Knot(i) = i) or open (where the first and last "order"
1463 * knots are equal). contr_isclosed determines knot kind - open contour means
1464 * open knot vector, and closed contour selects float knot vector.
1465 * Note the knot vector is not exist and this routine simulates it existance
1466 * Also note the indexes for the knot vector starts from 0.
1469 fetch_knot(TBOOLEAN contr_isclosed, int num_of_points, int order, int i)
1471 if(! contr_isclosed) {
1474 else if (i <= num_of_points)
1475 return (double) (i - order);
1477 return (double) (num_of_points - order);