diff options
Diffstat (limited to 'extern/carve/lib/triangulator.cpp')
| -rw-r--r-- | extern/carve/lib/triangulator.cpp | 1200 |
1 files changed, 0 insertions, 1200 deletions
diff --git a/extern/carve/lib/triangulator.cpp b/extern/carve/lib/triangulator.cpp deleted file mode 100644 index a48a9ccfa45..00000000000 --- a/extern/carve/lib/triangulator.cpp +++ /dev/null @@ -1,1200 +0,0 @@ -// Begin License: -// Copyright (C) 2006-2014 Tobias Sargeant (tobias.sargeant@gmail.com). -// All rights reserved. -// -// This file is part of the Carve CSG Library (http://carve-csg.com/) -// -// This file may be used under the terms of either the GNU General -// Public License version 2 or 3 (at your option) as published by the -// Free Software Foundation and appearing in the files LICENSE.GPL2 -// and LICENSE.GPL3 included in the packaging of this file. -// -// This file is provided "AS IS" with NO WARRANTY OF ANY KIND, -// INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR -// A PARTICULAR PURPOSE. -// End: - - -#if defined(HAVE_CONFIG_H) -# include <carve_config.h> -#endif - -#include <carve/csg.hpp> -#include <carve/triangulator.hpp> - -#include <fstream> -#include <sstream> - -#include <algorithm> - -namespace { - // private code related to hole patching. - - class order_h_loops_2d { - order_h_loops_2d &operator=(const order_h_loops_2d &); - - const std::vector<std::vector<carve::geom2d::P2> > &poly; - int axis; - public: - - order_h_loops_2d(const std::vector<std::vector<carve::geom2d::P2> > &_poly, int _axis) : - poly(_poly), axis(_axis) { - } - - bool operator()(const std::pair<size_t, size_t> &a, - const std::pair<size_t, size_t> &b) const { - return carve::triangulate::detail::axisOrdering(poly[a.first][a.second], poly[b.first][b.second], axis); - } - }; - - class heap_ordering_2d { - heap_ordering_2d &operator=(const heap_ordering_2d &); - - const std::vector<std::vector<carve::geom2d::P2> > &poly; - const std::vector<std::pair<size_t, size_t> > &loop; - const carve::geom2d::P2 p; - int axis; - - public: - - heap_ordering_2d(const std::vector<std::vector<carve::geom2d::P2> > &_poly, - const std::vector<std::pair<size_t, size_t> > &_loop, - const carve::geom2d::P2 _p, - int _axis) : poly(_poly), loop(_loop), p(_p), axis(_axis) { - } - - bool operator()(size_t a, size_t b) const { - double da = carve::geom::distance2(p, poly[loop[a].first][loop[a].second]); - double db = carve::geom::distance2(p, poly[loop[b].first][loop[b].second]); - if (da > db) return true; - if (da < db) return false; - return carve::triangulate::detail::axisOrdering(poly[loop[a].first][loop[a].second], poly[loop[b].first][loop[b].second], axis); - } - }; - - static inline void patchHoleIntoPolygon_2d(std::vector<std::pair<size_t, size_t> > &f_loop, - size_t f_loop_attach, - size_t h_loop, - size_t h_loop_attach, - size_t h_loop_size) { - f_loop.insert(f_loop.begin() + f_loop_attach + 1, h_loop_size + 2, std::make_pair(h_loop, 0)); - size_t f = f_loop_attach + 1; - - for (size_t h = h_loop_attach; h != h_loop_size; ++h) { - f_loop[f++].second = h; - } - - for (size_t h = 0; h <= h_loop_attach; ++h) { - f_loop[f++].second = h; - } - - f_loop[f] = f_loop[f_loop_attach]; - } - - static inline const carve::geom2d::P2 &pvert(const std::vector<std::vector<carve::geom2d::P2> > &poly, const std::pair<size_t, size_t> &idx) { - return poly[idx.first][idx.second]; - } -} - - -namespace { - // private code related to triangulation. - - using carve::triangulate::detail::vertex_info; - - struct vertex_info_ordering { - bool operator()(const vertex_info *a, const vertex_info *b) const { - return a->score < b->score; - } - }; - - struct vertex_info_l2norm_inc_ordering { - const vertex_info *v; - vertex_info_l2norm_inc_ordering(const vertex_info *_v) : v(_v) { - } - bool operator()(const vertex_info *a, const vertex_info *b) const { - return carve::geom::distance2(v->p, a->p) > carve::geom::distance2(v->p, b->p); - } - }; - - class EarQueue { - std::vector<vertex_info *> queue; - - void checkheap() { -#if defined(HAVE_IS_HEAP) - CARVE_ASSERT(std::__is_heap(queue.begin(), queue.end(), vertex_info_ordering())); -#endif - } - - public: - EarQueue() { - } - - size_t size() const { - return queue.size(); - } - - void push(vertex_info *v) { -#if defined(CARVE_DEBUG) - checkheap(); -#endif - queue.push_back(v); - std::push_heap(queue.begin(), queue.end(), vertex_info_ordering()); - } - - vertex_info *pop() { -#if defined(CARVE_DEBUG) - checkheap(); -#endif - std::pop_heap(queue.begin(), queue.end(), vertex_info_ordering()); - vertex_info *v = queue.back(); - queue.pop_back(); - return v; - } - - void remove(vertex_info *v) { -#if defined(CARVE_DEBUG) - checkheap(); -#endif - CARVE_ASSERT(std::find(queue.begin(), queue.end(), v) != queue.end()); - double score = v->score; - if (v != queue[0]) { - v->score = queue[0]->score + 1; - std::make_heap(queue.begin(), queue.end(), vertex_info_ordering()); - } - CARVE_ASSERT(v == queue[0]); - std::pop_heap(queue.begin(), queue.end(), vertex_info_ordering()); - CARVE_ASSERT(queue.back() == v); - queue.pop_back(); - v->score = score; - } - - void changeScore(vertex_info *v, double score) { -#if defined(CARVE_DEBUG) - checkheap(); -#endif - CARVE_ASSERT(std::find(queue.begin(), queue.end(), v) != queue.end()); - if (v->score != score) { - v->score = score; - std::make_heap(queue.begin(), queue.end(), vertex_info_ordering()); - } - } - - // 39% of execution time - void updateVertex(vertex_info *v) { - double spre = v->score; - bool qpre = v->isCandidate(); - v->recompute(); - bool qpost = v->isCandidate(); - double spost = v->score; - - v->score = spre; - - if (qpre) { - if (qpost) { - if (v->score != spre) { - changeScore(v, spost); - } - } else { - remove(v); - } - } else { - if (qpost) { - push(v); - } - } - } - }; - - - - int windingNumber(vertex_info *begin, const carve::geom2d::P2 &point) { - int wn = 0; - - vertex_info *v = begin; - do { - if (v->p.y <= point.y) { - if (v->next->p.y > point.y && carve::geom2d::orient2d(v->p, v->next->p, point) > 0.0) { - ++wn; - } - } else { - if (v->next->p.y <= point.y && carve::geom2d::orient2d(v->p, v->next->p, point) < 0.0) { - --wn; - } - } - v = v->next; - } while (v != begin); - - return wn; - } - - - - bool internalToAngle(const vertex_info *a, - const vertex_info *b, - const vertex_info *c, - const carve::geom2d::P2 &p) { - return carve::geom2d::internalToAngle(a->p, b->p, c->p, p); - } - - - - bool findDiagonal(vertex_info *begin, vertex_info *&v1, vertex_info *&v2) { - vertex_info *t; - std::vector<vertex_info *> heap; - - v1 = begin; - do { - heap.clear(); - - for (v2 = v1->next->next; v2 != v1->prev; v2 = v2->next) { - if (!internalToAngle(v1->next, v1, v1->prev, v2->p) || - !internalToAngle(v2->next, v2, v2->prev, v1->p)) continue; - - heap.push_back(v2); - std::push_heap(heap.begin(), heap.end(), vertex_info_l2norm_inc_ordering(v1)); - } - - while (heap.size()) { - std::pop_heap(heap.begin(), heap.end(), vertex_info_l2norm_inc_ordering(v1)); - v2 = heap.back(); heap.pop_back(); - -#if defined(CARVE_DEBUG) - std::cerr << "testing: " << v1 << " - " << v2 << std::endl; - std::cerr << " length = " << (v2->p - v1->p).length() << std::endl; - std::cerr << " pos: " << v1->p << " - " << v2->p << std::endl; -#endif - // test whether v1-v2 is a valid diagonal. - double v_min_x = std::min(v1->p.x, v2->p.x); - double v_max_x = std::max(v1->p.x, v2->p.x); - - bool intersected = false; - - for (t = v1->next; !intersected && t != v1->prev; t = t->next) { - vertex_info *u = t->next; - if (t == v2 || u == v2) continue; - - double l1 = carve::geom2d::orient2d(v1->p, v2->p, t->p); - double l2 = carve::geom2d::orient2d(v1->p, v2->p, u->p); - - if ((l1 > 0.0 && l2 > 0.0) || (l1 < 0.0 && l2 < 0.0)) { - // both on the same side; no intersection - continue; - } - - double dx13 = v1->p.x - t->p.x; - double dy13 = v1->p.y - t->p.y; - double dx43 = u->p.x - t->p.x; - double dy43 = u->p.y - t->p.y; - double dx21 = v2->p.x - v1->p.x; - double dy21 = v2->p.y - v1->p.y; - double ua_n = dx43 * dy13 - dy43 * dx13; - double ub_n = dx21 * dy13 - dy21 * dx13; - double u_d = dy43 * dx21 - dx43 * dy21; - - if (carve::math::ZERO(u_d)) { - // parallel - if (carve::math::ZERO(ua_n)) { - // colinear - if (std::max(t->p.x, u->p.x) >= v_min_x && std::min(t->p.x, u->p.x) <= v_max_x) { - // colinear and intersecting - intersected = true; - } - } - } else { - // not parallel - double ua = ua_n / u_d; - double ub = ub_n / u_d; - - if (0.0 <= ua && ua <= 1.0 && 0.0 <= ub && ub <= 1.0) { - intersected = true; - } - } -#if defined(CARVE_DEBUG) - if (intersected) { - std::cerr << " failed on edge: " << t << " - " << u << std::endl; - std::cerr << " pos: " << t->p << " - " << u->p << std::endl; - } -#endif - } - - if (!intersected) { - // test whether midpoint winding == 1 - - carve::geom2d::P2 mid = (v1->p + v2->p) / 2; - if (windingNumber(begin, mid) == 1) { - // this diagonal is ok - return true; - } - } - } - - // couldn't find a diagonal from v1 that was ok. - v1 = v1->next; - } while (v1 != begin); - return false; - } - - - -#if defined(CARVE_DEBUG_WRITE_PLY_DATA) - void dumpPoly(const std::vector<carve::geom2d::P2> &points, - const std::vector<carve::triangulate::tri_idx> &result) { - static int step = 0; - std::ostringstream filename; - filename << "poly_" << step++ << ".svg"; - std::cerr << "dumping to " << filename.str() << std::endl; - std::ofstream out(filename.str().c_str()); - - double minx = points[0].x, maxx = points[0].x; - double miny = points[0].y, maxy = points[0].y; - - for (size_t i = 1; i < points.size(); ++i) { - minx = std::min(points[i].x, minx); maxx = std::max(points[i].x, maxx); - miny = std::min(points[i].y, miny); maxy = std::max(points[i].y, maxy); - } - double scale = 100 / std::max(maxx-minx, maxy-miny); - - maxx *= scale; minx *= scale; - maxy *= scale; miny *= scale; - - double width = maxx - minx + 10; - double height = maxy - miny + 10; - - out << "\ -<?xml version=\"1.0\"?>\n\ -<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n\ -<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\" width=\"" << width << "\" height=\"" << height << "\">\n"; - - out << "<polygon fill=\"rgb(0,0,0)\" stroke=\"blue\" stroke-width=\"0.1\" points=\""; - for (size_t i = 0; i < points.size(); ++i) { - if (i) out << ' '; - double x, y; - x = scale * (points[i].x) - minx + 5; - y = scale * (points[i].y) - miny + 5; - out << x << ',' << y; - } - out << "\" />" << std::endl; - - for (size_t i = 0; i < result.size(); ++i) { - out << "<polygon fill=\"rgb(255,255,255)\" stroke=\"black\" stroke-width=\"0.1\" points=\""; - double x, y; - x = scale * (points[result[i].a].x) - minx + 5; - y = scale * (points[result[i].a].y) - miny + 5; - out << x << ',' << y << ' '; - x = scale * (points[result[i].b].x) - minx + 5; - y = scale * (points[result[i].b].y) - miny + 5; - out << x << ',' << y << ' '; - x = scale * (points[result[i].c].x) - minx + 5; - y = scale * (points[result[i].c].y) - miny + 5; - out << x << ',' << y; - out << "\" />" << std::endl; - } - - out << "</svg>" << std::endl; - } -#endif -} - - - -double carve::triangulate::detail::vertex_info::triScore(const vertex_info *p, const vertex_info *v, const vertex_info *n) { - - // different scoring functions. -#if 0 - bool convex = isLeft(p, v, n); - if (!convex) return -1e-5; - - double a1 = carve::geom2d::atan2(p->p - v->p) - carve::geom2d::atan2(n->p - v->p); - double a2 = carve::geom2d::atan2(v->p - n->p) - carve::geom2d::atan2(p->p - n->p); - if (a1 < 0) a1 += M_PI * 2; - if (a2 < 0) a2 += M_PI * 2; - - return std::min(a1, std::min(a2, M_PI - a1 - a2)) / (M_PI / 3); -#endif - -#if 1 - // range: 0 - 1 - double a, b, c; - - bool convex = isLeft(p, v, n); - if (!convex) return -1e-5; - - a = (n->p - v->p).length(); - b = (p->p - n->p).length(); - c = (v->p - p->p).length(); - - if (a < 1e-10 || b < 1e-10 || c < 1e-10) return 0.0; - - return std::max(std::min((a+b)/c, std::min((a+c)/b, (b+c)/a)) - 1.0, 0.0); -#endif -} - - - -double carve::triangulate::detail::vertex_info::calcScore() const { - -#if 0 - // examine only this triangle. - double this_tri = triScore(prev, this, next); - return this_tri; -#endif - -#if 1 - // attempt to look ahead in the neighbourhood to attempt to clip ears that have good neighbours. - double this_tri = triScore(prev, this, next); - double next_tri = triScore(prev, next, next->next); - double prev_tri = triScore(prev->prev, prev, next); - - return this_tri + std::max(next_tri, prev_tri) * .2; -#endif - -#if 0 - // attempt to penalise ears that will require producing a sliver triangle. - double score = triScore(prev, this, next); - - double a1, a2; - a1 = carve::geom2d::atan2(prev->p - next->p); - a2 = carve::geom2d::atan2(next->next->p - next->p); - if (fabs(a1 - a2) < 1e-5) score -= .5; - - a1 = carve::geom2d::atan2(next->p - prev->p); - a2 = carve::geom2d::atan2(prev->prev->p - prev->p); - if (fabs(a1 - a2) < 1e-5) score -= .5; - - return score; -#endif -} - - - -bool carve::triangulate::detail::vertex_info::isClipable() const { - for (const vertex_info *v_test = next->next; v_test != prev; v_test = v_test->next) { - if (v_test->convex) { - continue; - } - - if (v_test->p == prev->p || - v_test->p == next->p) { - continue; - } - - if (v_test->p == p) { - if (v_test->next->p == prev->p && - v_test->prev->p == next->p) { - return false; - } - if (v_test->next->p == prev->p || - v_test->prev->p == next->p) { - continue; - } - } - - if (pointInTriangle(prev, this, next, v_test)) { - return false; - } - } - return true; -} - - - -size_t carve::triangulate::detail::removeDegeneracies(vertex_info *&begin, std::vector<carve::triangulate::tri_idx> &result) { - vertex_info *v; - vertex_info *n; - size_t count = 0; - size_t remain = 0; - - v = begin; - do { - v = v->next; - ++remain; - } while (v != begin); - - v = begin; - do { - if (remain < 4) break; - - bool remove = false; - if (v->p == v->next->p) { - remove = true; - } else if (v->p == v->next->next->p) { - if (v->next->p == v->next->next->next->p) { - // a 'z' in the loop: z (a) b a b c -> remove a-b-a -> z (a) a b c -> remove a-a-b (next loop) -> z a b c - // z --(a)-- b - // / - // / - // a -- b -- d - remove = true; - } else { - // a 'shard' in the loop: z (a) b a c d -> remove a-b-a -> z (a) a b c d -> remove a-a-b (next loop) -> z a b c d - // z --(a)-- b - // / - // / - // a -- c -- d - // n.b. can only do this if the shard is pointing out of the polygon. i.e. b is outside z-a-c - remove = !internalToAngle(v->next->next->next, v, v->prev, v->next->p); - } - } - - if (remove) { - result.push_back(carve::triangulate::tri_idx(v->idx, v->next->idx, v->next->next->idx)); - n = v->next; - if (n == begin) begin = n->next; - n->remove(); - count++; - remain--; - delete n; - } else { - v = v->next; - } - } while (v != begin); - return count; -} - - - -bool carve::triangulate::detail::splitAndResume(vertex_info *begin, std::vector<carve::triangulate::tri_idx> &result) { - vertex_info *v1, *v2; - -#if defined(CARVE_DEBUG_WRITE_PLY_DATA) - { - std::vector<carve::triangulate::tri_idx> dummy; - std::vector<carve::geom2d::P2> dummy_p; - vertex_info *v = begin; - do { - dummy_p.push_back(v->p); - v = v->next; - } while (v != begin); - std::cerr << "input to splitAndResume:" << std::endl; - dumpPoly(dummy_p, dummy); - } -#endif - - - if (!findDiagonal(begin, v1, v2)) return false; - - vertex_info *v1_copy = new vertex_info(*v1); - vertex_info *v2_copy = new vertex_info(*v2); - - v1->next = v2; - v2->prev = v1; - - v1_copy->next->prev = v1_copy; - v2_copy->prev->next = v2_copy; - - v1_copy->prev = v2_copy; - v2_copy->next = v1_copy; - - bool r1 = doTriangulate(v1, result); - bool r2 = doTriangulate(v1_copy, result); - return r1 && r2; -} - - - -bool carve::triangulate::detail::doTriangulate(vertex_info *begin, std::vector<carve::triangulate::tri_idx> &result) { -#if defined(CARVE_DEBUG) - std::cerr << "entering doTriangulate" << std::endl; -#endif - -#if defined(CARVE_DEBUG_WRITE_PLY_DATA) - { - std::vector<carve::triangulate::tri_idx> dummy; - std::vector<carve::geom2d::P2> dummy_p; - vertex_info *v = begin; - do { - dummy_p.push_back(v->p); - v = v->next; - } while (v != begin); - dumpPoly(dummy_p, dummy); - } -#endif - - EarQueue vq; - - vertex_info *v = begin; - size_t remain = 0; - do { - if (v->isCandidate()) vq.push(v); - v = v->next; - remain++; - } while (v != begin); - -#if defined(CARVE_DEBUG) - std::cerr << "remain = " << remain << std::endl; -#endif - - while (remain > 3 && vq.size()) { - vertex_info *v = vq.pop(); - if (!v->isClipable()) { - v->failed = true; - continue; - } - - continue_clipping: - vertex_info *n = v->next; - vertex_info *p = v->prev; - - result.push_back(carve::triangulate::tri_idx(v->prev->idx, v->idx, v->next->idx)); - -#if defined(CARVE_DEBUG) - { - std::vector<carve::geom2d::P2> temp; - temp.push_back(v->prev->p); - temp.push_back(v->p); - temp.push_back(v->next->p); - std::cerr << "clip " << v << " idx = " << v->idx << " score = " << v->score << " area = " << carve::geom2d::signedArea(temp) << " " << temp[0] << " " << temp[1] << " " << temp[2] << std::endl; - } -#endif - - v->remove(); - if (v == begin) begin = v->next; - delete v; - - if (--remain == 3) break; - - vq.updateVertex(n); - vq.updateVertex(p); - - if (n->score < p->score) { std::swap(n, p); } - - if (n->score > 0.25 && n->isCandidate() && n->isClipable()) { - vq.remove(n); - v = n; -#if defined(CARVE_DEBUG) - std::cerr << " continue clipping (n), score = " << n->score << std::endl; -#endif - goto continue_clipping; - } - - if (p->score > 0.25 && p->isCandidate() && p->isClipable()) { - vq.remove(p); - v = p; -#if defined(CARVE_DEBUG) - std::cerr << " continue clipping (p), score = " << n->score << std::endl; -#endif - goto continue_clipping; - } - -#if defined(CARVE_DEBUG) - std::cerr << "looking for new start point" << std::endl; - std::cerr << "remain = " << remain << std::endl; -#endif - } - -#if defined(CARVE_DEBUG) - std::cerr << "doTriangulate complete; remain=" << remain << std::endl; -#endif - - if (remain > 3) { -#if defined(CARVE_DEBUG) - std::cerr << "before removeDegeneracies: remain=" << remain << std::endl; -#endif - remain -= removeDegeneracies(begin, result); -#if defined(CARVE_DEBUG) - std::cerr << "after removeDegeneracies: remain=" << remain << std::endl; -#endif - - if (remain > 3) { - return splitAndResume(begin, result); - } - } - - if (remain == 3) { - result.push_back(carve::triangulate::tri_idx(begin->idx, begin->next->idx, begin->next->next->idx)); - } - - vertex_info *d = begin; - do { - vertex_info *n = d->next; - delete d; - d = n; - } while (d != begin); - - return true; -} - - - -static bool testCandidateAttachment(const std::vector<std::vector<carve::geom2d::P2> > &poly, - std::vector<std::pair<size_t, size_t> > ¤t_f_loop, - size_t curr, - carve::geom2d::P2 hole_min) { - const size_t SZ = current_f_loop.size(); - - if (!carve::geom2d::internalToAngle(pvert(poly, current_f_loop[(curr+1) % SZ]), - pvert(poly, current_f_loop[curr]), - pvert(poly, current_f_loop[(curr+SZ-1) % SZ]), - hole_min)) { - return false; - } - - if (hole_min == pvert(poly, current_f_loop[curr])) { - return true; - } - - carve::geom2d::LineSegment2 test(hole_min, pvert(poly, current_f_loop[curr])); - - size_t v1 = current_f_loop.size() - 1; - size_t v2 = 0; - double v1_side = carve::geom2d::orient2d(test.v1, test.v2, pvert(poly, current_f_loop[v1])); - double v2_side = 0; - - while (v2 != current_f_loop.size()) { - v2_side = carve::geom2d::orient2d(test.v1, test.v2, pvert(poly, current_f_loop[v2])); - - if (v1_side != v2_side) { - // XXX: need to test vertices, not indices, because they may - // be duplicated. - if (pvert(poly, current_f_loop[v1]) != pvert(poly, current_f_loop[curr]) && - pvert(poly, current_f_loop[v2]) != pvert(poly, current_f_loop[curr])) { - carve::geom2d::LineSegment2 test2(pvert(poly, current_f_loop[v1]), pvert(poly, current_f_loop[v2])); - if (carve::geom2d::lineSegmentIntersection_simple(test, test2)) { - // intersection; failed. - return false; - } - } - } - - v1 = v2; - v1_side = v2_side; - ++v2; - } - return true; -} - - - -void -carve::triangulate::incorporateHolesIntoPolygon( - const std::vector<std::vector<carve::geom2d::P2> > &poly, - std::vector<std::pair<size_t, size_t> > &result, - size_t poly_loop, - const std::vector<size_t> &hole_loops) { - typedef std::vector<carve::geom2d::P2> loop_t; - - size_t N = poly[poly_loop].size(); - - // work out how much space to reserve for the patched in holes. - for (size_t i = 0; i < hole_loops.size(); i++) { - N += 2 + poly[hole_loops[i]].size(); - } - - // this is the vector that we will build the result in. - result.clear(); - result.reserve(N); - - // this is a heap of result indices that defines the vertex test order. - std::vector<size_t> f_loop_heap; - f_loop_heap.reserve(N); - - // add the poly loop to result. - for (size_t i = 0; i < poly[poly_loop].size(); ++i) { - result.push_back(std::make_pair((size_t)poly_loop, i)); - } - - if (hole_loops.size() == 0) { - return; - } - - std::vector<std::pair<size_t, size_t> > h_loop_min_vertex; - - h_loop_min_vertex.reserve(hole_loops.size()); - - // find the major axis for the holes - this is the axis that we - // will sort on for finding vertices on the polygon to join - // holes up to. - // - // it might also be nice to also look for whether it is better - // to sort ascending or descending. - // - // another trick that could be used is to modify the projection - // by 90 degree rotations or flipping about an axis. just as - // long as we keep the carve::geom3d::Vector pointers for the - // real data in sync, everything should be ok. then we wouldn't - // need to accomodate axes or sort order in the main loop. - - // find the bounding box of all the holes. - carve::geom2d::P2 h_min, h_max; - h_min = h_max = poly[hole_loops[0]][0]; - for (size_t i = 0; i < hole_loops.size(); ++i) { - const loop_t &hole = poly[hole_loops[i]]; - for (size_t j = 0; j < hole.size(); ++j) { - assign_op(h_min, h_min, hole[j], carve::util::min_functor()); - assign_op(h_max, h_max, hole[j], carve::util::max_functor()); - } - } - // choose the axis for which the bbox is largest. - int axis = (h_max.x - h_min.x) > (h_max.y - h_min.y) ? 0 : 1; - - // for each hole, find the minimum vertex in the chosen axis. - for (size_t i = 0; i < hole_loops.size(); ++i) { - const loop_t &hole = poly[hole_loops[i]]; - size_t best, curr; - best = 0; - for (curr = 1; curr != hole.size(); ++curr) { - if (detail::axisOrdering(hole[curr], hole[best], axis)) { - best = curr; - } - } - h_loop_min_vertex.push_back(std::make_pair(hole_loops[i], best)); - } - - // sort the holes by the minimum vertex. - std::sort(h_loop_min_vertex.begin(), h_loop_min_vertex.end(), order_h_loops_2d(poly, axis)); - - // now, for each hole, find a vertex in the current polygon loop that it can be joined to. - for (unsigned i = 0; i < h_loop_min_vertex.size(); ++i) { - // the index of the vertex in the hole to connect. - size_t hole_i = h_loop_min_vertex[i].first; - size_t hole_i_connect = h_loop_min_vertex[i].second; - - carve::geom2d::P2 hole_min = poly[hole_i][hole_i_connect]; - - f_loop_heap.clear(); - // we order polygon loop vertices that may be able to be connected - // to the hole vertex by their distance to the hole vertex - heap_ordering_2d _heap_ordering(poly, result, hole_min, axis); - - const size_t SZ = result.size(); - for (size_t j = 0; j < SZ; ++j) { - // it is guaranteed that there exists a polygon vertex with - // coord < the min hole coord chosen, which can be joined to - // the min hole coord without crossing the polygon - // boundary. also, because we merge holes in ascending - // order, it is also true that this join can never cross - // another hole (and that doesn't need to be tested for). - if (pvert(poly, result[j]).v[axis] <= hole_min.v[axis]) { - f_loop_heap.push_back(j); - std::push_heap(f_loop_heap.begin(), f_loop_heap.end(), _heap_ordering); - } - } - - // we are going to test each potential (according to the - // previous test) polygon vertex as a candidate join. we order - // by closeness to the hole vertex, so that the join we make - // is as small as possible. to test, we need to check the - // joining line segment does not cross any other line segment - // in the current polygon loop (excluding those that have the - // vertex that we are attempting to join with as an endpoint). - size_t attachment_point = result.size(); - - while (f_loop_heap.size()) { - std::pop_heap(f_loop_heap.begin(), f_loop_heap.end(), _heap_ordering); - size_t curr = f_loop_heap.back(); - f_loop_heap.pop_back(); - // test the candidate join from result[curr] to hole_min - - if (!testCandidateAttachment(poly, result, curr, hole_min)) { - continue; - } - - attachment_point = curr; - break; - } - - if (attachment_point == result.size()) { - CARVE_FAIL("didn't manage to link up hole!"); - } - - patchHoleIntoPolygon_2d(result, attachment_point, hole_i, hole_i_connect, poly[hole_i].size()); - } -} - - - -std::vector<std::pair<size_t, size_t> > -carve::triangulate::incorporateHolesIntoPolygon(const std::vector<std::vector<carve::geom2d::P2> > &poly) { -#if 1 - std::vector<std::pair<size_t, size_t> > result; - std::vector<size_t> hole_indices; - hole_indices.reserve(poly.size() - 1); - for (size_t i = 1; i < poly.size(); ++i) { - hole_indices.push_back(i); - } - - incorporateHolesIntoPolygon(poly, result, 0, hole_indices); - - return result; - -#else - typedef std::vector<carve::geom2d::P2> loop_t; - size_t N = poly[0].size(); - // - // work out how much space to reserve for the patched in holes. - for (size_t i = 0; i < poly.size(); i++) { - N += 2 + poly[i].size(); - } - - // this is the vector that we will build the result in. - std::vector<std::pair<size_t, size_t> > current_f_loop; - current_f_loop.reserve(N); - - // this is a heap of current_f_loop indices that defines the vertex test order. - std::vector<size_t> f_loop_heap; - f_loop_heap.reserve(N); - - // add the poly loop to current_f_loop. - for (size_t i = 0; i < poly[0].size(); ++i) { - current_f_loop.push_back(std::make_pair((size_t)0, i)); - } - - if (poly.size() == 1) { - return current_f_loop; - } - - std::vector<std::pair<size_t, size_t> > h_loop_min_vertex; - - h_loop_min_vertex.reserve(poly.size() - 1); - - // find the major axis for the holes - this is the axis that we - // will sort on for finding vertices on the polygon to join - // holes up to. - // - // it might also be nice to also look for whether it is better - // to sort ascending or descending. - // - // another trick that could be used is to modify the projection - // by 90 degree rotations or flipping about an axis. just as - // long as we keep the carve::geom3d::Vector pointers for the - // real data in sync, everything should be ok. then we wouldn't - // need to accomodate axes or sort order in the main loop. - - // find the bounding box of all the holes. - double min_x, min_y, max_x, max_y; - min_x = max_x = poly[1][0].x; - min_y = max_y = poly[1][0].y; - for (size_t i = 1; i < poly.size(); ++i) { - const loop_t &hole = poly[i]; - for (size_t j = 0; j < hole.size(); ++j) { - min_x = std::min(min_x, hole[j].x); - min_y = std::min(min_y, hole[j].y); - max_x = std::max(max_x, hole[j].x); - max_y = std::max(max_y, hole[j].y); - } - } - - // choose the axis for which the bbox is largest. - int axis = (max_x - min_x) > (max_y - min_y) ? 0 : 1; - - // for each hole, find the minimum vertex in the chosen axis. - for (size_t i = 1; i < poly.size(); ++i) { - const loop_t &hole = poly[i]; - size_t best, curr; - best = 0; - for (curr = 1; curr != hole.size(); ++curr) { - if (detail::axisOrdering(hole[curr], hole[best], axis)) { - best = curr; - } - } - h_loop_min_vertex.push_back(std::make_pair(i, best)); - } - - // sort the holes by the minimum vertex. - std::sort(h_loop_min_vertex.begin(), h_loop_min_vertex.end(), order_h_loops_2d(poly, axis)); - - // now, for each hole, find a vertex in the current polygon loop that it can be joined to. - for (unsigned i = 0; i < h_loop_min_vertex.size(); ++i) { - // the index of the vertex in the hole to connect. - size_t hole_i = h_loop_min_vertex[i].first; - size_t hole_i_connect = h_loop_min_vertex[i].second; - - carve::geom2d::P2 hole_min = poly[hole_i][hole_i_connect]; - - f_loop_heap.clear(); - // we order polygon loop vertices that may be able to be connected - // to the hole vertex by their distance to the hole vertex - heap_ordering_2d _heap_ordering(poly, current_f_loop, hole_min, axis); - - const size_t SZ = current_f_loop.size(); - for (size_t j = 0; j < SZ; ++j) { - // it is guaranteed that there exists a polygon vertex with - // coord < the min hole coord chosen, which can be joined to - // the min hole coord without crossing the polygon - // boundary. also, because we merge holes in ascending - // order, it is also true that this join can never cross - // another hole (and that doesn't need to be tested for). - if (pvert(poly, current_f_loop[j]).v[axis] <= hole_min.v[axis]) { - f_loop_heap.push_back(j); - std::push_heap(f_loop_heap.begin(), f_loop_heap.end(), _heap_ordering); - } - } - - // we are going to test each potential (according to the - // previous test) polygon vertex as a candidate join. we order - // by closeness to the hole vertex, so that the join we make - // is as small as possible. to test, we need to check the - // joining line segment does not cross any other line segment - // in the current polygon loop (excluding those that have the - // vertex that we are attempting to join with as an endpoint). - size_t attachment_point = current_f_loop.size(); - - while (f_loop_heap.size()) { - std::pop_heap(f_loop_heap.begin(), f_loop_heap.end(), _heap_ordering); - size_t curr = f_loop_heap.back(); - f_loop_heap.pop_back(); - // test the candidate join from current_f_loop[curr] to hole_min - - if (!testCandidateAttachment(poly, current_f_loop, curr, hole_min)) { - continue; - } - - attachment_point = curr; - break; - } - - if (attachment_point == current_f_loop.size()) { - CARVE_FAIL("didn't manage to link up hole!"); - } - - patchHoleIntoPolygon_2d(current_f_loop, attachment_point, hole_i, hole_i_connect, poly[hole_i].size()); - } - - return current_f_loop; -#endif -} - - - -std::vector<std::vector<std::pair<size_t, size_t> > > -carve::triangulate::mergePolygonsAndHoles(const std::vector<std::vector<carve::geom2d::P2> > &poly) { - std::vector<size_t> poly_indices, hole_indices; - - poly_indices.reserve(poly.size()); - hole_indices.reserve(poly.size()); - - for (size_t i = 0; i < poly.size(); ++i) { - if (carve::geom2d::signedArea(poly[i]) < 0) { - poly_indices.push_back(i); - } else { - hole_indices.push_back(i); - } - } - - std::vector<std::vector<std::pair<size_t, size_t> > > result; - result.resize(poly_indices.size()); - - if (hole_indices.size() == 0) { - for (size_t i = 0; i < poly.size(); ++i) { - result[i].resize(poly[i].size()); - for (size_t j = 0; j < poly[i].size(); ++j) { - result[i].push_back(std::make_pair(i, j)); - } - } - return result; - } - - if (poly_indices.size() == 1) { - incorporateHolesIntoPolygon(poly, result[0], poly_indices[0], hole_indices); - - return result; - } - - throw carve::exception("not implemented"); -} - - - -void carve::triangulate::triangulate(const std::vector<carve::geom2d::P2> &poly, - std::vector<carve::triangulate::tri_idx> &result) { - std::vector<detail::vertex_info *> vinfo; - const size_t N = poly.size(); - -#if defined(CARVE_DEBUG) - std::cerr << "TRIANGULATION BEGINS" << std::endl; -#endif - -#if defined(CARVE_DEBUG_WRITE_PLY_DATA) - dumpPoly(poly, result); -#endif - - result.clear(); - if (N < 3) { - return; - } - - result.reserve(poly.size() - 2); - - if (N == 3) { - result.push_back(tri_idx(0, 1, 2)); - return; - } - - vinfo.resize(N); - - vinfo[0] = new detail::vertex_info(poly[0], 0); - for (size_t i = 1; i < N-1; ++i) { - vinfo[i] = new detail::vertex_info(poly[i], i); - vinfo[i]->prev = vinfo[i-1]; - vinfo[i-1]->next = vinfo[i]; - } - vinfo[N-1] = new detail::vertex_info(poly[N-1], N-1); - vinfo[N-1]->prev = vinfo[N-2]; - vinfo[N-1]->next = vinfo[0]; - vinfo[0]->prev = vinfo[N-1]; - vinfo[N-2]->next = vinfo[N-1]; - - for (size_t i = 0; i < N; ++i) { - vinfo[i]->recompute(); - } - - detail::vertex_info *begin = vinfo[0]; - - removeDegeneracies(begin, result); - doTriangulate(begin, result); - -#if defined(CARVE_DEBUG) - std::cerr << "TRIANGULATION ENDS" << std::endl; -#endif - -#if defined(CARVE_DEBUG_WRITE_PLY_DATA) - dumpPoly(poly, result); -#endif -} - - - -void carve::triangulate::detail::tri_pair_t::flip(vert_edge_t &old_edge, - vert_edge_t &new_edge, - vert_edge_t perim[4]) { - unsigned ai, bi; - unsigned cross_ai, cross_bi; - - findSharedEdge(ai, bi); - old_edge = ordered_vert_edge_t(a->v[ai], b->v[bi]); - - cross_ai = P(ai); - cross_bi = P(bi); - new_edge = ordered_vert_edge_t(a->v[cross_ai], b->v[cross_bi]); - - score = -score; - - a->v[N(ai)] = b->v[cross_bi]; - b->v[N(bi)] = a->v[cross_ai]; - - perim[0] = ordered_vert_edge_t(a->v[P(ai)], a->v[ai]); - perim[1] = ordered_vert_edge_t(a->v[N(ai)], a->v[ai]); // this edge was a b-edge - - perim[2] = ordered_vert_edge_t(b->v[P(bi)], b->v[bi]); - perim[3] = ordered_vert_edge_t(b->v[N(bi)], b->v[bi]); // this edge was an a-edge -} - - - -void carve::triangulate::detail::tri_pairs_t::insert(unsigned a, unsigned b, carve::triangulate::tri_idx *t) { - tri_pair_t *tp; - if (a < b) { - tp = storage[vert_edge_t(a,b)]; - if (!tp) { - tp = storage[vert_edge_t(a,b)] = new tri_pair_t; - } - tp->a = t; - } else { - tp = storage[vert_edge_t(b,a)]; - if (!tp) { - tp = storage[vert_edge_t(b,a)] = new tri_pair_t; - } - tp->b = t; - } -} |