code for isoline(C)

Well, Thanks to all you folks who replied, I now have a contour plotting
routine, that does its job fairly well. I've posted this code below. It 
is not nearly runnable as other Unix programs go, but you could use it
if you need C code to start off with to do contour plotting.

It doesn't handle stuff like labelling etc. which is somewhat tricky.

Basically almost all of the respondents suggested dividing up the set of 
discrete points into a set of triangles and then computing the intersection
of these rectangles with planes representing the contour levels. This way
we can process all the levels at once instead of the "crawler" type algorithms
which is what I had originally been trying to get to work.

To summarize their ideas: if P1=P[j],P2=P[i+1][j],P3=P[i+1][j+1]
and P4=P[j+1] are four points, compute an intermediate point using
an interpolation scheme (in the code below, I just use averaging, but
one could obviously use fancier schemes). If this point is P, then the
foll. four triangles are used to compute the intersections with
contour levels l_0, .. l_n: P1 P2 P, P2 P3 P, P3 P4 P, P4 P1 P.

The intersection computation is quite trivial. Once these line segments have been
computed you could smooth them in 2d (the code below doesn't).

(see Michael Aramini's earlier posting for a variation on this that uses 
rectangles directly -- He also mentions that you could use the diagonal edges
between two vertices instead of computing another intermediate point, but this
method has problems).

-Sundar
-----Sample Code Follows----------
/*
* If you want try using this, you probably would need to provide
* graph_set_point(), and graph_draw_line(), besides do things like
* initialize window systems and so on.
*
* Please see the code for a note about data format
*/

/* Definitions */
struct _array {
float *a_elements;
int a_rows;
int a_cols;
float *a_filler;
};

#define ELTS(a) (a->a_elements)
#define NOROWS(a) (a->a_rows)
#define NOCOLS(a) (a->a_cols)
#define NOELTS(a) (NOROWS(a) * NOCOLS(a))

#define COLUMN_VECTOR(a) (NOCOLS((a)) == 1)
#define ROW_VECTOR(a) (NOROWS((a)) == 1)

typedef struct _array *array;

typedef struct _graph *graph;
typedef struct _graph_port *graph_port;

struct _graph_port {
char *p_window; /* window system dependent private */
double p_xmin; /* object space min x-coordinate */
double p_xmax; /* object space max x-coordinate */
double p_ymin; /* object space min y-coordinate */
double p_ymax; /* object space max y-coordinate */
double p_xscale; /* this is just x-max - x-min */
double p_yscale; /* this is just y-max - y-min */
int p_xorigin; /* image space x-origin co-ordinate */
int p_yorigin; /* image space y-origin */
int p_width; /* image space width */
int p_height; /* image space height */
int p_realwidth; /* image space width */
int p_realheight; /* image space height */
int p_curx; /* where current point is in image */
int p_cury; /* y-coordinate of above */
int p_line_offset; /* line offset for labels */
int p_autoscale; /* in auto scaling on? */
int p_include_xzero; /* include x zero axis */
int p_include_yzero; /* include y zero axis */
int p_titlefont;
int p_ticfont;
int p_labelfont;
int p_margin;
int p_ticstyle; /* styles of tics (large, small or none) */
};

struct _graph {
int g_no; /* number associated with this graph */
int g_points; /* total number of points in graph */
array g_x; /* graph x points */
array g_y; /* graph y points */
array g_z; /* graph z points */
char *g_line_label; /* graph line label */
char *g_xlabel; /* graph x axis label */
char *g_ylabel; /* graph y axis label */
char *g_title; /* graph title */
int g_color; /* graph color or pattern */
int g_type; /* type of graph */
};


/*
* Contour Plotting algorithm!
*/

/*
* This routine does the level comparison. 
* It takes three points, a level, and a color.
* It intersects the triangle with a plane at the given level,
* and draws the line segments that denote these intersections.
*/

#define swap_val(x, y, temp) {temp = x; x = y; y = x; }

static
graph_segment(gp, x1, y1, z1, x2, y2, z2, x3, y3, z3, level, color)
graph_port gp;
double x1, y1, z1, x2, y2, z2, x3, y3, z3, level;
int color;
{
double x, y, factor, temp;
if(z1 >= z2) {
swap_val(x1, x2, temp);
swap_val(y1, y2, temp);
swap_val(z1, z2, temp);
}
if(z2 > z3) {
if(z3 < z1) {
swap_val(x1, x3, temp); 
swap_val(y1, y3, temp);
swap_val(z1, z3, temp);
}
swap_val(x2, x3, temp); 
swap_val(y2, y3, temp);
swap_val(z2, z3, temp);
}
/* z1 <= z2 <= z3 */
if(level < z1 || level > z3) {
return;
}
if((level == z1) && (z1 == z2) && (z2 != z3)) {
/* draw a line from x1, y1 to x2, y2 */
graph_set_point(gp, x1, y1, color);
graph_lineto(gp, x2, y2, color);
return 1;
}
if((level == z3) && (z3 == z2) && (z1 != z2)) {
/* draw a line from x2, y2 to x3, y3 */
graph_set_point(gp, x2, y2, color);
graph_lineto(gp, x3, y3, color);
return 1;
}

if(level == z2) {
/* we know that z2 is definitely between z1 and z3 
draw a line from P2 to point on segment between P1 and P3
*/
if(z3 == z1) return 0;
factor = (level - z1)/(z3-z1);
x = x1 + (x3-x1) * factor;
y = y1 + (y3-y1) * factor;
graph_set_point(gp, x2, y2, color);
graph_lineto(gp, x, y, color);
return 1;
}
if(level < z2) {
/* P2 and P3 are above, P1 below */
if(z2 == z1 || z3 == z1) return 0;
factor = (level - z1) / (z2 - z1);
x = x1 + (x2 - x1) * factor;
y = y1 + (y2 - y1) * factor;
graph_set_point(gp, x, y, color);
factor = (level - z1) / (z3 - z1); 
x = x1 + (x3 - x1) * factor;
y = y1 + (y3 - y1) * factor;
graph_lineto(gp, x, y, color);
return 1;
} else {
/* P1 and P2 are below, P3 above */
if(z3 == z1 || z3 == z2) return 0;
factor = (level - z1) / (z3 - z1);
x = x1 + (x3 - x1) * factor;
y = y1 + (y3 - y1) * factor;
graph_set_point(gp, x, y, color);
factor = (level - z2) / (z3 - z2);
x = x2 + (x3 - x2) * factor;
y = y2 + (y3 - y2) * factor;
graph_lineto(gp, x, y, color);
return 1;
}
return 0;
}

/* This is the main function */
/* Data formats;
g->g_x -> contains N x values.
g->g_y -> contains N y values.
g->g_z -> contains N^2 z values, which are interpreted as a 
matrix of z values, as
z[0][0], z[0][1] .... z[0][N],
z[1][0], ... upto z[N][N]. 
(horizontal sweep from X_min to X_max, inner loop
going from Y_min to Y_max);

This way we only store N^2+2N values instead of 3N^2 values.
*/

graph_plot_contour(gp, g, color)
graph_port gp;
graph g;
int color;
{
register int i, j;
float *x, *y, *z;
float minz, maxz;
double xavg, yavg, zavg;
double level, diff;
double zij, zipj, zijp, zipjp;

x = g->g_x->a_elements;
y = g->g_y->a_elements;
z = g->g_z->a_elements;

/* 
* We find the min. and max values of the Z heights and divide that
* by a number of levels that is pre-determined. 
* We should really do something more reasonable here, similar to the
* stuff we do for figuring out tick marks on 2-d plots 
*/
array_minmax(g->g_z, &minz, &maxz);

diff = (maxz - minz)/10.0;

for(i=0;i<NOROWS(g->g_x)-1;i++) {
xavg = (x + x[i+1]) / 2.0;

for(j=0;j<NOROWS(g->g_y)-1;j++) {
yavg = (y[j] + y[j+1]) / 2.0;

zij = z[(i*NOROWS(g->g_y))+j];
zipj = z[((i+1)*NOROWS(g->g_y))+j];
zijp = z[(i*NOROWS(g->g_y))+j+1];
zipjp = z[((i+1)*NOROWS(g->g_y))+j+1];
zavg = (zij + zipj + zijp + zipjp)/4.0;

/* For each level */
for(level = minz; level < maxz; level += diff) {
/* Triangle i,j avg and i,j+1 */
graph_segment(gp, x, y[j], zij, x, y[j+1], zijp, 
xavg, yavg, zavg, level, color);
graph_segment(gp, x, y[j], zij, x[i+1], y[j], zipj, 
xavg, yavg, zavg, level, color);
graph_segment(gp, x[i+1], y[j], zipj, x[i+1], y[j+1], zipjp, 
xavg, yavg, zavg, level, color);
graph_segment(gp, x, y[j+1], zijp, x[i+1], y[j+1], zipjp, 
xavg, yavg, zavg, level, color);
}
}
}
}