Annotation of OpenXM_contrib/gnuplot/contour.c, Revision 1.1
1.1 ! maekawa 1: #ifndef lint
! 2: static char *RCSid = "$Id: contour.c,v 1.31 1998/04/14 00:15:15 drd Exp $";
! 3: #endif
! 4:
! 5: /* GNUPLOT - contour.c */
! 6:
! 7: /*[
! 8: * Copyright 1986 - 1993, 1998 Thomas Williams, Colin Kelley
! 9: *
! 10: * Permission to use, copy, and distribute this software and its
! 11: * documentation for any purpose with or without fee is hereby granted,
! 12: * provided that the above copyright notice appear in all copies and
! 13: * that both that copyright notice and this permission notice appear
! 14: * in supporting documentation.
! 15: *
! 16: * Permission to modify the software is granted, but not the right to
! 17: * distribute the complete modified source code. Modifications are to
! 18: * be distributed as patches to the released version. Permission to
! 19: * distribute binaries produced by compiling modified sources is granted,
! 20: * provided you
! 21: * 1. distribute the corresponding source modifications from the
! 22: * released version in the form of a patch file along with the binaries,
! 23: * 2. add special version identification to distinguish your version
! 24: * in addition to the base release version number,
! 25: * 3. provide your name and address as the primary contact for the
! 26: * support of your modified version, and
! 27: * 4. retain our contact information in regard to use of the base
! 28: * software.
! 29: * Permission to distribute the released version of the source code along
! 30: * with corresponding source modifications in the form of a patch file is
! 31: * granted with same provisions 2 through 4 for binary distributions.
! 32: *
! 33: * This software is provided "as is" without express or implied warranty
! 34: * to the extent permitted by applicable law.
! 35: ]*/
! 36:
! 37:
! 38: /*
! 39: * AUTHORS
! 40: *
! 41: * Original Software:
! 42: * Gershon Elber
! 43: *
! 44: * Improvements to the numerical algorithms:
! 45: * Hans-Martin Keller, 1995,1997 (hkeller@gwdg.de)
! 46: *
! 47: */
! 48:
! 49: #include "plot.h"
! 50: #include "setshow.h"
! 51:
! 52: #define DEFAULT_NUM_APPROX_PTS 5
! 53: #define DEFAULT_BSPLINE_ORDER 3
! 54: #define MAX_NUM_APPROX_PTS 100
! 55: #define MAX_BSPLINE_ORDER 10 /* ?? Not used ?? */
! 56:
! 57: /* for some reason these symbols are also defined in plot.h under different */
! 58: /* names */
! 59: #define INTERP_NOTHING CONTOUR_KIND_LINEAR /* Kind of interpolations on contours. */
! 60: #define INTERP_CUBIC CONTOUR_KIND_CUBIC_SPL /* Cubic spline interp. */
! 61: #define APPROX_BSPLINE CONTOUR_KIND_BSPLINE /* Bspline interpolation. */
! 62:
! 63: #define ACTIVE 1 /* Status of edges at certain Z level. */
! 64: #define INACTIVE 2
! 65: #define INNER_MESH 1 /* position of edge in mesh */
! 66: #define BOUNDARY 2
! 67: #define DIAGONAL 3
! 68:
! 69: #define OPEN_CONTOUR 1 /* Contour kinds. */
! 70: #define CLOSED_CONTOUR 2
! 71:
! 72: #define EPSILON 1e-5 /* Used to decide if two float are equal. */
! 73:
! 74: #ifndef TRUE
! 75: #define TRUE -1
! 76: #define FALSE 0
! 77: #endif
! 78:
! 79: #define MAX_POINTS_PER_CNTR 100
! 80:
! 81: #define ABS(x) ((x) > 0 ? (x) : (-(x)))
! 82: #define SQR(x) ((x) * (x))
! 83:
! 84: /*
! 85: * struct vrtx_struct {
! 86: * double X, Y, Z;
! 87: * struct vrtx_struct *next;
! 88: * };
! 89: *
! 90: * replaced by 'struct coordinate GPHUGE ', see plot.h (HMK 1997)
! 91: */
! 92:
! 93: struct edge_struct {
! 94: struct poly_struct *poly[2]; /* Each edge belongs to up to 2 polygons */
! 95: struct coordinate GPHUGE *vertex[2]; /* The two extreme points of this edge. */
! 96: struct edge_struct *next; /* To chain lists */
! 97: int status, /* Status flag to mark edges in scanning at certain Z level. */
! 98: position; /* position in mesh: INNER_MESH, BOUNDARY or DIAGONNAL. */
! 99: };
! 100:
! 101: struct poly_struct {
! 102: struct edge_struct *edge[3]; /* As we do triangolation here... */
! 103: struct poly_struct *next; /* To chain lists. */
! 104: };
! 105:
! 106: struct cntr_struct { /* Contours are saved using this struct list. */
! 107: double X, Y; /* The coordinates of this vertex. */
! 108: struct cntr_struct *next; /* To chain lists. */
! 109: };
! 110:
! 111: static struct gnuplot_contours *contour_list = NULL;
! 112: static double crnt_cntr[MAX_POINTS_PER_CNTR * 2];
! 113: static int crnt_cntr_pt_index = 0;
! 114: static double contour_level = 0.0;
! 115: static int num_approx_pts = DEFAULT_NUM_APPROX_PTS; /* # pts per approx/inter. */
! 116: static int bspline_order = DEFAULT_BSPLINE_ORDER; /* Bspline order to use. */
! 117: static int interp_kind = INTERP_NOTHING; /* Linear, Cubic interp., Bspline. */
! 118: static double x_min, y_min, z_min; /* Minimum values of x, y, and z */
! 119: static double x_max, y_max, z_max; /* Maximum values of x, y, and z */
! 120:
! 121: static void add_cntr_point __PROTO((double x, double y));
! 122: static void end_crnt_cntr __PROTO((void));
! 123: static void gen_contours __PROTO((struct edge_struct * p_edges, double z_level,
! 124: double xx_min, double xx_max, double yy_min, double yy_max));
! 125: static int update_all_edges __PROTO((struct edge_struct * p_edges,
! 126: double z_level));
! 127: static struct cntr_struct *gen_one_contour __PROTO((
! 128: struct edge_struct * p_edges, double z_level, int *contr_kind,
! 129: int *num_active));
! 130: static struct cntr_struct *trace_contour __PROTO((
! 131: struct edge_struct * pe_start, double z_level, int *num_active,
! 132: int contr_kind));
! 133: static struct cntr_struct *update_cntr_pt __PROTO((struct edge_struct * p_edge,
! 134: double z_level));
! 135: static int fuzzy_equal __PROTO((struct cntr_struct * p_cntr1,
! 136: struct cntr_struct * p_cntr2));
! 137:
! 138:
! 139: static void gen_triangle __PROTO((int num_isolines,
! 140: struct iso_curve * iso_lines, struct poly_struct ** p_polys,
! 141: struct edge_struct ** p_edges));
! 142: static void calc_min_max __PROTO((int num_isolines,
! 143: struct iso_curve * iso_lines, double *xx_min, double *yy_min, double *zz_min,
! 144: double *xx_max, double *yy_max, double *zz_max));
! 145: static struct edge_struct *add_edge __PROTO((struct coordinate GPHUGE * point0,
! 146: struct coordinate GPHUGE * point1, struct edge_struct ** p_edge,
! 147: struct edge_struct ** pe_tail));
! 148: static struct poly_struct *add_poly __PROTO((struct edge_struct * edge0,
! 149: struct edge_struct * edge1, struct edge_struct * edge2,
! 150: struct poly_struct ** p_poly, struct poly_struct ** pp_tail));
! 151:
! 152:
! 153: static void put_contour __PROTO((struct cntr_struct * p_cntr, double z_level,
! 154: double xx_min, double xx_max, double yy_min, double yy_max,
! 155: int contr_kind));
! 156: static void put_contour_nothing __PROTO((struct cntr_struct * p_cntr));
! 157: static int chk_contour_kind __PROTO((struct cntr_struct * p_cntr,
! 158: int contr_kind));
! 159: static void put_contour_cubic __PROTO((struct cntr_struct * p_cntr,
! 160: double z_level, double xx_min, double xx_max, double yy_min, double yy_max,
! 161: int contr_kind));
! 162: static void put_contour_bspline __PROTO((struct cntr_struct * p_cntr,
! 163: double z_level, double xx_min, double xx_max, double yy_min, double yy_max,
! 164: int contr_kind));
! 165: static void free_contour __PROTO((struct cntr_struct * p_cntr));
! 166: static int count_contour __PROTO((struct cntr_struct * p_cntr));
! 167: static int gen_cubic_spline __PROTO((int num_pts, struct cntr_struct * p_cntr,
! 168: double d2x[], double d2y[], double delta_t[], int contr_kind,
! 169: double unit_x, double unit_y));
! 170: static void intp_cubic_spline __PROTO((int n, struct cntr_struct * p_cntr,
! 171: double d2x[], double d2y[], double delta_t[], int n_intpol));
! 172: static int solve_cubic_1 __PROTO((tri_diag m[], int n));
! 173: static void solve_cubic_2 __PROTO((tri_diag m[], double x[], int n));
! 174: /*
! 175: * static int solve_tri_diag __PROTO((tri_diag m[], double r[], double x[],
! 176: * int n)); see "protos.h"
! 177: */
! 178: static void gen_bspline_approx __PROTO((struct cntr_struct * p_cntr,
! 179: int num_of_points, int order, int contr_kind));
! 180: static void eval_bspline __PROTO((double t, struct cntr_struct * p_cntr,
! 181: int num_of_points, int order, int j, int contr_kind, double *x,
! 182: double *y));
! 183: static double fetch_knot __PROTO((int contr_kind, int num_of_points,
! 184: int order, int i));
! 185:
! 186: /*
! 187: * Entry routine to this whole set of contouring module.
! 188: */
! 189: struct gnuplot_contours *contour(num_isolines, iso_lines, ZLevels, approx_pts, int_kind, order1, contour_levels_kind, cont_levels_list)
! 190: int num_isolines;
! 191: struct iso_curve *iso_lines;
! 192: int ZLevels, approx_pts, int_kind, order1, contour_levels_kind;
! 193: double *cont_levels_list;
! 194: {
! 195: int i;
! 196: int num_of_z_levels; /* # Z contour levels. */
! 197: struct poly_struct *p_polys, *p_poly;
! 198: struct edge_struct *p_edges, *p_edge;
! 199: double z = 0, dz = 0;
! 200: struct gnuplot_contours *save_contour_list;
! 201:
! 202: num_of_z_levels = ZLevels;
! 203: num_approx_pts = approx_pts;
! 204: bspline_order = order1 - 1;
! 205: interp_kind = int_kind;
! 206:
! 207: contour_list = NULL;
! 208:
! 209: /*
! 210: * Calculate min/max values :
! 211: */
! 212: calc_min_max(num_isolines, iso_lines,
! 213: &x_min, &y_min, &z_min, &x_max, &y_max, &z_max);
! 214:
! 215: /*
! 216: * Generate list of edges (p_edges) and list of triangles (p_polys):
! 217: */
! 218: gen_triangle(num_isolines, iso_lines, &p_polys, &p_edges);
! 219: crnt_cntr_pt_index = 0;
! 220:
! 221: if (contour_levels_kind == LEVELS_AUTO) {
! 222: dz = fabs(z_max - z_min);
! 223: if (dz == 0)
! 224: return NULL; /* empty z range ? */
! 225: dz = set_tic(log10(dz), ((int) ZLevels + 1) * 2);
! 226: z = floor(z_min / dz) * dz;
! 227: num_of_z_levels = (int) floor((z_max - z) / dz);
! 228: }
! 229: for (i = 0; i < num_of_z_levels; i++) {
! 230: switch (contour_levels_kind) {
! 231: case LEVELS_AUTO:
! 232: z += dz;
! 233: break;
! 234: case LEVELS_INCREMENTAL:
! 235: z = cont_levels_list[0] + i * cont_levels_list[1];
! 236: break;
! 237: case LEVELS_DISCRETE:
! 238: z = is_log_z ? log(cont_levels_list[i]) / log_base_log_z : cont_levels_list[i];
! 239: break;
! 240: }
! 241: contour_level = z;
! 242: save_contour_list = contour_list;
! 243: gen_contours(p_edges, z, x_min, x_max, y_min, y_max);
! 244: if (contour_list != save_contour_list) {
! 245: contour_list->isNewLevel = 1;
! 246: sprintf(contour_list->label, contour_format, is_log_z ? pow(base_log_z, z) : z);
! 247: }
! 248: }
! 249:
! 250: /* Free all contouring related temporary data. */
! 251: while (p_polys) {
! 252: p_poly = p_polys->next;
! 253: free(p_polys);
! 254: p_polys = p_poly;
! 255: }
! 256: while (p_edges) {
! 257: p_edge = p_edges->next;
! 258: free(p_edges);
! 259: p_edges = p_edge;
! 260: }
! 261:
! 262: return contour_list;
! 263: }
! 264:
! 265: /*
! 266: * Adds another point to the currently build contour.
! 267: */
! 268: static void add_cntr_point(x, y)
! 269: double x, y;
! 270: {
! 271: int index;
! 272:
! 273: if (crnt_cntr_pt_index >= MAX_POINTS_PER_CNTR - 1) {
! 274: index = crnt_cntr_pt_index - 1;
! 275: end_crnt_cntr();
! 276: crnt_cntr[0] = crnt_cntr[index * 2];
! 277: crnt_cntr[1] = crnt_cntr[index * 2 + 1];
! 278: crnt_cntr_pt_index = 1; /* Keep the last point as first of this one. */
! 279: }
! 280: crnt_cntr[crnt_cntr_pt_index * 2] = x;
! 281: crnt_cntr[crnt_cntr_pt_index * 2 + 1] = y;
! 282: crnt_cntr_pt_index++;
! 283: }
! 284:
! 285: /*
! 286: * Done with current contour - create gnuplot data structure for it.
! 287: */
! 288: static void end_crnt_cntr()
! 289: {
! 290: int i;
! 291: struct gnuplot_contours *cntr = (struct gnuplot_contours *)
! 292: gp_alloc((unsigned long) sizeof(struct gnuplot_contours), "gnuplot_contour");
! 293: cntr->coords = (struct coordinate GPHUGE *)
! 294: gp_alloc((unsigned long) sizeof(struct coordinate)
! 295: * (unsigned long) crnt_cntr_pt_index, "contour coords");
! 296:
! 297: for (i = 0; i < crnt_cntr_pt_index; i++) {
! 298: cntr->coords[i].x = crnt_cntr[i * 2];
! 299: cntr->coords[i].y = crnt_cntr[i * 2 + 1];
! 300: cntr->coords[i].z = contour_level;
! 301: }
! 302: cntr->num_pts = crnt_cntr_pt_index;
! 303:
! 304: cntr->next = contour_list;
! 305: contour_list = cntr;
! 306: contour_list->isNewLevel = 0;
! 307:
! 308: crnt_cntr_pt_index = 0;
! 309: }
! 310:
! 311: /*
! 312: * Generates all contours by tracing the intersecting triangles.
! 313: */
! 314: static void gen_contours(p_edges, z_level, xx_min, xx_max, yy_min, yy_max)
! 315: struct edge_struct *p_edges;
! 316: double z_level, xx_min, xx_max, yy_min, yy_max;
! 317: {
! 318: int num_active, /* Number of edges marked ACTIVE. */
! 319: contr_kind; /* One of OPEN_CONTOUR, CLOSED_CONTOUR. */
! 320: struct cntr_struct *p_cntr;
! 321:
! 322: num_active = update_all_edges(p_edges, z_level); /* Do pass 1. */
! 323:
! 324: contr_kind = OPEN_CONTOUR; /* Start to look for contour on boundaries. */
! 325:
! 326: while (num_active > 0) { /* Do Pass 2. */
! 327: /* Generate One contour (and update MumActive as needed): */
! 328: p_cntr = gen_one_contour(p_edges, z_level, &contr_kind, &num_active);
! 329: /* Emit it in requested format: */
! 330: put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind);
! 331: }
! 332: }
! 333:
! 334: /*
! 335: * Does pass 1, or marks the edges which are active (crosses this z_level)
! 336: * as ACTIVE, and the others as INACTIVE:
! 337: * Returns number of active edges (marked ACTIVE).
! 338: */
! 339: static int update_all_edges(p_edges, z_level)
! 340: struct edge_struct *p_edges;
! 341: double z_level;
! 342: {
! 343: int count = 0;
! 344:
! 345: while (p_edges) {
! 346: /* use the same test at both vertices to avoid roundoff errors */
! 347: if ((p_edges->vertex[0]->z >= z_level) !=
! 348: (p_edges->vertex[1]->z >= z_level)) {
! 349: p_edges->status = ACTIVE;
! 350: count++;
! 351: } else
! 352: p_edges->status = INACTIVE;
! 353: p_edges = p_edges->next;
! 354: }
! 355:
! 356: return count;
! 357: }
! 358:
! 359: /*
! 360: * Does pass 2, or find one complete contour out of the triangulation
! 361: * data base:
! 362: * Returns a pointer to the contour (as linked list), contr_kind is set to
! 363: * one of OPEN_CONTOUR, CLOSED_CONTOUR, and num_active is updated.
! 364: */
! 365: static struct cntr_struct *gen_one_contour(p_edges, z_level, contr_kind, num_active)
! 366: struct edge_struct *p_edges; /* list of edges input */
! 367: double z_level; /* Z level of contour input */
! 368: int *contr_kind; /* OPEN_ or CLOESED_CONTOUR in/out */
! 369: int *num_active; /* number of active edges in/out */
! 370: {
! 371: struct edge_struct *pe_temp;
! 372:
! 373: if (*contr_kind == OPEN_CONTOUR) {
! 374: /* Look for something to start with on boundary: */
! 375: pe_temp = p_edges;
! 376: while (pe_temp) {
! 377: if ((pe_temp->status == ACTIVE) && (pe_temp->position == BOUNDARY))
! 378: break;
! 379: pe_temp = pe_temp->next;
! 380: }
! 381: if (!pe_temp)
! 382: *contr_kind = CLOSED_CONTOUR; /* No more contours on boundary. */
! 383: else {
! 384: return trace_contour(pe_temp, z_level, num_active, *contr_kind);
! 385: }
! 386: }
! 387: if (*contr_kind == CLOSED_CONTOUR) {
! 388: /* Look for something to start with inside: */
! 389: pe_temp = p_edges;
! 390: while (pe_temp) {
! 391: if ((pe_temp->status == ACTIVE) && (!(pe_temp->position == BOUNDARY)))
! 392: break;
! 393: pe_temp = pe_temp->next;
! 394: }
! 395: if (!pe_temp) {
! 396: *num_active = 0;
! 397: fprintf(stderr, "gen_one_contour: no contour found\n");
! 398: return NULL;
! 399: } else {
! 400: *contr_kind = CLOSED_CONTOUR;
! 401: return trace_contour(pe_temp, z_level, num_active, *contr_kind);
! 402: }
! 403: }
! 404: return NULL; /* We should never be here, but lint... */
! 405: }
! 406:
! 407: /*
! 408: * Search the data base along a contour starts at the edge pe_start until
! 409: * a boundary edge is detected or until we close the loop back to pe_start.
! 410: * Returns a linked list of all the points on the contour
! 411: * Also decreases num_active by the number of points on contour.
! 412: */
! 413: static struct cntr_struct *trace_contour(pe_start, z_level, num_active, contr_kind)
! 414: struct edge_struct *pe_start; /* edge to start contour input */
! 415: double z_level; /* Z level of contour input */
! 416: int *num_active; /* number of active edges in/out */
! 417: int contr_kind; /* OPEN_ or CLOESED_CONTOUR input */
! 418: {
! 419: struct cntr_struct *p_cntr, *pc_tail;
! 420: struct edge_struct *p_edge, *p_next_edge;
! 421: struct poly_struct *p_poly, *PLastpoly = NULL;
! 422: int i;
! 423:
! 424: p_edge = pe_start; /* first edge to start contour */
! 425:
! 426: /* Generate the header of the contour - the point on pe_start. */
! 427: if (contr_kind == OPEN_CONTOUR) {
! 428: pe_start->status = INACTIVE;
! 429: (*num_active)--;
! 430: }
! 431: if (p_edge->poly[0] || p_edge->poly[1]) { /* more than one point */
! 432:
! 433: p_cntr = pc_tail = update_cntr_pt(pe_start, z_level); /* first point */
! 434:
! 435: do {
! 436: /* Find polygon to continue (Not where we came from - PLastpoly): */
! 437: if (p_edge->poly[0] == PLastpoly)
! 438: p_poly = p_edge->poly[1];
! 439: else
! 440: p_poly = p_edge->poly[0];
! 441: p_next_edge = NULL; /* In case of error, remains NULL. */
! 442: for (i = 0; i < 3; i++) /* Test the 3 edges of the polygon: */
! 443: if (p_poly->edge[i] != p_edge)
! 444: if (p_poly->edge[i]->status == ACTIVE)
! 445: p_next_edge = p_poly->edge[i];
! 446: if (!p_next_edge) { /* Error exit */
! 447: pc_tail->next = NULL;
! 448: free_contour(p_cntr);
! 449: fprintf(stderr, "trace_contour: unexpected end of contour\n");
! 450: return NULL;
! 451: }
! 452: p_edge = p_next_edge;
! 453: PLastpoly = p_poly;
! 454: p_edge->status = INACTIVE;
! 455: (*num_active)--;
! 456:
! 457: /* Do not allocate contour points on diagonal edges */
! 458: if (p_edge->position != DIAGONAL) {
! 459:
! 460: pc_tail->next = update_cntr_pt(p_edge, z_level);
! 461:
! 462: /* Remove nearby points */
! 463: if (fuzzy_equal(pc_tail, pc_tail->next)) {
! 464:
! 465: free((char *) pc_tail->next);
! 466: } else
! 467: pc_tail = pc_tail->next;
! 468: }
! 469: } while ((p_edge != pe_start) && (p_edge->position != BOUNDARY));
! 470:
! 471: pc_tail->next = NULL;
! 472:
! 473: /* For CLOSED_CONTOUR the first and last point should be equal */
! 474: if (pe_start == p_edge) {
! 475: (p_cntr->X) = (pc_tail->X);
! 476: (p_cntr->Y) = (pc_tail->Y);
! 477: }
! 478: } else { /* only one point, forget it */
! 479: p_cntr = NULL;
! 480: }
! 481:
! 482: return p_cntr;
! 483: }
! 484:
! 485: /*
! 486: * Allocates one contour location and update it to to correct position
! 487: * according to z_level and edge p_edge.
! 488: */
! 489: static struct cntr_struct *update_cntr_pt(p_edge, z_level)
! 490: struct edge_struct *p_edge;
! 491: double z_level;
! 492: {
! 493: double t;
! 494: struct cntr_struct *p_cntr;
! 495:
! 496: t = (z_level - p_edge->vertex[0]->z) /
! 497: (p_edge->vertex[1]->z - p_edge->vertex[0]->z);
! 498:
! 499: /* test if t is out of interval [0:1] (should not happen but who knows ...) */
! 500: t = (t < 0.0 ? 0.0 : t);
! 501: t = (t > 1.0 ? 1.0 : t);
! 502:
! 503: p_cntr = (struct cntr_struct *)
! 504: gp_alloc((unsigned long) sizeof(struct cntr_struct), "contour cntr_struct");
! 505:
! 506: p_cntr->X = p_edge->vertex[1]->x * t +
! 507: p_edge->vertex[0]->x * (1 - t);
! 508: p_cntr->Y = p_edge->vertex[1]->y * t +
! 509: p_edge->vertex[0]->y * (1 - t);
! 510: return p_cntr;
! 511: }
! 512:
! 513: /*
! 514: * Simple routine to decide if two contour points are equal by
! 515: * calculating the relative error (< EPSILON).
! 516: */
! 517: static int fuzzy_equal(p_cntr1, p_cntr2)
! 518: struct cntr_struct *p_cntr1, *p_cntr2;
! 519: {
! 520: double unit_x, unit_y;
! 521: unit_x = ABS(x_max - x_min) + zero; /* reference */
! 522: unit_y = ABS(y_max - y_min) + zero;
! 523: return (
! 524: ABS(p_cntr1->X - p_cntr2->X) / unit_x < EPSILON &&
! 525: ABS(p_cntr1->Y - p_cntr2->Y) / unit_y < EPSILON);
! 526: }
! 527:
! 528: /*
! 529: * Generate the triangles.
! 530: * Returns the lists (edges & polys) via pointers to their heads.
! 531: */
! 532: static void gen_triangle(num_isolines, iso_lines, p_polys, p_edges)
! 533: int num_isolines; /* number of iso-lines input */
! 534: struct iso_curve *iso_lines; /* iso-lines input */
! 535: struct poly_struct **p_polys; /* list of polygons output */
! 536: struct edge_struct **p_edges; /* list of edges output */
! 537: {
! 538: int i, j, grid_x_max = iso_lines->p_count;
! 539: struct edge_struct *p_edge1, *p_edge2, *edge0, *edge1, *edge2, *pe_tail,
! 540: *pe_tail1, *pe_tail2, *pe_temp;
! 541: struct poly_struct *pp_tail, *lower_tri, *upper_tri;
! 542: struct coordinate GPHUGE *p_vrtx1, GPHUGE * p_vrtx2; /* HBB 980308: need to tag *each* of them as GPHUGE! */
! 543:
! 544: (*p_polys) = pp_tail = NULL; /* clear lists */
! 545: (*p_edges) = pe_tail = NULL;
! 546:
! 547: p_vrtx1 = iso_lines->points; /* first row of vertices */
! 548: p_edge1 = pe_tail1 = NULL; /* clear list of edges */
! 549:
! 550: /* Generate edges of first row */
! 551: for (j = 0; j < grid_x_max - 1; j++)
! 552: add_edge(p_vrtx1 + j, p_vrtx1 + j + 1, &p_edge1, &pe_tail1);
! 553:
! 554: (*p_edges) = p_edge1; /* update main list */
! 555: pe_tail = pe_tail1;
! 556:
! 557:
! 558: /*
! 559: * Combines vertices to edges and edges to triangles:
! 560: * ==================================================
! 561: * The edges are stored in the edge list, referenced by p_edges
! 562: * (pe_tail points on last edge).
! 563: *
! 564: * Temporary pointers:
! 565: * 1. p_edge2: Top horizontal edge list: ----------------------- 2
! 566: * 2. pe_tail: middle edge list: |\ |\ |\ |\ |\ |\ |
! 567: * | \| \| \| \| \| \|
! 568: * 3. p_edge1: Bottom horizontal edge list: ----------------------- 1
! 569: *
! 570: * The routine generates two triangle Lower Upper 1
! 571: * upper one and lower one: | \ ----
! 572: * (Nums. are edges order in polys) 0| \1 0\ |2
! 573: * The polygons are stored in the polygon ---- \ |
! 574: * list (*p_polys) (pp_tail points on 2
! 575: * last polygon).
! 576: * 1
! 577: * -----------
! 578: * In addition, the edge lists are updated - | \ 0 |
! 579: * each edge has two pointers on the two | \ |
! 580: * (one active if boundary) polygons which 0|1 0\1 0|1
! 581: * uses it. These two pointer to polygons | \ |
! 582: * are named: poly[0], poly[1]. The diagram | 1 \ |
! 583: * on the right show how they are used for the -----------
! 584: * upper and lower polygons (INNER_MESH polygons only). 0
! 585: */
! 586:
! 587: for (i = 1; i < num_isolines; i++) {
! 588: /* Read next column and gen. polys. */
! 589: iso_lines = iso_lines->next;
! 590:
! 591: p_vrtx2 = iso_lines->points; /* next row of vertices */
! 592: p_edge2 = pe_tail2 = NULL; /* clear top horizontal list */
! 593: pe_temp = p_edge1; /* pointer in bottom list */
! 594:
! 595: /*
! 596: * Generate edges and triagles for next row:
! 597: */
! 598:
! 599: /* generate first vertical edge */
! 600: edge2 = add_edge(p_vrtx1, p_vrtx2, p_edges, &pe_tail);
! 601:
! 602: for (j = 0; j < grid_x_max - 1; j++) {
! 603:
! 604: /* copy vertical edge for lower triangle */
! 605: edge0 = edge2;
! 606:
! 607: if (pe_temp && pe_temp->vertex[0] == p_vrtx1 + j) {
! 608: /* test lower edge */
! 609: edge2 = pe_temp;
! 610: pe_temp = pe_temp->next;
! 611: } else {
! 612: edge2 = NULL; /* edge is undefined */
! 613: }
! 614:
! 615: /* generate diagonal edge */
! 616: edge1 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j, p_edges, &pe_tail);
! 617: if (edge1)
! 618: edge1->position = DIAGONAL;
! 619:
! 620: /* generate lower triangle */
! 621: lower_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
! 622:
! 623: /* copy diagonal edge for upper triangle */
! 624: edge0 = edge1;
! 625:
! 626: /* generate upper edge */
! 627: edge1 = add_edge(p_vrtx2 + j, p_vrtx2 + j + 1, &p_edge2, &pe_tail2);
! 628:
! 629: /* generate vertical edge */
! 630: edge2 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j + 1, p_edges, &pe_tail);
! 631:
! 632: /* generate upper triangle */
! 633: upper_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
! 634: }
! 635:
! 636: if ((*p_edges)) { /* Chain new edges to main list. */
! 637: pe_tail->next = p_edge2;
! 638: pe_tail = pe_tail2;
! 639: } else {
! 640: (*p_edges) = p_edge2;
! 641: pe_tail = pe_tail2;
! 642: }
! 643:
! 644: p_edge1 = p_edge2;
! 645: p_vrtx1 = p_vrtx2;
! 646: }
! 647:
! 648: /* Update the boundary flag, saved in each edge, and update indexes: */
! 649:
! 650: pe_temp = (*p_edges);
! 651:
! 652: while (pe_temp) {
! 653: if ((!(pe_temp->poly[0])) || (!(pe_temp->poly[1])))
! 654: (pe_temp->position) = BOUNDARY;
! 655: pe_temp = pe_temp->next;
! 656: }
! 657: }
! 658:
! 659: /*
! 660: * Calculate minimum and maximum values
! 661: */
! 662: static void calc_min_max(num_isolines, iso_lines, xx_min, yy_min, zz_min, xx_max, yy_max, zz_max)
! 663: int num_isolines; /* number of iso-lines input */
! 664: struct iso_curve *iso_lines; /* iso-lines input */
! 665: double *xx_min, *yy_min, *zz_min, *xx_max, *yy_max, *zz_max; /* min/max values in/out */
! 666: {
! 667: int i, j, grid_x_max;
! 668: struct coordinate GPHUGE *vertex;
! 669:
! 670: grid_x_max = iso_lines->p_count; /* number of vertices per iso_line */
! 671:
! 672: (*xx_min) = (*yy_min) = (*zz_min) = VERYLARGE; /* clear min/max values */
! 673: (*xx_max) = (*yy_max) = (*zz_max) = -VERYLARGE;
! 674:
! 675: for (j = 0; j < num_isolines; j++) {
! 676:
! 677: vertex = iso_lines->points;
! 678:
! 679: for (i = 0; i < grid_x_max; i++) {
! 680: if (vertex[i].type != UNDEFINED) {
! 681: if (vertex[i].x > (*xx_max))
! 682: (*xx_max) = vertex[i].x;
! 683: if (vertex[i].y > (*yy_max))
! 684: (*yy_max) = vertex[i].y;
! 685: if (vertex[i].z > (*zz_max))
! 686: (*zz_max) = vertex[i].z;
! 687: if (vertex[i].x < (*xx_min))
! 688: (*xx_min) = vertex[i].x;
! 689: if (vertex[i].y < (*yy_min))
! 690: (*yy_min) = vertex[i].y;
! 691: if (vertex[i].z < (*zz_min))
! 692: (*zz_min) = vertex[i].z;
! 693: }
! 694: }
! 695: iso_lines = iso_lines->next;
! 696: }
! 697: /*
! 698: * fprintf(stderr," x: %g, %g\n", (*xx_min), (*xx_max));
! 699: * fprintf(stderr," y: %g, %g\n", (*yy_min), (*yy_max));
! 700: * fprintf(stderr," z: %g, %g\n", (*zz_min), (*zz_max));
! 701: */
! 702: }
! 703:
! 704: /*
! 705: * Generate new edge and append it to list, but only if both vertices are
! 706: * defined. The list is referenced by p_edge and pe_tail (p_edge points on
! 707: * first edge and pe_tail on last one).
! 708: * Note, the list may be empty (pe_edge==pe_tail==NULL) on entry and exit.
! 709: */
! 710: static struct edge_struct *add_edge(point0, point1, p_edge, pe_tail)
! 711: struct coordinate GPHUGE * point0; /* 2 vertices input */
! 712: struct coordinate GPHUGE * point1;
! 713: struct edge_struct **p_edge, **pe_tail; /* pointers to edge list in/out */
! 714: {
! 715: struct edge_struct *pe_temp = NULL;
! 716:
! 717: if (point0->type != UNDEFINED && point1->type != UNDEFINED) {
! 718:
! 719: pe_temp = (struct edge_struct *)
! 720: gp_alloc((unsigned long) sizeof(struct edge_struct), "contour edge");
! 721:
! 722: pe_temp->poly[0] = NULL; /* clear links */
! 723: pe_temp->poly[1] = NULL;
! 724: pe_temp->vertex[0] = point0; /* First vertex of edge. */
! 725: pe_temp->vertex[1] = point1; /* Second vertex of edge. */
! 726: pe_temp->next = NULL;
! 727: pe_temp->position = INNER_MESH; /* default position in mesh */
! 728:
! 729: if ((*pe_tail)) {
! 730: (*pe_tail)->next = pe_temp; /* Stick new record as last one. */
! 731: } else {
! 732: (*p_edge) = pe_temp; /* start new list if empty */
! 733: }
! 734: (*pe_tail) = pe_temp; /* continue to last record. */
! 735:
! 736: }
! 737: return pe_temp; /* returns NULL, if no edge allocated */
! 738: }
! 739:
! 740: /*
! 741: * Generate new triangle and append it to list, but only if all edges are defined.
! 742: * The list is referenced by p_poly and pp_tail (p_poly points on first ploygon
! 743: * and pp_tail on last one).
! 744: * Note, the list may be empty (pe_ploy==pp_tail==NULL) on entry and exit.
! 745: */
! 746: static struct poly_struct *add_poly(edge0, edge1, edge2, p_poly, pp_tail)
! 747: struct edge_struct *edge0, *edge1, *edge2; /* 3 edges input */
! 748: struct poly_struct **p_poly, **pp_tail; /* pointers to polygon list in/out */
! 749: {
! 750: struct poly_struct *pp_temp = NULL;
! 751:
! 752: if (edge0 && edge1 && edge2) {
! 753:
! 754: pp_temp = (struct poly_struct *)
! 755: gp_alloc((unsigned long) sizeof(struct poly_struct), "contour polygon");
! 756:
! 757: pp_temp->edge[0] = edge0; /* First edge of triangle */
! 758: pp_temp->edge[1] = edge1; /* Second one */
! 759: pp_temp->edge[2] = edge2; /* Third one */
! 760: pp_temp->next = NULL;
! 761:
! 762: if (edge0->poly[0]) /* update edge0 */
! 763: edge0->poly[1] = pp_temp;
! 764: else
! 765: edge0->poly[0] = pp_temp;
! 766:
! 767: if (edge1->poly[0]) /* update edge1 */
! 768: edge1->poly[1] = pp_temp;
! 769: else
! 770: edge1->poly[0] = pp_temp;
! 771:
! 772: if (edge2->poly[0]) /* update edge2 */
! 773: edge2->poly[1] = pp_temp;
! 774: else
! 775: edge2->poly[0] = pp_temp;
! 776:
! 777: if ((*pp_tail)) /* Stick new record as last one. */
! 778: (*pp_tail)->next = pp_temp;
! 779: else
! 780: (*p_poly) = pp_temp; /* start new list if empty */
! 781:
! 782: (*pp_tail) = pp_temp; /* continue to last record. */
! 783:
! 784: }
! 785: return pp_temp; /* returns NULL, if no edge allocated */
! 786: }
! 787:
! 788:
! 789:
! 790: /*
! 791: * Calls the (hopefully) desired interpolation/approximation routine.
! 792: */
! 793: static void put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind)
! 794: struct cntr_struct *p_cntr; /* contour structure input */
! 795: double z_level, /* Z level of contour input */
! 796: xx_min, xx_max, yy_min, yy_max; /* minimum/maximum values input */
! 797: int contr_kind; /* OPEN_ or CLOESED_CONTOUR input */
! 798: {
! 799:
! 800: if (!p_cntr)
! 801: return; /* Nothing to do if it is empty contour. */
! 802:
! 803: switch (interp_kind) {
! 804: case INTERP_NOTHING: /* No interpolation/approximation. */
! 805: put_contour_nothing(p_cntr);
! 806: break;
! 807: case INTERP_CUBIC: /* Cubic spline interpolation. */
! 808: put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max,
! 809: chk_contour_kind(p_cntr, contr_kind));
! 810:
! 811: break;
! 812: case APPROX_BSPLINE: /* Bspline approximation. */
! 813: put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max,
! 814: chk_contour_kind(p_cntr, contr_kind));
! 815: break;
! 816: }
! 817: free_contour(p_cntr);
! 818: }
! 819:
! 820: /*
! 821: * Simply puts contour coordinates in order with no interpolation or
! 822: * approximation.
! 823: */
! 824: static void put_contour_nothing(p_cntr)
! 825: struct cntr_struct *p_cntr;
! 826: {
! 827: while (p_cntr) {
! 828: add_cntr_point(p_cntr->X, p_cntr->Y);
! 829: p_cntr = p_cntr->next;
! 830: }
! 831: end_crnt_cntr();
! 832: }
! 833:
! 834: /*
! 835: * for some reason contours are never flagged as CLOSED_CONTOUR
! 836: * if first point == last point, set flag accordingly
! 837: *
! 838: */
! 839:
! 840: static int chk_contour_kind(p_cntr, contr_kind)
! 841: struct cntr_struct *p_cntr;
! 842: int contr_kind;
! 843: {
! 844: struct cntr_struct *pc_tail = NULL;
! 845: int current_contr_kind;
! 846:
! 847: FPRINTF((stderr, "check_contour_kind: current contr_kind value is %d\n", contr_kind));
! 848:
! 849: current_contr_kind = contr_kind;
! 850:
! 851: if (contr_kind != CLOSED_CONTOUR) {
! 852: pc_tail = p_cntr;
! 853: while (pc_tail->next)
! 854: pc_tail = pc_tail->next; /* Find last point. */
! 855:
! 856: /* test if first and last point are equal */
! 857: if (fuzzy_equal(pc_tail, p_cntr)) {
! 858: current_contr_kind = CLOSED_CONTOUR;
! 859: FPRINTF((stderr, "check_contour_kind: contr_kind changed to %d\n", current_contr_kind));
! 860: }
! 861: }
! 862: return (current_contr_kind);
! 863: }
! 864:
! 865: /*
! 866: * Generate a cubic spline curve through the points (x_i,y_i) which are
! 867: * stored in the linked list p_cntr.
! 868: * The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the
! 869: * parameter t is the length of the linear stroke.
! 870: */
! 871: static void put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind)
! 872: struct cntr_struct *p_cntr;
! 873: double z_level, xx_min, xx_max, yy_min, yy_max;
! 874: int contr_kind;
! 875: {
! 876: int num_pts, num_intpol;
! 877: double unit_x, unit_y; /* To define norm (x,y)-plane */
! 878: double *delta_t; /* Interval length t_{i+1}-t_i */
! 879: double *d2x, *d2y; /* Second derivatives x''(t_i), y''(t_i) */
! 880: struct cntr_struct *pc_tail;
! 881:
! 882: num_pts = count_contour(p_cntr); /* Number of points in contour. */
! 883:
! 884: pc_tail = p_cntr; /* Find last point. */
! 885: while (pc_tail->next)
! 886: pc_tail = pc_tail->next;
! 887:
! 888: if (contr_kind == CLOSED_CONTOUR) {
! 889: /* Test if first and last point are equal (should be) */
! 890: if (!fuzzy_equal(pc_tail, p_cntr)) {
! 891: pc_tail->next = p_cntr; /* Close contour list - make it circular. */
! 892: num_pts++;
! 893: }
! 894: }
! 895: delta_t = (double *)
! 896: gp_alloc((unsigned long) (sizeof(double) * num_pts), "contour delta_t");
! 897: d2x = (double *)
! 898: gp_alloc((unsigned long) (sizeof(double) * num_pts), "contour d2x");
! 899: d2y = (double *)
! 900: gp_alloc((unsigned long) (sizeof(double) * num_pts), "contour d2y");
! 901:
! 902: /* Width and hight of the grid is used at unit length (2d-norm) */
! 903: unit_x = xx_max - x_min;
! 904: unit_y = yy_max - y_min;
! 905: unit_x = (unit_x > zero ? unit_x : zero); /* should not be zero */
! 906: unit_y = (unit_y > zero ? unit_y : zero);
! 907:
! 908: if (num_pts > 2) {
! 909: /*
! 910: * Calculate second derivatives d2x[], d2y[] and interval lengths delta_t[]:
! 911: */
! 912: if (!gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t,
! 913: contr_kind, unit_x, unit_y)) {
! 914: free((char *) delta_t);
! 915: free((char *) d2x);
! 916: free((char *) d2y);
! 917: if (contr_kind == CLOSED_CONTOUR)
! 918: pc_tail->next = NULL; /* Un-circular list */
! 919: return;
! 920: }
! 921: }
! 922: /* If following (num_pts > 1) is TRUE then exactly 2 points in contour. */
! 923: else if (num_pts > 1) {
! 924: /* set all second derivatives to zero, interval length to 1 */
! 925: d2x[0] = 0.;
! 926: d2y[0] = 0.;
! 927: d2x[1] = 0.;
! 928: d2y[1] = 0.;
! 929: delta_t[0] = 1.;
! 930: } else { /* Only one point ( ?? ) - ignore it. */
! 931: free((char *) delta_t);
! 932: free((char *) d2x);
! 933: free((char *) d2y);
! 934: if (contr_kind == CLOSED_CONTOUR)
! 935: pc_tail->next = NULL; /* Un-circular list */
! 936: return;
! 937: }
! 938:
! 939: /* Calculate "num_intpol" interpolated values */
! 940: num_intpol = 1 + (num_pts - 1) * num_approx_pts; /* global: num_approx_pts */
! 941: intp_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, num_intpol);
! 942:
! 943: free((char *) delta_t);
! 944: free((char *) d2x);
! 945: free((char *) d2y);
! 946:
! 947: if (contr_kind == CLOSED_CONTOUR)
! 948: pc_tail->next = NULL; /* Un-circular list */
! 949:
! 950: end_crnt_cntr();
! 951: }
! 952:
! 953:
! 954: /*
! 955: * Find Bspline approximation for this data set.
! 956: * Uses global variable num_approx_pts to determine number of samples per
! 957: * interval, where the knot vector intervals are assumed to be uniform, and
! 958: * Global variable bspline_order for the order of Bspline to use.
! 959: */
! 960: static void put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind)
! 961: struct cntr_struct *p_cntr;
! 962: double z_level, xx_min, xx_max, yy_min, yy_max;
! 963: int contr_kind;
! 964: {
! 965: int num_pts, order = bspline_order;
! 966:
! 967: num_pts = count_contour(p_cntr); /* Number of points in contour. */
! 968: if (num_pts < 2)
! 969: return; /* Can't do nothing if empty or one points! */
! 970: /* Order must be less than number of points in curve - fix it if needed. */
! 971: if (order > num_pts - 1)
! 972: order = num_pts - 1;
! 973:
! 974: gen_bspline_approx(p_cntr, num_pts, order, contr_kind);
! 975: end_crnt_cntr();
! 976: }
! 977:
! 978: /*
! 979: * Free all elements in the contour list.
! 980: */
! 981: static void free_contour(p_cntr)
! 982: struct cntr_struct *p_cntr;
! 983: {
! 984: struct cntr_struct *pc_temp;
! 985:
! 986: while (p_cntr) {
! 987: pc_temp = p_cntr;
! 988: p_cntr = p_cntr->next;
! 989: free((char *) pc_temp);
! 990: }
! 991: }
! 992:
! 993: /*
! 994: * Counts number of points in contour.
! 995: */
! 996: static int count_contour(p_cntr)
! 997: struct cntr_struct *p_cntr;
! 998: {
! 999: int count = 0;
! 1000:
! 1001: while (p_cntr) {
! 1002: count++;
! 1003: p_cntr = p_cntr->next;
! 1004: }
! 1005: return count;
! 1006: }
! 1007:
! 1008: /*
! 1009: * Find second derivatives (x''(t_i),y''(t_i)) of cubic spline interpolation
! 1010: * through list of points (x_i,y_i). The parameter t is calculated as the
! 1011: * length of the linear stroke. The number of points must be at least 3.
! 1012: * Note: For CLOSED_CONTOURs the first and last point must be equal.
! 1013: */
! 1014: static int gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, contr_kind, unit_x, unit_y)
! 1015: int num_pts; /* Number of points (num_pts>=3), input */
! 1016: struct cntr_struct *p_cntr; /* List of points (x(t_i),y(t_i)), input */
! 1017: double d2x[], d2y[], /* Second derivatives (x''(t_i),y''(t_i)), output */
! 1018: delta_t[]; /* List of interval lengths t_{i+1}-t_{i}, output */
! 1019: int contr_kind; /* CLOSED_CONTOUR or OPEN_CONTOUR, input */
! 1020: double unit_x, unit_y; /* Unit length in x and y (norm=1), input */
! 1021: {
! 1022: int n, i;
! 1023: double norm;
! 1024: tri_diag *m; /* The tri-diagonal matrix is saved here. */
! 1025: struct cntr_struct *pc_temp;
! 1026:
! 1027: m = (tri_diag *)
! 1028: gp_alloc((unsigned long) (sizeof(tri_diag) * num_pts), "contour tridiag m");
! 1029:
! 1030: /*
! 1031: * Calculate first differences in (d2x[i], d2y[i]) and interval lengths
! 1032: * in delta_t[i]:
! 1033: */
! 1034: pc_temp = p_cntr;
! 1035: for (i = 0; i < num_pts - 1; i++) {
! 1036: d2x[i] = pc_temp->next->X - pc_temp->X;
! 1037: d2y[i] = pc_temp->next->Y - pc_temp->Y;
! 1038: /*
! 1039: * The Norm of a linear stroke is calculated in "normal coordinates"
! 1040: * and used as interval length:
! 1041: */
! 1042: delta_t[i] = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y));
! 1043:
! 1044: d2x[i] /= delta_t[i]; /* first difference, with unit norm: */
! 1045: d2y[i] /= delta_t[i]; /* || (d2x[i], d2y[i]) || = 1 */
! 1046:
! 1047: pc_temp = pc_temp->next;
! 1048: }
! 1049:
! 1050: /*
! 1051: * Setup linear System: M * x = b
! 1052: */
! 1053: n = num_pts - 2; /* Without first and last point */
! 1054: if (contr_kind == CLOSED_CONTOUR) {
! 1055: /* First and last points must be equal for CLOSED_CONTOURs */
! 1056: delta_t[num_pts - 1] = delta_t[0];
! 1057: d2x[num_pts - 1] = d2x[0];
! 1058: d2y[num_pts - 1] = d2y[0];
! 1059: n++; /* Add last point (= first point) */
! 1060: }
! 1061: for (i = 0; i < n; i++) {
! 1062: /* Matrix M, mainly tridiagonal with cyclic second index ("j = j+n mod n") */
! 1063: m[i][0] = delta_t[i]; /* Off-diagonal element M_{i,i-1} */
! 1064: m[i][1] = 2. * (delta_t[i] + delta_t[i + 1]); /* M_{i,i} */
! 1065: m[i][2] = delta_t[i + 1]; /* Off-diagonal element M_{i,i+1} */
! 1066:
! 1067: /* Right side b_x and b_y */
! 1068: d2x[i] = (d2x[i + 1] - d2x[i]) * 6.;
! 1069: d2y[i] = (d2y[i + 1] - d2y[i]) * 6.;
! 1070:
! 1071: /*
! 1072: * If the linear stroke shows a cusps of more than 90 degree, the right
! 1073: * side is reduced to avoid oscillations in the spline:
! 1074: */
! 1075: norm = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y)) / 8.5;
! 1076:
! 1077: if (norm > 1.) {
! 1078: d2x[i] /= norm;
! 1079: d2y[i] /= norm;
! 1080: /* The first derivative will not be continuous */
! 1081: }
! 1082: }
! 1083:
! 1084: if (contr_kind != CLOSED_CONTOUR) {
! 1085: /* Third derivative is set to zero at both ends */
! 1086: m[0][1] += m[0][0]; /* M_{0,0} */
! 1087: m[0][0] = 0.; /* M_{0,n-1} */
! 1088: m[n - 1][1] += m[n - 1][2]; /* M_{n-1,n-1} */
! 1089: m[n - 1][2] = 0.; /* M_{n-1,0} */
! 1090: }
! 1091: /* Solve linear systems for d2x[] and d2y[] */
! 1092:
! 1093:
! 1094: if (solve_cubic_1(m, n)) { /* Calculate Cholesky decomposition */
! 1095: solve_cubic_2(m, d2x, n); /* solve M * d2x = b_x */
! 1096: solve_cubic_2(m, d2y, n); /* solve M * d2y = b_y */
! 1097:
! 1098: } else { /* Should not happen, but who knows ... */
! 1099: free((char *) m);
! 1100: return FALSE;
! 1101: }
! 1102:
! 1103: /* Shift all second derivatives one place right and abdate end points */
! 1104: for (i = n; i > 0; i--) {
! 1105: d2x[i] = d2x[i - 1];
! 1106: d2y[i] = d2y[i - 1];
! 1107: }
! 1108: if (contr_kind == CLOSED_CONTOUR) {
! 1109: d2x[0] = d2x[n];
! 1110: d2y[0] = d2y[n];
! 1111: } else {
! 1112: d2x[0] = d2x[1]; /* Third derivative is zero in */
! 1113: d2y[0] = d2y[1]; /* first and last interval */
! 1114: d2x[n + 1] = d2x[n];
! 1115: d2y[n + 1] = d2y[n];
! 1116: }
! 1117:
! 1118: free((char *) m);
! 1119: return TRUE;
! 1120: }
! 1121:
! 1122: /*
! 1123: * Calculate interpolated values of the spline function (defined via p_cntr
! 1124: * and the second derivatives d2x[] and d2y[]). The number of tabulated
! 1125: * values is n. On an equidistant grid n_intpol values are calculated.
! 1126: */
! 1127: static void intp_cubic_spline(n, p_cntr, d2x, d2y, delta_t, n_intpol)
! 1128: int n;
! 1129: struct cntr_struct *p_cntr;
! 1130: double d2x[], d2y[], delta_t[];
! 1131: int n_intpol;
! 1132: {
! 1133: double t, t_skip, t_max;
! 1134: double x0, x1, x, y0, y1, y;
! 1135: double d, hx, dx0, dx01, hy, dy0, dy01;
! 1136: int i;
! 1137:
! 1138: /* The length of the total interval */
! 1139: t_max = 0.;
! 1140: for (i = 0; i < n - 1; i++)
! 1141: t_max += delta_t[i];
! 1142:
! 1143: /* The distance between interpolated points */
! 1144: t_skip = (1. - 1e-7) * t_max / (n_intpol - 1);
! 1145:
! 1146: t = 0.; /* Parameter value */
! 1147: x1 = p_cntr->X;
! 1148: y1 = p_cntr->Y;
! 1149: add_cntr_point(x1, y1); /* First point. */
! 1150: t += t_skip;
! 1151:
! 1152: for (i = 0; i < n - 1; i++) {
! 1153: p_cntr = p_cntr->next;
! 1154:
! 1155: d = delta_t[i]; /* Interval length */
! 1156: x0 = x1;
! 1157: y0 = y1;
! 1158: x1 = p_cntr->X;
! 1159: y1 = p_cntr->Y;
! 1160: hx = (x1 - x0) / d;
! 1161: hy = (y1 - y0) / d;
! 1162: dx0 = (d2x[i + 1] + 2 * d2x[i]) / 6.;
! 1163: dy0 = (d2y[i + 1] + 2 * d2y[i]) / 6.;
! 1164: dx01 = (d2x[i + 1] - d2x[i]) / (6. * d);
! 1165: dy01 = (d2y[i + 1] - d2y[i]) / (6. * d);
! 1166: while (t <= delta_t[i]) { /* t in current interval ? */
! 1167: x = x0 + t * (hx + (t - d) * (dx0 + t * dx01));
! 1168: y = y0 + t * (hy + (t - d) * (dy0 + t * dy01));
! 1169: add_cntr_point(x, y); /* next point. */
! 1170: t += t_skip;
! 1171: }
! 1172: t -= delta_t[i]; /* Parameter t relative to start of next interval */
! 1173: }
! 1174: }
! 1175:
! 1176: /*
! 1177: * The following two procedures solve the special linear system which arise
! 1178: * in cubic spline interpolation. If x is assumed cyclic ( x[i]=x[n+i] ) the
! 1179: * equations can be written as (i=0,1,...,n-1):
! 1180: * m[i][0] * x[i-1] + m[i][1] * x[i] + m[i][2] * x[i+1] = b[i] .
! 1181: * In matrix notation one gets M * x = b, where the matrix M is tridiagonal
! 1182: * with additional elements in the upper right and lower left position:
! 1183: * m[i][0] = M_{i,i-1} for i=1,2,...,n-1 and m[0][0] = M_{0,n-1} ,
! 1184: * m[i][1] = M_{i, i } for i=0,1,...,n-1
! 1185: * m[i][2] = M_{i,i+1} for i=0,1,...,n-2 and m[n-1][2] = M_{n-1,0}.
! 1186: * M should be symmetric (m[i+1][0]=m[i][2]) and positiv definite.
! 1187: * The size of the system is given in n (n>=1).
! 1188: *
! 1189: * In the first procedure the Cholesky decomposition M = C^T * D * C
! 1190: * (C is upper triangle with unit diagonal, D is diagonal) is calculated.
! 1191: * Return TRUE if decomposition exist.
! 1192: */
! 1193: static int solve_cubic_1(m, n)
! 1194: tri_diag m[];
! 1195: int n;
! 1196: {
! 1197: int i;
! 1198: double m_ij, m_n, m_nn, d;
! 1199:
! 1200: if (n < 1)
! 1201: return FALSE; /* Dimension should be at least 1 */
! 1202:
! 1203: d = m[0][1]; /* D_{0,0} = M_{0,0} */
! 1204: if (d <= 0.)
! 1205: return FALSE; /* M (or D) should be positiv definite */
! 1206: m_n = m[0][0]; /* M_{0,n-1} */
! 1207: m_nn = m[n - 1][1]; /* M_{n-1,n-1} */
! 1208: for (i = 0; i < n - 2; i++) {
! 1209: m_ij = m[i][2]; /* M_{i,1} */
! 1210: m[i][2] = m_ij / d; /* C_{i,i+1} */
! 1211: m[i][0] = m_n / d; /* C_{i,n-1} */
! 1212: m_nn -= m[i][0] * m_n; /* to get C_{n-1,n-1} */
! 1213: m_n = -m[i][2] * m_n; /* to get C_{i+1,n-1} */
! 1214: d = m[i + 1][1] - m[i][2] * m_ij; /* D_{i+1,i+1} */
! 1215: if (d <= 0.)
! 1216: return FALSE; /* Elements of D should be positiv */
! 1217: m[i + 1][1] = d;
! 1218: }
! 1219: if (n >= 2) { /* Complete last column */
! 1220: m_n += m[n - 2][2]; /* add M_{n-2,n-1} */
! 1221: m[n - 2][0] = m_n / d; /* C_{n-2,n-1} */
! 1222: m[n - 1][1] = d = m_nn - m[n - 2][0] * m_n; /* D_{n-1,n-1} */
! 1223: if (d <= 0.)
! 1224: return FALSE;
! 1225: }
! 1226: return TRUE;
! 1227: }
! 1228:
! 1229: /*
! 1230: * The second procedure solves the linear system, with the Choleky
! 1231: * decomposition calculated above (in m[][]) and the right side b given
! 1232: * in x[]. The solution x overwrites the right side in x[].
! 1233: */
! 1234: static void solve_cubic_2(m, x, n)
! 1235: tri_diag m[];
! 1236: double x[];
! 1237: int n;
! 1238: {
! 1239: int i;
! 1240: double x_n;
! 1241:
! 1242: /* Division by transpose of C : b = C^{-T} * b */
! 1243: x_n = x[n - 1];
! 1244: for (i = 0; i < n - 2; i++) {
! 1245: x[i + 1] -= m[i][2] * x[i]; /* C_{i,i+1} * x_{i} */
! 1246: x_n -= m[i][0] * x[i]; /* C_{i,n-1} * x_{i} */
! 1247: }
! 1248: if (n >= 2)
! 1249: x[n - 1] = x_n - m[n - 2][0] * x[n - 2]; /* C_{n-2,n-1} * x_{n-1} */
! 1250:
! 1251: /* Division by D: b = D^{-1} * b */
! 1252: for (i = 0; i < n; i++)
! 1253: x[i] /= m[i][1];
! 1254:
! 1255: /* Division by C: b = C^{-1} * b */
! 1256: x_n = x[n - 1];
! 1257: if (n >= 2)
! 1258: x[n - 2] -= m[n - 2][0] * x_n; /* C_{n-2,n-1} * x_{n-1} */
! 1259: for (i = n - 3; i >= 0; i--) {
! 1260: /* C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1} */
! 1261: x[i] -= m[i][2] * x[i + 1] + m[i][0] * x_n;
! 1262: }
! 1263: return;
! 1264: }
! 1265:
! 1266: /*
! 1267: * Solve tri diagonal linear system equation. The tri diagonal matrix is
! 1268: * defined via matrix M, right side is r, and solution X i.e. M * X = R.
! 1269: * Size of system given in n. Return TRUE if solution exist.
! 1270: */
! 1271: /* not used any more in "contour.c", but in "spline.c" (21. Dec. 1995) ! */
! 1272:
! 1273: int solve_tri_diag(m, r, x, n)
! 1274: tri_diag m[];
! 1275: double r[], x[];
! 1276: int n;
! 1277: {
! 1278: int i;
! 1279: double t;
! 1280:
! 1281: for (i = 1; i < n; i++) { /* Eliminate element m[i][i-1] (lower diagonal). */
! 1282: if (m[i - 1][1] == 0)
! 1283: return FALSE;
! 1284: t = m[i][0] / m[i - 1][1]; /* Find ratio between the two lines. */
! 1285: /* m[i][0] = m[i][0] - m[i-1][1] * t; */
! 1286: /* m[i][0] is not used any more (and set to 0 in the above line) */
! 1287: m[i][1] = m[i][1] - m[i - 1][2] * t;
! 1288: r[i] = r[i] - r[i - 1] * t;
! 1289: }
! 1290: /* Now do back subtitution - update the solution vector X: */
! 1291: if (m[n - 1][1] == 0)
! 1292: return FALSE;
! 1293: x[n - 1] = r[n - 1] / m[n - 1][1]; /* Find last element. */
! 1294: for (i = n - 2; i >= 0; i--) {
! 1295: if (m[i][1] == 0)
! 1296: return FALSE;
! 1297: x[i] = (r[i] - x[i + 1] * m[i][2]) / m[i][1];
! 1298: }
! 1299: return TRUE;
! 1300: }
! 1301:
! 1302: /*
! 1303: * Generate a Bspline curve defined by all the points given in linked list p:
! 1304: * Algorithm: using deBoor algorithm
! 1305: * Note: if Curvekind is OPEN_CONTOUR than Open end knot vector is assumed,
! 1306: * else (CLOSED_CONTOUR) Float end knot vector is assumed.
! 1307: * It is assumed that num_of_points is at least 2, and order of Bspline is less
! 1308: * than num_of_points!
! 1309: */
! 1310: static void gen_bspline_approx(p_cntr, num_of_points, order, contr_kind)
! 1311: struct cntr_struct *p_cntr;
! 1312: int num_of_points, order, contr_kind;
! 1313: {
! 1314: int knot_index = 0, pts_count = 1;
! 1315: double dt, t, next_t, t_min, t_max, x, y;
! 1316: struct cntr_struct *pc_temp = p_cntr, *pc_tail = NULL;
! 1317:
! 1318: /* If the contour is Closed one we must update few things:
! 1319: * 1. Make the list temporary circular, so we can close the contour.
! 1320: * 2. Update num_of_points - increase it by "order-1" so contour will be
! 1321: * closed. This will evaluate order more sections to close it!
! 1322: */
! 1323: if (contr_kind == CLOSED_CONTOUR) {
! 1324: pc_tail = p_cntr;
! 1325: while (pc_tail->next)
! 1326: pc_tail = pc_tail->next; /* Find last point. */
! 1327:
! 1328: /* test if first and last point are equal */
! 1329: if (fuzzy_equal(pc_tail, p_cntr)) {
! 1330: /* Close contour list - make it circular. */
! 1331: pc_tail->next = p_cntr->next;
! 1332: num_of_points += order - 1;
! 1333: } else {
! 1334: pc_tail->next = p_cntr;
! 1335: num_of_points += order;
! 1336: }
! 1337: }
! 1338: /* Find first (t_min) and last (t_max) t value to eval: */
! 1339: t = t_min = fetch_knot(contr_kind, num_of_points, order, order);
! 1340: t_max = fetch_knot(contr_kind, num_of_points, order, num_of_points);
! 1341: next_t = t_min + 1.0;
! 1342: knot_index = order;
! 1343: dt = 1.0 / num_approx_pts; /* Number of points per one section. */
! 1344:
! 1345:
! 1346: while (t < t_max) {
! 1347: if (t > next_t) {
! 1348: pc_temp = pc_temp->next; /* Next order ctrl. pt. to blend. */
! 1349: knot_index++;
! 1350: next_t += 1.0;
! 1351: }
! 1352: eval_bspline(t, pc_temp, num_of_points, order, knot_index,
! 1353: contr_kind, &x, &y); /* Next pt. */
! 1354: add_cntr_point(x, y);
! 1355: pts_count++;
! 1356: /* As we might have some real number round off problems we do */
! 1357: /* the last point outside the loop */
! 1358: if (pts_count == num_approx_pts * (num_of_points - order) + 1)
! 1359: break;
! 1360: t += dt;
! 1361: }
! 1362:
! 1363: /* Now do the last point */
! 1364: eval_bspline(t_max - EPSILON, pc_temp, num_of_points, order, knot_index,
! 1365: contr_kind, &x, &y);
! 1366: add_cntr_point(x, y); /* Complete the contour. */
! 1367:
! 1368: if (contr_kind == CLOSED_CONTOUR) /* Update list - un-circular it. */
! 1369: pc_tail->next = NULL;
! 1370: }
! 1371:
! 1372: /*
! 1373: * The routine to evaluate the B-spline value at point t using knot vector
! 1374: * from function fetch_knot(), and the control points p_cntr.
! 1375: * Returns (x, y) of approximated B-spline. Note that p_cntr points on the
! 1376: * first control point to blend with. The B-spline is of order order.
! 1377: */
! 1378: static void eval_bspline(t, p_cntr, num_of_points, order, j, contr_kind, x, y)
! 1379: double t;
! 1380: struct cntr_struct *p_cntr;
! 1381: int num_of_points, order, j, contr_kind;
! 1382: double *x, *y;
! 1383: {
! 1384: int i, p;
! 1385: double ti, tikp, *dx, *dy; /* Copy p_cntr into it to make it faster. */
! 1386:
! 1387: dx = (double *)
! 1388: gp_alloc((unsigned long) (sizeof(double) * (order + j)), "contour b_spline");
! 1389: dy = (double *)
! 1390: gp_alloc((unsigned long) (sizeof(double) * (order + j)), "contour b_spline");
! 1391:
! 1392: /* Set the dx/dy - [0] iteration step, control points (p==0 iterat.): */
! 1393: for (i = j - order; i <= j; i++) {
! 1394: dx[i] = p_cntr->X;
! 1395: dy[i] = p_cntr->Y;
! 1396: p_cntr = p_cntr->next;
! 1397: }
! 1398:
! 1399: for (p = 1; p <= order; p++) { /* Iteration (b-spline level) counter. */
! 1400: for (i = j; i >= j - order + p; i--) { /* Control points indexing. */
! 1401: ti = fetch_knot(contr_kind, num_of_points, order, i);
! 1402: tikp = fetch_knot(contr_kind, num_of_points, order, i + order + 1 - p);
! 1403: if (ti == tikp) { /* Should not be a problems but how knows... */
! 1404: } else {
! 1405: dx[i] = dx[i] * (t - ti) / (tikp - ti) + /* Calculate x. */
! 1406: dx[i - 1] * (tikp - t) / (tikp - ti);
! 1407: dy[i] = dy[i] * (t - ti) / (tikp - ti) + /* Calculate y. */
! 1408: dy[i - 1] * (tikp - t) / (tikp - ti);
! 1409: }
! 1410: }
! 1411: }
! 1412: *x = dx[j];
! 1413: *y = dy[j];
! 1414: free((char *) dx);
! 1415: free((char *) dy);
! 1416: }
! 1417:
! 1418: /*
! 1419: * Routine to get the i knot from uniform knot vector. The knot vector
! 1420: * might be float (Knot(i) = i) or open (where the first and last "order"
! 1421: * knots are equal). contr_kind determines knot kind - OPEN_CONTOUR means
! 1422: * open knot vector, and CLOSED_CONTOUR selects float knot vector.
! 1423: * Note the knot vector is not exist and this routine simulates it existance
! 1424: * Also note the indexes for the knot vector starts from 0.
! 1425: */
! 1426: static double fetch_knot(contr_kind, num_of_points, order, i)
! 1427: int contr_kind, num_of_points, order, i;
! 1428: {
! 1429: switch (contr_kind) {
! 1430: case OPEN_CONTOUR:
! 1431: if (i <= order)
! 1432: return 0.0;
! 1433: else if (i <= num_of_points)
! 1434: return (double) (i - order);
! 1435: else
! 1436: return (double) (num_of_points - order);
! 1437: case CLOSED_CONTOUR:
! 1438: return (double) i;
! 1439: default: /* Should never happen */
! 1440: return 1.0;
! 1441: }
! 1442: #ifdef sequent
! 1443: return 1.0; /* ???? */
! 1444: #endif
! 1445: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to contour.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gnuplot