/*****************************************************************************/ /* */ /* 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 libraries */ /* */ #include #include #include #include /* */ /* get my header files */ /* */ #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" #ifdef GRAPHICS #include "graphics.h" #endif /* */ /* function prototypes of functions provided in this file */ /* */ boolean Desperate(global_struct *all, heapdef *hp, list_ind ind, int i, boolean *splitted); boolean ExistsCrossOver(global_struct *all, list_ind ind, list_ind *ind1, int *i1, list_ind *ind2, int *i2, list_ind *ind3, int *i3, list_ind *ind4, int *i4); void HandleCrossOver(global_struct *all, heapdef *hp, list_ind ind1, int i1, list_ind ind2, int i2, list_ind ind3, int i3, list_ind ind4, int i4); boolean LetsHope(global_struct *all, heapdef *hp, list_ind ind); void HandleSplit(global_struct *all, list_ind ind1, int i1, list_ind ind3, int i3); boolean FoundSplit(global_struct *all, list_ind ind5, int i5, list_ind ind, list_ind ind1, int i1, int i3, int i4, list_ind *ind2, int *i2); boolean ExistsSplit(global_struct *all, list_ind ind, list_ind *ind1, int *i1, list_ind *ind2, int *i2); void FreeDistances(bridgedef *bridge); int WindingNumber(global_struct *all, list_ind ind, point p); /* */ /* the functions in this file try to ensure that we always end up with */ /* something that (topologically) is a triangulation. */ /* */ /* the more desperate we get, the more aggressive means we choose for making */ /* diagonals "valid". */ /* */ boolean Desperate(global_struct *all, heapdef *hp, list_ind ind, int i, boolean *splitted) { griddef *grid = &all->c_grid; int i1, i2, i3, i4; list_ind ind1, ind2, ind3, ind4; *splitted = false; /* */ /* check whether there exist consecutive vertices i1, i2, i3, i4 such */ /* that i1, i2 and i3, i4 intersect */ /* */ if (ExistsCrossOver(all, ind, &ind1, &i1, &ind2, &i2, &ind3, &i3, &ind4, &i4)) { /* */ /* insert two new diagonals around the cross-over without checking */ /* whether they are intersection-free */ /* */ HandleCrossOver(all, hp, ind1, i1, ind2, i2, ind3, i3, ind4, i4); return false; } if (grid->no_grid_yet) { BuildGrid(all, i, i+1); } /* */ /* check whether there exists a valid diagonal that splits the polygon */ /* into two parts */ /* */ if (ExistsSplit(all, ind, &ind1, &i1, &ind2, &i2)) { /* */ /* break up the polygon by inserting this diagonal (which can't be an */ /* ear -- otherwise, we would not have ended up in this part of the */ /* code). then, let's treat the two polygons separately. hopefully, */ /* this will help to handle self-overlapping polygons in the "correct" */ /* way. */ /* */ HandleSplit(all, ind1, i1, ind2, i2); *splitted = true; return false; } return true; } boolean ExistsCrossOver(global_struct *all, list_ind ind, list_ind *ind1, int *i1, list_ind *ind2, int *i2, list_ind *ind3, int *i3, list_ind *ind4, int *i4) { datadef *data = &all->c_data; listdef *list = &all->c_list; bounding_box bb1, bb2; *ind1 = ind; *i1 = GetIndex(list, *ind1); *ind2 = GetNextNode(list, *ind1); *i2 = GetIndex(list, *ind2); *ind3 = GetNextNode(list, *ind2); *i3 = GetIndex(list, *ind3); *ind4 = GetNextNode(list, *ind3); *i4 = GetIndex(list, *ind4); do { BBox(data, *i1, *i2, bb1); BBox(data, *i3, *i4, bb2); if (BBox_Overlap(bb1, bb2)) { if (SegIntersect(all, bb1.imin, bb1.imax, bb2.imin, bb2.imax, -1)) return true; } *ind1 = *ind2; *i1 = *i2; *ind2 = *ind3; *i2 = *i3; *ind3 = *ind4; *i3 = *i4; *ind4 = GetNextNode(list, *ind3); *i4 = GetIndex(list, *ind4); } while (*ind1 != ind); return false; } void HandleCrossOver(global_struct *all, heapdef *hp, list_ind ind1, int i1, list_ind ind2, int i2, list_ind ind3, int i3, list_ind ind4, int i4) { datadef *data = &all->c_data; listdef *list = &all->c_list; machine_double ratio1, ratio4; boolean first; int angle1, angle4; /* */ /* which pair of triangles shall I insert?? we can use either i1, i2, i3 */ /* and i1, i3, i4, or we can use i2, i3, i4 and i1, i2, i4... */ /* */ angle1 = GetAngle(list, ind1); angle4 = GetAngle(list, ind4); if (angle1 < angle4) first = true; else if (angle1 > angle4) first = false; else if (hp->ears_sorted) { ratio1 = GetRatio(data, i3, i4, i1); ratio4 = GetRatio(data, i1, i2, i4); if (ratio4 < ratio1) first = false; else first = true; } else { first = true; } if (first) { /* */ /* first clip i1, i2, i3, then clip i1, i3, i4 */ /* */ DeleteLinks(list, ind2); StoreTriangle(all, ind1, ind2, ind3,TriTriColor); SetAngle(list, ind3, 1); InsertIntoHeap(hp, 0.0, ind3, ind1, ind4); #ifdef GRAPHICS if (all->rt_opt.graphics) { AddTriToBuffer(&all->c_redraw, i1, i2, i3, TriColor, TriFillColor); DrawTri(data->points[i1], data->points[i2], data->points[i3], TriColor, TriFillColor); DrawPnt(data->points[i3], EarColor); } #endif } else { /* */ /* first clip i2, i3, i4, then clip i1, i2, i4 */ /* */ DeleteLinks(list, ind3); StoreTriangle(all, ind2, ind3, ind4, TriTriColor); SetAngle(list, ind2, 1); InsertIntoHeap(hp, 0.0, ind2, ind1, ind4); #ifdef GRAPHICS if (all->rt_opt.graphics) { AddTriToBuffer(&all->c_redraw, i2, i3, i4, TriColor, TriFillColor); DrawTri(data->points[i2], data->points[i3], data->points[i4], TriColor, TriFillColor); DrawPnt(data->points[i2], EarColor); } #endif } return; } boolean LetsHope(global_struct *all, heapdef *hp, list_ind ind) { listdef *list = &all->c_list; #ifdef GRAPHICS datadef *data = &all->c_data; int i1; #endif list_ind ind0, ind1, ind2; /* */ /* perhaps we find a 180-degree corner? such a corner is not the best of */ /* all ears. but, at least, we'd reduce the vertex count by one without */ /* causing any real harm... */ /* */ ind1 = ind; do { if (GetAngle(list, ind1) == 0) { SetAngle(list, ind1, 1); ind0 = GetPrevNode(list, ind1); ind2 = GetNextNode(list, ind1); InsertIntoHeap(hp, 0.0, ind1, ind0, ind2); #ifdef GRAPHICS if (all->rt_opt.graphics) { i1 = GetIndex(list, ind1); DrawPnt(data->points[i1], EarColor); } #endif return true; } ind1 = GetNextNode(list, ind1); } while (ind1 != ind); /* */ /* let's clip the first convex corner. of course, we know that this is no */ /* ear in an ideal world. but this polygon isn't ideal, either! */ /* */ ind1 = ind; do { if (GetAngle(list, ind1) > 0) { ind0 = GetPrevNode(list, ind1); ind2 = GetNextNode(list, ind1); InsertIntoHeap(hp, 0.0, ind1, ind0, ind2); #ifdef GRAPHICS if (all->rt_opt.graphics) { i1 = GetIndex(list, ind1); DrawPnt(data->points[i1], EarColor); } #endif return true; } ind1 = GetNextNode(list, ind1); } while (ind1 != ind); /* */ /* no convex and tangential corners? so, let's cheat! this code ain't */ /* stop without some triangulation... ;-) g-i-g-o? right! perhaps, */ /* this is what you call a robust code?! */ /* */ SetAngle(list, ind, 1); ind0 = GetPrevNode(list, ind); ind2 = GetNextNode(list, ind); InsertIntoHeap(hp, 0.0, ind, ind0, ind2); #ifdef GRAPHICS if (all->rt_opt.graphics) { i1 = GetIndex(list, ind); DrawPnt(data->points[i1], EarColor); } #endif return true; /* */ /* see, we never have to return "false"... ;-) */ /* */ /* return false; */ } boolean ExistsSplit(global_struct *all, list_ind ind, list_ind *ind1, int *i1, list_ind *ind2, int *i2) { listdef *list = &all->c_list; bridgedef *bridge = &all->c_bridge; list_ind ind3, ind4, ind5; int i3, i4, i5, number; number = GetNumList(list); if (number > bridge->max_num_dist) { bridge->distances = CORE_ReallocateArray(bridge->memptr, bridge->distances, bridge->max_num_dist, number, distance, "desperate:distances"); bridge->max_num_dist = number; } *ind1 = ind; *i1 = GetIndex(list, *ind1); ind4 = GetNextNode(list, *ind1); i4 = GetIndex(list, ind4); assert(*ind1 != ind4); ind5 = GetNextNode(list, ind4); i5 = GetIndex(list, ind5); assert(*ind1 != *ind2); ind3 = GetPrevNode(list, *ind1); i3 = GetIndex(list, ind3); assert(*ind2 != ind3); if (FoundSplit(all, ind5, i5, ind3, *ind1, *i1, i3, i4, ind2, i2)) return true; i3 = *i1; *ind1 = ind4; *i1 = i4; ind4 = ind5; i4 = i5; ind5 = GetNextNode(list, ind4); i5 = GetIndex(list, ind5); while (ind5 != ind) { if (FoundSplit(all, ind5, i5, ind, *ind1, *i1, i3, i4, ind2, i2)) return true; i3 = *i1; *ind1 = ind4; *i1 = i4; ind4 = ind5; i4 = i5; ind5 = GetNextNode(list, ind4); i5 = GetIndex(list, ind5); } return false; } /* */ /* this function computes the winding number of a polygon with respect to a */ /* point p. no care is taken to handle cases where p lies on the */ /* boundary of the polygon. (this is no issue in our application, as we will */ /* always compute the winding number with respect to the mid-point of a */ /* valid diagonal.) */ /* */ int WindingNumber(global_struct *all, list_ind ind, point p) { listdef *list = &all->c_list; datadef *data = &all->c_data; double angle; list_ind ind2; int i1, i2, number; i1 = GetIndex(list, ind); ind2 = GetNextNode(list, ind); i2 = GetIndex(list, ind2); angle = Angle(p, data->points[i1], data->points[i2]); while (ind2 != ind) { i1 = i2; ind2 = GetNextNode(list, ind2); i2 = GetIndex(list, ind2); angle += Angle(p, data->points[i1], data->points[i2]); } angle += M_PI; number = REAL_TO_INT(angle / M_2PI); return number; } boolean FoundSplit(global_struct *all, list_ind ind5, int i5, list_ind ind, list_ind ind1, int i1, int i3, int i4, list_ind *ind2, int *i2) { listdef *list = &all->c_list; bridgedef *bridge = &all->c_bridge; datadef *data = &all->c_data; tmp_data_def tmp_data; point base, center; int num_dist = 0; int j, i6, i7; list_ind ind6, ind7; bounding_box bb; boolean convex, cone_ok; int d_comp(const void *a, const void *); /* */ /* sort the points according to their distance from i1 */ /* */ do { assert(num_dist < bridge->max_num_dist); BaseLength(data->points[i1], data->points[i5], base, bridge->distances[num_dist].dist); bridge->distances[num_dist].ind = ind5; ++num_dist; ind5 = GetNextNode(list, ind5); i5 = GetIndex(list, ind5); } while (ind5 != ind); qsort(bridge->distances, num_dist, sizeof(distance), &d_comp); /* */ /* find a valid diagonal. */ /* */ for (j = 0; j < num_dist; ++j) { *ind2 = bridge->distances[j].ind; *i2 = GetIndex(list, *ind2); if (i1 != *i2) { ind6 = GetPrevNode(list, *ind2); i6 = GetIndex(list, ind6); ind7 = GetNextNode(list, *ind2); i7 = GetIndex(list, ind7); convex = GetAngle(list, *ind2) > 0; IsInCone(data, tmp_data, i6, *i2, i7, i1, convex, cone_ok); if (cone_ok) { convex = GetAngle(list, ind1) > 0; IsInCone(data, tmp_data, i3, i1, i4, *i2, convex, cone_ok); if (cone_ok) { BBox(data, i1, *i2, bb); if (!GridIntersectionExists(all, bb, -1, -1, ind1, -1)) { /* */ /* check whether this is a good diagonal; we do not want a */ /* diagonal that may create figure-8's! */ /* */ VectorAdd2D(data->points[i1], data->points[*i2], center); MultScalar2D(0.5, center); if (WindingNumber(all, ind, center) == 1) return true; } } } } } return false; } void HandleSplit(global_struct *all, list_ind ind1, int i1, list_ind ind3, int i3) { listdef *list = &all->c_list; datadef *data = &all->c_data; list_ind ind2, ind4, prev, next; int prv, nxt, angle; /* */ /* duplicate nodes in order to form end points of the new diagonal */ /* */ ind2 = MakeNode(list, i1); InsertAfter(list, ind1, ind2); SetOriginal(list, ind2, GetOriginal(list, ind1)); ind4 = MakeNode(list, i3); InsertAfter(list, ind3, ind4); SetOriginal(list, ind4, GetOriginal(list, ind3)); /* */ /* insert the diagonal into the boundary loop, thus splitting the loop */ /* into two loops */ /* */ SplitSplice(list, ind1, ind2, ind3, ind4); /* */ /* store pointers to the two new loops */ /* */ StoreChain(list, ind1); StoreChain(list, ind3); /* */ /* update the graphics */ /* */ #ifdef GRAPHICS if (all->rt_opt.graphics) { DrawSeg(data->points[i1], data->points[i3], SplitColor); AddEdgeToBuffer(&all->c_redraw, i1, i3, BridgeColor); } #endif /* */ /* reset the angles */ /* */ next = GetNextNode(list, ind1); nxt = GetIndex(list, next); prev = GetPrevNode(list, ind1); prv = GetIndex(list, prev); angle = IsConvexAngle(list, data, prv, i1, nxt, ind1); SetAngle(list, ind1, angle); #ifdef GRAPHICS if (all->rt_opt.graphics) { if (angle == 1) DrawPnt(data->points[i1], ConvexColor); else if (angle == 2) DrawPnt(data->points[i1], ZeroColor); else if (angle == 0) DrawPnt(data->points[i1], TangentColor); else DrawPnt(data->points[i1], PntsColor); } #endif next = GetNextNode(list, ind2); nxt = GetIndex(list, next); prev = GetPrevNode(list, ind2); prv = GetIndex(list, prev); angle = IsConvexAngle(list, data, prv, i1, nxt, ind2); SetAngle(list, ind2, angle); next = GetNextNode(list, ind3); nxt = GetIndex(list, next); prev = GetPrevNode(list, ind3); prv = GetIndex(list, prev); angle = IsConvexAngle(list, data, prv, i3, nxt, ind3); SetAngle(list, ind3, angle); #ifdef GRAPHICS if (all->rt_opt.graphics) { if (angle == 1) DrawPnt(data->points[i3], ConvexColor); else if (angle == 2) DrawPnt(data->points[i3], ZeroColor); else if (angle == 0) DrawPnt(data->points[i3], TangentColor); else DrawPnt(data->points[i1], PntsColor); } #endif next = GetNextNode(list, ind4); nxt = GetIndex(list, next); prev = GetPrevNode(list, ind4); prv = GetIndex(list, prev); angle = IsConvexAngle(list, data, prv, i3, nxt, ind4); SetAngle(list, ind4, angle); return; } void FreeDistances(bridgedef *bridge) { bridge->max_num_dist = 0; CORE_FreeMemory(bridge->memptr, &bridge->distances, "bridge:distances"); return; }