Egalware/fist
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
fist/grid.cpp
T
SaraP cbec90699f FIST 6.8 : 
- creato il progetto in Visual Studio per compilare come libreria statica
- modifiche al codice originale per integrarlo nelle nostre librerie.
2025-03-04 16:19:35 +01:00

1850 lines
58 KiB
C++
Raw Permalink Blame History

/*****************************************************************************/
/* */
/* 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 <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <float.h>
/* 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;
}
Reference in New Issue View Git Blame Copy Permalink