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

#ifndef lint
static char *RCSid() { return RCSid("$Id: contour.c,v 1.1.1.1 2001-12-19 18:01:18 perderos Exp $"); }
#endif

/* GNUPLOT - contour.c */

/*[
 * Copyright 1986 - 1993, 1998   Thomas Williams, Colin Kelley
 *
 * Permission to use, copy, and distribute this software and its
 * documentation for any purpose with or without fee is hereby granted,
 * provided that the above copyright notice appear in all copies and
 * that both that copyright notice and this permission notice appear
 * in supporting documentation.
 *
 * Permission to modify the software is granted, but not the right to
 * distribute the complete modified source code.  Modifications are to
 * be distributed as patches to the released version.  Permission to
 * distribute binaries produced by compiling modified sources is granted,
 * provided you
 *   1. distribute the corresponding source modifications from the
 *    released version in the form of a patch file along with the binaries,
 *   2. add special version identification to distinguish your version
 *    in addition to the base release version number,
 *   3. provide your name and address as the primary contact for the
 *    support of your modified version, and
 *   4. retain our contact information in regard to use of the base
 *    software.
 * Permission to distribute the released version of the source code along
 * with corresponding source modifications in the form of a patch file is
 * granted with same provisions 2 through 4 for binary distributions.
 *
 * This software is provided "as is" without express or implied warranty
 * to the extent permitted by applicable law.
]*/

 
/* 
 * AUTHORS
 * 
 *   Original Software:
 *       Gershon Elber
 * 
 *   Improvements to the numerical algorithms:
 *        Hans-Martin Keller, 1995,1997 (hkeller@gwdg.de)
 *
 */

#include "contour.h"

#include "alloc.h"
#include "myaxis.h"
/*  #include "setshow.h" */

/* exported variables (to be handled by the 'set' and friends): */

char contour_format[32] = "%8.3g";      /* format for contour key entries */
t_contour_kind contour_kind = CONTOUR_KIND_LINEAR;
t_contour_levels_kind contour_levels_kind = LEVELS_AUTO;
int contour_levels = DEFAULT_CONTOUR_LEVELS;
int contour_order = DEFAULT_CONTOUR_ORDER;
int contour_pts = DEFAULT_NUM_APPROX_PTS;

dynarray dyn_contour_levels_list;/* storage for z levels to draw contours at */
double * contour_levels_list=NULL;

/* position of edge in mesh */
typedef enum en_edge_position {
    INNER_MESH=1,
    BOUNDARY,
    DIAGONAL
} t_edge_position;


/* FIXME HBB 2000052: yet another local copy of 'epsilon'. Why? */
#define EPSILON  1e-5           /* Used to decide if two float are equal. */

#ifndef TRUE
#define TRUE     -1
#define FALSE    0
#endif

#define MAX_POINTS_PER_CNTR     100

#define SQR(x)  ((x) * (x))

/*
 * struct vrtx_struct {
 *      double X, Y, Z;
 *      struct vrtx_struct *next;
 * };
 * 
 * replaced by 'struct coordinate  ', see plot.h (HMK 1997) 
 */

struct edge_struct {
    struct poly_struct *poly[2]; /* Each edge belongs to up to 2 polygons */
    struct coordinate  *vertex[2]; /* The two extreme points of this edge. */
    struct edge_struct *next;   /* To chain lists */
    TBOOLEAN is_active;         /* is edge is 'active' at certain Z level? */
    t_edge_position position;   /* position of edge in mesh */
};

struct poly_struct {
    struct edge_struct *edge[3];        /* As we do triangolation here... */
    struct poly_struct *next;   /* To chain lists. */
};

struct cntr_struct {            /* Contours are saved using this struct list. */
    double X, Y;                /* The coordinates of this vertex. */
    struct cntr_struct *next;   /* To chain lists. */
};

static struct gnuplot_contours *contour_list = NULL;
static double crnt_cntr[MAX_POINTS_PER_CNTR * 2];
static int crnt_cntr_pt_index = 0;
static double contour_level = 0.0;

/* Linear, Cubic interp., Bspline: */
static t_contour_kind interp_kind = CONTOUR_KIND_LINEAR;

static double x_min, y_min, z_min;      /* Minimum values of x, y, and z */
static double x_max, y_max, z_max;      /* Maximum values of x, y, and z */

static void add_cntr_point (double x, double y);
static void end_crnt_cntr (void);
static void gen_contours(struct edge_struct * p_edges, double z_level,
                                  double xx_min, double xx_max, double yy_min, double yy_max);
static int update_all_edges(struct edge_struct * p_edges,
                                     double z_level);
static struct cntr_struct *gen_one_contour (
                                                       struct edge_struct * p_edges, double
                                                       z_level, TBOOLEAN *contr_isclosed,
                                                       int *num_active);
static struct cntr_struct *trace_contour (
                                                     struct edge_struct * pe_start, double
                                                     z_level, int *num_active,
                                                     TBOOLEAN contr_isclosed);
static struct cntr_struct *update_cntr_pt (struct edge_struct * p_edge,
                                                   double z_level);
static int fuzzy_equal (struct cntr_struct * p_cntr1,
                                struct cntr_struct * p_cntr2);


static void gen_triangle (int num_isolines,
                                  struct iso_curve * iso_lines, struct poly_struct ** p_polys,
                                  struct edge_struct ** p_edges);
static void calc_min_max (int num_isolines,
                                  struct iso_curve * iso_lines, double *xx_min, double *yy_min,
                                  double *zz_min,
                                  double *xx_max, double *yy_max, double *zz_max);
static struct edge_struct *add_edge (struct coordinate  * point0,
                                             struct coordinate  * point1, struct edge_struct
                                             ** p_edge,
                                             struct edge_struct ** pe_tail);
static struct poly_struct *add_poly (struct edge_struct * edge0,
                                             struct edge_struct * edge1, struct edge_struct * edge2,
                                             struct poly_struct ** p_poly, struct poly_struct ** pp_tail);


static void put_contour (struct cntr_struct * p_cntr, double z_level,
                                 double xx_min, double xx_max, double yy_min, double yy_max,
                                 TBOOLEAN contr_isclosed);
static void put_contour_nothing (struct cntr_struct * p_cntr);
static int chk_contour_kind (struct cntr_struct * p_cntr,
                                     TBOOLEAN contr_isclosed);
static void put_contour_cubic (struct cntr_struct * p_cntr,
                                       double z_level, double xx_min, double xx_max, double
                                       yy_min, double yy_max,
                                       TBOOLEAN contr_isclosed);
static void put_contour_bspline (struct cntr_struct * p_cntr,
                                         double z_level, double xx_min, double xx_max, double
                                         yy_min, double yy_max,
                                         TBOOLEAN contr_isclosed);
static void free_contour (struct cntr_struct * p_cntr);
static int count_contour (struct cntr_struct * p_cntr);
static int gen_cubic_spline (int num_pts, struct cntr_struct * p_cntr,
                                     double d2x[], double d2y[], double delta_t[], TBOOLEAN contr_isclosed,
                                     double unit_x, double unit_y);
static void intp_cubic_spline (int n, struct cntr_struct * p_cntr,
                                       double d2x[], double d2y[], double delta_t[], int n_intpol);
