/*****************************************************************************/ /* */ /* F I S T : Fast, Industrial-Strength Triangulation */ /* */ /*****************************************************************************/ /* */ /* (C) Martin Held */ /* (C) Universitaet Salzburg, Salzburg, Austria */ /* */ /* This code is not in the public domain. All rights reserved! Please make */ /* sure to read the full copyright statement contained in api_functions.cpp. */ /* */ /*****************************************************************************/ /* get standard i/o library */ #include #include #include #include #include /* get my header file */ #include "fpkernel.h" #include "martin.h" #include "defs.h" #include "list.h" #include "header.h" #include "basic.h" #include "numerics.h" #include "bv_tree.h" #include "grid.h" #define CheckCell(all, data, grid, tmp, resettmp, k, bb, flag, i5) \ { \ flag = false; \ \ tmp.ind_seg = grid->grid[k]; \ do { \ tmp.seg_i = grid->seg_list[tmp.ind_seg].seg; \ tmp.seg_addr = &(grid->segments[tmp.seg_i]); \ if (!(tmp.seg_addr->checked)) { \ AddToSet(grid, tmp.seg_i); \ if (BBox_Overlap(bb, *tmp.seg_addr)) { \ if (SegIntersection(all,bb.imin, bb.imax, \ tmp.seg_addr->imin, tmp.seg_addr->imax, i5)) { \ ResetSet(grid, resettmp); \ flag = true; \ tmp.ind_seg = NIL; \ } \ else { \ tmp.ind_seg = grid->seg_list[tmp.ind_seg].next; \ } \ } \ else { \ tmp.ind_seg = grid->seg_list[tmp.ind_seg].next; \ } \ } \ else { \ tmp.ind_seg = grid->seg_list[tmp.ind_seg].next; \ } \ } while (tmp.ind_seg != NIL); \ } //void CheckCell(all, data, grid, tmp, tmp2, k, bb, flag, i5) //{ // flag = false; // // tmp.ind_seg = grid->grid[k]; // do { // tmp.seg_i = grid->seg_list[tmp.ind_seg].seg; // tmp.seg_addr = &(grid->segments[tmp.seg_i]); // if (!(tmp.seg_addr->checked)) { // AddToSet(grid, tmp.seg_i); // if (BBox_Overlap(bb, *tmp.seg_addr)) { // if (SegIntersection(all,bb.imin, bb.imax, // tmp.seg_addr->imin, tmp.seg_addr->imax, i5)) { // ResetSet(grid,tmp2); // flag = true; // tmp.ind_seg = NIL; // } // else { // tmp.ind_seg = grid->seg_list[tmp.ind_seg].next; // } // } // else { // tmp.ind_seg = grid->seg_list[tmp.ind_seg].next; // } // } // else { // tmp.ind_seg = grid->seg_list[tmp.ind_seg].next; // } // } while (tmp.ind_seg != NIL); //} // /* */ /* function prototypes of functions provided in this file */ /* */ void InitGridDefaults(griddef *grid); boolean GridIntersectionExists(global_struct *all, bounding_box bb, int i1, int i2, list_ind ind5, int i5); void FreeGrid(griddef *grid); void BuildGrid(global_struct *all, int loop_min, int loop_max); void InsertSegmentIntoGrid(datadef *data, griddef *grid, bounding_box bb); machine_double TopQualityGrid(global_struct *all, int i1, int i2); void BuildPntsGrid(global_struct *all, int loop_min); void BuildBuckets(global_struct *all, list_ind ind); boolean BucketIntersectionExists(global_struct *all, int i1, list_ind ind1, int i2, int i3, bounding_box bb, int *i_close, boolean check_close); void SetReflexNumber(global_struct *all,int number); #ifdef GRAZING void ScanBridgeBuckets(global_struct *all, int i0, int grid_offset, distance *distances, int *num_dist, boolean *grid_exhausted, double x_start); #else void ScanBridgeBuckets(global_struct *all, int i0, int grid_offset, distance *distances, int *num_dist, boolean *grid_exhausted); #endif void ResetGrid(griddef *grid); void InitBuckets(global_struct *all); void InitGridDefaults(griddef *grid) { grid->no_grid_yet = true; grid->grid = (seg_ind*)NULL; grid->max_num_grid = 0; grid->grid2 = (pnt_ind*)NULL; grid->max_num_grid2 = 0; grid->grid3 = (pnt_ind*)NULL; grid->max_num_grid3 = 0; grid->raster_x = (double*)NULL; grid->raster_y = (double*)NULL; grid->raster_x2 = (double*)NULL; grid->raster_y2 = (double*)NULL; grid->raster_x3 = (double*)NULL; grid->raster_y3 = (double*)NULL; grid->max_num_raster_x = 0; grid->max_num_raster_y = 0; grid->max_num_raster_x2 = 0; grid->max_num_raster_y2 = 0; grid->max_num_raster_x3 = 0; grid->max_num_raster_y3 = 0; grid->segments = (segment*)NULL; grid->num_segments = 0; grid->max_num_segments = 0; grid->seg_list = (segment_node*)NULL; grid->num_seg_list = 0; grid->max_num_seg_list = 0; grid->pnt_list = (pnt_node*)NULL; grid->num_pnt_list = 0; grid->max_num_pnt_list = 0; grid->vtx_list = (pnt_node*)NULL; grid->num_vtx_list = 0; grid->max_num_vtx_list = 0; grid->max_num_set = 0; grid->num_set = 0; grid->num_reflex = 0; grid->num_original_reflex = 0; grid->buckets_initialized = false; grid->set = (int*)NULL; } void SetReflexNumber(global_struct *all, int number) { griddef *grid = &all->c_grid; grid->num_reflex = number; grid->num_original_reflex = (int) (number / 10.0); if (grid->num_original_reflex < 5000) grid->num_original_reflex = 0; return; } void FreeGrid(griddef *grid) { FreeMemory(grid->memptr, (void**) &grid->grid, "grid:grid"); CORE_FreeMemory(grid->memptr, &grid->raster_x, "grid:raster_x"); CORE_FreeMemory(grid->memptr, &grid->raster_y, "grid:raster_y"); FreeMemory(grid->memptr, (void**) &grid->grid2, "grid:grid2"); CORE_FreeMemory(grid->memptr, &grid->raster_x2, "grid:raster_x2"); CORE_FreeMemory(grid->memptr, &grid->raster_y2, "grid:raster_y2"); FreeMemory(grid->memptr, (void**) &grid->grid3, "grid:grid3"); CORE_FreeMemory(grid->memptr, &grid->raster_x3, "grid:raster_x3"); CORE_FreeMemory(grid->memptr, &grid->raster_y3, "grid:raster_y3"); FreeMemory(grid->memptr, (void**) &grid->seg_list, "grid:seg_list"); FreeMemory(grid->memptr, (void**) &grid->pnt_list, "grid:pnt_list"); FreeMemory(grid->memptr, (void**) &grid->vtx_list, "grid:vtx_list"); CORE_FreeMemory(grid->memptr, &grid->segments, "grid:segments"); FreeMemory(grid->memptr, (void**) &grid->set, "grid:set"); grid->N = 0; grid->N_x = 0; grid->N_y = 0; grid->N2 = 0; grid->N_x2 = 0; grid->N_y2 = 0; grid->N3 = 0; grid->N_x3 = 0; grid->N_y3 = 0; grid->num_seg_list = 0; grid->max_num_seg_list = 0; grid->num_pnt_list = 0; grid->max_num_pnt_list = 0; grid->num_vtx_list = 0; grid->max_num_vtx_list = 0; grid->num_segments = 0; grid->max_num_segments = 0; grid->max_num_grid = 0; grid->max_num_raster_x = 0; grid->max_num_raster_y = 0; grid->max_num_grid2 = 0; grid->max_num_raster_x2 = 0; grid->max_num_raster_y2 = 0; grid->max_num_grid3 = 0; grid->max_num_raster_x3 = 0; grid->max_num_raster_y3 = 0; grid->max_num_set = 0; grid->num_set = 0; grid->num_reflex = 0; return; } void InitGrid(griddef *grid) { int k; assert(grid->N > 0); assert(grid->N_x * grid->N_y <= grid->N); if (grid->N > grid->max_num_grid) { grid->max_num_grid = grid->N; grid->grid = (seg_ind*) ReallocateArray(grid->memptr, grid->grid, grid->max_num_grid, sizeof(seg_ind), "grid:grid"); } for (k = 0; k < grid->N; ++k) { grid->grid[k] = NIL; } return; } boolean InGrid(griddef *grid, int cell) { return ((0 <= cell) && (cell < grid->max_num_grid)); } void InitRaster(griddef *grid) { int i, j; if (grid->N_x >= grid->max_num_raster_x) { grid->raster_x = CORE_ReallocateArray(grid->memptr, grid->raster_x, grid->max_num_raster_x, grid->N_x + 1, double, "grid:raster_x"); grid->max_num_raster_x = grid->N_x + 1; } if (grid->N_y >= grid->max_num_raster_y) { grid->raster_y = CORE_ReallocateArray(grid->memptr, grid->raster_y, grid->max_num_raster_y, grid->N_y + 1, double, "grid:raster_y"); grid->max_num_raster_y = grid->N_y + 1; } for (i = 0; i < grid->N_x; ++i) grid->raster_x[i] = grid->grid_min_x + i * grid->grid_x; grid->raster_x[grid->N_x] = grid->grid_max_x; for (j = 0; j < grid->N_y; ++j) grid->raster_y[j] = grid->grid_min_y + j * grid->grid_y; grid->raster_y[grid->N_y] = grid->grid_max_y; return; } void InitGrid2(griddef *grid) { int k; assert(grid->N2 > 0); assert(grid->N_x2 * grid->N_y2 <= grid->N2); if (grid->N2 > grid->max_num_grid2) { grid->max_num_grid2 = grid->N2; grid->grid2 = (pnt_ind*) ReallocateArray(grid->memptr, grid->grid2, grid->max_num_grid2, sizeof(pnt_ind), "grid:grid2"); } for (k = 0; k < grid->N2; ++k) { grid->grid2[k] = NIL; } return; } boolean InGrid2(griddef *grid, int cell) { return ((0 <= cell) && (cell < grid->max_num_grid2)); } void InitRaster2(griddef *grid) { int i, j; if (grid->N_x2 >= grid->max_num_raster_x2) { grid->raster_x2 = CORE_ReallocateArray(grid->memptr, grid->raster_x2, grid->max_num_raster_x2, grid->N_x2 + 1, double, "grid:raster_x2"); grid->max_num_raster_x2 = grid->N_x2 + 1; } if (grid->N_y2 >= grid->max_num_raster_y2) { grid->raster_y2 = CORE_ReallocateArray(grid->memptr, grid->raster_y2, grid->max_num_raster_y2, grid->N_y2 + 1, double, "grid:raster_y2"); grid->max_num_raster_y2 = grid->N_y2 + 1; } for (i = 0; i < grid->N_x2; ++i) grid->raster_x2[i] = grid->grid_min_x2 + i * grid->grid_x2; grid->raster_x2[grid->N_x2] = grid->grid_max_x2; for (j = 0; j < grid->N_y2; ++j) grid->raster_y2[j] = grid->grid_min_y2 + j * grid->grid_y2; grid->raster_y2[grid->N_y2] = grid->grid_max_y2; grid->step_x = grid->grid_x2 / C_4_0; grid->step_y = grid->grid_y2 / C_4_0; return; } void InitGrid3(griddef *grid) { int k; assert(grid->N3 > 0); assert(grid->N_x3 * grid->N_y3 <= grid->N3); if (grid->N3 > grid->max_num_grid3) { grid->max_num_grid3 = grid->N3; grid->grid3 = (pnt_ind*) ReallocateArray(grid->memptr, grid->grid3, grid->max_num_grid3, sizeof(pnt_ind), "grid:grid3"); } for (k = 0; k < grid->N3; ++k) { grid->grid3[k] = NIL; } return; } int CountEmptyCells(griddef *grid) { int k; int sum = 0; for (k = 0; k < grid->N3; ++k) { if (grid->grid3[k] == NIL) ++sum; } return sum; } boolean InGrid3(griddef *grid, int cell) { return ((0 <= cell) && (cell < grid->max_num_grid3)); } void InitRaster3(griddef *grid) { int i, j; if (grid->N_x3 >= grid->max_num_raster_x3) { grid->raster_x3 = CORE_ReallocateArray(grid->memptr, grid->raster_x3, grid->max_num_raster_x3, grid->N_x3 + 1, double, "grid:raster_x3"); grid->max_num_raster_x3 = grid->N_x3 + 1; } if (grid->N_y3 >= grid->max_num_raster_y3) { grid->raster_y3 = CORE_ReallocateArray(grid->memptr, grid->raster_y3, grid->max_num_raster_y3, grid->N_y3 + 1, double, "grid:raster_y3"); grid->max_num_raster_y3 = grid->N_y3 + 1; } for (i = 0; i < grid->N_x3; ++i) grid->raster_x3[i] = grid->grid_min_x3 + i * grid->grid_x3; grid->raster_x3[grid->N_x3] = grid->grid_max_x3; for (j = 0; j < grid->N_y3; ++j) grid->raster_y3[j] = grid->grid_min_y3 + j * grid->grid_y3; grid->raster_y3[grid->N_y3] = grid->grid_max_y3; return; } void InitSegs(griddef *grid, int number) { if (number > grid->max_num_segments) { grid->segments = CORE_ReallocateArray(grid->memptr, grid->segments, grid->max_num_segments, number, segment, "grid:segments"); grid->max_num_segments = number; } grid->num_segments = 0; return; } boolean InSegments(griddef *grid, seg_ind seg) { return ((0 <= seg) && (seg < grid->num_segments)); } boolean InVtxList(griddef *grid, int vtx) { return ((0 <= vtx) && (vtx < grid->num_vtx_list)); } void InitSegList(griddef *grid, int number) { if (number > grid->max_num_seg_list) { grid->max_num_seg_list = number; grid->seg_list = (segment_node*) ReallocateArray(grid->memptr, grid->seg_list, grid->max_num_seg_list, sizeof(segment_node), "grid:seg_list"); } grid->num_seg_list = 0; return; } void InitPntList(griddef *grid, int number) { if (number > grid->max_num_pnt_list) { grid->max_num_pnt_list = number; grid->pnt_list = (pnt_node*) ReallocateArray(grid->memptr, grid->pnt_list, grid->max_num_pnt_list, sizeof(pnt_node), "grid:pnt_list"); } grid->num_pnt_list = 0; return; } void InitVtxList(griddef *grid, int number) { if (number > grid->max_num_vtx_list) { grid->max_num_vtx_list = number; grid->vtx_list = (pnt_node*) ReallocateArray(grid->memptr, grid->vtx_list, grid->max_num_vtx_list, sizeof(pnt_node), "grid:vtx_list"); } grid->num_vtx_list = 0; return; } void RasterizeSegment(datadef *data, griddef *grid, seg_ind ind) { int k; int cx, cy, cx1, cx2, cy1, cy2, cxx1, cxx2, cyy1, cyy2; double x1, x2, y1, y2, det; point corner, p, q; p = data->points[grid->segments[ind].imin]; q = data->points[grid->segments[ind].imax]; x1 = p.x; y1 = p.y; x2 = q.x; y2 = q.y; VectorSub2D(p, q, p); /* */ /* determine the cells which contain the start and the end of the segment */ /* */ DetermineCell(grid, x1, y1, cx1, cy1, cxx1, cyy1); if (cxx1 < cx1) { Convert(grid, cxx1, cy1, k); InsertAfterSeg(grid, k, ind); } DetermineCell(grid, x2, y2, cx2, cy2, cxx2, cyy2); if (cxx2 > cx2) { Convert(grid, cxx2, cy2, k); InsertAfterSeg(grid, k, ind); } if (y2 >= y1) { /* */ /* upwards */ /* */ if (cyy1 < cy1) { Convert(grid, cx1, cyy1, k); InsertAfterSeg(grid, k, ind); if (cxx1 < cx1) { Convert(grid, cxx1, cyy1, k); InsertAfterSeg(grid, k, ind); } } if (cyy2 > cy2) { Convert(grid, cx2, cyy2, k); InsertAfterSeg(grid, k, ind); if (cxx2 > cx2) { Convert(grid, cxx2, cyy2, k); InsertAfterSeg(grid, k, ind); } } /* */ /* march through the grid from the start to the end */ /* */ cx = cx1 + 1; cy = cy1 + 1; corner.x = grid->raster_x[cx]; corner.y = grid->raster_y[cy]; while ((cx1 < cx2) && (cy1 < cy2)) { Convert(grid, cx1, cy1, k); InsertAfterSeg(grid, k, ind); det = Deter2D(p, q, corner); if (det > C_0_0) { if (eq(det, EPS)) { Convert(grid, cx1, cy, k); InsertAfterSeg(grid, k, ind); } cx1 = cx; ++cx; corner.x = grid->raster_x[cx]; } else { if (eq(det, EPS)) { Convert(grid, cx, cy1, k); InsertAfterSeg(grid, k, ind); } cy1 = cy; ++cy; corner.y = grid->raster_y[cy]; } } if (cx1 == cx2) { for (cy = cy1; cy <= cy2; ++cy) { Convert(grid, cx1, cy, k); InsertAfterSeg(grid, k, ind); } } else if (cy1 == cy2) { for (cx = cx1; cx <= cx2; ++cx) { Convert(grid, cx, cy1, k); InsertAfterSeg(grid, k, ind); } } } else { /* */ /* downwards */ /* */ if (cyy1 > cy1) { Convert(grid, cx1, cyy1, k); InsertAfterSeg(grid, k, ind); if (cxx1 < cx1) { Convert(grid, cxx1, cyy1, k); InsertAfterSeg(grid, k, ind); } } if (cyy2 < cy2) { Convert(grid, cx2, cyy2, k); InsertAfterSeg(grid, k, ind); if (cxx2 > cx2) { Convert(grid, cxx2, cyy2, k); InsertAfterSeg(grid, k, ind); } } /* */ /* march through the grid from the start to the end */ /* */ cx = cx1 + 1; cy = cy1 - 1; corner.x = grid->raster_x[cx]; corner.y = grid->raster_y[cy1]; while ((cx1 < cx2) && (cy1 > cy2)) { Convert(grid, cx1, cy1, k); InsertAfterSeg(grid, k, ind); det = Deter2D(p, q, corner); if (det > C_0_0) { if (eq(det, EPS)) { Convert(grid, cx, cy1, k); InsertAfterSeg(grid, k, ind); } cy1 = cy; --cy; corner.y = grid->raster_y[cy1]; } else { if (eq(det, EPS)) { Convert(grid, cx1, cy, k); InsertAfterSeg(grid, k, ind); } cx1 = cx; ++cx; corner.x = grid->raster_x[cx]; } } if (cx1 == cx2) { for (cy = cy1; cy >= cy2; --cy) { Convert(grid, cx1, cy, k); InsertAfterSeg(grid, k, ind); } } else if (cy1 == cy2) { for (cx = cx1; cx <= cx2; ++cx) { Convert(grid, cx, cy1, k); InsertAfterSeg(grid, k, ind); } } } return; } /* */ /* this function inserts the vertices of a simple polygon into the grid. the */ /* polygon is referenced by loops[loop_min]. */ /* */ void BuildPntsGrid(global_struct *all, int loop_min) { griddef *grid = &all->c_grid; listdef *list = &all->c_list; datadef *data = &all->c_data; list_ind ind, ind1; int k, i1, cx1, cy1, cxx1, cyy1; double x1, y1; /* */ /* initialization */ /* */ if (!grid->buckets_initialized) InitBuckets(all); grid->grid_min_x2 = grid->grid_min_x3; grid->grid_max_x2 = grid->grid_max_x3; grid->grid_min_y2 = grid->grid_min_y3; grid->grid_max_y2 = grid->grid_max_y3; grid->N_x2 = grid->N_x3; grid->N_y2 = grid->N_y3; if (grid->N_x2 < 1) grid->N_x2 = 1; if (grid->N_y2 < 1) grid->N_y2 = 1; grid->N2 = grid->N_x2 * grid->N_y2; grid->grid_x2 = grid->delta_x / grid->N_x2; grid->grid_y2 = grid->delta_y / grid->N_y2; InitGrid2(grid); InitRaster2(grid); InitPntList(grid, (int) (CD_2_5 * data->num_pnts)); /* */ /* insert the vertices */ /* */ ind = list->loops[loop_min]; ind1 = ind; i1 = GetIndex(list,ind1); do { /* */ /* store the vertex in the grid */ /* */ x1 = data->points[i1].x; y1 = data->points[i1].y; DetermineCell2(grid, x1, y1, cx1, cy1, cxx1, cyy1); Convert2(grid, cx1, cy1, k); InsertAfterPnt(list, grid, data, k, i1); if (cxx1 != cx1) { Convert2(grid, cxx1, cy1, k); InsertAfterPnt(list, grid, data, k, i1); if (cyy1 != cy1) { Convert2(grid, cxx1, cyy1, k); InsertAfterPnt(list, grid, data, k, i1); Convert2(grid, cx1, cyy1, k); InsertAfterPnt(list, grid, data, k, i1); } } else if (cyy1 != cy1) { Convert2(grid, cx1, cyy1, k); InsertAfterPnt(list, grid, data, k, i1); } ind1 = GetNextNode(list, ind1); i1 = GetIndex(list, ind1); } while (ind1 != ind); return; } /* */ /* this function inserts the boundary segments of the loops referenced by */ /* loops[loop_min,...loop_max-1] into a grid */ /* */ void BuildGrid(global_struct *all, int loop_min, int loop_max) { griddef *grid = &all->c_grid; listdef *list = &all->c_list; datadef *data = &all->c_data; list_ind ind, ind2; int i, i1, i2; machine_double n_seg, a, b, fudge_x, fudge_y, md_delta_x, md_delta_y; gridtmp tmp; assert(loop_min >= 0); assert(loop_max <= list->num_loops); InitSegs(grid, data->num_pnts); /* */ /* extract the boundary segments from the boundary loops */ /* */ for (i = loop_min; i < loop_max; ++i) { ind = list->loops[i]; ind2 = ind; i1 = GetIndex(list,ind2); do { ind2 = GetNextNode(list,ind2); i2 = GetIndex(list,ind2); /* */ /* store the segment */ /* */ StoreSegment(data, grid, tmp, i1, i2); i1 = i2; } while (ind2 != ind); } /* */ /* determine the grid size and location */ /* */ md_delta_x = (TO_MDOUBLE(all->bb_max.x) - TO_MDOUBLE(all->bb_min.x)) * CD_0_01; md_delta_y = (TO_MDOUBLE(all->bb_max.y) - TO_MDOUBLE(all->bb_min.y)) * CD_0_01; fudge_x = Max(Abs(TO_MDOUBLE(all->bb_max.x)), Abs(TO_MDOUBLE(all->bb_min.x))) * FUDGE; fudge_y = Max(Abs(TO_MDOUBLE(all->bb_max.y)), Abs(TO_MDOUBLE(all->bb_min.y))) * FUDGE; fudge_x = Max(fudge_x, FUDGE); fudge_y = Max(fudge_y, FUDGE); if (md_delta_x < fudge_x) md_delta_x = fudge_x; if (md_delta_y < fudge_y) md_delta_y = fudge_y; grid->grid_min_x = all->bb_min.x - md_delta_x; grid->grid_max_x = all->bb_max.x + md_delta_x; grid->grid_min_y = all->bb_min.y - md_delta_y; grid->grid_max_y = all->bb_max.y + md_delta_y; grid->delta_x = grid->grid_max_x - grid->grid_min_x; grid->delta_y = grid->grid_max_y - grid->grid_min_y; md_delta_x = TO_MDOUBLE(grid->grid_max_x) - TO_MDOUBLE(grid->grid_min_x); md_delta_y = TO_MDOUBLE(grid->grid_max_y) - TO_MDOUBLE(grid->grid_min_y); assert(grid->delta_x > ZERO); assert(grid->delta_y > ZERO); /* */ /* initialize the grid raster; we adapt the grid somewhat according to */ /* the aspect ratio of the bounding box of the segments. roughly, the */ /* number of cells equals the number of segments. */ /* */ n_seg = sqrt((machine_double)grid->num_segments); a = sqrt(md_delta_x / md_delta_y); if (a < CD_0_001) { a = CD_0_001; b = CD_1000_0; } else if (a > CD_1000_0) { a = CD_1000_0; b = CD_0_001; } else { b = CD_1_0 / a; } grid->N_x = (int)(N_GRID * n_seg * a); grid->N_y = (int)(N_GRID * n_seg * b); if (grid->N_x < 1) grid->N_x = 1; if (grid->N_y < 1) grid->N_y = 1; grid->N = grid->N_x * grid->N_y; grid->grid_x = grid->delta_x / grid->N_x; grid->grid_y = grid->delta_y / grid->N_y; InitGrid(grid); InitRaster(grid); InitSegList(grid, 5 * grid->num_segments); /* */ /* insert the segments into the grid */ /* */ for (i = 0; i < grid->num_segments; ++i) { RasterizeSegment(data,grid,i); } grid->no_grid_yet = false; return; } void InsertSegmentIntoGrid(datadef *data, griddef *grid, bounding_box bb) { assert(!grid->no_grid_yet); gridtmp tmp; if (grid->num_segments >= grid->max_num_segments) { grid->max_num_segments += BLOCK_SIZE; grid->segments = CORE_ReallocateArray(grid->memptr, grid->segments, grid->max_num_segments - BLOCK_SIZE, grid->max_num_segments, segment, "grid:segments"); } tmp.seg_addr = &(grid->segments[grid->num_segments]); BBox_Copy(bb, *tmp.seg_addr); tmp.seg_addr->checked = false; ++grid->num_segments; RasterizeSegment(data,grid,grid->num_segments - 1); return; } boolean GridIntersectionExists(global_struct *all, bounding_box bb, int i1, int i2, list_ind ind5, int i5) { griddef *grid = &all->c_grid; datadef *data = &all->c_data; int k; int cx, cy, cx1, cx2, cy1, cy2; double x1, x2, y1, y2, det; point corner, p, q; boolean intersection; gridtmp tmp, resettmp; assert(!grid->no_grid_yet); grid->ident_cntr = 0; p = data->points[bb.imin]; q = data->points[bb.imax]; x1 = p.x; y1 = p.y; x2 = q.x; y2 = q.y; VectorSub2D(p, q, p); /* */ /* determine the cells which contain the start and the end of the segment */ /* */ DetermineCellStrict(grid, x1, y1, cx1, cy1); DetermineCellStrict(grid, x2, y2, cx2, cy2); if (y2 >= y1) { /* */ /* upwards */ /* */ /* march through the grid from the start to the end */ /* */ cx = cx1 + 1; cy = cy1 + 1; corner.x = grid->raster_x[cx]; corner.y = grid->raster_y[cy]; while ((cx1 < cx2) && (cy1 < cy2)) { Convert(grid, cx1, cy1, k); if (grid->grid[k] != NIL) { CheckCell(all, data, grid, tmp, resettmp, k, bb, intersection, i5); if (intersection) return true; } det = Deter2D(p, q, corner); if (det > C_0_0) { cx1 = cx; ++cx; corner.x = grid->raster_x[cx]; } else { cy1 = cy; ++cy; corner.y = grid->raster_y[cy]; } } if (cx1 == cx2) { for (cy = cy1; cy <= cy2; ++cy) { Convert(grid, cx1, cy, k); if (grid->grid[k] != NIL) { CheckCell(all, data, grid, tmp, resettmp, k, bb, intersection, i5); if (intersection) return true; } } } else if (cy1 == cy2) { for (cx = cx1; cx <= cx2; ++cx) { Convert(grid, cx, cy1, k); if (grid->grid[k] != NIL) { CheckCell(all, data, grid, tmp, resettmp, k, bb, intersection, i5); if (intersection) return true; } } } } else { /* */ /* downwards */ /* */ /* march through the grid from the start to the end */ /* */ cx = cx1 + 1; cy = cy1 - 1; corner.x = grid->raster_x[cx]; corner.y = grid->raster_y[cy1]; while ((cx1 < cx2) && (cy1 > cy2)) { Convert(grid, cx1, cy1, k); if (grid->grid[k] != NIL) { CheckCell(all, data, grid, tmp, resettmp, k, bb, intersection, i5); if (intersection) return true; } det = Deter2D(p, q, corner); if (det > C_0_0) { cy1 = cy; --cy; corner.y = grid->raster_y[cy1]; } else { cx1 = cx; ++cx; corner.x = grid->raster_x[cx]; } } if (cx1 == cx2) { for (cy = cy1; cy >= cy2; --cy) { Convert(grid, cx1, cy, k); if (grid->grid[k] != NIL) { CheckCell(all, data, grid, tmp, resettmp, k, bb, intersection, i5); if (intersection) return true; } } } else if (cy1 == cy2) { for (cx = cx1; cx <= cx2; ++cx) { Convert(grid, cx, cy1, k); if (grid->grid[k] != NIL) { CheckCell(all, data, grid, tmp, resettmp, k, bb, intersection, i5); if (intersection) return true; } } } } ResetSet(grid, resettmp); /* */ /* oops! this segment shares one endpoint with at least four other */ /* boundary segments! oh well, yet another degenerate situation... */ /* */ if (grid->ident_cntr >= 4) { if (CheckBottleNeck(all, i5, i1, i2, ind5)) return true; else return false; } return false; } void EvaluateCellRatio(global_struct *all, int k, int i1, int i2, machine_double a, machine_double *ratio) { griddef *grid = &all->c_grid; datadef *data = &all->c_data; tmp_data_def tmp_data; pnt_ind ind_pnt; int i3; double area; machine_double b, c; machine_double this_ratio; point p, q, r; md_point side; assert(InPointsList(data, i1)); assert(InPointsList(data, i2)); p = data->points[i1]; q = data->points[i2]; ind_pnt = grid->grid2[k]; do { i3 = grid->pnt_list[ind_pnt].pnt; assert(InPointsList(data, i3)); if ((i3 != i1) && (i3 != i2)) { StableDet2D(data, tmp_data, i1, i2, i3, area); if (area >= C_0_0) { r = data->points[i3]; PntPntDistSqd(p, r, side, b); PntPntDistSqd(r, q, side, c); ComputeRatio(a, b, c, this_ratio); if (this_ratio > *ratio) *ratio = this_ratio; } } ind_pnt = grid->pnt_list[ind_pnt].next; } while (ind_pnt != NIL); return; } machine_double TopQualityGrid(global_struct *all, int i1, int i2) { griddef *grid = &all->c_grid; datadef *data = &all->c_data; int k; int cx, cy, cx1, cx2, cy1, cy2, cxx1, cxx2, cyy1, cyy2; double x1, x2, y1, y2, det; machine_double ratio = 0.0, a; point p,q, corner, u, v; if (i1 < i2) { p = data->points[i1]; q = data->points[i2]; x1 = p.x; y1 = p.y; x2 = q.x; y2 = q.y; VectorSub2D(p, q, u); v = q; } else if (i1 > i2) { p = data->points[i1]; q = data->points[i2]; x1 = q.x; y1 = q.y; x2 = p.x; y2 = p.y; VectorSub2D(q, p, u); v = p; } else { return CD_0_0; } a = DotProduct2DTMD(u, u); /* */ /* determine the cells which contain the start and the end of the segment */ /* */ DetermineCell2(grid, x1, y1, cx1, cy1, cxx1, cyy1); DetermineCell2(grid, x2, y2, cx2, cy2, cxx2, cyy2); if (y2 >= y1) { /* */ /* upwards */ /* */ /* march through the grid from the start to the end */ /* */ cx = cx1 + 1; cy = cy1 + 1; corner.x = grid->raster_x2[cx]; corner.y = grid->raster_y2[cy]; while ((cx1 < cx2) && (cy1 < cy2)) { Convert2(grid, cx1, cy1, k); if (grid->grid2[k] != NIL) { EvaluateCellRatio(all, k, i1, i2, a, &ratio); } det = Deter2D(u, v, corner); if (det > C_0_0) { cx1 = cx; ++cx; corner.x = grid->raster_x2[cx]; } else { cy1 = cy; ++cy; corner.y = grid->raster_y2[cy]; } } if (cx1 == cx2) { for (cy = cy1; cy <= cy2; ++cy) { Convert2(grid, cx1, cy, k); if (grid->grid2[k] != NIL) { EvaluateCellRatio(all, k, i1, i2, a, &ratio); } } } else if (cy1 == cy2) { for (cx = cx1; cx <= cx2; ++cx) { Convert2(grid, cx, cy1, k); if (grid->grid2[k] != NIL) { EvaluateCellRatio(all, k, i1, i2, a, &ratio); } } } } else { /* */ /* downwards */ /* */ /* march through the grid from the start to the end */ /* */ cx = cx1 + 1; cy = cy1 - 1; corner.x = grid->raster_x2[cx]; corner.y = grid->raster_y2[cy1]; while ((cx1 < cx2) && (cy1 > cy2)) { Convert2(grid, cx1, cy1, k); if (grid->grid2[k] != NIL) { EvaluateCellRatio(all, k, i1, i2, a, &ratio); } det = Deter2D(u, v, corner); if (det > C_0_0) { cy1 = cy; --cy; corner.y = grid->raster_y2[cy1]; } else { cx1 = cx; ++cx; corner.x = grid->raster_x2[cx]; } } if (cx1 == cx2) { for (cy = cy1; cy >= cy2; --cy) { Convert2(grid, cx1, cy, k); if (grid->grid2[k] != NIL) { EvaluateCellRatio(all, k, i1, i2, a, &ratio); } } } else if (cy1 == cy2) { for (cx = cx1; cx <= cx2; ++cx) { Convert2(grid, cx, cy1, k); if (grid->grid2[k] != NIL) { EvaluateCellRatio(all, k, i1, i2, a, &ratio); } } } } assert(gt(data->bbox_diagonal_sqd, ZERO)); ratio *= (1.0 + sqrt(sqrt(a / data->bbox_diagonal_sqd))); return ratio; } void ResetGrid(griddef *grid) { grid->buckets_initialized = false; return; } void InitBuckets(global_struct *all) { machine_double n_seg, a, b, fudge_x, fudge_y, md_delta_x, md_delta_y; griddef *grid = &all->c_grid; n_seg = sqrt((machine_double)grid->num_reflex); /* */ /* determine the grid size and location */ /* */ md_delta_x = (TO_MDOUBLE(all->bb_max.x) - TO_MDOUBLE(all->bb_min.x)) * CD_0_01; md_delta_y = (TO_MDOUBLE(all->bb_max.y) - TO_MDOUBLE(all->bb_min.y)) * CD_0_01; fudge_x = Max(Abs(TO_MDOUBLE(all->bb_max.x)), Abs(TO_MDOUBLE(all->bb_min.x))) * FUDGE; fudge_y = Max(Abs(TO_MDOUBLE(all->bb_max.y)), Abs(TO_MDOUBLE(all->bb_min.y))) * FUDGE; fudge_x = Max(fudge_x, FUDGE); fudge_y = Max(fudge_y, FUDGE); if (md_delta_x < fudge_x) md_delta_x = fudge_x; if (md_delta_y < fudge_y) md_delta_y = fudge_y; grid->grid_min_x3 = all->bb_min.x - md_delta_x; grid->grid_max_x3 = all->bb_max.x + md_delta_x; grid->grid_min_y3 = all->bb_min.y - md_delta_y; grid->grid_max_y3 = all->bb_max.y + md_delta_y; grid->delta_x = grid->grid_max_x3 - grid->grid_min_x3; grid->delta_y = grid->grid_max_y3 - grid->grid_min_y3; md_delta_x = TO_MDOUBLE(grid->grid_max_x3) - TO_MDOUBLE(grid->grid_min_x3); md_delta_y = TO_MDOUBLE(grid->grid_max_y3) - TO_MDOUBLE(grid->grid_min_y3); assert(grid->delta_x > ZERO); assert(grid->delta_y > ZERO); /* */ /* initialize the grid raster; we adapt the grid somewhat according to */ /* the aspect ratio of the bounding box of the segments. roughly, the */ /* number of cells equals the number of segments. */ /* */ a = sqrt(md_delta_x / md_delta_y); if (a < CD_0_05) { a = CD_0_05; b = CD_20_0; } else if (a > CD_20_0) { a = CD_20_0; b = CD_0_05; } else { b = CD_1_0 / a; } grid->N_x3 = (int)(CD_2_0 * n_seg * a); grid->N_y3 = (int)(CD_2_0 * n_seg * b); if (grid->N_x3 < 1) grid->N_x3 = 1; if (grid->N_y3 < 1) grid->N_y3 = 1; grid->buckets_initialized = true; return; } void BuildBuckets(global_struct *all, list_ind ind) { griddef *grid = &all->c_grid; listdef *list = &all->c_list; datadef *data = &all->c_data; list_ind ind1; int i, j, k, ii, jj, i1, N_x3_old, N_y3_old; double x, y, factor; int full_cells; boolean keep_going = false; boolean refine_grid = true; InitBuckets(all); grid->num_reflex = 0; full_cells = 0; do { grid->N3 = grid->N_x3 * grid->N_y3; grid->grid_x3 = grid->delta_x / grid->N_x3; grid->grid_y3 = grid->delta_y / grid->N_y3; InitGrid3(grid); InitRaster3(grid); InitVtxList(grid, data->num_pnts); /* */ /* insert reflex and tangential vertices into the grid */ /* */ ind1 = ind; grid->num_reflex = 0; do { if (GetAngle(list,ind1) <= 0) { i1 = GetIndex(list, ind1); x = data->points[i1].x; y = data->points[i1].y; DetermineCell3(grid, x, y, i, j, ii, jj); Convert3(grid, i, j, k); assert(InGrid3(grid, k)); InsertAfterVtxCells(list, grid, k, ind1, full_cells); if (i != ii) { Convert3(grid, ii, j, k); assert(InGrid3(grid, k)); InsertAfterVtxCells(list, grid, k, ind1, full_cells); if (j != jj) { Convert3(grid, i, jj, k); assert(InGrid3(grid, k)); InsertAfterVtxCells(list, grid, k, ind1, full_cells); Convert3(grid, ii, jj, k); assert(InGrid3(grid, k)); InsertAfterVtxCells(list, grid, k, ind1, full_cells); } } else if (j != jj) { Convert3(grid, i, jj, k); assert(InGrid3(grid, k)); InsertAfterVtxCells(list, grid, k, ind1, full_cells); } } ind1 = GetNextNode(list, ind1); } while (ind1 != ind); all->is_convex_face = (grid->num_reflex == 0); if (refine_grid && !all->is_convex_face) { assert(full_cells > 0); if (((double) grid->num_reflex) / ((double) (full_cells)) > N_REFLEX) { factor = ((double) grid->num_reflex) / ((double) (full_cells) * N_REFLEX); refine_grid = false; factor = sqrt(factor); N_x3_old = grid->N_x3; N_y3_old = grid->N_y3; grid->N_x3 = REAL_TO_INT((double) grid->N_x3 * C_1_1 * factor); grid->N_y3 = REAL_TO_INT((double) grid->N_y3 * C_1_1 * factor); if ((grid->N_x3 * grid->N_y3) > (grid->num_reflex * 1000.0)) { keep_going = false; grid->N_x3 = N_x3_old; grid->N_y3 = N_y3_old; } else { keep_going = true; } } else { keep_going = false; } } else { keep_going = false; } } while (keep_going); /* */ /* export the convexity status */ /* */ SetReflexNumber(all,grid->num_reflex); SetConvexityStatus(&all->c_ear, all->is_convex_face); return; } /* */ /* this function checks whether the cell with index k contains any reflex */ /* vertex that lies inside the triangle i1, i2, i3. it is assumed that the */ /* new diagonal is given by i2, i3, and that the triangle is oriented CCW. */ /* */ boolean CheckBucket(global_struct *all, int k, int i1, list_ind ind1, int i2, int i3, double *area, int *i_close, boolean check_close) { griddef *grid = &all->c_grid; listdef *list = &all->c_list; datadef *data = &all->c_data; tmp_grid_def tmp; int i4, pnt_on_edge; list_ind ind5; double area1; point p, q, r, v, w; tmp.ind_pnt = grid->grid3[k]; do { assert(InVtxList(grid, tmp.ind_pnt)); tmp.ind_vtx = grid->vtx_list[tmp.ind_pnt].pnt; assert(InPolyList(list, tmp.ind_vtx)); i4 = GetIndex(list, tmp.ind_vtx); ind5 = GetNextNode(list, tmp.ind_vtx); if ((tmp.ind_vtx != ind1) && (tmp.ind_vtx != ind5)) { /* */ /* only if this node isn't i1, and if it still belongs to the */ /* polygon */ /* */ if (i4 == i1) { if (!GetBridgeNode(list,tmp.ind_vtx)) { if (HandleDegeneracies(all, i1, ind1, i2, i3, i4, tmp.ind_vtx)) { return true; } } } else if ((i4 != i2) && (i4 != i3)) { if (PntInTriClass(data, i1, i2, i3, i4, &area1, &pnt_on_edge)) { if (pnt_on_edge <= 1) { return true; } else if (GetAngle(list, tmp.ind_vtx) < 0) { if (HandleDegeneracies(all, i1, ind1, i2, i3, i4, tmp.ind_vtx)) { return true; } } } else if (check_close) { p = data->points[i2]; q = data->points[i3]; r = data->points[i4]; VectorSub2D(q, p, v); VectorSub2D(r, p, w); if (ge(DotProduct2D(v, w), ZERO)) { VectorSub2D(r, q, w); if (le(DotProduct2D(v, w), ZERO)) { if (area1 > *area) { *area = area1; *i_close = i4; } } } } } } tmp.ind_pnt = grid->vtx_list[tmp.ind_pnt].next; } while (tmp.ind_pnt != NIL); return false; } boolean BucketIntersectionExists(global_struct *all, int i1, list_ind ind1, int i2, int i3, bounding_box bb, int *i_close, boolean check_close) { griddef *grid = &all->c_grid; listdef *list = &all->c_list; datadef *data = &all->c_data; double x, y; int i, j, k, cx1, cx2, cy1, cy2; double area = -fpkernel_DBL_MAX; assert(InPointsList(data, i1)); assert(InPointsList(data, i2)); assert(InPointsList(data, i3)); if (grid->num_reflex <= 0) { /* */ /* tell ear clipping that no reflex vertices are left */ /* */ SetConvexityStatus(&all->c_ear, true); return false; } else if (grid->num_reflex < grid->num_original_reflex) { #ifdef PARTITION_FIST #pragma omp barrier #pragma omp single #endif { ComputeBoundingBox(data, list, ind1, &all->bb_min, &all->bb_max); BuildBuckets(all, ind1); } } /* */ /* determine those cells of the grid that are intersected by the bounding */ /* box bb of the triangle i1, i2, i3. first, let's extend the BBox of */ /* the line segment i2, i3 to a BBox of the entire triangle. */ /* */ if (i1 < bb.imin) bb.imin = i1; else if (i1 > bb.imax) bb.imax = i1; y = data->points[i1].y; if (y < bb.ymin) bb.ymin = y; else if (y > bb.ymax) bb.ymax = y; x = data->points[bb.imin].x; y = bb.ymin; DetermineCell4(grid, x, y, cx1, cy1); x = data->points[bb.imax].x; y = bb.ymax; DetermineCell4(grid, x, y, cx2, cy2); assert(cx1 <= cx2); assert(cy1 <= cy2); /* */ /* check whether the triangle i1, i2, i3 contains any reflex vertex; we */ /* assume that i2, i3 is the new diagonal */ /* */ for (i = cx1; i <= cx2; ++i) { for (j = cy1; j <= cy2; ++j) { Convert3(grid, i, j, k); assert(InGrid3(grid, k)); if (grid->grid3[k] != NIL) { /* */ /* this cell contains at least one reflex vertex */ /* */ if (CheckBucket(all, k, i1, ind1, i2, i3, &area, i_close, check_close)) return true; } } } return false; } void DeleteReflexVertex(global_struct *all, list_ind ind) { griddef *grid = &all->c_grid; listdef *list = &all->c_list; datadef *data = &all->c_data; tmp_grid_def tmp; int i1, i, ii, j, jj, k; double x, y; /* */ /* determine the cell(s) which contain(s) this vertex */ /* */ assert(InPolyList(list, ind)); i1 = GetIndex(list, ind); x = data->points[i1].x; y = data->points[i1].y; DetermineCell3(grid, x, y, i, j, ii, jj); /* */ /* delete it from all the cells */ /* */ Convert3(grid, i, j, k); assert(InGrid3(grid, k)); if (grid->grid3[k] != NIL) DeleteFromCell(grid, tmp, k, ind); if (i != ii) { Convert3(grid, ii, j, k); assert(InGrid3(grid,k)); if (grid->grid3[k] != NIL) DeleteFromCell(grid, tmp, k, ind); if (j != jj) { Convert3(grid, i, jj, k); assert(InGrid3(grid,k)); if (grid->grid3[k] != NIL) DeleteFromCell(grid, tmp, k, ind); Convert3(grid, ii, jj, k); assert(InGrid3(grid,k)); if (grid->grid3[k] != NIL) DeleteFromCell(grid, tmp, k, ind); } } else if (j != jj) { Convert3(grid, i, jj, k); assert(InGrid3(grid,k)); if (grid->grid3[k] != NIL) DeleteFromCell(grid, tmp, k, ind); } if (grid->num_reflex <= 0) { /* */ /* tell ear clipping that no reflex vertices are left */ /* */ SetConvexityStatus(&all->c_ear, true); } return; } void BuildBridgeBuckets(global_struct *all, list_ind ind0, boolean init_grid) { griddef *grid = &all->c_grid; listdef *list = &all->c_list; datadef *data = &all->c_data; list_ind ind1; int i, j, k, i1; double x, y; machine_double n_seg, a, b, fudge_x, fudge_y, md_delta_x, md_delta_y; if (init_grid) { n_seg = sqrt((machine_double)data->num_pnts); /* */ /* determine the grid size and location */ /* */ md_delta_x = (TO_MDOUBLE(all->bb_max.x) - TO_MDOUBLE(all->bb_min.x)) * CD_0_01; md_delta_y = (TO_MDOUBLE(all->bb_max.y) - TO_MDOUBLE(all->bb_min.y)) * CD_0_01; fudge_x = Max(Abs(TO_MDOUBLE(all->bb_max.x)), Abs(TO_MDOUBLE(all->bb_min.x))) * FUDGE; fudge_y = Max(Abs(TO_MDOUBLE(all->bb_max.y)), Abs(TO_MDOUBLE(all->bb_min.y))) * FUDGE; fudge_x = Max(fudge_x, FUDGE); fudge_y = Max(fudge_y, FUDGE); if (md_delta_x < fudge_x) md_delta_x = fudge_x; if (md_delta_y < fudge_y) md_delta_y = fudge_y; grid->grid_min_x3 = all->bb_min.x - md_delta_x; grid->grid_max_x3 = all->bb_max.x + md_delta_x; grid->grid_min_y3 = all->bb_min.y - md_delta_y; grid->grid_max_y3 = all->bb_max.y + md_delta_y; grid->delta_x = grid->grid_max_x3 - grid->grid_min_x3; grid->delta_y = grid->grid_max_y3 - grid->grid_min_y3; md_delta_x = TO_MDOUBLE(grid->grid_max_x3) - TO_MDOUBLE(grid->grid_min_x3); md_delta_y = TO_MDOUBLE(grid->grid_max_y3) - TO_MDOUBLE(grid->grid_min_y3); assert(grid->delta_x > ZERO); assert(grid->delta_y > ZERO); /* */ /* initialize the grid raster; we adapt the grid somewhat according to */ /* the aspect ratio of the bounding box of the segments. roughly, the */ /* number of cells equals the number of segments. */ /* */ a = sqrt(md_delta_x / md_delta_y); if (a < CD_0_05) { a = CD_0_05; b = CD_20_0; } else if (a > CD_20_0) { a = CD_20_0; b = CD_0_05; } else { b = CD_1_0 / a; } grid->N_x3 = (int)(N_BUCKET * n_seg * a); grid->N_y3 = (int)(N_BUCKET * n_seg * b); if (grid->N_x3 < 1) grid->N_x3 = 1; if (grid->N_y3 < 1) grid->N_y3 = 1; grid->N3 = grid->N_x3 * grid->N_y3; grid->grid_x3 = grid->delta_x / grid->N_x3; grid->grid_y3 = grid->delta_y / grid->N_y3; InitGrid3(grid); InitRaster3(grid); InitVtxList(grid, data->num_pnts); } /* */ /* insert the vertices of loop loop_id into the grid */ /* */ ind1 = ind0; i1 = GetIndex(list,ind1); do { x = data->points[i1].x; y = data->points[i1].y; DetermineCell4(grid, x, y, i, j); Convert3(grid, i, j, k); assert(InGrid3(grid, k)); InsertAfterVtx(list, grid, k, ind1); ind1 = GetNextNode(list, ind1); i1 = GetIndex(list, ind1); } while (ind1 != ind0); return; } void InsertPntIntoBridgeBuckets(global_struct *all, list_ind ind1) { griddef *grid = &all->c_grid; listdef *list = &all->c_list; datadef *data = &all->c_data; int i1, i, j, k; double x, y; i1 = GetIndex(list,ind1); x = data->points[i1].x; y = data->points[i1].y; DetermineCell4(grid, x, y, i, j); Convert3(grid, i, j, k); assert(InGrid3(grid, k)); InsertAfterVtx(list, grid, k, ind1); return; } #ifdef GRAZING void ScanBridgeBuckets(global_struct *all, int i0, int grid_offset, distance *distances, int *num_dist, boolean *grid_exhausted, double x_start) #else void ScanBridgeBuckets(global_struct *all, int i0, int grid_offset, distance *distances, int *num_dist, boolean *grid_exhausted) #endif { griddef *grid = &all->c_grid; listdef *list = &all->c_list; datadef *data = &all->c_data; int i1, i, j, k, ii, jj, jj_min, jj_max, delta; double x, y; list_ind ind1; pnt_ind ind_pnt; md_point base; *num_dist = 0; x = data->points[i0].x; y = data->points[i0].y; DetermineCell4(grid, x, y, i, j); if (grid_offset == 0) { assert((0 <= i) && (i < grid->N_x)); assert((0 <= j) && (j < grid->N_y)); Convert3(grid, i, j, k); assert(InGrid3(grid, k)); #ifdef GRAZING ScanBridgeBucketCell(list, data, grid, k, i0, i1, ind1, ind_pnt, distances, num_dist, base, x_start); #else ScanBridgeBucketCell(list, data, grid, k, i0, i1, ind1, ind_pnt, distances, num_dist, base); #endif *grid_exhausted = false; } else { *grid_exhausted = true; jj_min = j - grid_offset; if (jj_min >= 0) { *grid_exhausted = false; ii = i; delta = 0; while ((delta <= grid_offset) && (ii >= 0)) { Convert3(grid, ii, jj_min, k); assert(InGrid3(grid, k)); #ifdef GRAZING ScanBridgeBucketCell(list, data, grid, k, i0, i1, ind1, ind_pnt, distances, num_dist, base, x_start); #else ScanBridgeBucketCell(list, data, grid, k, i0, i1, ind1, ind_pnt, distances,num_dist, base); #endif --ii; ++delta; } ++jj_min; } else { jj_min = 0; } jj_max = j + grid_offset; if (jj_max < grid->N_y3) { *grid_exhausted = false; ii = i; delta = 0; while ((delta <= grid_offset) && (ii >= 0)) { Convert3(grid, ii, jj_max, k); assert(InGrid3(grid, k)); #ifdef GRAZING ScanBridgeBucketCell(list, data, grid, k, i0, i1, ind1, ind_pnt, distances, num_dist, base, x_start); #else ScanBridgeBucketCell(list, data, grid, k, i0, i1, ind1, ind_pnt, distances, num_dist, base); #endif --ii; ++delta; } --jj_max; } else { jj_max = grid->N_y3 - 1; } ii = i - grid_offset; if (ii >= 0) { *grid_exhausted = false; for (jj = jj_min; jj <= jj_max; ++jj) { Convert3(grid, ii, jj, k); assert(InGrid3(grid, k)); #ifdef GRAZING ScanBridgeBucketCell(list, data, grid, k, i0, i1, ind1, ind_pnt, distances, num_dist, base, x_start); #else ScanBridgeBucketCell(list, data, grid, k, i0, i1, ind1, ind_pnt, distances, num_dist, base); #endif } } } return; }