Annotation of OpenXM_contrib/gnuplot/contour.c, Revision 1.1.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