static int solve_cubic_1 (tri_diag m[], int n);
static void solve_cubic_2 (tri_diag m[], double x[], int n);
static void gen_bspline_approx (struct cntr_struct * p_cntr,
                                        int num_of_points, int order, TBOOLEAN contr_isclosed);
static void eval_bspline (double t, struct cntr_struct * p_cntr,
                                  int num_of_points, int order, int j, TBOOLEAN contr_isclosed, double *x,
                                  double *y);
static double fetch_knot (TBOOLEAN contr_isclosed, int num_of_points,
                                  int order, int i);


static int num_of_z_levels;/*_____ OP ___________*/     /* # Z contour levels. */

int Get_Num_Of_Z_Levels(){
return num_of_z_levels;
} /* OP       __________ */

/*
 * Entry routine to this whole set of contouring module.
 */
struct gnuplot_contours *
contour(num_isolines, iso_lines)
int num_isolines;
struct iso_curve *iso_lines;
{
    int i;
    /*OP    int num_of_z_levels;*/      /* # Z contour levels. */
    struct poly_struct *p_polys, *p_poly;
    struct edge_struct *p_edges, *p_edge;
    double z = 0, dz = 0;
    struct gnuplot_contours *save_contour_list;

    num_of_z_levels = contour_levels;
    interp_kind = contour_kind;

    contour_list = NULL;

    /* 
     * Calculate min/max values :
     */
    calc_min_max(num_isolines, iso_lines,
                 &x_min, &y_min, &z_min, &x_max, &y_max, &z_max);

    dz = fabs(z_max - z_min);
    /* 
     *  printf(" contour z_max %g z_min %g dz=%g kind %d   \n",z_max, z_min , dz,contour_levels_kind);
     * Generate list of edges (p_edges) and list of triangles (p_polys):
     */

    gen_triangle(num_isolines, iso_lines, &p_polys, &p_edges);
    crnt_cntr_pt_index = 0;
    /*AJOUT OP */
    if (contour_levels_kind == LEVELS_NUM) {
        dz = fabs(z_max - z_min)/(num_of_z_levels);
        z = z_min - dz/2.;
    }

    if (contour_levels_kind == LEVELS_AUTO) {
        dz = fabs(z_max - z_min);
        /*printf(" contour z_max %g z_min %g dz=%g\n",z_max, z_min , dz);*/
        if (dz == 0)
            return NULL;        /* empty z range ? */
        /* what is the deeper sense of this ? (joze) */
        dz = set_tic(log10(dz), ((int) contour_levels + 1) * 2);
        z = floor(z_min / dz) * dz;
        num_of_z_levels = (int) floor((z_max - z) / dz);
        /*printf("contour() : num_of_z_levels %d\n",num_of_z_levels);*/
    }
    for (i = 0; i < num_of_z_levels; i++) {
        switch (contour_levels_kind) {
        case LEVELS_AUTO:
        case LEVELS_NUM:
            z += dz;
            break;
        case LEVELS_INCREMENTAL:
            z = contour_levels_list[0] + i * contour_levels_list[1];
            break;
        case LEVELS_DISCRETE:
          /*printf("????? contour z=%f\n",contour_levels_list[i] );*/
            /*z = AXIS_LOG_VALUE(FIRST_Z_AXIS, contour_levels_list[i]); PAS BESOIN ? OP */
            z = contour_levels_list[i];
            break;
        }
        contour_level = z;
        /*printf(" contour z=%f\n",z);*/
        save_contour_list = contour_list;
        gen_contours(p_edges, z, x_min, x_max, y_min, y_max);
        if (contour_list != save_contour_list) {
            contour_list->isNewLevel = 1;
            sprintf(contour_list->label, contour_format, AXIS_DE_LOG_VALUE(FIRST_Z_AXIS,z));
#ifdef PM3D
            contour_list->z = AXIS_DE_LOG_VALUE(FIRST_Z_AXIS, z);
#endif
        }
    }

    /* Free all contouring related temporary data. */ 
    while (p_polys) {
        p_poly = p_polys->next;
        free(p_polys);
        p_polys = p_poly; 
    }
    while (p_edges) {
        p_edge = p_edges->next;
        free(p_edges);
        p_edges = p_edge;
    }

    return contour_list;
}

/*
 * Adds another point to the currently build contour.
 */
static void
add_cntr_point(x, y)
double x, y;
{
    int index;

    if (crnt_cntr_pt_index >= MAX_POINTS_PER_CNTR - 1) {
        index = crnt_cntr_pt_index - 1;
        end_crnt_cntr();
        crnt_cntr[0] = crnt_cntr[index * 2];
        crnt_cntr[1] = crnt_cntr[index * 2 + 1];
        crnt_cntr_pt_index = 1; /* Keep the last point as first of this one. */
    }
    crnt_cntr[crnt_cntr_pt_index * 2] = x;
    crnt_cntr[crnt_cntr_pt_index * 2 + 1] = y;
    crnt_cntr_pt_index++;
}
 
/*
 * Done with current contour - create gnuplot data structure for it.
 */
static void
end_crnt_cntr()
{
    int i;
    struct gnuplot_contours *cntr = (struct gnuplot_contours *)
    gp_alloc(sizeof(struct gnuplot_contours), "gnuplot_contour");
    cntr->coords = (struct coordinate  *)
        gp_alloc(sizeof(struct coordinate) * crnt_cntr_pt_index,
                 "contour coords");

    for (i = 0; i < crnt_cntr_pt_index; i++) {
        cntr->coords[i].x = crnt_cntr[i * 2];
        cntr->coords[i].y = crnt_cntr[i * 2 + 1];
        cntr->coords[i].z = contour_level;
    }
    cntr->num_pts = crnt_cntr_pt_index;

    cntr->next = contour_list;
    contour_list = cntr;
    contour_list->isNewLevel = 0;

    crnt_cntr_pt_index = 0;
}

/*
 * Generates all contours by tracing the intersecting triangles.
 */
static void
gen_contours(p_edges, z_level, xx_min, xx_max, yy_min, yy_max)
struct edge_struct *p_edges;
double z_level, xx_min, xx_max, yy_min, yy_max;
{
    int num_active;             /* Number of edges marked ACTIVE. */
    TBOOLEAN contr_isclosed;    /* Is this contour a closed line? */
    struct cntr_struct *p_cntr;

    num_active = update_all_edges(p_edges, z_level);    /* Do pass 1. */

    contr_isclosed = FALSE;     /* Start to look for contour on boundaries. */
    /*printf("<gen_contour> z=%g num_active %d \n",z_level,num_active);*/
    while (num_active > 0) {    /* Do Pass 2. */
        /* Generate One contour (and update MumActive as needed): */
        p_cntr = gen_one_contour(p_edges, z_level, &contr_isclosed, &num_active);
        if (p_cntr ==NULL) printf("<gen_contour> gen_one_contour retourne NULL \n");
        /* Emit it in requested format: */
        put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_isclosed);
    }
}

/*
 * Does pass 1, or marks the edges which are active (crosses this z_level)
 * Returns number of active edges (marked ACTIVE).
 */
static int
update_all_edges(p_edges, z_level)
struct edge_struct *p_edges;
double z_level;
{
    int count = 0;

    while (p_edges) {
        /* use the same test at both vertices to avoid roundoff errors */
        
        if ((p_edges->vertex[0]->z >= z_level) !=
            (p_edges->vertex[1]->z >= z_level)) {
            p_edges->is_active = TRUE;
            count++;
        } else
            p_edges->is_active = FALSE;
        p_edges = p_edges->next;
    }

    return count;
}

