From sundar@wheaties.ai.mit.edu Wed Mar 8 10:29:02 1989 Path: leah!csd4.milw.wisc.edu!mailrus!tut.cis.ohio-state.edu!ucbvax!bloom-beacon!ai-lab!sundar From: sundar@wheaties.ai.mit.edu (Sundar Narasimhan) Newsgroups: comp.graphics Subject: contour plotting algorithms (Summary of responses + code) Keywords: contour plotting, contours Message-ID: <789@special-k.ai.mit.edu> Date: 8 Mar 89 15:29:02 GMT Distribution: na Organization: MIT AI Lab, Cambridge, MA Lines: 276 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[i][j],P2=P[i+1][j],P3=P[i+1][j+1] and P4=P[i][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). References: A Decomposable Algorithm for Contour Surface Display Generation Michael, J. Zyda in ACM Transactions on Graphics, Vol 7, No 2, April 1988, pp 129-148. Contouring over rectangular and skewed rectangular grids - An introduction D. C. Sutcliffe in Mathematical Methods in Computer Graphics and Design edited by K. W. Brodlie, Academic Press, New York, 1980. Macsyma contour algorithm: Computer Journal, Vol 15 No 4, pp 382, 1972 Macsyma 3d hidden line: CACM algorithm 420, Vol 15. 1972 A Contouring Subroutine, Paul D. Bourke, June 1987, pp 143-150 Michael Aramini's master's thesis (refer earlier article) TOMS ACM Transactions on Mathematical Software Algorithm #531 Thanks again. -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;ig_x)-1;i++) { xavg = (x[i] + x[i+1]) / 2.0; for(j=0;jg_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[i], y[j], zij, x[i], y[j+1], zijp, xavg, yavg, zavg, level, color); graph_segment(gp, x[i], 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[i], y[j+1], zijp, x[i+1], y[j+1], zipjp, xavg, yavg, zavg, level, color); } } } }