/* * $Id: areas.c,v 1.6 2006/07/13 05:06:37 cho Exp $ * * Cell-file area extraction * Boundary tracing algorithm * * Jean Ezell * US Army Corps of Engineers * Construction Engineering Research Lab * Modelling and Simulation Team * Champaign, IL 61820 * January 1988 * * Fixed 3/2001 by Andrea Aime * fixed area label problem with r.poly (it was creating more * labels than required if there were islands) * ************************* * The algorithm allocates COOR structures as necessary to mark endpoints of * and bends in boundaries between areas. Each COOR structure contains the * row and column coordinates of a point which is either an endpoint or a * bend. If the point represents a bend, fptr and bptr point to the adjacent * endpoint(s) or bend(s); if an endpoint, one of fptr or bptr points to the * point itself and the other to the adjacent bend or endpoint. While a * boundary is under construction, fptr of the "free" end is NULL. * The right and left fields contain the area numbers to the right and * left of the line under construction. When two lines are joined, it * will be necessary to make the area numbers agree. This is done by mapping * one of the area numbers to another. The mapping is done in such a way * that the smallest number in each equivalence class is chosen to represent * the class. We also keep track of the longest horizontal strip in each * area. The left end of this strip will be used as the position of the * label in the dlg label file. * * Global variables: * v_list pointer to allocated array of pointers to endpoints of * vertical lines currently under construction * h_ptr pointer to endpoint of horizontal line currently * under construction * col, row column and row of current position in file * buffer[] pointers to the two allocated areas which hold the * two rows currently being processed * top, bottom which row of buffer[] is "top" line from file and * which is "bottom" * scan_length length of a row from the data file * tl, tr, bl, br top left and right, bottom left and right elements * in current 2 by 2 data window; one-time calculation * prevents multiple indexing and indirection * computations * a_list pointer to allocated array of correspondences between * areas and categories * n_areas current length of a_list * area_num lowest area number available for use * a_list_new pointer to next a_list entry to be used * a_list_old pointer to a_list entry just filled * e_list pointer to allocated array which is used to form * mapping of all equivalent area numbers onto one * number from the class * n_equiv current length of e_list * * Entry points: * extract_areas driver for boundary extraction, area labelling * algorithm * alloc_bufs allocate buffers for cell file data and storage of * information obtained in the extraction process (v_list, * a_list, e_list, buffer[0], buffer[1]) * ********************************************************************/ #include #include #include #include #include #include "global.h" static int col, row, top, bottom; static double tl, tr, bl, br; static struct COOR **v_list; static struct COOR *h_ptr; static void *buffer[2]; static int scan_length; static int n_areas, area_num, n_equiv, tl_area; static struct area_table *a_list, *a_list_new, *a_list_old; static struct equiv_table *e_list; /* function prototypes */ static int update_list(int); static int end_vline(); static int end_hline(); static int start_vline(); static int start_hline(); static struct COOR *get_ptr(); static int read_next(); static int equiv_areas(int,int); static int map_area(int,int); static int add_to_list(int,int); static int assign_area(double,int); static int more_areas(); static int update_width(struct area_table *,int); static int nabors (void); #define get_raster_value(ptr, col) \ G_get_raster_value_d(G_incr_void_ptr(ptr, (col)*data_size), data_type) /* extract_areas - trace boundaries of polygons in file */ int extract_areas (int quiet) { double nullVal; row = col = top = 0; /* get started for read of first */ bottom = 1; /* line from cell file */ area_num = 0; tl_area = 0; G_set_d_null_value(&nullVal, 1); /* represents the "outside", the external null values */ assign_area(nullVal,0); if (!quiet) fprintf(stderr, _("Extracting areas ... ")); scan_length = read_next(); while (read_next()) /* read rest of file, one row at */ { /* a time */ if (!quiet) G_percent(row, n_rows, 2); for (col = 0; col < scan_length - 1; col++) { tl = get_raster_value(buffer[top], col); /* top left in window */ tr = get_raster_value(buffer[top], col+1); /* top right */ bl = get_raster_value(buffer[bottom], col); /* bottom left */ br = get_raster_value(buffer[bottom], col+1); /* bottom right */ update_list(nabors()); } if (h_ptr != NULPTR) /* if we have a loose end, */ end_hline(); /* tie it down */ row++; } if (!quiet) G_percent(row, n_rows, 2); write_area(a_list,e_list,area_num,n_equiv); G_free(a_list); G_free(e_list); return 0; } /* extract_areas */ /* update_list - maintains linked list of COOR structures which resprsent */ /* bends in and endpoints of lines separating areas in input file; */ /* compiles a list of area to category number correspondences; */ /* for pictures of what each case in the switch represents, see comments */ /* before nabors() */ static int update_list (int i) { struct COOR *new_ptr, *new_ptr1, *new_ptr2, *new_ptr3; double right, left; switch (i) { case 0: /* Vertical line - Just update width information */ /* update_width(a_list + v_list[col]->left,0);*/ tl_area = v_list[col]->left; break; case 1: /* Bottom right corner - Point in middle of new line */ /* (growing) <- ptr2 -><- ptr1 -><- ptr3 -> (growing) */ /* (?, col) (row, col) (row, ?) */ new_ptr1 = get_ptr(); /* corner point */ new_ptr2 = get_ptr(); /* downward-growing point */ new_ptr3 = get_ptr(); /* right-growing point */ new_ptr1->bptr = new_ptr2; new_ptr1->fptr = new_ptr3; new_ptr2->bptr = new_ptr3->bptr = new_ptr1; /* if(G_is_c_null_value(&tl_area)) { new_ptr1->left = new_ptr2->right = new_ptr3->left = 0; assign_area(tl,1); } else { */ new_ptr1->left = new_ptr2->right = new_ptr3->left = tl_area; new_ptr1->right = new_ptr2->left = new_ptr3->right = area_num; assign_area(br,1); update_width(a_list_old,1); v_list[col] = new_ptr2; h_ptr = new_ptr3; break; case 3: /* Bottom left corner - Add point to line already under construction */ /* (fixed) -><- original h_ptr -><- new_ptr -> (growing) */ /* (row, col) (?, col) */ tl_area = h_ptr->left; new_ptr = get_ptr(); /* downward-growing point */ h_ptr->col = col; h_ptr->fptr = new_ptr; new_ptr->bptr = h_ptr; new_ptr->left = h_ptr->left; new_ptr->right = h_ptr->right; /* update_width(a_list + new_ptr->left,3); */ v_list[col] = new_ptr; h_ptr = NULPTR; break; case 4: /* Top left corner - Join two lines already under construction */ /* (fixed) -><- original v_list -><- (fixed) */ /* (row, col) */ tl_area = v_list[col]->left; equiv_areas(h_ptr->left,v_list[col]->right); equiv_areas(h_ptr->right,v_list[col]->left); v_list[col]->row = row; /* keep downward-growing point */ v_list[col]->fptr = h_ptr->bptr; /* and join it to predecessor */ h_ptr->bptr->fptr = v_list[col]; /* of right-growing point */ G_free(h_ptr); /* right-growing point disappears */ h_ptr = NULPTR; /* turn loose of pointers */ write_boundary(v_list[col]); /* try to write line */ v_list[col] = NULPTR; /* turn loose of pointers */ break; case 5: /* Top right corner - Add point to line already under construction */ /* (fixed) -><- original v_list -><- new_ptr -> (growing) */ /* (row, col) (row, ?) */ new_ptr = get_ptr(); /* right-growing point */ v_list[col]->row = row; new_ptr->bptr = v_list[col]; new_ptr->left = v_list[col]->left; new_ptr->right = v_list[col]->right; v_list[col]->fptr = new_ptr; h_ptr = new_ptr; v_list[col] = NULPTR; break; case 6: /* T upward - End one vertical and one horizontal line */ /* Start horizontal line */ v_list[col]->node = h_ptr->node = 1; left = v_list[col]->left; right = h_ptr->right; end_vline(); end_hline(); start_hline(); h_ptr->bptr->node = 1; /* where we came from is a node */ h_ptr->left = h_ptr->bptr->left = left; h_ptr->right = h_ptr->bptr->right = right; break; case 7: /* T downward - End horizontal line */ /* Start one vertical and one horizontal line */ h_ptr->node = 1; right = h_ptr->right; left = h_ptr->left; end_hline(); start_hline(); start_vline(); h_ptr->bptr->node = v_list[col]->bptr->node = 1; h_ptr->left = h_ptr->bptr->left = left; h_ptr->right = h_ptr->bptr->right = v_list[col]->left = v_list[col]->bptr->left = area_num; assign_area(br,7); update_width(a_list_old,7); v_list[col]->right = v_list[col]->bptr->right = right; break; case 8: /* T left - End one vertical and one horizontal line */ /* Start one vertical line */ tl_area = v_list[col]->left; h_ptr->node = v_list[col]->node = 1; right = h_ptr->right; left = v_list[col]->left; end_vline(); end_hline(); start_vline(); v_list[col]->bptr->node = 1; /* where we came from is a node */ v_list[col]->left = v_list[col]->bptr->left = left; v_list[col]->right = v_list[col]->bptr->right = right; /* update_width(a_list + v_list[col]->left,8); */ break; case 9: /* T right - End one vertical line */ /* Start one vertical and one horizontal line */ v_list[col]->node = 1; right = v_list[col]->right; left = v_list[col]->left; end_vline(); start_vline(); start_hline(); v_list[col]->bptr->node = h_ptr->bptr->node = 1; h_ptr->left = h_ptr->bptr->left = left; h_ptr->right = h_ptr->bptr->right = v_list[col]->left = v_list[col]->bptr->left = area_num; assign_area(br,9); update_width(a_list_old,9); v_list[col]->right = v_list[col]->bptr->right = right; break; case 10: /* Cross - End one vertical and one horizontal line */ /* Start one vertical and one horizontal line */ v_list[col]->node = h_ptr->node = 1; left = v_list[col]->left; right = h_ptr->right; end_vline(); end_hline(); start_vline(); start_hline(); v_list[col]->bptr->node = h_ptr->bptr->node = 1; h_ptr->left = h_ptr->bptr->left = left; v_list[col]->left = v_list[col]->bptr->left = h_ptr->right = h_ptr->bptr->right = area_num; assign_area(br,10); update_width(a_list_old,10); v_list[col]->right = v_list[col]->bptr->right = right; break; } /* switch */ return 0; } /* update_list */ /* end_vline, end_hline - end vertical or horizontal line and try */ /* to write it out */ static int end_vline (void) { v_list[col]->row = row; v_list[col]->fptr = v_list[col]; write_boundary(v_list[col]); v_list[col] = NULPTR; return 0; } static int end_hline (void) { h_ptr->col = col; h_ptr->fptr = h_ptr; write_boundary(h_ptr); h_ptr = NULPTR; return 0; } /* start_vline, start_hline - begin line in vertical or horizontal */ /* direction */ static int start_vline (void) { struct COOR *new_ptr1, *new_ptr2; new_ptr1 = get_ptr(); new_ptr2 = get_ptr(); new_ptr1->fptr = new_ptr2; new_ptr2->bptr = new_ptr1->bptr = new_ptr1; new_ptr2->fptr = NULPTR; v_list[col] = new_ptr2; return 0; } static int start_hline (void) { struct COOR *new_ptr1, *new_ptr2; new_ptr1 = get_ptr(); new_ptr2 = get_ptr(); new_ptr1->bptr = new_ptr2->bptr = new_ptr1; new_ptr1->fptr = new_ptr2; new_ptr2->fptr = NULPTR; h_ptr = new_ptr2; return 0; } /* get_ptr - allocate storage for yet another COOR structure */ static struct COOR *get_ptr (void) { static struct COOR *ptr; ptr = (struct COOR *) G_malloc(sizeof (struct COOR)); ptr->row = row; ptr->col = col; ptr->fptr = ptr->bptr = NULPTR; ptr->node = ptr->left = ptr->right = 0; return(ptr); } /* nabors - check 2 x 2 matrix and return case from table below */ /* *--*--* *--*--* *--*--* *--*--* */ /* | | | | | | | | | */ /* * | * * *--* *-----* *--* * */ /* | | | | | | | | | | | */ /* *--*--* *--*--* *--*--* *--*--* */ /* 0 1 2 3 */ /* *--*--* *--*--* *--*--* *--*--* */ /* | | | | | | | | | | | */ /* *--* * * *--* *--*--* *--*--* */ /* | | | | | | | | | */ /* *--*--* *--*--* *--*--* *--*--* */ /* 4 5 6 7 */ /* *--*--* *--*--* *--*--* *--*--* */ /* | | | | | | | | | | | */ /* *--* * * *--* *--*--* * * */ /* | | | | | | | | | | | */ /* *--*--* *--*--* *--*--* *--*--* */ /* */ /* 8 9 10 11 */ static int nabors (void) { int tl_null = G_is_d_null_value(&tl); int tr_null = G_is_d_null_value(&tr); int bl_null = G_is_d_null_value(&bl); int br_null = G_is_d_null_value(&br); /* if both a and b are NULLs, thery are equal */ #define cmp(a, b) (a##_null+b##_null==1 || (a##_null+b##_null==0 && a != b)) if (cmp(tl, tr) != 0) /* 0, 4, 5, 6, 8, 9, 10 */ { if (cmp(tl, bl) != 0) /* 4, 6, 8, 10 */ { if (cmp(bl, br) != 0) /* 8, 10 */ { if (cmp(tr, br) != 0) return(10); else return(8); } else /* 4, 6 */ { if (cmp(tr, br) != 0) return(6); else return(4); } } else /* 0, 5, 9 */ { if (cmp(bl, br) != 0) /* 0, 9 */ { if (cmp(tr, br) != 0) return(9); else return(0); } else return(5); } } else /* 1, 2, 3, 7, 11 */ { if (cmp(tl, bl) != 0) /* 2, 3, 7 */ { if (cmp(bl, br) != 0) /* 3, 7 */ { if (cmp(tr, br) != 0) return(7); else return(3); } else return(2); } else /* 1, 11 */ { if (cmp(bl, br) != 0) return(1); else return(11); } } return 0; } /* read_next - read another line from input file */ static int read_next (void) { int n; top = bottom++; /* switch top and bottom, */ bottom = 1 & bottom; /* which are always 0 or 1 */ n = read_row(buffer[bottom]); return(n); } /* alloc_bufs - allocate buffers we will need for storing cell file */ /* data, pointers to extracted lines, area number information */ int alloc_areas_bufs (int size) { int i; buffer[0] = (void *) G_malloc(size * data_size); buffer[1] = (void *) G_malloc(size * data_size); v_list = (struct COOR **) G_malloc(size * sizeof(*v_list)); n_areas = n_equiv = 500; /* guess at number of areas, equivs */ a_list = (struct area_table *) G_malloc(n_areas * sizeof(struct area_table)); for (i = 0; i < n_areas; i++) { (a_list + i)->width = (a_list + i)->row = (a_list + i)->col = 0; (a_list + i)->free = 1; } a_list_new = a_list_old = a_list; e_list = (struct equiv_table *) G_malloc(n_equiv * sizeof(struct equiv_table)); for (i = 0; i < n_equiv; i++) { (e_list + i)->mapped = (e_list + i)->count = 0; (e_list + i)->ptr = NULL; } return 0; } /* equiv_areas - force two areas to be equivalent and generate */ /* mapping information */ static int equiv_areas (int a1, int a2) { int small, large, small_obj, large_obj; if(a1 == -1 || a2 == -2) return 0; if (a1 == a2) return(0); if (a1 < a2) { small = a1; large = a2; } else { small = a2; large = a1; } while (large >= n_equiv) /* make sure our equivalence tables */ more_equivs(); /* are large enough */ if ((e_list + large)->mapped) { if ((e_list + small)->mapped) /* small mapped, large mapped */ { large_obj = (e_list + large)->where; small_obj = (e_list + small)->where; if (large_obj == small_obj) /* both mapped to same place */ return(0); if (small_obj < large_obj) /* map where large goes to where */ map_area(large_obj,small_obj); /* small goes */ else /* map where small goes to where */ map_area(small_obj,large_obj); /* large goes */ } else /* small not mapped, large mapped */ { large_obj = (e_list + large)->where; if (small == large_obj) /* large already mapped to small */ return(0); if (small < large_obj) /* map where large goes to small */ map_area(large_obj,small); else /* map small to where large goes */ map_area(small,large_obj); } } else { if ((e_list + small)->mapped) /* small mapped, large not mapped */ map_area(large,(e_list + small)->where); else /* small not mapped, large not mapped */ map_area(large,small); } return (0); } /* map_area - establish a mapping from one area to another */ /* the area number mapping looks like the following: */ /* area numbers index into e_list to get equiv_table structures */ /* for each of these structures, */ /* .mapped indicates whether area number is mapped */ /* if area number is mapped */ /* .where tells to what other area number */ /* if area number is not mapped */ /* .count tells how many areas are mapped to this one */ /* if count != 0 */ /* .ptr gives a pointer to the list of area numbers */ static int map_area ( int x, int y /* map x to y */ ) { int n, i, *p; (e_list + x)->mapped = 1; (e_list + x)->where = y; if ((a_list + x)->width > (a_list + y)->width) { (a_list + y)->width = (a_list + x)->width; (a_list + y)->row = (a_list + x)->row; (a_list + y)->col = (a_list + x)->col; } if (add_to_list(x,y)) /* if x is not already in y's list */ { n = (e_list + x)->count; p = (e_list + x)->ptr; for (i = 0; i < n; i++) /* map everything that is currently */ { /* mapped onto x onto y; because */ (e_list + *p)->where = y; /* of this reshuffle, only one */ add_to_list(*p++,y); /* level of mapping is ever needed */ } } return 0; } /* add_to_list - add another area number to an equivalence list */ static int add_to_list (int x, int y) { int n, i; struct equiv_table *e_list_y; e_list_y = e_list + y; n = e_list_y->count; if (n == 0) /* first time through--start list */ { e_list_y->length = 20; /* initial guess at storage needed */ e_list_y->ptr = (int *) G_malloc(e_list_y->length * sizeof(int)); *(e_list_y->ptr) = x; (e_list_y->count)++; } else /* list already started */ { for (i = 0; i < n; i++) /* check whether x is already */ { /* in list */ if (*(e_list_y->ptr + i) == x) return(0); /* if so, we don't need to add it */ } if (n + 1 >= e_list_y->length) /* add more space for storage */ { /* if necessary */ e_list_y->length += 10; e_list_y->ptr = (int *) G_realloc(e_list_y->ptr,e_list_y->length * sizeof(int)); } *(e_list_y->ptr + n) = x; /* add x to list */ (e_list_y->count)++; } return(1); /* indicate addition made */ } /* assign_area - make current area number correspond to the passed */ /* category number and allocate more space to store areas if necessary */ static int assign_area (double cat, int kase) { a_list_new->free = 0; a_list_new->cat = cat; area_num++; if (area_num >= n_areas) more_areas(); a_list_old = a_list + area_num - 1; a_list_new = a_list + area_num; return 0; } /* more_areas - allocate larger space to store area correspondences */ static int more_areas (void) { int old_n, i; old_n = n_areas; n_areas += 250; a_list = (struct area_table *) G_realloc(a_list,n_areas * sizeof(struct area_table)); for (i = old_n; i < n_areas; i++) { (a_list + i)->width = -1; (a_list + i)->free = 1; } return 0; } /* more_equivs - allocate more space to construct equivalence information */ int more_equivs (void) { int old_n, i; old_n = n_equiv; n_equiv += 250; e_list = (struct equiv_table *) G_realloc(e_list,n_equiv * sizeof(struct equiv_table)); for (i = old_n; i < n_equiv; i++) { (e_list + i)->mapped = (e_list + i)->count = 0; (e_list + i)->ptr = NULL; } return 0; } /* update_width - update position of longest horizontal strip in an area */ static int update_width (struct area_table *ptr, int kase) { int w, j, a; struct equiv_table *ep; a = (ptr - a_list); for (j = col + 1, w = 0; j < scan_length && get_raster_value(buffer[bottom], j) == br; j++, w++); if(a == 0) G_message(_("Area 0, %d \t%d \t%d \t%d \t%d"), kase, row, col, ptr->width, w); if (a < n_equiv) { ep = e_list + a; if (ep->mapped) ptr = a_list + ep->where; } if (w > ptr->width) { ptr->width = w; ptr->row = row; ptr->col = col; } return 0; }