/*
 * Does pass 2, or find one complete contour out of the triangulation
 * data base:
 *
 * Returns a pointer to the contour (as linked list), contr_isclosed
 * tells if the contour is a closed line or not, and num_active is
 * updated.  
 */
static struct cntr_struct *
gen_one_contour(p_edges, z_level, contr_isclosed, num_active)
struct edge_struct *p_edges;    /* list of edges input */
double z_level;                 /* Z level of contour input */
TBOOLEAN *contr_isclosed;       /* open or closed contour, in/out */
int *num_active;                /* number of active edges     in/out */
{
    struct edge_struct *pe_temp;
    
    if (! *contr_isclosed) {
      /*printf("<gen_one_contour> contr_isclosed FALSE \n");*/
        /* Look for something to start with on boundary: */
        pe_temp = p_edges;
        while (pe_temp) {
            if (pe_temp->is_active && (pe_temp->position == BOUNDARY))
                break;
            pe_temp = pe_temp->next;
        }
        if (!pe_temp)
            *contr_isclosed = TRUE;     /* No more contours on boundary. */
        else {
          /*printf("<gen_one_contour> contr_isclosed FALSE return \n");*/
            return trace_contour(pe_temp, z_level, num_active, *contr_isclosed);
        }
    }
    if (*contr_isclosed) {
        /* Look for something to start with inside: */
      /*printf("<gen_one_contour> contr_isclosed TRUE \n");*/
        pe_temp = p_edges;
        while (pe_temp) {
            if (pe_temp->is_active && (pe_temp->position != BOUNDARY))
                break;
            pe_temp = pe_temp->next;
        }
        if (!pe_temp) {
            *num_active = 0;
            fprintf(stderr, "gen_one_contour: no contour found\n");
            return NULL;
        } else {
            *contr_isclosed = TRUE;
            return trace_contour(pe_temp, z_level, num_active, *contr_isclosed);
        }
    }
    printf("<gen_one_contour> We should never be here, but lint... \n"); 
    return NULL;                /* We should never be here, but lint... */
}

/*
 * Search the data base along a contour starts at the edge pe_start until
 * a boundary edge is detected or until we close the loop back to pe_start.
 * Returns a linked list of all the points on the contour
 * Also decreases num_active by the number of points on contour.
 */
static struct cntr_struct *
trace_contour(pe_start, z_level, num_active, contr_isclosed)
    struct edge_struct *pe_start; /* edge to start contour input */
    double z_level;             /* Z level of contour input */
    int *num_active;            /* number of active edges in/out */
    TBOOLEAN contr_isclosed;    /* open or closed contour line (input) */
{
    struct cntr_struct *p_cntr, *pc_tail;
    struct edge_struct *p_edge, *p_next_edge;
    struct poly_struct *p_poly, *PLastpoly = NULL;
    int i;

    p_edge = pe_start;          /* first edge to start contour */

    /* Generate the header of the contour - the point on pe_start. */
    if (! contr_isclosed) {
        pe_start->is_active = FALSE;
        (*num_active)--;
    }
    if (p_edge->poly[0] || p_edge->poly[1]) {   /* more than one point */

        p_cntr = pc_tail = update_cntr_pt(pe_start, z_level);   /* first point */

        do {
            /* Find polygon to continue (Not where we came from - PLastpoly): */
            if (p_edge->poly[0] == PLastpoly)
                p_poly = p_edge->poly[1];
            else
                p_poly = p_edge->poly[0];
            p_next_edge = NULL; /* In case of error, remains NULL. */
            for (i = 0; i < 3; i++)     /* Test the 3 edges of the polygon: */
                if (p_poly->edge[i] != p_edge)
                    if (p_poly->edge[i]->is_active)
                        p_next_edge = p_poly->edge[i];
            if (!p_next_edge) { /* Error exit */
                pc_tail->next = NULL;
                free_contour(p_cntr);
                fprintf(stderr, "trace_contour: unexpected end of contour\n");
                return NULL;
            }
            p_edge = p_next_edge;
            PLastpoly = p_poly;
            p_edge->is_active = FALSE;
            (*num_active)--;

            /* Do not allocate contour points on diagonal edges */
            if (p_edge->position != DIAGONAL) {

                pc_tail->next = update_cntr_pt(p_edge, z_level);

                /* Remove nearby points */
                if (fuzzy_equal(pc_tail, pc_tail->next)) {

                    free((char *) pc_tail->next);
                } else
                    pc_tail = pc_tail->next;
            }
        } while ((p_edge != pe_start) && (p_edge->position != BOUNDARY));

        pc_tail->next = NULL;

        /* For closed contour the first and last point should be equal */
        if (pe_start == p_edge) {
            (p_cntr->X) = (pc_tail->X);
            (p_cntr->Y) = (pc_tail->Y);
        }
    } else {                    /* only one point, forget it */
        p_cntr = NULL;
    }

    return p_cntr;
}

/*
 * Allocates one contour location and update it to to correct position
 * according to z_level and edge p_edge.
 */
static struct cntr_struct *
update_cntr_pt(p_edge, z_level)
struct edge_struct *p_edge;
double z_level;
{
    double t;
    struct cntr_struct *p_cntr;

    t = (z_level - p_edge->vertex[0]->z) /
        (p_edge->vertex[1]->z - p_edge->vertex[0]->z);

    /* test if t is out of interval [0:1] (should not happen but who knows ...) */
    /*if(t>1) printf(" <update_cntr_pt> t >1 !\n");*/
    /*if(t<0) printf(" <update_cntr_pt> t negatif !\n");*/
    t = (t < 0.0 ? 0.0 : t);
    t = (t > 1.0 ? 1.0 : t);
    /*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);*/
    /*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);*/
    p_cntr = (struct cntr_struct *)
        gp_alloc(sizeof(struct cntr_struct), "contour cntr_struct");

    p_cntr->X = p_edge->vertex[1]->x * t +
        p_edge->vertex[0]->x * (1 - t);
    p_cntr->Y = p_edge->vertex[1]->y * t +
        p_edge->vertex[0]->y * (1 - t);
    /*printf(" <update_cntr_pt> p_cntr X %g Y %g \n",p_cntr->X,p_cntr->Y);*/
    return p_cntr;
}

/* Simple routine to decide if two contour points are equal by
 * calculating the relative error (< EPSILON).  */
/* HBB 20010121: don't use absolute value 'zero' to compare to data
 * values. */
static int
fuzzy_equal(p_cntr1, p_cntr2)
    struct cntr_struct *p_cntr1, *p_cntr2;
{
    double unit_x, unit_y;
    unit_x = fabs(x_max - x_min);               /* reference */
    unit_y = fabs(y_max - y_min);
    return ((fabs(p_cntr1->X - p_cntr2->X) < unit_x * EPSILON)
            && (fabs(p_cntr1->Y - p_cntr2->Y) < unit_y * EPSILON));
}

/*
 * Generate the triangles.
 * Returns the lists (edges & polys) via pointers to their heads.
 */
