/* File: drac.c Authors: Andrew W. Moore, Mary Soon Lee Created: Fri Jan 21 09:56:04 EST 1994 Description: DRAwing Contours Last update: 7 Mar 94 Copyright (C) 1994, Andrew W. Moore, Mary Soon Lee */ #include #include #include "drac.h" /* Drawing contours of functions and 2d arrays */ #define RADIUS 10.0 #define NOT_NEIGH 10 #define NUM_CONTOURS 10 #define NUM_COLUMNS 25 #define NUM_ROWS 25 #define DEFAULT_BL_WINDOW_X 0.0 #define DEFAULT_BL_WINDOW_Y 0.0 #define DEFAULT_TR_WINDOW_X 512.0 #define DEFAULT_TR_WINDOW_Y 512.0 #define FRAME_X_BOTTOM 60.0 #define FRAME_Y_BOTTOM 60.0 #define FRAME_X_TOP 30.0 #define FRAME_Y_TOP 30.0 #define DEFAULT_BL_DOMAIN_X 0.0 #define DEFAULT_BL_DOMAIN_Y 0.0 #define DEFAULT_TR_DOMAIN_X 10.0 #define DEFAULT_TR_DOMAIN_Y 10.0 typedef tgpoint *tgpoint_ptr; typedef tgpoint_ptr *tgpoint_ptr_ptr; float default_h_fn(char *data,float x, float y) { float result; float tem1 = (x - 2.0) * (x - 2.0) + (y - 4.0) * (y - 4.0); float tem2 = (2 * x - 8.0) * (x - 8.0) + (y - 6.0) * (y - 6.0); if ( data != NULL ) my_error("default_h_fn expected NULL data"); if (tem1 < tem2) result = tem1; else result = tem2; return(result); } tgrid *malloc_tgrid(int cols, int rows) { tgrid *tg = AM_MALLOC(tgrid); /* see ../amut/amma.h for definition of AM_MALLOC macro */ int c,r; tg -> num_cols = cols; tg -> num_rows = rows; tg -> tgpoints = AM_MALLOC_ARRAY(tgpoint_ptr_ptr,rows); for ( r = 0 ; r < rows ; r++ ) { tg -> tgpoints[r] = AM_MALLOC_ARRAY(tgpoint_ptr,cols); for ( c = 0 ; c < cols ; c++ ) tg -> tgpoints[r][c] = AM_MALLOC(tgpoint); } return(tg); } void free_tgrid(tgrid *tg) { int cols = tg->num_cols; int rows = tg->num_rows; int c,r; for ( r = 0 ; r < rows ; r++ ) { for ( c = 0 ; c < cols ; c++ ) am_free((char *)tg -> tgpoints[r][c],sizeof(tgpoint)); am_free((char *)tg -> tgpoints[r],sizeof(tgpoint_ptr) * cols); } am_free((char *)tg->tgpoints,sizeof(tgpoint_ptr_ptr) * rows); am_free((char *)tg,sizeof(tgrid)); } tgpoint *tg_ref(tgrid *tg, int c, int r) { if ( c < 0 || c >= tg -> num_cols || r < 0 || r >= tg -> num_rows ) /* my_error("pwicnbap"); */ { printf("ERROR"); wait_for_key(); my_error("blob"); } return(tg->tgpoints[r][c]); } /* Initializes a stage */ void init_stage(stage *st, drac_params *par) { (st->bl_window).x = DEFAULT_BL_WINDOW_X + FRAME_X_BOTTOM; (st->bl_window).y = DEFAULT_BL_WINDOW_Y + FRAME_Y_BOTTOM; (st->tr_window).x = DEFAULT_TR_WINDOW_X - FRAME_X_TOP; (st->tr_window).y = DEFAULT_TR_WINDOW_Y - FRAME_Y_TOP; (st->bl_domain).x = par->x_low; (st->bl_domain).y = par->y_low; (st->tr_domain).x = par->x_high; (st->tr_domain).y = par->y_high; } /* Draws in the boundary of the stage */ void frame_stage(stage *st) { float x_min = (st->bl_window).x; float y_min = (st->bl_window).y; float x_max = (st->tr_window).x; float y_max = (st->tr_window).y; int i; int size; float lo,hi,delta; ag_rectangle(x_min,y_min,x_max,y_max); sensible_limits(st->bl_domain.x,st->tr_domain.x,&lo,&hi,&delta); size = (int) floor(0.5 + (hi - lo)/delta); if ( Verbosity > 15.55 ) { printf("lo = %g\n",lo); printf("hi = %g\n",hi); printf("delta = %g\n",delta); printf("size = %d\n",size); } for ( i = 0 ; i <= size ; i++ ) { float d = lo + delta * i; float normal = (d - st->bl_domain.x) / (st->tr_domain.x - st->bl_domain.x); if ( normal >= 0.0 && normal <= 1.0 ) { float s = st->bl_window.x * (1 - normal) + st->tr_window.x * normal; char buff[100]; ag_line(s,y_min,s,y_min - 4.0); sprintf(buff,"%g",d); ag_print(s,y_min - 15.0,buff); } } sensible_limits(st->bl_domain.y,st->tr_domain.y,&lo,&hi,&delta); size = (int) floor(0.5 + (hi - lo)/delta); if ( Verbosity > 15.55 ) { printf("lo = %g\n",lo); printf("hi = %g\n",hi); printf("delta = %g\n",delta); printf("size = %d\n",size); } for ( i = 0 ; i <= size ; i++ ) { float d = lo + delta * i; float normal = (d - st->bl_domain.y) / (st->tr_domain.y - st->bl_domain.y); if ( normal >= 0.0 && normal <= 1.0 ) { float s = st->bl_window.y * (1 - normal) + st->tr_window.y * normal; char buff[100]; ag_line(x_min,s,x_min - 4.0,s); sprintf(buff,"%g",d); ag_print(x_min - FRAME_X_BOTTOM + 5.0,s,buff); } } } /* Draws in the boundary of the stage and titles it */ void frame_and_title_stage(stage *st, char *title) { frame_stage(st); ag_print((st->bl_window).x, (st->tr_window).y + FRAME_Y_TOP - 10.0, title); } /* Given an x in h_fn coordinates, this returns the corresponding x in the graphics window */ float x_tog(stage *st, float x) { float dom_xmin = (st->bl_domain).x; float dom_width = (st->tr_domain).x - dom_xmin; float win_xmin = (st->bl_window).x; float win_width = (st->tr_window).x - win_xmin; float result = win_xmin + (((x - dom_xmin) * win_width) / dom_width); if ((x < dom_xmin) || (x > (st->tr_domain).x)) my_error("x out of bounds"); return(result); } /* Given y in h_fn coordinates, this returns the corresponding y in the graphics window */ float y_tog(stage *st, float y) { float dom_ymin = (st->bl_domain).y; float dom_height = (st->tr_domain).y - dom_ymin; float win_ymin = (st->bl_window).y; float win_height = (st->tr_window).y - win_ymin; float result = win_ymin + (((y - dom_ymin) * win_height) / dom_height); if ((y < dom_ymin) || (y > (st->tr_domain).y)) my_error("y out of bounds"); return(result); } /* Given an x in graphics coordinates, this returns the corresponding x in h_fn coordinates */ float x_fromg(stage *st, float x) { float dom_xmin = (st->bl_domain).x; float dom_width = (st->tr_domain).x - dom_xmin; float win_xmin = (st->bl_window).x; float win_width = (st->tr_window).x - win_xmin; float result = dom_xmin + (((x - win_xmin) * dom_width) / win_width); if ((x < win_xmin) || (x > (st->tr_window).x)) my_error("x_g out of bounds"); return(result); } /* Given y in graphics coordinates, this returns the corresponding y in h_fn coordinates */ float y_fromg(stage *st, float y) { float dom_ymin = (st->bl_domain).y; float dom_height = (st->tr_domain).y - dom_ymin; float win_ymin = (st->bl_window).y; float win_height = (st->tr_window).y - win_ymin; float result = dom_ymin + (((y - win_ymin) * dom_height) / win_height); if ((y < win_ymin) || (y > (st->tr_window).y)) my_error("y_g out of bounds"); return(result); } /* Given a point p in h_fn coordinates, this draws it in the right place in the graphics window. */ void drac_dot(stage *st, point *p) { ag_dot(x_tog(st,p->x), y_tog(st,p->y)); } /* Given a point p in h_fn coordinates, this highlights the corresponding place in the graphics window with a circle. */ void drac_highlight(stage *st, point *p) { ag_circle(x_tog(st,p->x), y_tog(st,p->y), RADIUS); } /* Given points p1, p2 in h_fn coordinates, this draws a line between the corresponding points in the graphics window. */ void drac_line(stage *st, point *p1, point *p2) { ag_line(x_tog(st,p1->x), y_tog(st,p1->y), x_tog(st,p2->x), y_tog(st,p2->y)); } /* Draws the point represented by tp */ void draw_tgpoint(stage *st, tgpoint *tp) { drac_dot(st, &(tp->point)); } /* Highlights the point represented by tp with a circle */ void highlight_tgpoint(stage *st, tgpoint *tp) { drac_highlight(st, &(tp->point)); } /* Initializes the tg_i, tg_j, index, point, and height fields in the tg_points of a tg_grid. */ void init_tg_basics( drac_params *par, tgrid *tg, float (*h_fn)(char *data,float x,float y), char *data ) { int cols = tg->num_cols; int rows = tg->num_rows; int i, j; float tpx, tpy; tgpoint *tp; if ((cols <= 2) || (rows <= 2)) my_error("Too few columns or rows for a tg_grid"); for (i = 0; i < cols; i++) for (j = 0; j < rows; j++) { tp = tg_ref(tg, i, j); tp->tg_i = i; tp->tg_j = j; tp->index = i + j * cols; tpx = par->x_low + (i / ((float) (cols - 1))) * (par->x_high - par->x_low); tpy = par->y_low + (j / ((float) (rows - 1))) * (par->y_high - par->y_low); (tp->point).x = tpx; (tp->point).y = tpy; tp->height = h_fn(data,tpx, tpy); } } /* Initializes the neighbor fields in the tg_points of a tg_grid. This is an ugly function, split into lots of special cases. */ void init_tg_neighs(tgrid *tg) { int cols = tg->num_cols; int rows = tg->num_rows; int i, j; tgpoint *tp; tgpoint **neighs; /* first set loop to be FALSE for all edge-points, TRUE otherwise */ for (i = 0; i < cols; i++) for (j = 0; j < rows; j++) { tp = tg_ref(tg, i, j); if ((i == 0) || (j == 0) || (i == (cols - 1)) || (j == (rows - 1))) tp->loop = FALSE; else tp->loop = TRUE; } /* the 4 corner points */ tp = tg_ref(tg, 0, 0); tp->num_neighs = 2; neighs = tp->neighs; neighs[0] = tg_ref(tg, 0, 1); neighs[1] = tg_ref(tg, 1, 0); tp = tg_ref(tg, 0, rows - 1); tp->num_neighs = 3; neighs = tp->neighs; neighs[0] = tg_ref(tg, 1, rows - 1); neighs[1] = tg_ref(tg, 1, rows - 2); neighs[2] = tg_ref(tg, 0, rows - 2); tp = tg_ref(tg, cols - 1, rows - 1); tp->num_neighs = 2; neighs = tp->neighs; neighs[0] = tg_ref(tg, cols - 1, rows - 2); neighs[1] = tg_ref(tg, cols - 2, rows - 1); tp = tg_ref(tg, cols - 1, 0); tp->num_neighs = 3; neighs = tp->neighs; neighs[0] = tg_ref(tg, cols - 2, 0); neighs[1] = tg_ref(tg, cols - 2, 1); neighs[2] = tg_ref(tg, cols - 1, 1); /* the center points */ for (i = 1; i < cols - 1; i++) for (j = 1; j < rows - 1; j++) { tp = tg_ref(tg, i, j); tp->num_neighs = 6; neighs = tp->neighs; neighs[0] = tg_ref(tg, i, j + 1); neighs[1] = tg_ref(tg, i + 1, j); neighs[2] = tg_ref(tg, i + 1, j - 1); neighs[3] = tg_ref(tg, i, j - 1); neighs[4] = tg_ref(tg, i - 1, j); neighs[5] = tg_ref(tg, i - 1, j + 1); } /* the 4 edges */ for (j = 1; j < rows - 1; j++) /* left edge */ { tp = tg_ref(tg, 0, j); tp->num_neighs = 4; neighs = tp->neighs; neighs[0] = tg_ref(tg, 0, j + 1); neighs[1] = tg_ref(tg, 1, j); neighs[2] = tg_ref(tg, 1, j - 1); neighs[3] = tg_ref(tg, 0, j - 1); } for (i = 1; i < cols - 1; i++) /* top edge */ { tp = tg_ref(tg, i, rows - 1); tp->num_neighs = 4; neighs = tp->neighs; neighs[0] = tg_ref(tg, i + 1, rows - 1); neighs[1] = tg_ref(tg, i + 1, rows - 2); neighs[2] = tg_ref(tg, i, rows - 2); neighs[3] = tg_ref(tg, i - 1, rows - 1); } for (j = 1; j < rows - 1; j++) /* right edge */ { tp = tg_ref(tg, cols - 1, j); tp->num_neighs = 4; neighs = tp->neighs; neighs[0] = tg_ref(tg, cols - 1, j - 1); neighs[1] = tg_ref(tg, cols - 2, j); neighs[2] = tg_ref(tg, cols - 2, j + 1); neighs[3] = tg_ref(tg, cols - 1, j + 1); } for (i = 1; i < cols - 1; i++) /* bottom edge */ { tp = tg_ref(tg, i, 0); tp->num_neighs = 4; neighs = tp->neighs; neighs[0] = tg_ref(tg, i - 1, 0); neighs[1] = tg_ref(tg, i - 1, 1); neighs[2] = tg_ref(tg, i, 1); neighs[3] = tg_ref(tg, i + 1, 0); } } /* Initializes the tg_points of a tg_grid. */ void init_tg( drac_params *par, tgrid *tg, float (*h_fn)(char *data,float x,float y), char *data ) { init_tg_basics(par, tg, h_fn, data); init_tg_neighs(tg); } /* Draws in all the points on a tgrid. */ void draw_tgrid(stage *st, tgrid *tg) { int cols = tg->num_cols; int rows = tg->num_rows; int i, j; for (i = 0; i < cols; i++) for (j = 0; j < rows; j++) { draw_tgpoint(st, tg_ref(tg, i, j)); } } /* Returns the nearest tgpoint to p (p in h_fn coordinates). "nearest" means the tgpoint whose x coordinate is closest to p's x coordinate, and similarly for y. */ tgpoint *nearest_tgpoint(stage *st, tgrid *tg, point *p) { int cols = tg->num_cols; int rows = tg->num_rows; float dom_xmin = st->bl_domain.x; float dom_ymin = st->bl_domain.y; float dom_xmax = st->tr_domain.x; float dom_ymax = st->tr_domain.y; float dom_width = dom_xmax - dom_xmin; float dom_height = dom_ymax - dom_ymin; int i_l; /* the tg_i of the column of tgpoints just to the left of p */ int j_b; /* the tg_j of the row of tgpoints just below p */ tgpoint *p1, *p2; int i_needed, j_needed; if ((p->x <= dom_xmin) || (p->x >= dom_xmax) || (p->y <= dom_ymin) || (p->y >= dom_ymax)) my_error("p is outside the region being considered"); i_l = floor(((p->x - dom_xmin) * ((float) (cols - 1))) / dom_width); j_b = floor(((p->y - dom_ymin) * ((float)( rows - 1))) / dom_height); /* the i we need is either i_l or i_l + 1, similarly for j.... */ p1 = tg_ref(tg, i_l, j_b); p2 = tg_ref(tg, i_l + 1, j_b + 1); if ((p->x - p1->point.x) < (p2->point.x - p->x)) i_needed = i_l; else i_needed = i_l + 1; if ((p->y - p1->point.y) < (p2->point.y - p->y)) j_needed = j_b; else j_needed = j_b + 1; return(tg_ref(tg, i_needed, j_needed)); } /* Waits for the mouse to be clicked within the stage; then returns the click's location IN H_FN COORDINATES in p. */ void drac_use_mouse(stage *st, point *p) { float x, y; bool not_yet_found = TRUE; float win_xmin = st->bl_window.x; float win_ymin = st->bl_window.y; float win_xmax = st->tr_window.x; float win_ymax = st->tr_window.y; while (not_yet_found) { ag_get_xy(&x, &y); if ((x <= win_xmin) || (x >= win_xmax) || (y <= win_ymin) || (y >= win_ymax)) printf("Clicked outside stage.... Try again.\n"); else not_yet_found = FALSE; } p->x = x_fromg(st, x); p->y = y_fromg(st, y); } /* Draws lines connecting all the neighbors of tgp. */ void display_loop(stage *st, tgpoint *tgp) { bool is_loop = tgp->loop; int n; int max_n = tgp->num_neighs; tgpoint **neighs = tgp->neighs; tgpoint *tgp1, *tgp2; if (is_loop) { for (n = 0; n < max_n; n++) { tgp1 = neighs[n]; if (n == (max_n - 1)) tgp2 = neighs[0]; else tgp2 = neighs[n + 1]; drac_line(st, &(tgp1->point), &(tgp2->point)); } } else { for (n = 0; n < (max_n - 1); n++) { tgp1 = neighs[n]; tgp2 = neighs[n + 1]; drac_line(st, &(tgp1->point), &(tgp2->point)); } } } /* If tgp1 and tgp2 are neighbors, this returns i s.t. tgp1->neighs[i] = tgp2 else it returns NOT_NEIGH */ int neigh_num(tgpoint *tgp1, tgpoint *tgp2) { int i; int result = NOT_NEIGH; for (i = 0; i < tgp1->num_neighs; i++) { if (tgp1->neighs[i] == tgp2) result = i; } return(result); } /* This draws an egde */ void draw_tg_edge(stage *st, edge *e) { drac_line(st, &((e->start)->point), &((e->end)->point)); } typedef edge *edge_ptr; typedef edge_ptr *edge_ptr_ptr; /* This creates the egrid given the number of columns and rows in the triangulation grid. */ egrid *malloc_egrid(int cols, int rows) { egrid *eg = AM_MALLOC(egrid); /* see ../amut/amma.h for definition of AM_MALLOC macro */ int i, tgpoint_index; int num_tgpoints = cols * rows; eg -> num_cols = cols; eg -> num_rows = rows; eg -> edges = AM_MALLOC_ARRAY(edge_ptr_ptr, num_tgpoints); for (tgpoint_index = 0; tgpoint_index < num_tgpoints; tgpoint_index++) { eg -> edges[tgpoint_index] = AM_MALLOC_ARRAY(edge_ptr, MAX_NEIGHS); for (i = 0; i < MAX_NEIGHS; i++) eg -> edges[tgpoint_index][i] = AM_MALLOC(edge); } return(eg); } void free_egrid(egrid *eg) { int rows = eg -> num_rows; int cols = eg -> num_cols; int i, tgpoint_index; int num_tgpoints = cols * rows; for (tgpoint_index = 0; tgpoint_index < num_tgpoints; tgpoint_index++) { for (i = 0; i < MAX_NEIGHS; i++) am_free((char *)eg -> edges[tgpoint_index][i],sizeof(edge)); am_free((char *)eg -> edges[tgpoint_index],MAX_NEIGHS * sizeof(edge_ptr)); } am_free((char *)eg->edges,num_tgpoints * sizeof(edge_ptr_ptr)); am_free((char *)eg,sizeof(egrid)); } /* Given two neighboring tgpoints, this returns the edge between them */ edge *e_ref(egrid *eg, tgpoint *tgp1, tgpoint *tgp2) { int i = neigh_num(tgp1, tgp2); if (i == NOT_NEIGH) my_error("Trying to find edge between non-neighboring points"); return(eg->edges[tgp1->index][i]); } /* Given a tgpoint and the neigh_num of its neighbors this returns the edge between them with no checks */ edge *fast_e_ref(egrid *eg, tgpoint *tgp1, int i) { return(eg->edges[tgp1->index][i]); } /* Given the index for a tgpoint and the neigh_num of its neighbors this returns the edge between them with no checks */ edge *vfast_e_ref(egrid *eg, int index, int i) { return(eg->edges[index][i]); } /* Initializes edge to run between tgp1 and tgp2, and fills in the is_edge and visited fields */ void init_edge_basic(edge *e, tgpoint *tgp1, tgpoint *tgp2) { e->start = tgp1; e->end = tgp2; e->visited = FALSE; e->is_edge = TRUE; } /* Initializes the start, end, visited and is_edge fields of the edges in an egrid */ void init_eg_basics(tgrid *tg, egrid *eg) { int cols = tg->num_cols; int rows = tg->num_rows; int i, j, k; tgpoint *tgp; edge *e; for (j = 0; j < rows ; j++) for (i = 0; i < cols; i++) { tgp = tg_ref(tg, i, j); for (k = 0; k < tgp->num_neighs; k++) { e = fast_e_ref(eg, tgp, k); init_edge_basic(e, tgp, tgp->neighs[k]); } for (k = tgp->num_neighs; k < MAX_NEIGHS; k++) { e = fast_e_ref(eg, tgp, k); e->is_edge = FALSE; } } } /* Initializes the remaining fields in the edges of an egrid, i.e. border_edge, tr1_e1, tr1_e2, tr2_e1, tr2_e2 */ void init_eg_rest(stage *st, tgrid *tg, egrid *eg) { int num_tgpoints = tg->num_cols * tg->num_rows; int i, j, k, l; edge *e; tgpoint *tgp1, *tgp2; bool found_neigh1, found_neigh2; for (i = 0; i < num_tgpoints; i++) for (j = 0; j < MAX_NEIGHS; j++) { e = vfast_e_ref(eg, i, j); if (e->is_edge == TRUE) { tgp1 = e->start; tgp2 = e->end; e->border_edge = TRUE; /* if e is not on the border, this should later get set to FALSE */ /* Find the common neighbors of tgp1 and tgp2; there should be at least one and at most two. Once a neighbor is found, fill in the appropriate fields in e */ found_neigh1 = FALSE; found_neigh2 = FALSE; for (k = 0; k < tgp1->num_neighs; k++) for (l = 0; l < tgp2->num_neighs; l++) { if (tgp1->neighs[k] == tgp2->neighs[l]) { if (found_neigh1 == FALSE) { found_neigh1 = TRUE; e->tr1_e1 = e_ref(eg, tgp1, tgp1->neighs[k]); e->tr1_e2 = e_ref(eg, tgp2, tgp1->neighs[k]); } else if (found_neigh2 == FALSE) { found_neigh2 = TRUE; e->border_edge = FALSE; e->tr2_e1 = e_ref(eg, tgp1, tgp1->neighs[k]); e->tr2_e2 = e_ref(eg, tgp2, tgp1->neighs[k]); } else my_error("Edge points should not have three common neighbors"); } } if (found_neigh1 == FALSE) my_error("Failed to find a common neighbor during init_eg_rest"); } } } /* Fully initializes an egrid */ void init_eg(stage *st, tgrid *tg, egrid *eg) { init_eg_basics(tg, eg); init_eg_rest(st, tg, eg); } /* Draws the triangles in the triangulation grid that have e as an edge */ void draw_triangles(stage *st, edge *e) { draw_tg_edge(st, e->tr1_e1); draw_tg_edge(st, e->tr1_e2); if (e->border_edge == FALSE) { draw_tg_edge(st, e->tr2_e1); draw_tg_edge(st, e->tr2_e2); } } /* Returns TRUE iff edge e crosses height h, i.e. one end-point is strictly higher than h, the other less than or equal to h. */ bool crosses_h(edge *e, float h) { float h1 = (e->start)->height; float h2 = (e->end)->height; if (((h1 <= h) && (h2 > h)) || ((h1 > h) && (h2 <= h))) return(TRUE); else return(FALSE); } /* Given an edge, e, whose end tgpoints are on opposite sides of height h, this sets cross_pt to be the (x, y) point in h_fn coordinates that lies on the edge, e, and that linear interpolation gives as having height h */ void set_cross_pt(edge *e, float h, point *cross_pt) { tgpoint *tgp1 = e->start; tgpoint *tgp2 = e->end; float ratio = (tgp1->height - h) / (tgp1->height - tgp2->height); /* the ratio d1/d2, where d1 is the distance from tgp1 to the estimated crossing point of the h contour on e, and d2 is the distance from tgp2 to tgp1 */ cross_pt->x = (tgp1->point).x + ((tgp2->point).x - (tgp1->point).x) * ratio; cross_pt->y = (tgp1->point).y + ((tgp2->point).y - (tgp1->point).y) * ratio; } /* Given an edge, e, that crosses height h, this first performs linear interpolation to estimate where the h contour crosses the edge. Then for each triangle that e is in (at most two), it finds the edge that also crosses the h contour, draws in the estimate of the h contour across the triangle, and then iterates to continue tracing the contour around. The algorithm stops trying to extend the contour when all the relevant edges have either already been visited (as will eventually happen in a loop) or lead off the boundary of the grid. It updates p to be the (x, y) point in h_fn coordinates where the contour crosses e. */ void draw_contour(stage *st, edge *e, float h, point *cross_pt) { edge *e1, *e2; point p; int i; int num_triangles_to_try; e->visited = TRUE; set_cross_pt(e, h, cross_pt); if (e->border_edge == TRUE) num_triangles_to_try = 1; else num_triangles_to_try = 2; for (i = 0; i < num_triangles_to_try; i++) { if (i == 0) { /* will try first triangle */ e1 = e->tr1_e1; e2 = e->tr1_e2; } else { /* will try second triangle */ e1 = e->tr2_e1; e2 = e->tr2_e2; } if (crosses_h(e1, h) == TRUE) { if (e1->visited == FALSE) draw_contour(st, e1, h, &p); else /* even if the edge has been visited, we need to find its crossing point, p, and draw the line segment across this triangle, otherwise loops would not be completed */ set_cross_pt(e1, h, &p); drac_line(st, cross_pt, &p); } else if (crosses_h(e2, h) == TRUE) { if (e2->visited == FALSE) draw_contour(st, e2, h, &p); else set_cross_pt(e2, h, &p); drac_line(st, cross_pt, &p); } else my_error("h contour should cross one of e1, e2"); } } /* This marks all the edges in the grid as unvisited */ void clear_edges(egrid *eg) { int num_tgpoints = eg->num_cols * eg->num_rows; int i, j; edge *e; for (i = 0; i < num_tgpoints; i++) for (j = 0; j < MAX_NEIGHS; j++) { e = vfast_e_ref(eg, i, j); e->visited = FALSE; } } /* This draws in the contours at height h. We progress through the edges in the egrid. Whenever we find an unvisited edge whose endpoints are opposite sides of h (so that height h lies in between), we look at its neighboring edges and trace the contour round. Linear interpolation is used to decide where to draw the contour. */ void draw_h_contours(stage *st, egrid *eg, float h) { int num_tgpoints = eg->num_cols * eg->num_rows; int i, j; edge *e; point p; clear_edges(eg); for (i = 0; i < num_tgpoints; i++) for (j = 0; j < MAX_NEIGHS; j++) { e = vfast_e_ref(eg, i, j); if ((e->is_edge == TRUE) && (e->visited == FALSE)) { if (crosses_h(e, h) == TRUE) draw_contour(st, e, h, &p); } } } float tgrid_min(tgrid *tg) { int i,j; float result = 0.0; for (i = 0; i < tg->num_cols; i++) for(j = 0; j < tg->num_rows; j++) { float cur_h = tg_ref(tg, i, j)->height; if ( (i==0 && j==0) || cur_h < result ) result = cur_h; } return(result); } float tgrid_max(tgrid *tg) { int i,j; float result = 0.0; for (i = 0; i < tg->num_cols; i++) for(j = 0; j < tg->num_rows; j++) { float cur_h = tg_ref(tg, i, j)->height; if ( (i==0 && j==0) || cur_h > result ) result = cur_h; } return(result); } /* Autonomously draws in n contours for the egrid. It starts by finding the lowest and highest heights on the tgrid--low_h and high_h, then draws in contours at low_h + d, low_h + 2d, ... , low_h + n * d where d = (high_h - low_h)/(n + 1). N.B. There is little point in wasting effort attempting to draw contours at low_h and high_h. */ void draw_contours(stage *st, tgrid *tg, egrid *eg, int n) { int num_cols = tg->num_cols; int num_rows = tg->num_rows; float low_height; float high_height; float d; int i; char height_info[100]; /* first some checks */ if (n < 1) my_error("Asked to draw too few contours"); if ((num_cols != eg->num_cols) || (num_rows != eg->num_rows)) my_error("Egrid and tgrid have different dimensions"); /* now find high and low values */ low_height = tgrid_min(tg); high_height = tgrid_max(tg); /* now print the height information and draw the contours */ d = (high_height - low_height) / (n + 1); sprintf(height_info, "Contours from %.3g to %.3g in increments of %.3g", low_height + d, low_height + n * d, d); ag_print((st->bl_window).x, (st->bl_window).y - FRAME_Y_BOTTOM + 10.0, height_info ); for (i = 1; i <= n; i++) draw_h_contours(st, eg, low_height + i * d); } /* Reads the parameters from the command line arguments and puts them into par. */ void drac_params_from_args(drac_params *par, int argc, char **argv) { par->num_cols = int_from_args("-cols",argc,argv,NUM_COLUMNS); par->num_rows = int_from_args("-rows",argc,argv,NUM_ROWS); par->num_contours = int_from_args("-nc",argc,argv, NUM_CONTOURS); sprintf(par->rfile, string_from_args("-res",argc, argv,"contour.ps")); par->x_low = float_from_args("-xlow",argc,argv,DEFAULT_BL_DOMAIN_X); par->y_low = float_from_args("-ylow",argc,argv,DEFAULT_BL_DOMAIN_Y); par->x_high = float_from_args("-xhigh",argc,argv,DEFAULT_TR_DOMAIN_X); par->y_high = float_from_args("-yhigh",argc,argv,DEFAULT_TR_DOMAIN_Y); sprintf(par->title, string_from_args("-title",argc, argv,"CONTOURS")); sprintf(par->points, string_from_args("-points",argc, argv,"")); } /* Given the name of a file, f, this scans through what should be sets of (x, y, z) triples to determine the number of rows and columns in the grid. The points must occur in the following order: scanning along rows from the minimum x to the maximum x, starting with the bottom row (minimum y) and rising to the top row. The points must lie on a rectangular grid, but there is no check for this. */ void size_of_grid_in_file(char *f, int *cols, int *rows) { FILE *s = fopen(f, "r"); float cur_x, cur_y, cur_z; float min_x = 0.0; /* Just to prevent unitialized warning */ int i = 0; int j = 0; int n = 0; /* number of triples scanned */ int q; bool NOTHING_READ = TRUE; bool STOP = FALSE; if (s == NULL) my_error("Failed when trying to open a file to read in a grid"); while (STOP == FALSE) { q = fscanf(s, "%f %f %f", &cur_x, &cur_y, &cur_z); if (q == EOF) /* test for end of file */ STOP = TRUE; else if (q == 3) /* found a new triple */ { n++; if (NOTHING_READ == TRUE) /* this was the 1st triple; set min_x */ { min_x = cur_x; NOTHING_READ = FALSE; i = 1; } else if (min_x == cur_x) /* at the start of a new row */ { i = 1; j++; } else i++; } else my_error("Error in format of grid file"); } if (n != i * (j + 1)) my_error("Wrong number of points in grid file"); *cols = i; *rows = j + 1; fclose(s); } /* Initializes a tgrid from a file. The dimensions of the tgrid (i.e. num_cols, num_rows) should already be correctly set. */ void init_tg_from_file(tgrid *tg, char *f, drac_params *par) { int cols = tg->num_cols; int rows = tg->num_rows; int i, j; tgpoint *tp; FILE *s = fopen(f, "r"); float cur_x, cur_y, cur_z; if (s == NULL) my_error("Failed when trying to re-open a file to read in a grid"); if ((cols <= 2) || (rows <= 2)) my_error("Too few columns or rows for a tg_grid"); for (j = 0; j < rows; j++) for (i = 0; i < cols; i++) { fscanf(s, "%f %f %f", &cur_x, &cur_y, &cur_z); tp = tg_ref(tg, i, j); tp->tg_i = i; tp->tg_j = j; tp->index = i + j * cols; (tp->point).x = cur_x; (tp->point).y = cur_y; tp->height = cur_z; } init_tg_neighs(tg); /* Now set the bl_domain and tr_domain points correctly in par */ par->x_low = tg_ref(tg, 0, 0)->point.x; par->y_low = tg_ref(tg, 0, 0)->point.y; par->x_high = tg_ref(tg, cols - 1, 0)->point.x; par->y_high = tg_ref(tg, 0, rows - 1)->point.y; } /* Initializes the structures from the parameters. Note that if the points are being read in from a grid in a file, then parameters such as num_cols are deduced from that file, not taken from the command line. */ void init_structs_from_drac_params( drac_params *par, tgrid **tg, egrid **eg, stage *st, float (*h_fn)(char *data,float x,float y), char *data ) { bool read_grid_from_file; if (strcmp(par->points, "") != 0) read_grid_from_file = TRUE; else read_grid_from_file = FALSE; if (read_grid_from_file == TRUE) size_of_grid_in_file(par->points, &par->num_cols, &par->num_rows); *tg = malloc_tgrid(par->num_cols, par->num_rows); *eg = malloc_egrid(par->num_cols, par->num_rows); if (read_grid_from_file == TRUE) init_tg_from_file(*tg, par->points, par); else init_tg(par, *tg, h_fn,data); init_stage(st, par); init_eg(st, *tg, *eg); } /* Gives the user help on using the program */ void drac_help(void) { printf("\n This is a contour drawing program. There are two ways to\n"); printf("use it. Either you can define the height function by changing\n"); printf("h_fn(x, y) in the file ~awm/w/drac/drac.c, or else you can specify\n"); printf("a file that contains the (x,y,z) coordinates of the points on a\n"); printf("rectangular grid.\n"); printf(" The program constructs a triangulation grid, finds the\n"); printf("contours, displays them on the screen, and saves the diagram\n"); printf("to a postscript file.\n"); printf(" To run the program, type drac followed by any arguments:-\n"); printf("-cols the number of columns in the grid (default = %d)\n", NUM_COLUMNS); printf("-rows the number of rows in the grid (default = %d)\n", NUM_ROWS); printf("-nc the number of contours to draw (default = %d)\n", NUM_CONTOURS); printf("-res the file to which to save the postscript (default = contour.ps)\n"); printf("-title the title for the diagram (default = CONTOURS)\n"); printf("-xlow the minimum x value to be mapped (default = %g)\n", DEFAULT_BL_DOMAIN_X); printf("-ylow the minimum y value to be mapped (default = %g)\n", DEFAULT_BL_DOMAIN_Y); printf("-xhigh the maximum x value to be mapped (default = %g)\n", DEFAULT_TR_DOMAIN_X); printf("-yhigh the maximum y value to be mapped (default = %g)\n", DEFAULT_TR_DOMAIN_Y); printf("E.g. 'drac -nc 4' would produce a diagram with 4 contours\n"); printf(" If you use a file to give the (x, y, z) coordinates, the\n"); printf("points must be given in the following order: scan along the\n"); printf("rows from the bottom row (minimum y) to the top row, going\n"); printf("from the minimum x to the maximum x along each row. For each\n"); printf("point, give the x, y, z values in turn. The points must lie\n"); printf("on a rectangular grid. Values for the number of columns, the\n"); printf("minimum x value to be mapped, etc, are then deduced from the\n"); printf("file and override any values that were specified in the\n"); printf("command line.\n"); } void drac_main(int argc, char **argv) { drac_params par; tgrid *tg; egrid *eg; stage st; if ((argc > 1) && (strcmp(argv[1], "help") == 0)) drac_help(); else { drac_params_from_args(&par, argc, argv); init_structs_from_drac_params(&par, &tg, &eg, &st, default_h_fn,NULL); ag_on(par.rfile); frame_and_title_stage(&st, par.title); draw_contours(&st, tg, eg, par.num_contours); printf("Hit return to quit\n"); wait_for_key(); ag_off(); } } void fprintf_dproj(FILE *s, char *m1, dproj *dp, char *m2) { char buff[100]; fprintf(s,"%s -> screen_x_min = %g%s",m1,dp->screen_x_min,m2); fprintf(s,"%s -> screen_x_max = %g%s",m1,dp->screen_x_max,m2); fprintf(s,"%s -> screen_y_min = %g%s",m1,dp->screen_y_min,m2); fprintf(s,"%s -> screen_y_max = %g%s",m1,dp->screen_y_max,m2); fprintf(s,"%s -> indim = %d%s",m1,dp->indim,m2); sprintf(buff,"%s -> in_min",m1); fprintf_floats(s,buff,dp->in_min,dp->indim,m2); sprintf(buff,"%s -> in_max",m1); fprintf_floats(s,buff,dp->in_max,dp->indim,m2); fprintf(s,"%s -> height_min = %g%s",m1,dp->height_min,m2); fprintf(s,"%s -> height_max = %g%s",m1,dp->height_max,m2); } float ag_value( float x, float x_min, float x_max, float screen_min, float screen_max ) { return(screen_min + (screen_max - screen_min) * (x - x_min) / (x_max - x_min) ); } float dproj_to_ag_x(dproj *dp,float *farr,float height) { if ( dp->indim < 1 || dp -> indim > 2 ) my_error("drac.c :: dprooj_to_ag_x: only do this for 1 or 2 dimensions"); return(ag_value(farr[0],dp->in_min[0],dp->in_max[0], dp->screen_x_min,dp->screen_x_max ) ); } float dproj_to_ag_y(dproj *dp,float *farr,float height) { float y,y_min,y_max; float result; if ( dp->indim < 1 || dp -> indim > 2 ) my_error("drac.c :: dprooj_to_ag_y: only do this for 1 or 2 dimensions"); if ( dp->indim == 1 ) { y = height; y_min = dp->height_min; y_max = dp->height_max; } else { y = farr[1]; y_min = dp->in_min[1]; y_max = dp->in_max[1]; } result = ag_value(y,y_min,y_max,dp->screen_y_min,dp->screen_y_max); return(result); } void dproj_dot(dproj *dp,float *farr,float hgt) { /* fprintf_floats(stdout,"dproj_dot, farr = ",farr,dp->indim,"\n"); */ ag_dot(dproj_to_ag_x(dp,farr,hgt),dproj_to_ag_y(dp,farr,hgt)); } void dproj_line(dproj *dp,float *farr1,float hgt1,float *farr2,float hgt2) { ag_line(dproj_to_ag_x(dp,farr1,hgt1),dproj_to_ag_y(dp,farr1,hgt1), dproj_to_ag_x(dp,farr2,hgt2),dproj_to_ag_y(dp,farr2,hgt2) ); } void dproj_string(dproj *dp,float *farr,float hgt,char *string) { ag_print(dproj_to_ag_x(dp,farr,hgt),dproj_to_ag_y(dp,farr,hgt), string ); } void dproj_small_cross(dproj *dp,float *farr,float hgt) { float x = dproj_to_ag_x(dp,farr,hgt); float y = dproj_to_ag_y(dp,farr,hgt); float r = 5.0; ag_line(x-r,y-r,x+r,y+r); ag_line(x-r,y+r,x+r,y-r); } void dproj_small_circle(dproj *dp,float *farr,float hgt) { float x = dproj_to_ag_x(dp,farr,hgt); float y = dproj_to_ag_y(dp,farr,hgt); float r = 5.0; ag_circle(x,y,r); } void dproj_circle(dproj *dp,float *farr,float hgt,float radius) { float x = dproj_to_ag_x(dp,farr,hgt); float y = dproj_to_ag_y(dp,farr,hgt); float farr2[2]; float r; farr2[0] = x + radius; farr2[1] = 0.0; r = dproj_to_ag_x(dp,farr2,hgt) - x; ag_circle(x,y,r); } void dproj_dyv_dot(dproj *dp,dyv *d,float hgt) { float *farr = mk_farr_from_dyv(d); dproj_dot(dp,farr,hgt); am_free_floats(farr,dyv_size(d)); } void dproj_dyv_line(dproj *dp,dyv *d1,float hgt1,dyv *d2,float hgt2) { float *farr1 = mk_farr_from_dyv(d1); float *farr2 = mk_farr_from_dyv(d2); dproj_line(dp,farr1,hgt1,farr2,hgt2); am_free_floats(farr1,dyv_size(d1)); am_free_floats(farr2,dyv_size(d2)); } void dproj_dyv_string(dproj *dp,dyv *d,float hgt,char *string) { float *farr = mk_farr_from_dyv(d); dproj_string(dp,farr,hgt,string); am_free_floats(farr,dyv_size(d)); } void dproj_dyv_small_cross(dproj *dp,dyv *dv,float hgt) { float *farr = mk_farr_from_dyv(dv); dproj_small_cross(dp,farr,hgt); am_free_floats(farr,dyv_size(dv)); } void dproj_dyv_small_circle(dproj *dp,dyv *dv,float hgt) { float *farr = mk_farr_from_dyv(dv); dproj_small_circle(dp,farr,hgt); am_free_floats(farr,dyv_size(dv)); } void dproj_dyv_circle(dproj *dp,dyv *dv,float hgt,float radius) { float *farr = mk_farr_from_dyv(dv); dproj_circle(dp,farr,hgt,radius); am_free_floats(farr,dyv_size(dv)); } void draw_a_height( dproj *dp, float (*height_function)(char *data, float x, float y), char *data, float *farr ) { char buff[100]; float height = height_function(data,farr[0],farr[1]); if ( fabs(height) < 1e-5 ) height = 0.0; sprintf(buff,"%.3g",height); /* dproj_dot(dp,farr,height); */ dproj_string(dp,farr,height,buff); } void draw_nine_heights( dproj *dp, float (*height_function)(char *data, float x, float y), char *data ) { int i,j; for ( i = -1 ; i <= 1 ; i++ ) for ( j = -1 ; j <= 1 ; j++ ) { float farr[2]; farr[0] = (dp->in_min[0] + dp->in_max[0]) / 2.0 + 0.45 * i * (dp->in_max[0] - dp->in_min[0]); farr[1] = (dp->in_min[1] + dp->in_max[1]) / 2.0 + 0.45 * j * (dp->in_max[1] - dp->in_min[1]); draw_a_height(dp,height_function,data,farr); } } void basic_draw_2d_function( bool draw_borders_only, float (*height_function)(char *data, float x, float y), char *data, int num_contours, int grid_size, char *title, float x_low, float y_low, float x_high, float y_high, dproj *result_dproj ) { int fake_argc = 0; char **fake_argv = NULL; drac_params pa; tgrid *tg; egrid *eg; stage st; drac_params_from_args(&pa,fake_argc,fake_argv); pa.num_cols = grid_size; /* Number of columns in the grids */ pa.num_rows = grid_size; /* Number of rows in the grids */ pa.num_contours = num_contours; /* The number of contours to draw */ pa.x_low = x_low; pa.x_high = x_high; pa.y_low = y_low; pa.y_high = y_high; /* Domains for x and y */ init_structs_from_drac_params(&pa,&tg,&eg,&st,height_function,data); frame_and_title_stage(&st, (title == NULL) ? "Contours" : title); result_dproj -> screen_x_min = st.bl_window.x; result_dproj -> screen_y_min = st.bl_window.y; result_dproj -> screen_x_max = st.tr_window.x; result_dproj -> screen_y_max = st.tr_window.y; result_dproj -> indim = 2; result_dproj -> in_min[0] = st.bl_domain.x; result_dproj -> in_max[0] = st.tr_domain.x; result_dproj -> in_min[1] = st.bl_domain.y; result_dproj -> in_max[1] = st.tr_domain.y; result_dproj -> height_min = tgrid_min(tg); result_dproj -> height_max = tgrid_max(tg); if ( Verbosity > 15.55 ) fprintf_dproj(stdout,"dproj",result_dproj,"\n"); if ( !draw_borders_only ) { draw_contours(&st, tg, eg, pa.num_contours); draw_nine_heights(result_dproj,height_function,data); } free_egrid(eg); free_tgrid(tg); } void draw_2d_function( float (*height_function)(char *data, float x, float y), char *data, int num_contours, int grid_size, char *title, float x_low, float y_low, float x_high, float y_high, dproj *result_dproj ) { basic_draw_2d_function(FALSE, height_function, data, num_contours, grid_size, title, x_low, y_low, x_high, y_high, result_dproj ); } float dummy_height_fn(char *data,float x,float y) { return(0.0); } void draw_2d_border(dproj *dp,float xlo,float ylo,float xhi,float yhi) { basic_draw_2d_function(TRUE, dummy_height_fn, (char *)NULL, 3, 3, "", xlo, ylo, xhi, yhi, dp ); } void simple_draw_2d_function( float (*height_function)(char *data, float x, float y), char *data, float x_low, float y_low, float x_high, float y_high ) { dproj dummy[1]; draw_2d_function(height_function,data, 10,21,(char *)NULL, x_low,y_low,x_high,y_high,dummy); } void draw_1d_function( float (*height_function)(char *data, float x), char *data, int grid_size, char *xlabel, char *ylabel, float x_low, float x_high, dproj *result_dproj ) { int num_points = grid_size+1; int i; ongr on[1]; frame fr[1]; float *x_arr = am_malloc_floats(num_points); float *y_arr = am_malloc_floats(num_points); for ( i = 0 ; i < num_points ; i++ ) { float z = i / (float) (num_points-1); float x = x_low + (x_high - x_low) * z; float y = height_function(data,x); x_arr[i] = x; y_arr[i] = y; } full_screen_frame(fr); clear_ongr(on); compute_axis_limits(x_arr,num_points,&on->x_axis); compute_axis_limits(y_arr,num_points,&on->y_axis); compute_axes_details(on,fr); draw_axes(fr,on); plot_in_frame(fr,on,x_arr,y_arr,num_points,"LN"); result_dproj -> screen_x_min = on->x_axis.start.x; result_dproj -> screen_y_min = on->y_axis.start.y; result_dproj -> screen_x_max = on->x_axis.end.x; result_dproj -> screen_y_max = on->y_axis.end.y; result_dproj -> indim = 1; result_dproj -> in_min[0] = on->x_axis.lo; result_dproj -> in_max[0] = on->x_axis.hi; result_dproj -> height_min = on->y_axis.lo; result_dproj -> height_max = on->y_axis.hi; am_free_floats(x_arr,num_points); am_free_floats(y_arr,num_points); } void simple_draw_1d_function( float (*height_function)(char *data, float x), char *data, float x_low, float x_high ) { dproj dummy[1]; draw_1d_function(height_function,data, 200,"x","f(x)",x_low,x_high,dummy ); }