source: Sophya/trunk/SophyaPI/PIGcont/contour.c@ 1837

Last change on this file since 1837 was 1829, checked in by perderos, 24 years ago

Creation du module PIGcont, Drawer de trace de contour / a partir
de l'adaptation du code de trace de contour de GNUplot, interfacage
avec la librairie PI .

Olivier Perdereau 19/12/2001

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