static void
gen_triangle(num_isolines, iso_lines, p_polys, p_edges)
int num_isolines;               /* number of iso-lines input */
struct iso_curve *iso_lines;    /* iso-lines input */
struct poly_struct **p_polys;   /* list of polygons output */
struct edge_struct **p_edges;   /* list of edges output */
{
    int i, j, grid_x_max = iso_lines->p_count;
    struct edge_struct *p_edge1, *p_edge2, *edge0, *edge1, *edge2, *pe_tail,
           *pe_tail2, *pe_temp;
    struct poly_struct *pp_tail, *lower_tri, *upper_tri;
    /* HBB 980308: need to tag *each* of them as ! */
    struct coordinate  *p_vrtx1,  * p_vrtx2;

    (*p_polys) = pp_tail = NULL;        /* clear lists */
    (*p_edges) = pe_tail = NULL;

    p_vrtx1 = iso_lines->points;        /* first row of vertices */
    p_edge1 = pe_tail = NULL;   /* clear list of edges */

    /* Generate edges of first row */
    for (j = 0; j < grid_x_max - 1; j++)
        add_edge(p_vrtx1 + j, p_vrtx1 + j + 1, &p_edge1, &pe_tail);

    (*p_edges) = p_edge1;       /* update main list */


    /*
     * Combines vertices to edges and edges to triangles:
     * ==================================================
     * The edges are stored in the edge list, referenced by p_edges
     * (pe_tail points on last edge).
     *
     * Temporary pointers:
     * 1. p_edge2: Top horizontal edge list:      +-----------------------+ 2
     * 2. p_tail : end of middle edge list:       |\  |\  |\  |\  |\  |\  |
     *                                            |  \|  \|  \|  \|  \|  \|
     * 3. p_edge1: Bottom horizontal edge list:   +-----------------------+ 1
     *
     * pe_tail2  : end of list beginning at p_edge2
     * pe_temp   : position inside list beginning at p_edge1
     * p_edges   : head of the master edge list (part of our output)
     * p_vrtx1   : start of lower row of input vertices
     * p_vrtx2   : start of higher row of input vertices
     *
     * The routine generates two triangle            Lower      Upper 1  
     * upper one and lower one:                     | \           ----   
     * (Nums. are edges order in polys)            0|   \1       0\   |2 
     * The polygons are stored in the polygon        ----           \ |  
     * list (*p_polys) (pp_tail points on             2                  
     * last polygon).
     *                                                        1
     *                                                   -----------
     * In addition, the edge lists are updated -        | \   0     |
     * each edge has two pointers on the two            |   \       |
     * (one active if boundary) polygons which         0|1   0\1   0|1
     * uses it. These two pointer to polygons           |       \   |
     * are named: poly[0], poly[1]. The diagram         |    1    \ |
     * on the right show how they are used for the       -----------
     * upper and lower polygons (INNER_MESH polygons only).  0
     */

    for (i = 1; i < num_isolines; i++) {
        /* Read next column and gen. polys. */
        iso_lines = iso_lines->next;

        p_vrtx2 = iso_lines->points;    /* next row of vertices */
        p_edge2 = pe_tail2 = NULL;      /* clear top horizontal list */
        pe_temp = p_edge1;      /* pointer in bottom list */

        /*
         * Generate edges and triagles for next row:
         */

        /* generate first vertical edge */
        edge2 = add_edge(p_vrtx1, p_vrtx2, p_edges, &pe_tail);

        for (j = 0; j < grid_x_max - 1; j++) {

            /* copy vertical edge for lower triangle */
            edge0 = edge2;

            if (pe_temp && pe_temp->vertex[0] == p_vrtx1 + j) {
                /* test lower edge */
                edge2 = pe_temp;
                pe_temp = pe_temp->next;
            } else {
                edge2 = NULL;   /* edge is undefined */
            }

            /* generate diagonal edge */
            edge1 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j, p_edges, &pe_tail);
            if (edge1)
                edge1->position = DIAGONAL;

            /* generate lower triangle */
            lower_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);

            /* copy diagonal edge for upper triangle */
            edge0 = edge1;

            /* generate upper edge */
            edge1 = add_edge(p_vrtx2 + j, p_vrtx2 + j + 1, &p_edge2, &pe_tail2);

            /* generate vertical edge */
            edge2 = add_edge(p_vrtx1 + j + 1, p_vrtx2 + j + 1, p_edges, &pe_tail);

            /* generate upper triangle */
            upper_tri = add_poly(edge0, edge1, edge2, p_polys, &pp_tail);
        }

        if (p_edge2) {
            /* HBB 19991130 bugfix: if p_edge2 list is empty,
             * don't change p_edges list! Crashes by access
             * to NULL pointer pe_tail, the second time through,
             * otherwise */
            if ((*p_edges)) {   /* Chain new edges to main list. */
                pe_tail->next = p_edge2;
                pe_tail = pe_tail2;
            } else {
                (*p_edges) = p_edge2;
                pe_tail = pe_tail2;
            }
        }

        /* this row finished, move list heads up one row: */
        p_edge1 = p_edge2;
        p_vrtx1 = p_vrtx2;
    }

    /* Update the boundary flag, saved in each edge, and update indexes: */

    pe_temp = (*p_edges);

    while (pe_temp) {
        if ((!(pe_temp->poly[0])) || (!(pe_temp->poly[1])))
            (pe_temp->position) = BOUNDARY;
        pe_temp = pe_temp->next;
    }
}

/*
 * Calculate minimum and maximum values
 */
static void
calc_min_max(num_isolines, iso_lines, xx_min, yy_min, zz_min, xx_max, yy_max, zz_max)
    int num_isolines;           /* number of iso-lines input */
    struct iso_curve *iso_lines; /* iso-lines input */
    double *xx_min, *yy_min, *zz_min, *xx_max, *yy_max, *zz_max; /* min/max values in/out */
{
    int i, j, grid_x_max;
    struct coordinate  *vertex;
    /*printf("<calc_min_max> :  iso_lines->p_count %ld \n",iso_lines->p_count);*/
    grid_x_max = iso_lines->p_count;    /* number of vertices per iso_line */

    (*xx_min) = (*yy_min) = (*zz_min) = VERYLARGE;      /* clear min/max values */
    (*xx_max) = (*yy_max) = (*zz_max) = -VERYLARGE;
    /*printf(" <calc_min_max> %d \n",num_isolines);*/
    for (j = 0; j < num_isolines; j++) {
      /*printf(" <calc_min_max> iso_lines %lx %lx %d min %g max %g\n",      */
      /*       iso_lines,iso_lines->points,grid_x_max,(*zz_min),(*zz_max)); */
        vertex = iso_lines->points;

        for (i = 0; i < grid_x_max; i++) {
            if (vertex[i].type != UNDEFINED) {
                if (vertex[i].x > (*xx_max))
                    (*xx_max) = vertex[i].x;
                if (vertex[i].y > (*yy_max))
                    (*yy_max) = vertex[i].y;
                if (vertex[i].z > (*zz_max))
                    (*zz_max) = vertex[i].z;
                if (vertex[i].x < (*xx_min))
                    (*xx_min) = vertex[i].x;
                if (vertex[i].y < (*yy_min))
                    (*yy_min) = vertex[i].y;
                if (vertex[i].z < (*zz_min))
                    (*zz_min) = vertex[i].z;
                
            }
        }
        iso_lines = iso_lines->next;
        /*printf(" End of loop calc_min_max %d \n", j);*/
    }
    /* HBB 20000426: this code didn't take into account that axes might
     * be logscaled... */
#if 0
    /* HBB 20001220: DON'T. The values are actually already stored
     * logarithmized, as should be! */
    axis_unlog_interval(FIRST_X_AXIS, xx_min, xx_max, 0);
    axis_unlog_interval(FIRST_Y_AXIS, yy_min, yy_max, 0);
    axis_unlog_interval(FIRST_Z_AXIS, zz_min, zz_max, 0);
#endif

    /* 
     * fprintf(stderr," x: %g, %g\n", (*xx_min), (*xx_max));
     * fprintf(stderr," y: %g, %g\n", (*yy_min), (*yy_max));
     * fprintf(stderr," z: %g, %g\n", (*zz_min), (*zz_max));
     */
}

