Context Navigation

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

Visit:

Last change on this file since 1838 was 1838, checked in by ansari, 24 years ago
Nettoyage des .h en trop, compil avec g++ sous LinuxPPC (rename alloc.h gp_alloc.h) - Reza 21/12/2001
File size: 48.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: