Annotation of OpenXM_contrib/gnuplot/contour.c, Revision 1.1.1.3
1.1 maekawa 1: #ifndef lint
1.1.1.3 ! ohara 2: static char *RCSid = "$Id: contour.c,v 1.9.2.2 2002/01/31 21:18:22 lhecking Exp $";
1.1 maekawa 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 *)
1.1.1.3 ! ohara 292: gp_alloc(sizeof(struct gnuplot_contours), "gnuplot_contour");
1.1 maekawa 293: cntr->coords = (struct coordinate GPHUGE *)
1.1.1.3 ! ohara 294: gp_alloc(sizeof(struct coordinate) * crnt_cntr_pt_index,
! 295: "contour coords");
1.1 maekawa 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 *)
1.1.1.3 ! ohara 504: gp_alloc(sizeof(struct cntr_struct), "contour cntr_struct");
1.1 maekawa 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;
1.1.1.3 ! ohara 539: struct edge_struct *p_edge1, *p_edge2, *edge0, *edge1, *edge2,
! 540: *pe_tail, *pe_tail2, *pe_temp;
1.1 maekawa 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 */
1.1.1.3 ! ohara 548: p_edge1 = pe_tail = NULL; /* clear list of edges */
1.1 maekawa 549:
550: /* Generate edges of first row */
1.1.1.3 ! ohara 551: /* HBB 19991130: removed effectively unused variable 'pe_tail1' */
1.1 maekawa 552: for (j = 0; j < grid_x_max - 1; j++)
1.1.1.3 ! ohara 553: add_edge(p_vrtx1 + j, p_vrtx1 + j + 1, &p_edge1, &pe_tail);
1.1 maekawa 554:
555: (*p_edges) = p_edge1; /* update main list */
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:
1.1.1.3 ! ohara 565: * 1. p_edge2: Top horizontal edge list: +-----------------------+ 2
! 566: * 2. p_tail : end of middle edge list: |\ |\ |\ |\ |\ |\ |
1.1 maekawa 567: * | \| \| \| \| \| \|
1.1.1.3 ! ohara 568: * 3. p_edge1: Bottom horizontal edge list: +-----------------------+ 1
! 569: *
! 570: * pe_tail2 : end of list beginning at p_edge2
! 571: * pe_temp : position inside list beginning at p_edge1
! 572: * p_edges : head of the master edge list (part of our output)
! 573: * p_vrtx1 : start of lower row of input vertices
! 574: * p_vrtx2 : start of higher row of input vertices
1.1 maekawa 575: *
576: * The routine generates two triangle Lower Upper 1
577: * upper one and lower one: | \ ----
578: * (Nums. are edges order in polys) 0| \1 0\ |2
579: * The polygons are stored in the polygon ---- \ |
580: * list (*p_polys) (pp_tail points on 2
581: * last polygon).
582: * 1
583: * -----------
584: * In addition, the edge lists are updated - | \ 0 |
585: * each edge has two pointers on the two | \ |
586: * (one active if boundary) polygons which 0|1 0\1 0|1
587: * uses it. These two pointer to polygons | \ |
588: * are named: poly[0], poly[1]. The diagram | 1 \ |
589: * on the right show how they are used for the -----------
590: * upper and lower polygons (INNER_MESH polygons only). 0
591: */
592:
593: for (i = 1; i < num_isolines; i++) {
594: /* Read next column and gen. polys. */
595: iso_lines = iso_lines->next;
596:
597: p_vrtx2 = iso_lines->points; /* next row of vertices */
598: p_edge2 = pe_tail2 = NULL; /* clear top horizontal list */
599: pe_temp = p_edge1; /* pointer in bottom list */
600:
601: /*
602: * Generate edges and triagles for next row:
603: */
604:
605: /* generate first vertical edge */
606: edge2 = add_edge(p_vrtx1, p_vrtx2, p_edges, &pe_tail);
607:
608: for (j = 0; j < grid_x_max - 1; j++) {
609:
610: /* copy vertical edge for lower triangle */
611: edge0 = edge2;
612:
613: if (pe_temp && pe_temp->vertex[0] == p_vrtx1 + j) {
614: /* test lower edge */
615: edge2 = pe_temp;
616: pe_temp = pe_temp->next;
617: } else {
618: edge2 = NULL; /* edge is undefined */
619: }
620:
621: /* generate diagonal edge */
622: edge1 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j, p_edges, &pe_tail);
623: if (edge1)
624: edge1->position = DIAGONAL;
625:
626: /* generate lower triangle */
627: lower_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
628:
629: /* copy diagonal edge for upper triangle */
630: edge0 = edge1;
631:
632: /* generate upper edge */
633: edge1 = add_edge(p_vrtx2 + j, p_vrtx2 + j + 1, &p_edge2, &pe_tail2);
634:
635: /* generate vertical edge */
636: edge2 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j + 1, p_edges, &pe_tail);
637:
638: /* generate upper triangle */
639: upper_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
640: }
641:
1.1.1.3 ! ohara 642: if (p_edge2) {
! 643: /* HBB 19991130 bugfix: if p_edge2 list is empty,
! 644: * don't change p_edges list! Crashes by access
! 645: * to NULL pointer pe_tail, the second time through,
! 646: * otherwise */
! 647: if ((*p_edges)) { /* Chain new edges to main list. */
! 648: pe_tail->next = p_edge2;
! 649: pe_tail = pe_tail2;
! 650: } else {
! 651: (*p_edges) = p_edge2;
! 652: pe_tail = pe_tail2;
! 653: }
1.1 maekawa 654: }
655:
1.1.1.3 ! ohara 656: /* this row finished, move list heads up one row: */
1.1 maekawa 657: p_edge1 = p_edge2;
658: p_vrtx1 = p_vrtx2;
659: }
660:
661: /* Update the boundary flag, saved in each edge, and update indexes: */
662:
663: pe_temp = (*p_edges);
664:
665: while (pe_temp) {
666: if ((!(pe_temp->poly[0])) || (!(pe_temp->poly[1])))
667: (pe_temp->position) = BOUNDARY;
668: pe_temp = pe_temp->next;
669: }
670: }
671:
672: /*
673: * Calculate minimum and maximum values
674: */
675: static void calc_min_max(num_isolines, iso_lines, xx_min, yy_min, zz_min, xx_max, yy_max, zz_max)
676: int num_isolines; /* number of iso-lines input */
677: struct iso_curve *iso_lines; /* iso-lines input */
678: double *xx_min, *yy_min, *zz_min, *xx_max, *yy_max, *zz_max; /* min/max values in/out */
679: {
680: int i, j, grid_x_max;
681: struct coordinate GPHUGE *vertex;
682:
683: grid_x_max = iso_lines->p_count; /* number of vertices per iso_line */
684:
685: (*xx_min) = (*yy_min) = (*zz_min) = VERYLARGE; /* clear min/max values */
686: (*xx_max) = (*yy_max) = (*zz_max) = -VERYLARGE;
687:
688: for (j = 0; j < num_isolines; j++) {
689:
690: vertex = iso_lines->points;
691:
692: for (i = 0; i < grid_x_max; i++) {
693: if (vertex[i].type != UNDEFINED) {
694: if (vertex[i].x > (*xx_max))
695: (*xx_max) = vertex[i].x;
696: if (vertex[i].y > (*yy_max))
697: (*yy_max) = vertex[i].y;
698: if (vertex[i].z > (*zz_max))
699: (*zz_max) = vertex[i].z;
700: if (vertex[i].x < (*xx_min))
701: (*xx_min) = vertex[i].x;
702: if (vertex[i].y < (*yy_min))
703: (*yy_min) = vertex[i].y;
704: if (vertex[i].z < (*zz_min))
705: (*zz_min) = vertex[i].z;
706: }
707: }
708: iso_lines = iso_lines->next;
709: }
710: /*
711: * fprintf(stderr," x: %g, %g\n", (*xx_min), (*xx_max));
712: * fprintf(stderr," y: %g, %g\n", (*yy_min), (*yy_max));
713: * fprintf(stderr," z: %g, %g\n", (*zz_min), (*zz_max));
714: */
715: }
716:
717: /*
718: * Generate new edge and append it to list, but only if both vertices are
719: * defined. The list is referenced by p_edge and pe_tail (p_edge points on
720: * first edge and pe_tail on last one).
721: * Note, the list may be empty (pe_edge==pe_tail==NULL) on entry and exit.
722: */
723: static struct edge_struct *add_edge(point0, point1, p_edge, pe_tail)
724: struct coordinate GPHUGE * point0; /* 2 vertices input */
725: struct coordinate GPHUGE * point1;
726: struct edge_struct **p_edge, **pe_tail; /* pointers to edge list in/out */
727: {
728: struct edge_struct *pe_temp = NULL;
729:
1.1.1.3 ! ohara 730: #if 1
! 731: if (point0->type == INRANGE && point1->type == INRANGE) {
! 732: #else
1.1 maekawa 733: if (point0->type != UNDEFINED && point1->type != UNDEFINED) {
1.1.1.3 ! ohara 734: #endif
1.1 maekawa 735:
736: pe_temp = (struct edge_struct *)
1.1.1.3 ! ohara 737: gp_alloc(sizeof(struct edge_struct), "contour edge");
1.1 maekawa 738:
739: pe_temp->poly[0] = NULL; /* clear links */
740: pe_temp->poly[1] = NULL;
741: pe_temp->vertex[0] = point0; /* First vertex of edge. */
742: pe_temp->vertex[1] = point1; /* Second vertex of edge. */
743: pe_temp->next = NULL;
744: pe_temp->position = INNER_MESH; /* default position in mesh */
745:
746: if ((*pe_tail)) {
747: (*pe_tail)->next = pe_temp; /* Stick new record as last one. */
748: } else {
749: (*p_edge) = pe_temp; /* start new list if empty */
750: }
751: (*pe_tail) = pe_temp; /* continue to last record. */
752:
753: }
754: return pe_temp; /* returns NULL, if no edge allocated */
755: }
756:
757: /*
758: * Generate new triangle and append it to list, but only if all edges are defined.
759: * The list is referenced by p_poly and pp_tail (p_poly points on first ploygon
760: * and pp_tail on last one).
761: * Note, the list may be empty (pe_ploy==pp_tail==NULL) on entry and exit.
762: */
763: static struct poly_struct *add_poly(edge0, edge1, edge2, p_poly, pp_tail)
764: struct edge_struct *edge0, *edge1, *edge2; /* 3 edges input */
765: struct poly_struct **p_poly, **pp_tail; /* pointers to polygon list in/out */
766: {
767: struct poly_struct *pp_temp = NULL;
768:
769: if (edge0 && edge1 && edge2) {
770:
771: pp_temp = (struct poly_struct *)
1.1.1.3 ! ohara 772: gp_alloc(sizeof(struct poly_struct), "contour polygon");
1.1 maekawa 773:
774: pp_temp->edge[0] = edge0; /* First edge of triangle */
775: pp_temp->edge[1] = edge1; /* Second one */
776: pp_temp->edge[2] = edge2; /* Third one */
777: pp_temp->next = NULL;
778:
779: if (edge0->poly[0]) /* update edge0 */
780: edge0->poly[1] = pp_temp;
781: else
782: edge0->poly[0] = pp_temp;
783:
784: if (edge1->poly[0]) /* update edge1 */
785: edge1->poly[1] = pp_temp;
786: else
787: edge1->poly[0] = pp_temp;
788:
789: if (edge2->poly[0]) /* update edge2 */
790: edge2->poly[1] = pp_temp;
791: else
792: edge2->poly[0] = pp_temp;
793:
794: if ((*pp_tail)) /* Stick new record as last one. */
795: (*pp_tail)->next = pp_temp;
796: else
797: (*p_poly) = pp_temp; /* start new list if empty */
798:
799: (*pp_tail) = pp_temp; /* continue to last record. */
800:
801: }
802: return pp_temp; /* returns NULL, if no edge allocated */
803: }
804:
805:
806:
807: /*
808: * Calls the (hopefully) desired interpolation/approximation routine.
809: */
810: static void put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind)
811: struct cntr_struct *p_cntr; /* contour structure input */
812: double z_level, /* Z level of contour input */
813: xx_min, xx_max, yy_min, yy_max; /* minimum/maximum values input */
814: int contr_kind; /* OPEN_ or CLOESED_CONTOUR input */
815: {
816:
817: if (!p_cntr)
818: return; /* Nothing to do if it is empty contour. */
819:
820: switch (interp_kind) {
821: case INTERP_NOTHING: /* No interpolation/approximation. */
822: put_contour_nothing(p_cntr);
823: break;
824: case INTERP_CUBIC: /* Cubic spline interpolation. */
825: put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max,
826: chk_contour_kind(p_cntr, contr_kind));
827:
828: break;
829: case APPROX_BSPLINE: /* Bspline approximation. */
830: put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max,
831: chk_contour_kind(p_cntr, contr_kind));
832: break;
833: }
834: free_contour(p_cntr);
835: }
836:
837: /*
838: * Simply puts contour coordinates in order with no interpolation or
839: * approximation.
840: */
841: static void put_contour_nothing(p_cntr)
842: struct cntr_struct *p_cntr;
843: {
844: while (p_cntr) {
845: add_cntr_point(p_cntr->X, p_cntr->Y);
846: p_cntr = p_cntr->next;
847: }
848: end_crnt_cntr();
849: }
850:
851: /*
852: * for some reason contours are never flagged as CLOSED_CONTOUR
853: * if first point == last point, set flag accordingly
854: *
855: */
856:
857: static int chk_contour_kind(p_cntr, contr_kind)
858: struct cntr_struct *p_cntr;
859: int contr_kind;
860: {
861: struct cntr_struct *pc_tail = NULL;
862: int current_contr_kind;
863:
864: FPRINTF((stderr, "check_contour_kind: current contr_kind value is %d\n", contr_kind));
865:
866: current_contr_kind = contr_kind;
867:
868: if (contr_kind != CLOSED_CONTOUR) {
869: pc_tail = p_cntr;
870: while (pc_tail->next)
871: pc_tail = pc_tail->next; /* Find last point. */
872:
873: /* test if first and last point are equal */
874: if (fuzzy_equal(pc_tail, p_cntr)) {
875: current_contr_kind = CLOSED_CONTOUR;
876: FPRINTF((stderr, "check_contour_kind: contr_kind changed to %d\n", current_contr_kind));
877: }
878: }
879: return (current_contr_kind);
880: }
881:
882: /*
883: * Generate a cubic spline curve through the points (x_i,y_i) which are
884: * stored in the linked list p_cntr.
885: * The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the
886: * parameter t is the length of the linear stroke.
887: */
888: static void put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind)
889: struct cntr_struct *p_cntr;
890: double z_level, xx_min, xx_max, yy_min, yy_max;
891: int contr_kind;
892: {
893: int num_pts, num_intpol;
894: double unit_x, unit_y; /* To define norm (x,y)-plane */
895: double *delta_t; /* Interval length t_{i+1}-t_i */
896: double *d2x, *d2y; /* Second derivatives x''(t_i), y''(t_i) */
897: struct cntr_struct *pc_tail;
898:
899: num_pts = count_contour(p_cntr); /* Number of points in contour. */
900:
901: pc_tail = p_cntr; /* Find last point. */
902: while (pc_tail->next)
903: pc_tail = pc_tail->next;
904:
905: if (contr_kind == CLOSED_CONTOUR) {
906: /* Test if first and last point are equal (should be) */
907: if (!fuzzy_equal(pc_tail, p_cntr)) {
908: pc_tail->next = p_cntr; /* Close contour list - make it circular. */
909: num_pts++;
910: }
911: }
912: delta_t = (double *)
1.1.1.3 ! ohara 913: gp_alloc((sizeof(double) * num_pts), "contour delta_t");
1.1 maekawa 914: d2x = (double *)
1.1.1.3 ! ohara 915: gp_alloc((sizeof(double) * num_pts), "contour d2x");
1.1 maekawa 916: d2y = (double *)
1.1.1.3 ! ohara 917: gp_alloc((sizeof(double) * num_pts), "contour d2y");
1.1 maekawa 918:
919: /* Width and hight of the grid is used at unit length (2d-norm) */
920: unit_x = xx_max - x_min;
921: unit_y = yy_max - y_min;
922: unit_x = (unit_x > zero ? unit_x : zero); /* should not be zero */
923: unit_y = (unit_y > zero ? unit_y : zero);
924:
925: if (num_pts > 2) {
926: /*
927: * Calculate second derivatives d2x[], d2y[] and interval lengths delta_t[]:
928: */
929: if (!gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t,
930: contr_kind, unit_x, unit_y)) {
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: /* If following (num_pts > 1) is TRUE then exactly 2 points in contour. */
940: else if (num_pts > 1) {
941: /* set all second derivatives to zero, interval length to 1 */
942: d2x[0] = 0.;
943: d2y[0] = 0.;
944: d2x[1] = 0.;
945: d2y[1] = 0.;
946: delta_t[0] = 1.;
947: } else { /* Only one point ( ?? ) - ignore it. */
948: free((char *) delta_t);
949: free((char *) d2x);
950: free((char *) d2y);
951: if (contr_kind == CLOSED_CONTOUR)
952: pc_tail->next = NULL; /* Un-circular list */
953: return;
954: }
955:
956: /* Calculate "num_intpol" interpolated values */
957: num_intpol = 1 + (num_pts - 1) * num_approx_pts; /* global: num_approx_pts */
958: intp_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, num_intpol);
959:
960: free((char *) delta_t);
961: free((char *) d2x);
962: free((char *) d2y);
963:
964: if (contr_kind == CLOSED_CONTOUR)
965: pc_tail->next = NULL; /* Un-circular list */
966:
967: end_crnt_cntr();
968: }
969:
970:
971: /*
972: * Find Bspline approximation for this data set.
973: * Uses global variable num_approx_pts to determine number of samples per
974: * interval, where the knot vector intervals are assumed to be uniform, and
975: * Global variable bspline_order for the order of Bspline to use.
976: */
977: static void put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_kind)
978: struct cntr_struct *p_cntr;
979: double z_level, xx_min, xx_max, yy_min, yy_max;
980: int contr_kind;
981: {
982: int num_pts, order = bspline_order;
983:
984: num_pts = count_contour(p_cntr); /* Number of points in contour. */
985: if (num_pts < 2)
986: return; /* Can't do nothing if empty or one points! */
987: /* Order must be less than number of points in curve - fix it if needed. */
988: if (order > num_pts - 1)
989: order = num_pts - 1;
990:
991: gen_bspline_approx(p_cntr, num_pts, order, contr_kind);
992: end_crnt_cntr();
993: }
994:
995: /*
996: * Free all elements in the contour list.
997: */
998: static void free_contour(p_cntr)
999: struct cntr_struct *p_cntr;
1000: {
1001: struct cntr_struct *pc_temp;
1002:
1003: while (p_cntr) {
1004: pc_temp = p_cntr;
1005: p_cntr = p_cntr->next;
1006: free((char *) pc_temp);
1007: }
1008: }
1009:
1010: /*
1011: * Counts number of points in contour.
1012: */
1013: static int count_contour(p_cntr)
1014: struct cntr_struct *p_cntr;
1015: {
1016: int count = 0;
1017:
1018: while (p_cntr) {
1019: count++;
1020: p_cntr = p_cntr->next;
1021: }
1022: return count;
1023: }
1024:
1025: /*
1026: * Find second derivatives (x''(t_i),y''(t_i)) of cubic spline interpolation
1027: * through list of points (x_i,y_i). The parameter t is calculated as the
1028: * length of the linear stroke. The number of points must be at least 3.
1029: * Note: For CLOSED_CONTOURs the first and last point must be equal.
1030: */
1031: static int gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, contr_kind, unit_x, unit_y)
1032: int num_pts; /* Number of points (num_pts>=3), input */
1033: struct cntr_struct *p_cntr; /* List of points (x(t_i),y(t_i)), input */
1034: double d2x[], d2y[], /* Second derivatives (x''(t_i),y''(t_i)), output */
1035: delta_t[]; /* List of interval lengths t_{i+1}-t_{i}, output */
1036: int contr_kind; /* CLOSED_CONTOUR or OPEN_CONTOUR, input */
1037: double unit_x, unit_y; /* Unit length in x and y (norm=1), input */
1038: {
1039: int n, i;
1040: double norm;
1041: tri_diag *m; /* The tri-diagonal matrix is saved here. */
1042: struct cntr_struct *pc_temp;
1043:
1044: m = (tri_diag *)
1.1.1.3 ! ohara 1045: gp_alloc((sizeof(tri_diag) * num_pts), "contour tridiag m");
1.1 maekawa 1046:
1047: /*
1048: * Calculate first differences in (d2x[i], d2y[i]) and interval lengths
1049: * in delta_t[i]:
1050: */
1051: pc_temp = p_cntr;
1052: for (i = 0; i < num_pts - 1; i++) {
1053: d2x[i] = pc_temp->next->X - pc_temp->X;
1054: d2y[i] = pc_temp->next->Y - pc_temp->Y;
1055: /*
1056: * The Norm of a linear stroke is calculated in "normal coordinates"
1057: * and used as interval length:
1058: */
1059: delta_t[i] = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y));
1060:
1061: d2x[i] /= delta_t[i]; /* first difference, with unit norm: */
1062: d2y[i] /= delta_t[i]; /* || (d2x[i], d2y[i]) || = 1 */
1063:
1064: pc_temp = pc_temp->next;
1065: }
1066:
1067: /*
1068: * Setup linear System: M * x = b
1069: */
1070: n = num_pts - 2; /* Without first and last point */
1071: if (contr_kind == CLOSED_CONTOUR) {
1072: /* First and last points must be equal for CLOSED_CONTOURs */
1073: delta_t[num_pts - 1] = delta_t[0];
1074: d2x[num_pts - 1] = d2x[0];
1075: d2y[num_pts - 1] = d2y[0];
1076: n++; /* Add last point (= first point) */
1077: }
1078: for (i = 0; i < n; i++) {
1079: /* Matrix M, mainly tridiagonal with cyclic second index ("j = j+n mod n") */
1080: m[i][0] = delta_t[i]; /* Off-diagonal element M_{i,i-1} */
1081: m[i][1] = 2. * (delta_t[i] + delta_t[i + 1]); /* M_{i,i} */
1082: m[i][2] = delta_t[i + 1]; /* Off-diagonal element M_{i,i+1} */
1083:
1084: /* Right side b_x and b_y */
1085: d2x[i] = (d2x[i + 1] - d2x[i]) * 6.;
1086: d2y[i] = (d2y[i + 1] - d2y[i]) * 6.;
1087:
1088: /*
1089: * If the linear stroke shows a cusps of more than 90 degree, the right
1090: * side is reduced to avoid oscillations in the spline:
1091: */
1092: norm = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y)) / 8.5;
1093:
1094: if (norm > 1.) {
1095: d2x[i] /= norm;
1096: d2y[i] /= norm;
1097: /* The first derivative will not be continuous */
1098: }
1099: }
1100:
1101: if (contr_kind != CLOSED_CONTOUR) {
1102: /* Third derivative is set to zero at both ends */
1103: m[0][1] += m[0][0]; /* M_{0,0} */
1104: m[0][0] = 0.; /* M_{0,n-1} */
1105: m[n - 1][1] += m[n - 1][2]; /* M_{n-1,n-1} */
1106: m[n - 1][2] = 0.; /* M_{n-1,0} */
1107: }
1108: /* Solve linear systems for d2x[] and d2y[] */
1109:
1110:
1111: if (solve_cubic_1(m, n)) { /* Calculate Cholesky decomposition */
1112: solve_cubic_2(m, d2x, n); /* solve M * d2x = b_x */
1113: solve_cubic_2(m, d2y, n); /* solve M * d2y = b_y */
1114:
1115: } else { /* Should not happen, but who knows ... */
1116: free((char *) m);
1117: return FALSE;
1118: }
1119:
1120: /* Shift all second derivatives one place right and abdate end points */
1121: for (i = n; i > 0; i--) {
1122: d2x[i] = d2x[i - 1];
1123: d2y[i] = d2y[i - 1];
1124: }
1125: if (contr_kind == CLOSED_CONTOUR) {
1126: d2x[0] = d2x[n];
1127: d2y[0] = d2y[n];
1128: } else {
1129: d2x[0] = d2x[1]; /* Third derivative is zero in */
1130: d2y[0] = d2y[1]; /* first and last interval */
1131: d2x[n + 1] = d2x[n];
1132: d2y[n + 1] = d2y[n];
1133: }
1134:
1135: free((char *) m);
1136: return TRUE;
1137: }
1138:
1139: /*
1140: * Calculate interpolated values of the spline function (defined via p_cntr
1141: * and the second derivatives d2x[] and d2y[]). The number of tabulated
1142: * values is n. On an equidistant grid n_intpol values are calculated.
1143: */
1144: static void intp_cubic_spline(n, p_cntr, d2x, d2y, delta_t, n_intpol)
1145: int n;
1146: struct cntr_struct *p_cntr;
1147: double d2x[], d2y[], delta_t[];
1148: int n_intpol;
1149: {
1150: double t, t_skip, t_max;
1151: double x0, x1, x, y0, y1, y;
1152: double d, hx, dx0, dx01, hy, dy0, dy01;
1153: int i;
1154:
1155: /* The length of the total interval */
1156: t_max = 0.;
1157: for (i = 0; i < n - 1; i++)
1158: t_max += delta_t[i];
1159:
1160: /* The distance between interpolated points */
1161: t_skip = (1. - 1e-7) * t_max / (n_intpol - 1);
1162:
1163: t = 0.; /* Parameter value */
1164: x1 = p_cntr->X;
1165: y1 = p_cntr->Y;
1166: add_cntr_point(x1, y1); /* First point. */
1167: t += t_skip;
1168:
1169: for (i = 0; i < n - 1; i++) {
1170: p_cntr = p_cntr->next;
1171:
1172: d = delta_t[i]; /* Interval length */
1173: x0 = x1;
1174: y0 = y1;
1175: x1 = p_cntr->X;
1176: y1 = p_cntr->Y;
1177: hx = (x1 - x0) / d;
1178: hy = (y1 - y0) / d;
1179: dx0 = (d2x[i + 1] + 2 * d2x[i]) / 6.;
1180: dy0 = (d2y[i + 1] + 2 * d2y[i]) / 6.;
1181: dx01 = (d2x[i + 1] - d2x[i]) / (6. * d);
1182: dy01 = (d2y[i + 1] - d2y[i]) / (6. * d);
1183: while (t <= delta_t[i]) { /* t in current interval ? */
1184: x = x0 + t * (hx + (t - d) * (dx0 + t * dx01));
1185: y = y0 + t * (hy + (t - d) * (dy0 + t * dy01));
1186: add_cntr_point(x, y); /* next point. */
1187: t += t_skip;
1188: }
1189: t -= delta_t[i]; /* Parameter t relative to start of next interval */
1190: }
1191: }
1192:
1193: /*
1194: * The following two procedures solve the special linear system which arise
1195: * in cubic spline interpolation. If x is assumed cyclic ( x[i]=x[n+i] ) the
1196: * equations can be written as (i=0,1,...,n-1):
1197: * m[i][0] * x[i-1] + m[i][1] * x[i] + m[i][2] * x[i+1] = b[i] .
1198: * In matrix notation one gets M * x = b, where the matrix M is tridiagonal
1199: * with additional elements in the upper right and lower left position:
1200: * m[i][0] = M_{i,i-1} for i=1,2,...,n-1 and m[0][0] = M_{0,n-1} ,
1201: * m[i][1] = M_{i, i } for i=0,1,...,n-1
1202: * m[i][2] = M_{i,i+1} for i=0,1,...,n-2 and m[n-1][2] = M_{n-1,0}.
1203: * M should be symmetric (m[i+1][0]=m[i][2]) and positiv definite.
1204: * The size of the system is given in n (n>=1).
1205: *
1206: * In the first procedure the Cholesky decomposition M = C^T * D * C
1207: * (C is upper triangle with unit diagonal, D is diagonal) is calculated.
1208: * Return TRUE if decomposition exist.
1209: */
1210: static int solve_cubic_1(m, n)
1211: tri_diag m[];
1212: int n;
1213: {
1214: int i;
1215: double m_ij, m_n, m_nn, d;
1216:
1217: if (n < 1)
1218: return FALSE; /* Dimension should be at least 1 */
1219:
1220: d = m[0][1]; /* D_{0,0} = M_{0,0} */
1221: if (d <= 0.)
1222: return FALSE; /* M (or D) should be positiv definite */
1223: m_n = m[0][0]; /* M_{0,n-1} */
1224: m_nn = m[n - 1][1]; /* M_{n-1,n-1} */
1225: for (i = 0; i < n - 2; i++) {
1226: m_ij = m[i][2]; /* M_{i,1} */
1227: m[i][2] = m_ij / d; /* C_{i,i+1} */
1228: m[i][0] = m_n / d; /* C_{i,n-1} */
1229: m_nn -= m[i][0] * m_n; /* to get C_{n-1,n-1} */
1230: m_n = -m[i][2] * m_n; /* to get C_{i+1,n-1} */
1231: d = m[i + 1][1] - m[i][2] * m_ij; /* D_{i+1,i+1} */
1232: if (d <= 0.)
1233: return FALSE; /* Elements of D should be positiv */
1234: m[i + 1][1] = d;
1235: }
1236: if (n >= 2) { /* Complete last column */
1237: m_n += m[n - 2][2]; /* add M_{n-2,n-1} */
1238: m[n - 2][0] = m_n / d; /* C_{n-2,n-1} */
1239: m[n - 1][1] = d = m_nn - m[n - 2][0] * m_n; /* D_{n-1,n-1} */
1240: if (d <= 0.)
1241: return FALSE;
1242: }
1243: return TRUE;
1244: }
1245:
1246: /*
1247: * The second procedure solves the linear system, with the Choleky
1248: * decomposition calculated above (in m[][]) and the right side b given
1249: * in x[]. The solution x overwrites the right side in x[].
1250: */
1251: static void solve_cubic_2(m, x, n)
1252: tri_diag m[];
1253: double x[];
1254: int n;
1255: {
1256: int i;
1257: double x_n;
1258:
1259: /* Division by transpose of C : b = C^{-T} * b */
1260: x_n = x[n - 1];
1261: for (i = 0; i < n - 2; i++) {
1262: x[i + 1] -= m[i][2] * x[i]; /* C_{i,i+1} * x_{i} */
1263: x_n -= m[i][0] * x[i]; /* C_{i,n-1} * x_{i} */
1264: }
1265: if (n >= 2)
1266: x[n - 1] = x_n - m[n - 2][0] * x[n - 2]; /* C_{n-2,n-1} * x_{n-1} */
1267:
1268: /* Division by D: b = D^{-1} * b */
1269: for (i = 0; i < n; i++)
1270: x[i] /= m[i][1];
1271:
1272: /* Division by C: b = C^{-1} * b */
1273: x_n = x[n - 1];
1274: if (n >= 2)
1275: x[n - 2] -= m[n - 2][0] * x_n; /* C_{n-2,n-1} * x_{n-1} */
1276: for (i = n - 3; i >= 0; i--) {
1277: /* C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1} */
1278: x[i] -= m[i][2] * x[i + 1] + m[i][0] * x_n;
1279: }
1280: return;
1281: }
1282:
1283: /*
1284: * Solve tri diagonal linear system equation. The tri diagonal matrix is
1285: * defined via matrix M, right side is r, and solution X i.e. M * X = R.
1286: * Size of system given in n. Return TRUE if solution exist.
1287: */
1288: /* not used any more in "contour.c", but in "spline.c" (21. Dec. 1995) ! */
1289:
1290: int solve_tri_diag(m, r, x, n)
1291: tri_diag m[];
1292: double r[], x[];
1293: int n;
1294: {
1295: int i;
1296: double t;
1297:
1298: for (i = 1; i < n; i++) { /* Eliminate element m[i][i-1] (lower diagonal). */
1299: if (m[i - 1][1] == 0)
1300: return FALSE;
1301: t = m[i][0] / m[i - 1][1]; /* Find ratio between the two lines. */
1302: /* m[i][0] = m[i][0] - m[i-1][1] * t; */
1303: /* m[i][0] is not used any more (and set to 0 in the above line) */
1304: m[i][1] = m[i][1] - m[i - 1][2] * t;
1305: r[i] = r[i] - r[i - 1] * t;
1306: }
1307: /* Now do back subtitution - update the solution vector X: */
1308: if (m[n - 1][1] == 0)
1309: return FALSE;
1310: x[n - 1] = r[n - 1] / m[n - 1][1]; /* Find last element. */
1311: for (i = n - 2; i >= 0; i--) {
1312: if (m[i][1] == 0)
1313: return FALSE;
1314: x[i] = (r[i] - x[i + 1] * m[i][2]) / m[i][1];
1315: }
1316: return TRUE;
1317: }
1318:
1319: /*
1320: * Generate a Bspline curve defined by all the points given in linked list p:
1321: * Algorithm: using deBoor algorithm
1322: * Note: if Curvekind is OPEN_CONTOUR than Open end knot vector is assumed,
1323: * else (CLOSED_CONTOUR) Float end knot vector is assumed.
1324: * It is assumed that num_of_points is at least 2, and order of Bspline is less
1325: * than num_of_points!
1326: */
1327: static void gen_bspline_approx(p_cntr, num_of_points, order, contr_kind)
1328: struct cntr_struct *p_cntr;
1329: int num_of_points, order, contr_kind;
1330: {
1331: int knot_index = 0, pts_count = 1;
1332: double dt, t, next_t, t_min, t_max, x, y;
1333: struct cntr_struct *pc_temp = p_cntr, *pc_tail = NULL;
1334:
1335: /* If the contour is Closed one we must update few things:
1336: * 1. Make the list temporary circular, so we can close the contour.
1337: * 2. Update num_of_points - increase it by "order-1" so contour will be
1338: * closed. This will evaluate order more sections to close it!
1339: */
1340: if (contr_kind == CLOSED_CONTOUR) {
1341: pc_tail = p_cntr;
1342: while (pc_tail->next)
1343: pc_tail = pc_tail->next; /* Find last point. */
1344:
1345: /* test if first and last point are equal */
1346: if (fuzzy_equal(pc_tail, p_cntr)) {
1347: /* Close contour list - make it circular. */
1348: pc_tail->next = p_cntr->next;
1349: num_of_points += order - 1;
1350: } else {
1351: pc_tail->next = p_cntr;
1352: num_of_points += order;
1353: }
1354: }
1355: /* Find first (t_min) and last (t_max) t value to eval: */
1356: t = t_min = fetch_knot(contr_kind, num_of_points, order, order);
1357: t_max = fetch_knot(contr_kind, num_of_points, order, num_of_points);
1358: next_t = t_min + 1.0;
1359: knot_index = order;
1360: dt = 1.0 / num_approx_pts; /* Number of points per one section. */
1361:
1362:
1363: while (t < t_max) {
1364: if (t > next_t) {
1365: pc_temp = pc_temp->next; /* Next order ctrl. pt. to blend. */
1366: knot_index++;
1367: next_t += 1.0;
1368: }
1369: eval_bspline(t, pc_temp, num_of_points, order, knot_index,
1370: contr_kind, &x, &y); /* Next pt. */
1371: add_cntr_point(x, y);
1372: pts_count++;
1373: /* As we might have some real number round off problems we do */
1374: /* the last point outside the loop */
1375: if (pts_count == num_approx_pts * (num_of_points - order) + 1)
1376: break;
1377: t += dt;
1378: }
1379:
1380: /* Now do the last point */
1381: eval_bspline(t_max - EPSILON, pc_temp, num_of_points, order, knot_index,
1382: contr_kind, &x, &y);
1383: add_cntr_point(x, y); /* Complete the contour. */
1384:
1385: if (contr_kind == CLOSED_CONTOUR) /* Update list - un-circular it. */
1386: pc_tail->next = NULL;
1387: }
1388:
1389: /*
1390: * The routine to evaluate the B-spline value at point t using knot vector
1391: * from function fetch_knot(), and the control points p_cntr.
1392: * Returns (x, y) of approximated B-spline. Note that p_cntr points on the
1393: * first control point to blend with. The B-spline is of order order.
1394: */
1395: static void eval_bspline(t, p_cntr, num_of_points, order, j, contr_kind, x, y)
1396: double t;
1397: struct cntr_struct *p_cntr;
1398: int num_of_points, order, j, contr_kind;
1399: double *x, *y;
1400: {
1401: int i, p;
1402: double ti, tikp, *dx, *dy; /* Copy p_cntr into it to make it faster. */
1403:
1404: dx = (double *)
1.1.1.3 ! ohara 1405: gp_alloc((sizeof(double) * (order + j)), "contour b_spline");
1.1 maekawa 1406: dy = (double *)
1.1.1.3 ! ohara 1407: gp_alloc((sizeof(double) * (order + j)), "contour b_spline");
1.1 maekawa 1408:
1409: /* Set the dx/dy - [0] iteration step, control points (p==0 iterat.): */
1410: for (i = j - order; i <= j; i++) {
1411: dx[i] = p_cntr->X;
1412: dy[i] = p_cntr->Y;
1413: p_cntr = p_cntr->next;
1414: }
1415:
1416: for (p = 1; p <= order; p++) { /* Iteration (b-spline level) counter. */
1417: for (i = j; i >= j - order + p; i--) { /* Control points indexing. */
1418: ti = fetch_knot(contr_kind, num_of_points, order, i);
1419: tikp = fetch_knot(contr_kind, num_of_points, order, i + order + 1 - p);
1420: if (ti == tikp) { /* Should not be a problems but how knows... */
1421: } else {
1422: dx[i] = dx[i] * (t - ti) / (tikp - ti) + /* Calculate x. */
1423: dx[i - 1] * (tikp - t) / (tikp - ti);
1424: dy[i] = dy[i] * (t - ti) / (tikp - ti) + /* Calculate y. */
1425: dy[i - 1] * (tikp - t) / (tikp - ti);
1426: }
1427: }
1428: }
1429: *x = dx[j];
1430: *y = dy[j];
1431: free((char *) dx);
1432: free((char *) dy);
1433: }
1434:
1435: /*
1436: * Routine to get the i knot from uniform knot vector. The knot vector
1437: * might be float (Knot(i) = i) or open (where the first and last "order"
1438: * knots are equal). contr_kind determines knot kind - OPEN_CONTOUR means
1439: * open knot vector, and CLOSED_CONTOUR selects float knot vector.
1440: * Note the knot vector is not exist and this routine simulates it existance
1441: * Also note the indexes for the knot vector starts from 0.
1442: */
1443: static double fetch_knot(contr_kind, num_of_points, order, i)
1444: int contr_kind, num_of_points, order, i;
1445: {
1446: switch (contr_kind) {
1447: case OPEN_CONTOUR:
1448: if (i <= order)
1449: return 0.0;
1450: else if (i <= num_of_points)
1451: return (double) (i - order);
1452: else
1453: return (double) (num_of_points - order);
1454: case CLOSED_CONTOUR:
1455: return (double) i;
1456: default: /* Should never happen */
1457: return 1.0;
1458: }
1459: #ifdef sequent
1460: return 1.0; /* ???? */
1461: #endif
1462: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>