/*
 * Generate new edge and append it to list, but only if both vertices are 
 * defined. The list is referenced by p_edge and pe_tail (p_edge points on 
 * first edge and pe_tail on last one).
 * Note, the list may be empty (pe_edge==pe_tail==NULL) on entry and exit.
 */
static struct edge_struct *
add_edge(point0, point1, p_edge, pe_tail)
struct coordinate  *point0;     /* 2 vertices input */
struct coordinate  *point1;
struct edge_struct **p_edge, **pe_tail;         /* pointers to edge list in/out */
{
    struct edge_struct *pe_temp = NULL;

#if 1
    if (point0->type == INRANGE && point1->type == INRANGE) {
#else
    if (point0->type != UNDEFINED && point1->type != UNDEFINED) {
#endif

        pe_temp = (struct edge_struct *)
            gp_alloc(sizeof(struct edge_struct), "contour edge");

        pe_temp->poly[0] = NULL;        /* clear links           */
        pe_temp->poly[1] = NULL;
        pe_temp->vertex[0] = point0;    /* First vertex of edge. */
        pe_temp->vertex[1] = point1;    /* Second vertex of edge. */
        pe_temp->next = NULL;
        pe_temp->position = INNER_MESH;         /* default position in mesh */

        if ((*pe_tail)) {
            (*pe_tail)->next = pe_temp;         /* Stick new record as last one. */
        } else {
            (*p_edge) = pe_temp;        /* start new list if empty */
        }
        (*pe_tail) = pe_temp;   /* continue to last record. */

    }
    return pe_temp;             /* returns NULL, if no edge allocated */
}

/*
 * Generate new triangle and append it to list, but only if all edges are defined.
 * The list is referenced by p_poly and pp_tail (p_poly points on first ploygon
 * and pp_tail on last one).
 * Note, the list may be empty (pe_ploy==pp_tail==NULL) on entry and exit.
 */
static struct poly_struct *
add_poly(edge0, edge1, edge2, p_poly, pp_tail)
    struct edge_struct *edge0, *edge1, *edge2;  /* 3 edges input */
    struct poly_struct **p_poly, **pp_tail;             /* pointers to polygon list in/out */
{
    struct poly_struct *pp_temp = NULL;

    if (edge0 && edge1 && edge2) {

        pp_temp = (struct poly_struct *)
            gp_alloc(sizeof(struct poly_struct), "contour polygon");

        pp_temp->edge[0] = edge0;       /* First edge of triangle */
        pp_temp->edge[1] = edge1;       /* Second one             */
        pp_temp->edge[2] = edge2;       /* Third one              */
        pp_temp->next = NULL;

        if (edge0->poly[0])     /* update edge0 */
            edge0->poly[1] = pp_temp;
        else
            edge0->poly[0] = pp_temp;

        if (edge1->poly[0])     /* update edge1 */
            edge1->poly[1] = pp_temp;
        else
            edge1->poly[0] = pp_temp;

        if (edge2->poly[0])     /* update edge2 */
            edge2->poly[1] = pp_temp;
        else
            edge2->poly[0] = pp_temp;

        if ((*pp_tail))         /* Stick new record as last one. */
            (*pp_tail)->next = pp_temp;
        else
            (*p_poly) = pp_temp;        /* start new list if empty */

        (*pp_tail) = pp_temp;   /* continue to last record. */

    }
    return pp_temp;             /* returns NULL, if no edge allocated */
}



/*
 * Calls the (hopefully) desired interpolation/approximation routine.
 */
static void
put_contour(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_isclosed)
    struct cntr_struct *p_cntr; /* contour structure input */
    double z_level;             /* Z level of contour input */
    double xx_min, xx_max, yy_min, yy_max; /* minimum/maximum values input */
    TBOOLEAN contr_isclosed;            /* contour line closed? (input) */
{

    if (!p_cntr)
        return;                 /* Nothing to do if it is empty contour. */

    switch (interp_kind) {
    case CONTOUR_KIND_LINEAR:   /* No interpolation/approximation. */
        put_contour_nothing(p_cntr);
        break;
    case CONTOUR_KIND_CUBIC_SPL: /* Cubic spline interpolation. */
        put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max,
                          chk_contour_kind(p_cntr, contr_isclosed));

        break;
    case CONTOUR_KIND_BSPLINE:  /* Bspline approximation. */
        put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max,
                            chk_contour_kind(p_cntr, contr_isclosed));
        break;
    }
    free_contour(p_cntr);
}

/*
 * Simply puts contour coordinates in order with no interpolation or
 * approximation.
 */
static void
put_contour_nothing(p_cntr)
struct cntr_struct *p_cntr;
{
    while (p_cntr) {
        add_cntr_point(p_cntr->X, p_cntr->Y);
        p_cntr = p_cntr->next;
    }
    end_crnt_cntr();
}

/*
 * for some reason contours are never flagged as 'isclosed'
 * if first point == last point, set flag accordingly
 *
 */

static int
chk_contour_kind(p_cntr, contr_isclosed)
    struct cntr_struct *p_cntr;
    TBOOLEAN contr_isclosed;
{
    struct cntr_struct *pc_tail = NULL;
    TBOOLEAN current_contr_isclosed;

    /*fprintf(stderr, "check_contour_kind: current contour_kind value is %d\n", contour_kind);*/

    current_contr_isclosed = contr_isclosed;

    if (! contr_isclosed) {
        pc_tail = p_cntr;
        while (pc_tail->next)
            pc_tail = pc_tail->next;    /* Find last point. */

        /* test if first and last point are equal */
        if (fuzzy_equal(pc_tail, p_cntr)) {
            current_contr_isclosed = TRUE;
            fprintf(stderr, "check_contour_kind: contr_isclosed changed to %d\n", current_contr_isclosed);
        }
    }
    return (current_contr_isclosed);
}

/*
 * Generate a cubic spline curve through the points (x_i,y_i) which are
 * stored in the linked list p_cntr.
 * The spline is defined as a 2d-function s(t) = (x(t),y(t)), where the
 * parameter t is the length of the linear stroke.
 */
