#ifndef lint static char *RCSid = "$Id: contour.c,v 1.9 1998/12/03 22:23:46 lhecking Exp $"; #endif /* GNUPLOT - contour.c */ /*[ * Copyright 1986 - 1993, 1998 Thomas Williams, Colin Kelley * * Permission to use, copy, and distribute this software and its * documentation for any purpose with or without fee is hereby granted, * provided that the above copyright notice appear in all copies and * that both that copyright notice and this permission notice appear * in supporting documentation. * * Permission to modify the software is granted, but not the right to * distribute the complete modified source code. Modifications are to * be distributed as patches to the released version. Permission to * distribute binaries produced by compiling modified sources is granted, * provided you * 1. distribute the corresponding source modifications from the * released version in the form of a patch file along with the binaries, * 2. add special version identification to distinguish your version * in addition to the base release version number, * 3. provide your name and address as the primary contact for the * support of your modified version, and * 4. retain our contact information in regard to use of the base * software. * Permission to distribute the released version of the source code along * with corresponding source modifications in the form of a patch file is * granted with same provisions 2 through 4 for binary distributions. * * This software is provided "as is" without express or implied warranty * to the extent permitted by applicable law. ]*/ /* * AUTHORS * * Original Software: * Gershon Elber * * Improvements to the numerical algorithms: * Hans-Martin Keller, 1995,1997 (hkeller@gwdg.de) * */ #include "plot.h" #include "setshow.h" #define DEFAULT_NUM_APPROX_PTS 5 #define DEFAULT_BSPLINE_ORDER 3 #define MAX_NUM_APPROX_PTS 100 #define MAX_BSPLINE_ORDER 10 /* ?? Not used ?? */ /* for some reason these symbols are also defined in plot.h under different */ /* names */ #define INTERP_NOTHING CONTOUR_KIND_LINEAR /* Kind of interpolations on contours. */ #define INTERP_CUBIC CONTOUR_KIND_CUBIC_SPL /* Cubic spline interp. */ #define APPROX_BSPLINE CONTOUR_KIND_BSPLINE /* Bspline interpolation. */ #define ACTIVE 1 /* Status of edges at certain Z level. */ #define INACTIVE 2 #define INNER_MESH 1 /* position of edge in mesh */ #define BOUNDARY 2 #define DIAGONAL 3 #define OPEN_CONTOUR 1 /* Contour kinds. */ #define CLOSED_CONTOUR 2 #define EPSILON 1e-5 /* Used to decide if two float are equal. */ #ifndef TRUE #define TRUE -1 #define FALSE 0 #endif #define MAX_POINTS_PER_CNTR 100 #define ABS(x) ((x) > 0 ? (x) : (-(x))) #define SQR(x) ((x) * (x)) /* * struct vrtx_struct { * double X, Y, Z; * struct vrtx_struct *next; * }; * * replaced by 'struct coordinate GPHUGE ', see plot.h (HMK 1997) */ struct edge_struct { struct poly_struct *poly[2]; /* Each edge belongs to up to 2 polygons */ struct coordinate GPHUGE *vertex[2]; /* The two extreme points of this edge. */ struct edge_struct *next; /* To chain lists */ int status, /* Status flag to mark edges in scanning at certain Z level. */ position; /* position in mesh: INNER_MESH, BOUNDARY or DIAGONNAL. */ }; struct poly_struct { struct edge_struct *edge[3]; /* As we do triangolation here... */ struct poly_struct *next; /* To chain lists. */ }; struct cntr_struct { /* Contours are saved using this struct list. */ double X, Y; /* The coordinates of this vertex. */ struct cntr_struct *next; /* To chain lists. */ }; static struct gnuplot_contours *contour_list = NULL; static double crnt_cntr[MAX_POINTS_PER_CNTR * 2]; static int crnt_cntr_pt_index = 0; static double contour_level = 0.0; static int num_approx_pts = DEFAULT_NUM_APPROX_PTS; /* # pts per approx/inter. */ static int bspline_order = DEFAULT_BSPLINE_ORDER; /* Bspline order to use. */ static int interp_kind = INTERP_NOTHING; /* Linear, Cubic interp., Bspline. */ static double x_min, y_min, z_min; /* Minimum values of x, y, and z */ static double x_max, y_max, z_max; /* Maximum values of x, y, and z */ static void add_cntr_point __PROTO((double x, double y)); static void end_crnt_cntr __PROTO((void)); static void gen_contours __PROTO((struct edge_struct * p_edges, double z_level, double xx_min, double xx_max, double yy_min, double yy_max)); static int update_all_edges __PROTO((struct edge_struct * p_edges, double z_level)); static struct cntr_struct *gen_one_contour __PROTO(( struct edge_struct * p_edges, double z_level, int *contr_kind, int *num_active)); static struct cntr_struct *trace_contour __PROTO(( struct edge_struct * pe_start, double z_level, int *num_active, int contr_kind)); static struct cntr_struct *update_cntr_pt __PROTO((struct edge_struct * p_edge, double z_level)); static int fuzzy_equal __PROTO((struct cntr_struct * p_cntr1, struct cntr_struct * p_cntr2)); static void gen_triangle __PROTO((int num_isolines, struct iso_curve * iso_lines, struct poly_struct ** p_polys, struct edge_struct ** p_edges)); static void calc_min_max __PROTO((int num_isolines, struct iso_curve * iso_lines, double *xx_min, double *yy_min, double *zz_min, double *xx_max, double *yy_max, double *zz_max)); static struct edge_struct *add_edge __PROTO((struct coordinate GPHUGE * point0, struct coordinate GPHUGE * point1, struct edge_struct ** p_edge, struct edge_struct ** pe_tail)); static struct poly_struct *add_poly __PROTO((struct edge_struct * edge0, struct edge_struct * edge1, struct edge_struct * edge2, struct poly_struct ** p_poly, struct poly_struct ** pp_tail)); static void put_contour __PROTO((struct cntr_struct * p_cntr, double z_level, double xx_min, double xx_max, double yy_min, double yy_max, int contr_kind)); static void put_contour_nothing __PROTO((struct cntr_struct * p_cntr)); static int chk_contour_kind __PROTO((struct cntr_struct * p_cntr, int contr_kind)); static void put_contour_cubic __PROTO((struct cntr_struct * p_cntr, double z_level, double xx_min, double xx_max, double yy_min, double yy_max, int contr_kind)); static void put_contour_bspline __PROTO((struct cntr_struct * p_cntr, double z_level, double xx_min, double xx_max, double yy_min, double yy_max, int contr_kind)); static void free_contour __PROTO((struct cntr_struct * p_cntr)); static int count_contour __PROTO((struct cntr_struct * p_cntr)); static int gen_cubic_spline __PROTO((int num_pts, struct cntr_struct * p_cntr, double d2x[], double d2y[], double delta_t[], int contr_kind, double unit_x, double unit_y)); static void intp_cubic_spline __PROTO((int n, struct cntr_struct * p_cntr, double d2x[], double d2y[], double delta_t[], int n_intpol)); static int solve_cubic_1 __PROTO((tri_diag m[], int n)); static void solve_cubic_2 __PROTO((tri_diag m[], double x[], int n)); /* * static int solve_tri_diag __PROTO((tri_diag m[], double r[], double x[], * int n)); see "protos.h" */ static void gen_bspline_approx __PROTO((struct cntr_struct * p_cntr, int num_of_points, int order, int contr_kind)); static void eval_bspline __PROTO((double t, struct cntr_struct * p_cntr, int num_of_points, int order, int j, int contr_kind, double *x, double *y)); static double fetch_knot __PROTO((int contr_kind, int num_of_points, int order, int i)); /* * Entry routine to this whole set of contouring module. */ struct gnuplot_contours *contour(num_isolines, iso_lines, ZLevels, approx_pts, int_kind, order1, contour_levels_kind, cont_levels_list) int num_isolines; struct iso_curve *iso_lines; int ZLevels, approx_pts, int_kind, order1, contour_levels_kind; double *cont_levels_list; { int i; int num_of_z_levels; /* # Z contour levels. */ struct poly_struct *p_polys, *p_poly; struct edge_struct *p_edges, *p_edge; double z = 0, dz = 0; struct gnuplot_contours *save_contour_list; num_of_z_levels = ZLevels; num_approx_pts = approx_pts; bspline_order = order1 - 1; interp_kind = int_kind; contour_list = NULL; /* * Calculate min/max values : */ calc_min_max(num_isolines, iso_lines, &x_min, &y_min, &z_min, &x_max, &y_max, &z_max); /* * Generate list of edges (p_edges) and list of triangles (p_polys): */ gen_triangle(num_isolines, iso_lines, &p_polys, &p_edges); crnt_cntr_pt_index = 0; if (contour_levels_kind == LEVELS_AUTO) { dz = fabs(z_max - z_min); if (dz == 0) return NULL; /* empty z range ? */ dz = set_tic(log10(dz), ((int) ZLevels + 1) * 2); z = floor(z_min / dz) * dz; num_of_z_levels = (int) floor((z_max - z) / dz); } for (i = 0; i < num_of_z_levels; i++) { switch (contour_levels_kind) { case LEVELS_AUTO: z += dz; break; case LEVELS_INCREMENTAL: z = cont_levels_list[0] + i * cont_levels_list[1]; break; case LEVELS_DISCRETE: z = is_log_z ? log(cont_levels_list[i]) / log_base_log_z : cont_levels_list[i]; break; } contour_level = z; save_contour_list = contour_list; gen_contours(p_edges, z, x_min, x_max, y_min, y_max); if (contour_list != save_contour_list) { contour_list->isNewLevel = 1; sprintf(contour_list->label, contour_format, is_log_z ? pow(base_log_z, z) : z); } } /* Free all contouring related temporary data. */ while (p_polys) { p_poly = p_polys->next; free(p_polys); p_polys = p_poly; } while (p_edges) { p_edge = p_edges->next; free(p_edges); p_edges = p_edge; } return contour_list; } /* * Adds another point to the currently build contour. */ static void add_cntr_point(x, y) double x, y; { int index; if (crnt_cntr_pt_index >= MAX_POINTS_PER_CNTR - 1) { index = crnt_cntr_pt_index - 1; end_crnt_cntr(); crnt_cntr[0] = crnt_cntr[index * 2]; crnt_cntr[1] = crnt_cntr[index * 2 + 1]; crnt_cntr_pt_index = 1; /* Keep the last point as first of this one. */ } crnt_cntr[crnt_cntr_pt_index * 2] = x; crnt_cntr[crnt_cntr_pt_index * 2 + 1] = y; crnt_cntr_pt_index++; } /* * Done with current contour - create gnuplot data structure for it. */ static void end_crnt_cntr() { int i; struct gnuplot_contours *cntr = (struct gnuplot_contours *) gp_alloc((unsigned long) sizeof(struct gnuplot_contours), "gnuplot_contour"); cntr->coords = (struct coordinate GPHUGE *) gp_alloc((unsigned long) sizeof(struct coordinate) * (unsigned long) crnt_cntr_pt_index, "contour coords"); for (i = 0; i < crnt_cntr_pt_index; i++) { cntr->coords[i].x = crnt_cntr[i * 2]; cntr->coords[i].y = crnt_cntr[i * 2 + 1]; cntr->coords[i].z = contour_level; } cntr->num_pts = crnt_cntr_pt_index; cntr->next = contour_list; contour_list = cntr; contour_list->isNewLevel = 0; crnt_cntr_pt_index = 0; } /* * Generates all contours by tracing the intersecting triangles. */ static void gen_contours(p_edges, z_level, xx_min, xx_max, yy_min, yy_max) struct edge_struct *p_edges; double z_level, xx_min, xx_max, yy_min, yy_max; { int num_active, /* Number of edges marked ACTIVE. */ contr_kind; /* One of OPEN_CONTOUR, CLOSED_CONTOUR. */ struct cntr_struct *p_cntr; num_active = update_all_edges(p_edges, z_level); /* Do pass 1. */ contr_kind = OPEN_CONTOUR; /* Start to look for contour on boundaries. */ while (num_active > 0) { /* Do Pass 2. */ /* Generate One contour (and update MumActive as needed): */ p_cntr = gen_one_contour(p_edges, z_level, &contr_kind, &num_active); /* Emit it in requested format: */ put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind); } } /* * Does pass 1, or marks the edges which are active (crosses this z_level) * as ACTIVE, and the others as INACTIVE: * Returns number of active edges (marked ACTIVE). */ static int update_all_edges(p_edges, z_level) struct edge_struct *p_edges; double z_level; { int count = 0; while (p_edges) { /* use the same test at both vertices to avoid roundoff errors */ if ((p_edges->vertex[0]->z >= z_level) != (p_edges->vertex[1]->z >= z_level)) { p_edges->status = ACTIVE; count++; } else p_edges->status = INACTIVE; p_edges = p_edges->next; } return count; } /* * Does pass 2, or find one complete contour out of the triangulation * data base: * Returns a pointer to the contour (as linked list), contr_kind is set to * one of OPEN_CONTOUR, CLOSED_CONTOUR, and num_active is updated. */ static struct cntr_struct *gen_one_contour(p_edges, z_level, contr_kind, num_active) struct edge_struct *p_edges; /* list of edges input */ double z_level; /* Z level of contour input */ int *contr_kind; /* OPEN_ or CLOESED_CONTOUR in/out */ int *num_active; /* number of active edges in/out */ { struct edge_struct *pe_temp; if (*contr_kind == OPEN_CONTOUR) { /* Look for something to start with on boundary: */ pe_temp = p_edges; while (pe_temp) { if ((pe_temp->status == ACTIVE) && (pe_temp->position == BOUNDARY)) break; pe_temp = pe_temp->next; } if (!pe_temp) *contr_kind = CLOSED_CONTOUR; /* No more contours on boundary. */ else { return trace_contour(pe_temp, z_level, num_active, *contr_kind); } } if (*contr_kind == CLOSED_CONTOUR) { /* Look for something to start with inside: */ pe_temp = p_edges; while (pe_temp) { if ((pe_temp->status == ACTIVE) && (!(pe_temp->position == BOUNDARY))) break; pe_temp = pe_temp->next; } if (!pe_temp) { *num_active = 0; fprintf(stderr, "gen_one_contour: no contour found\n"); return NULL; } else { *contr_kind = CLOSED_CONTOUR; return trace_contour(pe_temp, z_level, num_active, *contr_kind); } } return NULL; /* We should never be here, but lint... */ } /* * Search the data base along a contour starts at the edge pe_start until * a boundary edge is detected or until we close the loop back to pe_start. * Returns a linked list of all the points on the contour * Also decreases num_active by the number of points on contour. */ static struct cntr_struct *trace_contour(pe_start, z_level, num_active, contr_kind) struct edge_struct *pe_start; /* edge to start contour input */ double z_level; /* Z level of contour input */ int *num_active; /* number of active edges in/out */ int contr_kind; /* OPEN_ or CLOESED_CONTOUR input */ { struct cntr_struct *p_cntr, *pc_tail; struct edge_struct *p_edge, *p_next_edge; struct poly_struct *p_poly, *PLastpoly = NULL; int i; p_edge = pe_start; /* first edge to start contour */ /* Generate the header of the contour - the point on pe_start. */ if (contr_kind == OPEN_CONTOUR) { pe_start->status = INACTIVE; (*num_active)--; } if (p_edge->poly[0] || p_edge->poly[1]) { /* more than one point */ p_cntr = pc_tail = update_cntr_pt(pe_start, z_level); /* first point */ do { /* Find polygon to continue (Not where we came from - PLastpoly): */ if (p_edge->poly[0] == PLastpoly) p_poly = p_edge->poly[1]; else p_poly = p_edge->poly[0]; p_next_edge = NULL; /* In case of error, remains NULL. */ for (i = 0; i < 3; i++) /* Test the 3 edges of the polygon: */ if (p_poly->edge[i] != p_edge) if (p_poly->edge[i]->status == ACTIVE) p_next_edge = p_poly->edge[i]; if (!p_next_edge) { /* Error exit */ pc_tail->next = NULL; free_contour(p_cntr); fprintf(stderr, "trace_contour: unexpected end of contour\n"); return NULL; } p_edge = p_next_edge; PLastpoly = p_poly; p_edge->status = INACTIVE; (*num_active)--; /* Do not allocate contour points on diagonal edges */ if (p_edge->position != DIAGONAL) { pc_tail->next = update_cntr_pt(p_edge, z_level); /* Remove nearby points */ if (fuzzy_equal(pc_tail, pc_tail->next)) { free((char *) pc_tail->next); } else pc_tail = pc_tail->next; } } while ((p_edge != pe_start) && (p_edge->position != BOUNDARY)); pc_tail->next = NULL; /* For CLOSED_CONTOUR the first and last point should be equal */ if (pe_start == p_edge) { (p_cntr->X) = (pc_tail->X); (p_cntr->Y) = (pc_tail->Y); } } else { /* only one point, forget it */ p_cntr = NULL; } return p_cntr; } /* * Allocates one contour location and update it to to correct position * according to z_level and edge p_edge. */ static struct cntr_struct *update_cntr_pt(p_edge, z_level) struct edge_struct *p_edge; double z_level; { double t; struct cntr_struct *p_cntr; t = (z_level - p_edge->vertex[0]->z) / (p_edge->vertex[1]->z - p_edge->vertex[0]->z); /* test if t is out of interval [0:1] (should not happen but who knows ...) */ t = (t < 0.0 ? 0.0 : t); t = (t > 1.0 ? 1.0 : t); p_cntr = (struct cntr_struct *) gp_alloc((unsigned long) sizeof(struct cntr_struct), "contour cntr_struct"); p_cntr->X = p_edge->vertex[1]->x * t + p_edge->vertex[0]->x * (1 - t); p_cntr->Y = p_edge->vertex[1]->y * t + p_edge->vertex[0]->y * (1 - t); return p_cntr; } /* * Simple routine to decide if two contour points are equal by * calculating the relative error (< EPSILON). */ static int fuzzy_equal(p_cntr1, p_cntr2) struct cntr_struct *p_cntr1, *p_cntr2; { double unit_x, unit_y; unit_x = ABS(x_max - x_min) + zero; /* reference */ unit_y = ABS(y_max - y_min) + zero; return ( ABS(p_cntr1->X - p_cntr2->X) / unit_x < EPSILON && ABS(p_cntr1->Y - p_cntr2->Y) / unit_y < EPSILON); } /* * Generate the triangles. * Returns the lists (edges & polys) via pointers to their heads. */ static void gen_triangle(num_isolines, iso_lines, p_polys, p_edges) int num_isolines; /* number of iso-lines input */ struct iso_curve *iso_lines; /* iso-lines input */ struct poly_struct **p_polys; /* list of polygons output */ struct edge_struct **p_edges; /* list of edges output */ { int i, j, grid_x_max = iso_lines->p_count; struct edge_struct *p_edge1, *p_edge2, *edge0, *edge1, *edge2, *pe_tail, *pe_tail1, *pe_tail2, *pe_temp; struct poly_struct *pp_tail, *lower_tri, *upper_tri; struct coordinate GPHUGE *p_vrtx1, GPHUGE * p_vrtx2; /* HBB 980308: need to tag *each* of them as GPHUGE! */ (*p_polys) = pp_tail = NULL; /* clear lists */ (*p_edges) = pe_tail = NULL; p_vrtx1 = iso_lines->points; /* first row of vertices */ p_edge1 = pe_tail1 = NULL; /* clear list of edges */ /* Generate edges of first row */ for (j = 0; j < grid_x_max - 1; j++) add_edge(p_vrtx1 + j, p_vrtx1 + j + 1, &p_edge1, &pe_tail1); (*p_edges) = p_edge1; /* update main list */ pe_tail = pe_tail1; /* * Combines vertices to edges and edges to triangles: * ================================================== * The edges are stored in the edge list, referenced by p_edges * (pe_tail points on last edge). * * Temporary pointers: * 1. p_edge2: Top horizontal edge list: ----------------------- 2 * 2. pe_tail: middle edge list: |\ |\ |\ |\ |\ |\ | * | \| \| \| \| \| \| * 3. p_edge1: Bottom horizontal edge list: ----------------------- 1 * * The routine generates two triangle Lower Upper 1 * upper one and lower one: | \ ---- * (Nums. are edges order in polys) 0| \1 0\ |2 * The polygons are stored in the polygon ---- \ | * list (*p_polys) (pp_tail points on 2 * last polygon). * 1 * ----------- * In addition, the edge lists are updated - | \ 0 | * each edge has two pointers on the two | \ | * (one active if boundary) polygons which 0|1 0\1 0|1 * uses it. These two pointer to polygons | \ | * are named: poly[0], poly[1]. The diagram | 1 \ | * on the right show how they are used for the ----------- * upper and lower polygons (INNER_MESH polygons only). 0 */ for (i = 1; i < num_isolines; i++) { /* Read next column and gen. polys. */ iso_lines = iso_lines->next; p_vrtx2 = iso_lines->points; /* next row of vertices */ p_edge2 = pe_tail2 = NULL; /* clear top horizontal list */ pe_temp = p_edge1; /* pointer in bottom list */ /* * Generate edges and triagles for next row: */ /* generate first vertical edge */ edge2 = add_edge(p_vrtx1, p_vrtx2, p_edges, &pe_tail); for (j = 0; j < grid_x_max - 1; j++) { /* copy vertical edge for lower triangle */ edge0 = edge2; if (pe_temp && pe_temp->vertex[0] == p_vrtx1 + j) { /* test lower edge */ edge2 = pe_temp; pe_temp = pe_temp->next; } else { edge2 = NULL; /* edge is undefined */ } /* generate diagonal edge */ edge1 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j, p_edges, &pe_tail); if (edge1) edge1->position = DIAGONAL; /* generate lower triangle */ lower_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail); /* copy diagonal edge for upper triangle */ edge0 = edge1; /* generate upper edge */ edge1 = add_edge(p_vrtx2 + j, p_vrtx2 + j + 1, &p_edge2, &pe_tail2); /* generate vertical edge */ edge2 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j + 1, p_edges, &pe_tail); /* generate upper triangle */ upper_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail); } if ((*p_edges)) { /* Chain new edges to main list. */ pe_tail->next = p_edge2; pe_tail = pe_tail2; } else { (*p_edges) = p_edge2; pe_tail = pe_tail2; } p_edge1 = p_edge2; p_vrtx1 = p_vrtx2; } /* Update the boundary flag, saved in each edge, and update indexes: */ pe_temp = (*p_edges); while (pe_temp) { if ((!(pe_temp->poly[0])) || (!(pe_temp->poly[1]))) (pe_temp->position) = BOUNDARY; pe_temp = pe_temp->next; } } /* * Calculate minimum and maximum values */ static void calc_min_max(num_isolines, iso_lines, xx_min, yy_min, zz_min, xx_max, yy_max, zz_max) int num_isolines; /* number of iso-lines input */ struct iso_curve *iso_lines; /* iso-lines input */ double *xx_min, *yy_min, *zz_min, *xx_max, *yy_max, *zz_max; /* min/max values in/out */ { int i, j, grid_x_max; struct coordinate GPHUGE *vertex; grid_x_max = iso_lines->p_count; /* number of vertices per iso_line */ (*xx_min) = (*yy_min) = (*zz_min) = VERYLARGE; /* clear min/max values */ (*xx_max) = (*yy_max) = (*zz_max) = -VERYLARGE; for (j = 0; j < num_isolines; j++) { vertex = iso_lines->points; for (i = 0; i < grid_x_max; i++) { if (vertex[i].type != UNDEFINED) { if (vertex[i].x > (*xx_max)) (*xx_max) = vertex[i].x; if (vertex[i].y > (*yy_max)) (*yy_max) = vertex[i].y; if (vertex[i].z > (*zz_max)) (*zz_max) = vertex[i].z; if (vertex[i].x < (*xx_min)) (*xx_min) = vertex[i].x; if (vertex[i].y < (*yy_min)) (*yy_min) = vertex[i].y; if (vertex[i].z < (*zz_min)) (*zz_min) = vertex[i].z; } } iso_lines = iso_lines->next; } /* * fprintf(stderr," x: %g, %g\n", (*xx_min), (*xx_max)); * fprintf(stderr," y: %g, %g\n", (*yy_min), (*yy_max)); * fprintf(stderr," z: %g, %g\n", (*zz_min), (*zz_max)); */ } /* * Generate new edge and append it to list, but only if both vertices are * defined. The list is referenced by p_edge and pe_tail (p_edge points on * first edge and pe_tail on last one). * Note, the list may be empty (pe_edge==pe_tail==NULL) on entry and exit. */ static struct edge_struct *add_edge(point0, point1, p_edge, pe_tail) struct coordinate GPHUGE * point0; /* 2 vertices input */ struct coordinate GPHUGE * point1; struct edge_struct **p_edge, **pe_tail; /* pointers to edge list in/out */ { struct edge_struct *pe_temp = NULL; if (point0->type != UNDEFINED && point1->type != UNDEFINED) { pe_temp = (struct edge_struct *) gp_alloc((unsigned long) sizeof(struct edge_struct), "contour edge"); pe_temp->poly[0] = NULL; /* clear links */ pe_temp->poly[1] = NULL; pe_temp->vertex[0] = point0; /* First vertex of edge. */ pe_temp->vertex[1] = point1; /* Second vertex of edge. */ pe_temp->next = NULL; pe_temp->position = INNER_MESH; /* default position in mesh */ if ((*pe_tail)) { (*pe_tail)->next = pe_temp; /* Stick new record as last one. */ } else { (*p_edge) = pe_temp; /* start new list if empty */ } (*pe_tail) = pe_temp; /* continue to last record. */ } return pe_temp; /* returns NULL, if no edge allocated */ } /* * Generate new triangle and append it to list, but only if all edges are defined. * The list is referenced by p_poly and pp_tail (p_poly points on first ploygon * and pp_tail on last one). * Note, the list may be empty (pe_ploy==pp_tail==NULL) on entry and exit. */ static struct poly_struct *add_poly(edge0, edge1, edge2, p_poly, pp_tail) struct edge_struct *edge0, *edge1, *edge2; /* 3 edges input */ struct poly_struct **p_poly, **pp_tail; /* pointers to polygon list in/out */ { struct poly_struct *pp_temp = NULL; if (edge0 && edge1 && edge2) { pp_temp = (struct poly_struct *) gp_alloc((unsigned long) sizeof(struct poly_struct), "contour polygon"); pp_temp->edge[0] = edge0; /* First edge of triangle */ pp_temp->edge[1] = edge1; /* Second one */ pp_temp->edge[2] = edge2; /* Third one */ pp_temp->next = NULL; if (edge0->poly[0]) /* update edge0 */ edge0->poly[1] = pp_temp; else edge0->poly[0] = pp_temp; if (edge1->poly[0]) /* update edge1 */ edge1->poly[1] = pp_temp; else edge1->poly[0] = pp_temp; if (edge2->poly[0]) /* update edge2 */ edge2->poly[1] = pp_temp; else edge2->poly[0] = pp_temp; if ((*pp_tail)) /* Stick new record as last one. */ (*pp_tail)->next = pp_temp; else (*p_poly) = pp_temp; /* start new list if empty */ (*pp_tail) = pp_temp; /* continue to last record. */ } return pp_temp; /* returns NULL, if no edge allocated */ } /* * Calls the (hopefully) desired interpolation/approximation routine. */ static void put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind) struct cntr_struct *p_cntr; /* contour structure input */ double z_level, /* Z level of contour input */ xx_min, xx_max, yy_min, yy_max; /* minimum/maximum values input */ int contr_kind; /* OPEN_ or CLOESED_CONTOUR input */ { if (!p_cntr) return; /* Nothing to do if it is empty contour. */ switch (interp_kind) { case INTERP_NOTHING: /* No interpolation/approximation. */ put_contour_nothing(p_cntr); break; case INTERP_CUBIC: /* Cubic spline interpolation. */ put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, chk_contour_kind(p_cntr, contr_kind)); break; case APPROX_BSPLINE: /* Bspline approximation. */ put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, chk_contour_kind(p_cntr, contr_kind)); break; } free_contour(p_cntr); } /* * Simply puts contour coordinates in order with no interpolation or * approximation. */ static void put_contour_nothing(p_cntr) struct cntr_struct *p_cntr; { while (p_cntr) { add_cntr_point(p_cntr->X, p_cntr->Y); p_cntr = p_cntr->next; } end_crnt_cntr(); } /* * for some reason contours are never flagged as CLOSED_CONTOUR * if first point == last point, set flag accordingly * */ static int chk_contour_kind(p_cntr, contr_kind) struct cntr_struct *p_cntr; int contr_kind; { struct cntr_struct *pc_tail = NULL; int current_contr_kind; FPRINTF((stderr, "check_contour_kind: current contr_kind value is %d\n", contr_kind)); current_contr_kind = contr_kind; if (contr_kind != CLOSED_CONTOUR) { pc_tail = p_cntr; while (pc_tail->next) pc_tail = pc_tail->next; /* Find last point. */ /* test if first and last point are equal */ if (fuzzy_equal(pc_tail, p_cntr)) { current_contr_kind = CLOSED_CONTOUR; FPRINTF((stderr, "check_contour_kind: contr_kind changed to %d\n", current_contr_kind)); } } return (current_contr_kind); } /* * Generate a cubic spline curve through the points (x_i,y_i) which are * stored in the linked list p_cntr. * The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the * parameter t is the length of the linear stroke. */ static void put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind) struct cntr_struct *p_cntr; double z_level, xx_min, xx_max, yy_min, yy_max; int contr_kind; { int num_pts, num_intpol; double unit_x, unit_y; /* To define norm (x,y)-plane */ double *delta_t; /* Interval length t_{i+1}-t_i */ double *d2x, *d2y; /* Second derivatives x''(t_i), y''(t_i) */ struct cntr_struct *pc_tail; num_pts = count_contour(p_cntr); /* Number of points in contour. */ pc_tail = p_cntr; /* Find last point. */ while (pc_tail->next) pc_tail = pc_tail->next; if (contr_kind == CLOSED_CONTOUR) { /* Test if first and last point are equal (should be) */ if (!fuzzy_equal(pc_tail, p_cntr)) { pc_tail->next = p_cntr; /* Close contour list - make it circular. */ num_pts++; } } delta_t = (double *) gp_alloc((unsigned long) (sizeof(double) * num_pts), "contour delta_t"); d2x = (double *) gp_alloc((unsigned long) (sizeof(double) * num_pts), "contour d2x"); d2y = (double *) gp_alloc((unsigned long) (sizeof(double) * num_pts), "contour d2y"); /* Width and hight of the grid is used at unit length (2d-norm) */ unit_x = xx_max - x_min; unit_y = yy_max - y_min; unit_x = (unit_x > zero ? unit_x : zero); /* should not be zero */ unit_y = (unit_y > zero ? unit_y : zero); if (num_pts > 2) { /* * Calculate second derivatives d2x[], d2y[] and interval lengths delta_t[]: */ if (!gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, contr_kind, unit_x, unit_y)) { free((char *) delta_t); free((char *) d2x); free((char *) d2y); if (contr_kind == CLOSED_CONTOUR) pc_tail->next = NULL; /* Un-circular list */ return; } } /* If following (num_pts > 1) is TRUE then exactly 2 points in contour. */ else if (num_pts > 1) { /* set all second derivatives to zero, interval length to 1 */ d2x[0] = 0.; d2y[0] = 0.; d2x[1] = 0.; d2y[1] = 0.; delta_t[0] = 1.; } else { /* Only one point ( ?? ) - ignore it. */ free((char *) delta_t); free((char *) d2x); free((char *) d2y); if (contr_kind == CLOSED_CONTOUR) pc_tail->next = NULL; /* Un-circular list */ return; } /* Calculate "num_intpol" interpolated values */ num_intpol = 1 + (num_pts - 1) * num_approx_pts; /* global: num_approx_pts */ intp_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, num_intpol); free((char *) delta_t); free((char *) d2x); free((char *) d2y); if (contr_kind == CLOSED_CONTOUR) pc_tail->next = NULL; /* Un-circular list */ end_crnt_cntr(); } /* * Find Bspline approximation for this data set. * Uses global variable num_approx_pts to determine number of samples per * interval, where the knot vector intervals are assumed to be uniform, and * Global variable bspline_order for the order of Bspline to use. */ static void put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind) struct cntr_struct *p_cntr; double z_level, xx_min, xx_max, yy_min, yy_max; int contr_kind; { int num_pts, order = bspline_order; num_pts = count_contour(p_cntr); /* Number of points in contour. */ if (num_pts < 2) return; /* Can't do nothing if empty or one points! */ /* Order must be less than number of points in curve - fix it if needed. */ if (order > num_pts - 1) order = num_pts - 1; gen_bspline_approx(p_cntr, num_pts, order, contr_kind); end_crnt_cntr(); } /* * Free all elements in the contour list. */ static void free_contour(p_cntr) struct cntr_struct *p_cntr; { struct cntr_struct *pc_temp; while (p_cntr) { pc_temp = p_cntr; p_cntr = p_cntr->next; free((char *) pc_temp); } } /* * Counts number of points in contour. */ static int count_contour(p_cntr) struct cntr_struct *p_cntr; { int count = 0; while (p_cntr) { count++; p_cntr = p_cntr->next; } return count; } /* * Find second derivatives (x''(t_i),y''(t_i)) of cubic spline interpolation * through list of points (x_i,y_i). The parameter t is calculated as the * length of the linear stroke. The number of points must be at least 3. * Note: For CLOSED_CONTOURs the first and last point must be equal. */ static int gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, contr_kind, unit_x, unit_y) int num_pts; /* Number of points (num_pts>=3), input */ struct cntr_struct *p_cntr; /* List of points (x(t_i),y(t_i)), input */ double d2x[], d2y[], /* Second derivatives (x''(t_i),y''(t_i)), output */ delta_t[]; /* List of interval lengths t_{i+1}-t_{i}, output */ int contr_kind; /* CLOSED_CONTOUR or OPEN_CONTOUR, input */ double unit_x, unit_y; /* Unit length in x and y (norm=1), input */ { int n, i; double norm; tri_diag *m; /* The tri-diagonal matrix is saved here. */ struct cntr_struct *pc_temp; m = (tri_diag *) gp_alloc((unsigned long) (sizeof(tri_diag) * num_pts), "contour tridiag m"); /* * Calculate first differences in (d2x[i], d2y[i]) and interval lengths * in delta_t[i]: */ pc_temp = p_cntr; for (i = 0; i < num_pts - 1; i++) { d2x[i] = pc_temp->next->X - pc_temp->X; d2y[i] = pc_temp->next->Y - pc_temp->Y; /* * The Norm of a linear stroke is calculated in "normal coordinates" * and used as interval length: */ delta_t[i] = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y)); d2x[i] /= delta_t[i]; /* first difference, with unit norm: */ d2y[i] /= delta_t[i]; /* || (d2x[i], d2y[i]) || = 1 */ pc_temp = pc_temp->next; } /* * Setup linear System: M * x = b */ n = num_pts - 2; /* Without first and last point */ if (contr_kind == CLOSED_CONTOUR) { /* First and last points must be equal for CLOSED_CONTOURs */ delta_t[num_pts - 1] = delta_t[0]; d2x[num_pts - 1] = d2x[0]; d2y[num_pts - 1] = d2y[0]; n++; /* Add last point (= first point) */ } for (i = 0; i < n; i++) { /* Matrix M, mainly tridiagonal with cyclic second index ("j = j+n mod n") */ m[i][0] = delta_t[i]; /* Off-diagonal element M_{i,i-1} */ m[i][1] = 2. * (delta_t[i] + delta_t[i + 1]); /* M_{i,i} */ m[i][2] = delta_t[i + 1]; /* Off-diagonal element M_{i,i+1} */ /* Right side b_x and b_y */ d2x[i] = (d2x[i + 1] - d2x[i]) * 6.; d2y[i] = (d2y[i + 1] - d2y[i]) * 6.; /* * If the linear stroke shows a cusps of more than 90 degree, the right * side is reduced to avoid oscillations in the spline: */ norm = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y)) / 8.5; if (norm > 1.) { d2x[i] /= norm; d2y[i] /= norm; /* The first derivative will not be continuous */ } } if (contr_kind != CLOSED_CONTOUR) { /* Third derivative is set to zero at both ends */ m[0][1] += m[0][0]; /* M_{0,0} */ m[0][0] = 0.; /* M_{0,n-1} */ m[n - 1][1] += m[n - 1][2]; /* M_{n-1,n-1} */ m[n - 1][2] = 0.; /* M_{n-1,0} */ } /* Solve linear systems for d2x[] and d2y[] */ if (solve_cubic_1(m, n)) { /* Calculate Cholesky decomposition */ solve_cubic_2(m, d2x, n); /* solve M * d2x = b_x */ solve_cubic_2(m, d2y, n); /* solve M * d2y = b_y */ } else { /* Should not happen, but who knows ... */ free((char *) m); return FALSE; } /* Shift all second derivatives one place right and abdate end points */ for (i = n; i > 0; i--) { d2x[i] = d2x[i - 1]; d2y[i] = d2y[i - 1]; } if (contr_kind == CLOSED_CONTOUR) { d2x[0] = d2x[n]; d2y[0] = d2y[n]; } else { d2x[0] = d2x[1]; /* Third derivative is zero in */ d2y[0] = d2y[1]; /* first and last interval */ d2x[n + 1] = d2x[n]; d2y[n + 1] = d2y[n]; } free((char *) m); return TRUE; } /* * Calculate interpolated values of the spline function (defined via p_cntr * and the second derivatives d2x[] and d2y[]). The number of tabulated * values is n. On an equidistant grid n_intpol values are calculated. */ static void intp_cubic_spline(n, p_cntr, d2x, d2y, delta_t, n_intpol) int n; struct cntr_struct *p_cntr; double d2x[], d2y[], delta_t[]; int n_intpol; { double t, t_skip, t_max; double x0, x1, x, y0, y1, y; double d, hx, dx0, dx01, hy, dy0, dy01; int i; /* The length of the total interval */ t_max = 0.; for (i = 0; i < n - 1; i++) t_max += delta_t[i]; /* The distance between interpolated points */ t_skip = (1. - 1e-7) * t_max / (n_intpol - 1); t = 0.; /* Parameter value */ x1 = p_cntr->X; y1 = p_cntr->Y; add_cntr_point(x1, y1); /* First point. */ t += t_skip; for (i = 0; i < n - 1; i++) { p_cntr = p_cntr->next; d = delta_t[i]; /* Interval length */ x0 = x1; y0 = y1; x1 = p_cntr->X; y1 = p_cntr->Y; hx = (x1 - x0) / d; hy = (y1 - y0) / d; dx0 = (d2x[i + 1] + 2 * d2x[i]) / 6.; dy0 = (d2y[i + 1] + 2 * d2y[i]) / 6.; dx01 = (d2x[i + 1] - d2x[i]) / (6. * d); dy01 = (d2y[i + 1] - d2y[i]) / (6. * d); while (t <= delta_t[i]) { /* t in current interval ? */ x = x0 + t * (hx + (t - d) * (dx0 + t * dx01)); y = y0 + t * (hy + (t - d) * (dy0 + t * dy01)); add_cntr_point(x, y); /* next point. */ t += t_skip; } t -= delta_t[i]; /* Parameter t relative to start of next interval */ } } /* * The following two procedures solve the special linear system which arise * in cubic spline interpolation. If x is assumed cyclic ( x[i]=x[n+i] ) the * equations can be written as (i=0,1,...,n-1): * m[i][0] * x[i-1] + m[i][1] * x[i] + m[i][2] * x[i+1] = b[i] . * In matrix notation one gets M * x = b, where the matrix M is tridiagonal * with additional elements in the upper right and lower left position: * m[i][0] = M_{i,i-1} for i=1,2,...,n-1 and m[0][0] = M_{0,n-1} , * m[i][1] = M_{i, i } for i=0,1,...,n-1 * m[i][2] = M_{i,i+1} for i=0,1,...,n-2 and m[n-1][2] = M_{n-1,0}. * M should be symmetric (m[i+1][0]=m[i][2]) and positiv definite. * The size of the system is given in n (n>=1). * * In the first procedure the Cholesky decomposition M = C^T * D * C * (C is upper triangle with unit diagonal, D is diagonal) is calculated. * Return TRUE if decomposition exist. */ static int solve_cubic_1(m, n) tri_diag m[]; int n; { int i; double m_ij, m_n, m_nn, d; if (n < 1) return FALSE; /* Dimension should be at least 1 */ d = m[0][1]; /* D_{0,0} = M_{0,0} */ if (d <= 0.) return FALSE; /* M (or D) should be positiv definite */ m_n = m[0][0]; /* M_{0,n-1} */ m_nn = m[n - 1][1]; /* M_{n-1,n-1} */ for (i = 0; i < n - 2; i++) { m_ij = m[i][2]; /* M_{i,1} */ m[i][2] = m_ij / d; /* C_{i,i+1} */ m[i][0] = m_n / d; /* C_{i,n-1} */ m_nn -= m[i][0] * m_n; /* to get C_{n-1,n-1} */ m_n = -m[i][2] * m_n; /* to get C_{i+1,n-1} */ d = m[i + 1][1] - m[i][2] * m_ij; /* D_{i+1,i+1} */ if (d <= 0.) return FALSE; /* Elements of D should be positiv */ m[i + 1][1] = d; } if (n >= 2) { /* Complete last column */ m_n += m[n - 2][2]; /* add M_{n-2,n-1} */ m[n - 2][0] = m_n / d; /* C_{n-2,n-1} */ m[n - 1][1] = d = m_nn - m[n - 2][0] * m_n; /* D_{n-1,n-1} */ if (d <= 0.) return FALSE; } return TRUE; } /* * The second procedure solves the linear system, with the Choleky * decomposition calculated above (in m[][]) and the right side b given * in x[]. The solution x overwrites the right side in x[]. */ static void solve_cubic_2(m, x, n) tri_diag m[]; double x[]; int n; { int i; double x_n; /* Division by transpose of C : b = C^{-T} * b */ x_n = x[n - 1]; for (i = 0; i < n - 2; i++) { x[i + 1] -= m[i][2] * x[i]; /* C_{i,i+1} * x_{i} */ x_n -= m[i][0] * x[i]; /* C_{i,n-1} * x_{i} */ } if (n >= 2) x[n - 1] = x_n - m[n - 2][0] * x[n - 2]; /* C_{n-2,n-1} * x_{n-1} */ /* Division by D: b = D^{-1} * b */ for (i = 0; i < n; i++) x[i] /= m[i][1]; /* Division by C: b = C^{-1} * b */ x_n = x[n - 1]; if (n >= 2) x[n - 2] -= m[n - 2][0] * x_n; /* C_{n-2,n-1} * x_{n-1} */ for (i = n - 3; i >= 0; i--) { /* C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1} */ x[i] -= m[i][2] * x[i + 1] + m[i][0] * x_n; } return; } /* * Solve tri diagonal linear system equation. The tri diagonal matrix is * defined via matrix M, right side is r, and solution X i.e. M * X = R. * Size of system given in n. Return TRUE if solution exist. */ /* not used any more in "contour.c", but in "spline.c" (21. Dec. 1995) ! */ int solve_tri_diag(m, r, x, n) tri_diag m[]; double r[], x[]; int n; { int i; double t; for (i = 1; i < n; i++) { /* Eliminate element m[i][i-1] (lower diagonal). */ if (m[i - 1][1] == 0) return FALSE; t = m[i][0] / m[i - 1][1]; /* Find ratio between the two lines. */ /* m[i][0] = m[i][0] - m[i-1][1] * t; */ /* m[i][0] is not used any more (and set to 0 in the above line) */ m[i][1] = m[i][1] - m[i - 1][2] * t; r[i] = r[i] - r[i - 1] * t; } /* Now do back subtitution - update the solution vector X: */ if (m[n - 1][1] == 0) return FALSE; x[n - 1] = r[n - 1] / m[n - 1][1]; /* Find last element. */ for (i = n - 2; i >= 0; i--) { if (m[i][1] == 0) return FALSE; x[i] = (r[i] - x[i + 1] * m[i][2]) / m[i][1]; } return TRUE; } /* * Generate a Bspline curve defined by all the points given in linked list p: * Algorithm: using deBoor algorithm * Note: if Curvekind is OPEN_CONTOUR than Open end knot vector is assumed, * else (CLOSED_CONTOUR) Float end knot vector is assumed. * It is assumed that num_of_points is at least 2, and order of Bspline is less * than num_of_points! */ static void gen_bspline_approx(p_cntr, num_of_points, order, contr_kind) struct cntr_struct *p_cntr; int num_of_points, order, contr_kind; { int knot_index = 0, pts_count = 1; double dt, t, next_t, t_min, t_max, x, y; struct cntr_struct *pc_temp = p_cntr, *pc_tail = NULL; /* If the contour is Closed one we must update few things: * 1. Make the list temporary circular, so we can close the contour. * 2. Update num_of_points - increase it by "order-1" so contour will be * closed. This will evaluate order more sections to close it! */ if (contr_kind == CLOSED_CONTOUR) { pc_tail = p_cntr; while (pc_tail->next) pc_tail = pc_tail->next; /* Find last point. */ /* test if first and last point are equal */ if (fuzzy_equal(pc_tail, p_cntr)) { /* Close contour list - make it circular. */ pc_tail->next = p_cntr->next; num_of_points += order - 1; } else { pc_tail->next = p_cntr; num_of_points += order; } } /* Find first (t_min) and last (t_max) t value to eval: */ t = t_min = fetch_knot(contr_kind, num_of_points, order, order); t_max = fetch_knot(contr_kind, num_of_points, order, num_of_points); next_t = t_min + 1.0; knot_index = order; dt = 1.0 / num_approx_pts; /* Number of points per one section. */ while (t < t_max) { if (t > next_t) { pc_temp = pc_temp->next; /* Next order ctrl. pt. to blend. */ knot_index++; next_t += 1.0; } eval_bspline(t, pc_temp, num_of_points, order, knot_index, contr_kind, &x, &y); /* Next pt. */ add_cntr_point(x, y); pts_count++; /* As we might have some real number round off problems we do */ /* the last point outside the loop */ if (pts_count == num_approx_pts * (num_of_points - order) + 1) break; t += dt; } /* Now do the last point */ eval_bspline(t_max - EPSILON, pc_temp, num_of_points, order, knot_index, contr_kind, &x, &y); add_cntr_point(x, y); /* Complete the contour. */ if (contr_kind == CLOSED_CONTOUR) /* Update list - un-circular it. */ pc_tail->next = NULL; } /* * The routine to evaluate the B-spline value at point t using knot vector * from function fetch_knot(), and the control points p_cntr. * Returns (x, y) of approximated B-spline. Note that p_cntr points on the * first control point to blend with. The B-spline is of order order. */ static void eval_bspline(t, p_cntr, num_of_points, order, j, contr_kind, x, y) double t; struct cntr_struct *p_cntr; int num_of_points, order, j, contr_kind; double *x, *y; { int i, p; double ti, tikp, *dx, *dy; /* Copy p_cntr into it to make it faster. */ dx = (double *) gp_alloc((unsigned long) (sizeof(double) * (order + j)), "contour b_spline"); dy = (double *) gp_alloc((unsigned long) (sizeof(double) * (order + j)), "contour b_spline"); /* Set the dx/dy - [0] iteration step, control points (p==0 iterat.): */ for (i = j - order; i <= j; i++) { dx[i] = p_cntr->X; dy[i] = p_cntr->Y; p_cntr = p_cntr->next; } for (p = 1; p <= order; p++) { /* Iteration (b-spline level) counter. */ for (i = j; i >= j - order + p; i--) { /* Control points indexing. */ ti = fetch_knot(contr_kind, num_of_points, order, i); tikp = fetch_knot(contr_kind, num_of_points, order, i + order + 1 - p); if (ti == tikp) { /* Should not be a problems but how knows... */ } else { dx[i] = dx[i] * (t - ti) / (tikp - ti) + /* Calculate x. */ dx[i - 1] * (tikp - t) / (tikp - ti); dy[i] = dy[i] * (t - ti) / (tikp - ti) + /* Calculate y. */ dy[i - 1] * (tikp - t) / (tikp - ti); } } } *x = dx[j]; *y = dy[j]; free((char *) dx); free((char *) dy); } /* * Routine to get the i knot from uniform knot vector. The knot vector * might be float (Knot(i) = i) or open (where the first and last "order" * knots are equal). contr_kind determines knot kind - OPEN_CONTOUR means * open knot vector, and CLOSED_CONTOUR selects float knot vector. * Note the knot vector is not exist and this routine simulates it existance * Also note the indexes for the knot vector starts from 0. */ static double fetch_knot(contr_kind, num_of_points, order, i) int contr_kind, num_of_points, order, i; { switch (contr_kind) { case OPEN_CONTOUR: if (i <= order) return 0.0; else if (i <= num_of_points) return (double) (i - order); else return (double) (num_of_points - order); case CLOSED_CONTOUR: return (double) i; default: /* Should never happen */ return 1.0; } #ifdef sequent return 1.0; /* ???? */ #endif }