static void
put_contour_cubic(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_isclosed)
    struct cntr_struct *p_cntr;
    double z_level, xx_min, xx_max, yy_min, yy_max;
    TBOOLEAN contr_isclosed;
{
    int num_pts, num_intpol;
    double unit_x, unit_y;      /* To define norm (x,y)-plane */
    double *delta_t;            /* Interval length t_{i+1}-t_i */
    double *d2x, *d2y;          /* Second derivatives x''(t_i), y''(t_i) */
    struct cntr_struct *pc_tail;

    num_pts = count_contour(p_cntr);    /* Number of points in contour. */

    pc_tail = p_cntr;           /* Find last point. */
    while (pc_tail->next)
        pc_tail = pc_tail->next;

    if (contr_isclosed) {
        /* Test if first and last point are equal (should be) */
        if (!fuzzy_equal(pc_tail, p_cntr)) {
            pc_tail->next = p_cntr;     /* Close contour list - make it circular. */
            num_pts++;
        }
    }
    delta_t = (double *) gp_alloc(num_pts * sizeof(double), "contour delta_t");
    d2x = (double *) gp_alloc(num_pts * sizeof(double), "contour d2x");
    d2y = (double *) gp_alloc(num_pts * sizeof(double), "contour d2y");

    /* Width and height of the grid is used as a unit length (2d-norm) */
    unit_x = xx_max - x_min;
    unit_y = yy_max - y_min;
    /* FIXME HBB 20010121: 'zero' should not be used as an absolute
     * figure to compare to data */
    unit_x = (unit_x > zero ? unit_x : zero);   /* should not be zero */
    unit_y = (unit_y > zero ? unit_y : zero);

    if (num_pts > 2) {
        /*
         * Calculate second derivatives d2x[], d2y[] and interval lengths delta_t[]:
         */
        if (!gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t,
                              contr_isclosed, unit_x, unit_y)) {
            free((char *) delta_t);
            free((char *) d2x);
            free((char *) d2y);
            if (contr_isclosed)
                pc_tail->next = NULL;   /* Un-circular list */
            return;
        }
    }
    /* If following (num_pts > 1) is TRUE then exactly 2 points in contour.  */
    else if (num_pts > 1) {
        /* set all second derivatives to zero, interval length to 1 */
        d2x[0] = 0.;
        d2y[0] = 0.;
        d2x[1] = 0.;
        d2y[1] = 0.;
        delta_t[0] = 1.;
    } else {                    /* Only one point ( ?? ) - ignore it. */
        free((char *) delta_t);
        free((char *) d2x);
        free((char *) d2y);
        if (contr_isclosed)
            pc_tail->next = NULL;       /* Un-circular list */
        return;
    }

    /* Calculate "num_intpol" interpolated values */
    num_intpol = 1 + (num_pts - 1) * contour_pts;       /* global: contour_pts */
    intp_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, num_intpol);

    free((char *) delta_t);
    free((char *) d2x);
    free((char *) d2y);

    if (contr_isclosed)
        pc_tail->next = NULL;   /* Un-circular list */

    end_crnt_cntr();
}


/*
 * Find Bspline approximation for this data set.
 * Uses global variable contour_pts to determine number of samples per
 * interval, where the knot vector intervals are assumed to be uniform, and
 * global variable contour_order for the order of Bspline to use.
 */
static void
put_contour_bspline(p_cntr, z_level, xx_min, xx_max, yy_min, yy_max, contr_isclosed)
struct cntr_struct *p_cntr;
double z_level, xx_min, xx_max, yy_min, yy_max;
TBOOLEAN contr_isclosed;
{
    int num_pts;
    int order = contour_order - 1;

    num_pts = count_contour(p_cntr);    /* Number of points in contour. */
    if (num_pts < 2)
        return;                 /* Can't do nothing if empty or one points! */
    /* Order must be less than number of points in curve - fix it if needed. */
    if (order > num_pts - 1)
        order = num_pts - 1;

    gen_bspline_approx(p_cntr, num_pts, order, contr_isclosed);
    end_crnt_cntr();
}

/*
 * Free all elements in the contour list.
 */
static void
free_contour(p_cntr)
struct cntr_struct *p_cntr;
{
    struct cntr_struct *pc_temp;

    while (p_cntr) {
        pc_temp = p_cntr;
        p_cntr = p_cntr->next;
        free((char *) pc_temp);
    }
}

/*
 * Counts number of points in contour.
 */
static int
count_contour(p_cntr)
struct cntr_struct *p_cntr;
{
    int count = 0;

    while (p_cntr) {
        count++;
        p_cntr = p_cntr->next;
    }
    return count;
}

/*
 * Find second derivatives (x''(t_i),y''(t_i)) of cubic spline interpolation
 * through list of points (x_i,y_i). The parameter t is calculated as the
 * length of the linear stroke. The number of points must be at least 3.
 * Note: For closed contours the first and last point must be equal.
 */
static int
gen_cubic_spline(num_pts, p_cntr, d2x, d2y, delta_t, contr_isclosed, unit_x, unit_y)
int num_pts;                    /* Number of points (num_pts>=3), input */
struct cntr_struct *p_cntr;     /* List of points (x(t_i),y(t_i)), input */
double d2x[], d2y[],            /* Second derivatives (x''(t_i),y''(t_i)), output */
 delta_t[];                     /* List of interval lengths t_{i+1}-t_{i}, output */
TBOOLEAN contr_isclosed;        /* Closed or open contour?, input  */
double unit_x, unit_y;          /* Unit length in x and y (norm=1), input */
{
    int n, i;
    double norm;
    tri_diag *m;                /* The tri-diagonal matrix is saved here. */
    struct cntr_struct *pc_temp;

    m = (tri_diag *) gp_alloc(num_pts * sizeof(tri_diag), "contour tridiag m");

    /*
     * Calculate first differences in (d2x[i], d2y[i]) and interval lengths
     * in delta_t[i]:
     */
    pc_temp = p_cntr;
    for (i = 0; i < num_pts - 1; i++) {
        d2x[i] = pc_temp->next->X - pc_temp->X;
        d2y[i] = pc_temp->next->Y - pc_temp->Y;
        /*
         * The norm of a linear stroke is calculated in "normal coordinates"
         * and used as interval length:
         */
        delta_t[i] = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y));

        d2x[i] /= delta_t[i];   /* first difference, with unit norm: */
        d2y[i] /= delta_t[i];   /*   || (d2x[i], d2y[i]) || = 1      */

        pc_temp = pc_temp->next;
    }

    /*
     * Setup linear system:  m * x = b
     */
    n = num_pts - 2;            /* Without first and last point */
    if (contr_isclosed) {
        /* First and last points must be equal for closed contours */
        delta_t[num_pts - 1] = delta_t[0];
        d2x[num_pts - 1] = d2x[0];
        d2y[num_pts - 1] = d2y[0];
        n++;                    /* Add last point (= first point) */
    }
    for (i = 0; i < n; i++) {
        /* Matrix M, mainly tridiagonal with cyclic second index ("j = j+n mod n") */
        m[i][0] = delta_t[i];   /* Off-diagonal element M_{i,i-1} */
        m[i][1] = 2. * (delta_t[i] + delta_t[i + 1]);   /* M_{i,i} */
        m[i][2] = delta_t[i + 1];       /* Off-diagonal element M_{i,i+1} */

        /* Right side b_x and b_y */
        d2x[i] = (d2x[i + 1] - d2x[i]) * 6.;
        d2y[i] = (d2y[i + 1] - d2y[i]) * 6.;

        /*
         * If the linear stroke shows a cusps of more than 90 degree, the right
         * side is reduced to avoid oscillations in the spline:
         */
        norm = sqrt(SQR(d2x[i] / unit_x) + SQR(d2y[i] / unit_y)) / 8.5;

        if (norm > 1.) {
            d2x[i] /= norm;
            d2y[i] /= norm;
            /* The first derivative will not be continuous */
        }
    }

    if (!contr_isclosed) {
        /* Third derivative is set to zero at both ends */
        m[0][1] += m[0][0];     /* M_{0,0}     */
        m[0][0] = 0.;           /* M_{0,n-1}   */
        m[n - 1][1] += m[n - 1][2];     /* M_{n-1,n-1} */
        m[n - 1][2] = 0.;       /* M_{n-1,0}   */
    }
    /* Solve linear systems for d2x[] and d2y[] */


    if (solve_cubic_1(m, n)) {  /* Calculate Cholesky decomposition */
        solve_cubic_2(m, d2x, n);       /* solve M * d2x = b_x */
        solve_cubic_2(m, d2y, n);       /* solve M * d2y = b_y */

    } else {                    /* Should not happen, but who knows ... */
        free((char *) m);
        return FALSE;
    }

    /* Shift all second derivatives one place right and abdate end points */
    for (i = n; i > 0; i--) {
        d2x[i] = d2x[i - 1];
        d2y[i] = d2y[i - 1];
    }
    if (contr_isclosed) {
        d2x[0] = d2x[n];
        d2y[0] = d2y[n];
    } else {
        d2x[0] = d2x[1];        /* Third derivative is zero in */
        d2y[0] = d2y[1];        /*     first and last interval */
        d2x[n + 1] = d2x[n];
        d2y[n + 1] = d2y[n];
    }

    free((char *) m);
    return TRUE;
}

/*
 * Calculate interpolated values of the spline function (defined via p_cntr
 * and the second derivatives d2x[] and d2y[]). The number of tabulated
 * values is n. On an equidistant grid n_intpol values are calculated.
 */
static void
intp_cubic_spline(n, p_cntr, d2x, d2y, delta_t, n_intpol)
int n;
struct cntr_struct *p_cntr;
double d2x[], d2y[], delta_t[];
int n_intpol;
{
    double t, t_skip, t_max;
    double x0, x1, x, y0, y1, y;
    double d, hx, dx0, dx01, hy, dy0, dy01;
    int i;

    /* The length of the total interval */
    t_max = 0.;
    for (i = 0; i < n - 1; i++)
        t_max += delta_t[i];

    /* The distance between interpolated points */
    t_skip = (1. - 1e-7) * t_max / (n_intpol - 1);

    t = 0.;                     /* Parameter value */
    x1 = p_cntr->X;
    y1 = p_cntr->Y;
    add_cntr_point(x1, y1);     /* First point. */
    t += t_skip;

    for (i = 0; i < n - 1; i++) {
        p_cntr = p_cntr->next;

        d = delta_t[i];         /* Interval length */
        x0 = x1;
        y0 = y1;
        x1 = p_cntr->X;
        y1 = p_cntr->Y;
        hx = (x1 - x0) / d;
        hy = (y1 - y0) / d;
        dx0 = (d2x[i + 1] + 2 * d2x[i]) / 6.;
        dy0 = (d2y[i + 1] + 2 * d2y[i]) / 6.;
        dx01 = (d2x[i + 1] - d2x[i]) / (6. * d);
        dy01 = (d2y[i + 1] - d2y[i]) / (6. * d);
        while (t <= delta_t[i]) {       /* t in current interval ? */
            x = x0 + t * (hx + (t - d) * (dx0 + t * dx01));
            y = y0 + t * (hy + (t - d) * (dy0 + t * dy01));
            add_cntr_point(x, y);       /* next point. */
            t += t_skip;
        }
        t -= delta_t[i];        /* Parameter t relative to start of next interval */
    }
}

/*
 * The following two procedures solve the special linear system which arise
 * in cubic spline interpolation. If x is assumed cyclic ( x[i]=x[n+i] ) the
 * equations can be written as (i=0,1,...,n-1):
 *     m[i][0] * x[i-1] + m[i][1] * x[i] + m[i][2] * x[i+1] = b[i] .
 * In matrix notation one gets M * x = b, where the matrix M is tridiagonal
 * with additional elements in the upper right and lower left position:
 *   m[i][0] = M_{i,i-1}  for i=1,2,...,n-1    and    m[0][0] = M_{0,n-1} ,
 *   m[i][1] = M_{i, i }  for i=0,1,...,n-1
 *   m[i][2] = M_{i,i+1}  for i=0,1,...,n-2    and    m[n-1][2] = M_{n-1,0}.
 * M should be symmetric (m[i+1][0]=m[i][2]) and positiv definite.
 * The size of the system is given in n (n>=1).
 *
 * In the first procedure the Cholesky decomposition M = C^T * D * C
 * (C is upper triangle with unit diagonal, D is diagonal) is calculated.
 * Return TRUE if decomposition exist.
 */
static int
solve_cubic_1(m, n)
tri_diag m[];
int n;
{
    int i;
    double m_ij, m_n, m_nn, d;

    if (n < 1)
        return FALSE;           /* Dimension should be at least 1 */

    d = m[0][1];                /* D_{0,0} = M_{0,0} */
    if (d <= 0.)
        return FALSE;           /* M (or D) should be positiv definite */
    m_n = m[0][0];              /*  M_{0,n-1}  */
    m_nn = m[n - 1][1];         /* M_{n-1,n-1} */
    for (i = 0; i < n - 2; i++) {
        m_ij = m[i][2];         /*  M_{i,1}  */
        m[i][2] = m_ij / d;     /* C_{i,i+1} */
        m[i][0] = m_n / d;      /* C_{i,n-1} */
        m_nn -= m[i][0] * m_n;  /* to get C_{n-1,n-1} */
        m_n = -m[i][2] * m_n;   /* to get C_{i+1,n-1} */
        d = m[i + 1][1] - m[i][2] * m_ij;       /* D_{i+1,i+1} */
        if (d <= 0.)
            return FALSE;       /* Elements of D should be positiv */
        m[i + 1][1] = d;
    }
    if (n >= 2) {               /* Complete last column */
        m_n += m[n - 2][2];     /* add M_{n-2,n-1} */
        m[n - 2][0] = m_n / d;  /* C_{n-2,n-1} */
        m[n - 1][1] = d = m_nn - m[n - 2][0] * m_n;     /* D_{n-1,n-1} */
        if (d <= 0.)
            return FALSE;
    }
    return TRUE;
}

/*
 * The second procedure solves the linear system, with the Choleky
 * decomposition calculated above (in m[][]) and the right side b given
 * in x[]. The solution x overwrites the right side in x[].
 */
static void
solve_cubic_2(m, x, n)
tri_diag m[];
double x[];
int n;
{
    int i;
    double x_n;

    /* Division by transpose of C : b = C^{-T} * b */
    x_n = x[n - 1];
    for (i = 0; i < n - 2; i++) {
        x[i + 1] -= m[i][2] * x[i];     /* C_{i,i+1} * x_{i} */
        x_n -= m[i][0] * x[i];  /* C_{i,n-1} * x_{i} */
    }
    if (n >= 2)
        x[n - 1] = x_n - m[n - 2][0] * x[n - 2];        /* C_{n-2,n-1} * x_{n-1} */

    /* Division by D: b = D^{-1} * b */
    for (i = 0; i < n; i++)
        x[i] /= m[i][1];

    /* Division by C: b = C^{-1} * b */
    x_n = x[n - 1];
    if (n >= 2)
        x[n - 2] -= m[n - 2][0] * x_n;  /* C_{n-2,n-1} * x_{n-1} */
    for (i = n - 3; i >= 0; i--) {
        /*      C_{i,i+1} * x_{i+1} + C_{i,n-1} * x_{n-1} */
        x[i] -= m[i][2] * x[i + 1] + m[i][0] * x_n;
    }
    return;
}

/*
 * Solve tri diagonal linear system equation. The tri diagonal matrix is
 * defined via matrix M, right side is r, and solution X i.e. M * X = R.
 * Size of system given in n. Return TRUE if solution exist.
 */
/* not used any more in "contour.c", but in "spline.c" (21. Dec. 1995) ! */

int
solve_tri_diag(m, r, x, n)
tri_diag m[];
double r[], x[];
int n;
{
    int i;
    double t;

    for (i = 1; i < n; i++) {   /* Eliminate element m[i][i-1] (lower diagonal). */
        if (m[i - 1][1] == 0)
            return FALSE;
        t = m[i][0] / m[i - 1][1];      /* Find ratio between the two lines. */
/*      m[i][0] = m[i][0] - m[i-1][1] * t; */
/* m[i][0] is not used any more (and set to 0 in the above line) */
        m[i][1] = m[i][1] - m[i - 1][2] * t;
        r[i] = r[i] - r[i - 1] * t;
    }
    /* Now do back subtitution - update the solution vector X: */
    if (m[n - 1][1] == 0)
        return FALSE;
    x[n - 1] = r[n - 1] / m[n - 1][1];  /* Find last element. */
    for (i = n - 2; i >= 0; i--) {
        if (m[i][1] == 0)
            return FALSE;
        x[i] = (r[i] - x[i + 1] * m[i][2]) / m[i][1];
    }
    return TRUE;
}

/*
 * Generate a Bspline curve defined by all the points given in linked list p:
 * Algorithm: using deBoor algorithm
 * Note: if Curvekind is open contour than Open end knot vector is assumed,
 *       else (closed contour) Float end knot vector is assumed.
 * It is assumed that num_of_points is at least 2, and order of Bspline is less
 * than num_of_points!
 */
static void
gen_bspline_approx(p_cntr, num_of_points, order, contr_isclosed)
struct cntr_struct *p_cntr;
int num_of_points, order;
TBOOLEAN contr_isclosed;
{
    int knot_index = 0, pts_count = 1;
    double dt, t, next_t, t_min, t_max, x, y;
    struct cntr_struct *pc_temp = p_cntr, *pc_tail = NULL;

    /* If the contour is Closed one we must update few things:
     * 1. Make the list temporary circular, so we can close the contour.
     * 2. Update num_of_points - increase it by "order-1" so contour will be
     *    closed. This will evaluate order more sections to close it!
     */
    if (contr_isclosed) {
        pc_tail = p_cntr;
        while (pc_tail->next)
            pc_tail = pc_tail->next;    /* Find last point. */

        /* test if first and last point are equal */
        if (fuzzy_equal(pc_tail, p_cntr)) {
            /* Close contour list - make it circular. */
            pc_tail->next = p_cntr->next;
            num_of_points += order - 1;
        } else {
            pc_tail->next = p_cntr;
            num_of_points += order;
        }
    }
    /* Find first (t_min) and last (t_max) t value to eval: */
    t = t_min = fetch_knot(contr_isclosed, num_of_points, order, order);
    t_max = fetch_knot(contr_isclosed, num_of_points, order, num_of_points);
    next_t = t_min + 1.0;
    knot_index = order;
    dt = 1.0 / contour_pts;     /* Number of points per one section. */


    while (t < t_max) {
        if (t > next_t) {
            pc_temp = pc_temp->next;    /* Next order ctrl. pt. to blend. */
            knot_index++;
            next_t += 1.0;
        }
        eval_bspline(t, pc_temp, num_of_points, order, knot_index,
                     contr_isclosed, &x, &y);   /* Next pt. */
        add_cntr_point(x, y);
        pts_count++;
        /* As we might have some real number round off problems we do      */
        /* the last point outside the loop                                 */
        if (pts_count == contour_pts * (num_of_points - order) + 1)
            break;
        t += dt;
    }

    /* Now do the last point */
    eval_bspline(t_max - EPSILON, pc_temp, num_of_points, order, knot_index,
                 contr_isclosed, &x, &y);
    add_cntr_point(x, y);       /* Complete the contour. */

    if (contr_isclosed) /* Update list - un-circular it. */
        pc_tail->next = NULL;
}

/*
 * The routine to evaluate the B-spline value at point t using knot vector
 * from function fetch_knot(), and the control points p_cntr.
 * Returns (x, y) of approximated B-spline. Note that p_cntr points on the
 * first control point to blend with. The B-spline is of order order.
 */
static void
eval_bspline(t, p_cntr, num_of_points, order, j, contr_isclosed, x, y)
double t;
struct cntr_struct *p_cntr;
int num_of_points, order, j;
TBOOLEAN contr_isclosed;
double *x, *y;
{
    int i, p;
    double ti, tikp, *dx, *dy;  /* Copy p_cntr into it to make it faster. */

    dx = (double *) gp_alloc((order + j) * sizeof(double), "contour b_spline");
    dy = (double *) gp_alloc((order + j) * sizeof(double), "contour b_spline");

    /* Set the dx/dy - [0] iteration step, control points (p==0 iterat.): */
    for (i = j - order; i <= j; i++) {
        dx[i] = p_cntr->X;
        dy[i] = p_cntr->Y;
        p_cntr = p_cntr->next;
    }

    for (p = 1; p <= order; p++) {      /* Iteration (b-spline level) counter. */
        for (i = j; i >= j - order + p; i--) {  /* Control points indexing. */
            ti = fetch_knot(contr_isclosed, num_of_points, order, i);
            tikp = fetch_knot(contr_isclosed, num_of_points, order, i + order + 1 - p);
            if (ti == tikp) {   /* Should not be a problems but how knows... */
            } else {
                dx[i] = dx[i] * (t - ti) / (tikp - ti) +        /* Calculate x. */
                    dx[i - 1] * (tikp - t) / (tikp - ti);
                dy[i] = dy[i] * (t - ti) / (tikp - ti) +        /* Calculate y. */
                    dy[i - 1] * (tikp - t) / (tikp - ti);
            }
        }
    }
    *x = dx[j];
    *y = dy[j];
    free((char *) dx);
    free((char *) dy);
}

/*
 * Routine to get the i knot from uniform knot vector. The knot vector
 * might be float (Knot(i) = i) or open (where the first and last "order"
 * knots are equal). contr_isclosed determines knot kind - open contour means
 * open knot vector, and closed contour selects float knot vector.
 * Note the knot vector is not exist and this routine simulates it existance
 * Also note the indexes for the knot vector starts from 0.
 */
static double
fetch_knot(contr_isclosed, num_of_points, order, i)
    TBOOLEAN contr_isclosed;
    int num_of_points, order, i;
{
    if(! contr_isclosed) {
        if (i <= order)
            return 0.0;
        else if (i <= num_of_points)
            return (double) (i - order);
        else
            return (double) (num_of_points - order);
    } else {
        return (double) i;
    }
}