blender/extern/carve/include/carve/triangulator_impl.hpp
Sergey Sharybin e81f2853c8 Carve booleans library integration
==================================

Merging Carve library integration project into the trunk.

This commit switches Boolean modifier to another library which handles
mesh boolean operations in much stable and faster way, resolving old
well-known limitations of intern boolop library.

Carve is integrating as alternative interface for boolop library and
which makes it totally transparent for blender sources to switch between
old-fashioned boolop and new Carve backends.

Detailed changes in this commit:

- Integrated needed subset of Carve library sources into extern/
  Added script for re-bundling it (currently works only if repo
  was cloned by git-svn).
- Added BOP_CarveInterface for boolop library which can be used by
  Boolean modifier.
- Carve backend is enabled by default, can be disabled by WITH_BF_CARVE
  SCons option and WITH_CARVE CMake option.
- If Boost library is found in build environment it'll be used for
  unordered collections. If Boost isn't found, it'll fallback to TR1
  implementation for GCC compilers. Boost is obligatory if MSVC is used.

Tested on Linux 64bit and Windows 7 64bit.

NOTE: behavior of flat objects was changed. E.g. Plane-Sphere now gives
      plane with circle hole, not plane with semisphere. Don't think
      it's really issue because it's not actually defined behavior in
      such situations and both of ways might be useful. Since it's
      only known "regression" think it's OK to deal with it.

Details are there http://wiki.blender.org/index.php/User:Nazg-gul/CarveBooleans

Special thanks to:

- Ken Hughes: author of original carve integration patch.
- Campbell Barton: help in project development, review tests.
- Tobias Sargeant: author of Carve library, help in resolving some
                   merge stoppers, bug fixing.
2012-01-16 16:46:00 +00:00

852 lines
28 KiB
C++

// Begin License:
// Copyright (C) 2006-2011 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 the GNU General Public
// License version 2.0 as published by the Free Software Foundation
// and appearing in the file LICENSE.GPL2 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:
#pragma once
#include <carve/geom2d.hpp>
#if defined(CARVE_DEBUG)
# include <iostream>
#endif
namespace carve {
namespace triangulate {
namespace detail {
static inline bool axisOrdering(const carve::geom2d::P2 &a,
const carve::geom2d::P2 &b,
int axis) {
return a.v[axis] < b.v[axis] || (a.v[axis] == b.v[axis] && a.v[1-axis] < b.v[1-axis]);
}
/**
* \class order_h_loops
* \brief Provides an ordering of hole loops based upon a single
* projected axis.
*
* @tparam project_t A functor which converts vertices to a 2d
* projection.
* @tparam hole_t A collection of vertices.
*/
template<typename project_t, typename vert_t>
class order_h_loops {
const project_t &project;
int axis;
public:
/**
*
* @param _project The projection functor.
* @param _axis The axis of the 2d projection upon which hole
* loops are ordered.
*/
order_h_loops(const project_t &_project, int _axis) : project(_project), axis(_axis) { }
bool operator()(const vert_t &a,
const vert_t &b) const {
return axisOrdering(project(a), project(b), axis);
}
bool operator()(
const std::pair<const typename std::vector<vert_t> *, typename std::vector<vert_t>::const_iterator> &a,
const std::pair<const typename std::vector<vert_t> *, typename std::vector<vert_t>::const_iterator> &b) {
return axisOrdering(project(*(a.second)), project(*(b.second)), axis);
}
};
/**
* \class heap_ordering
* \brief Provides an ordering of vertex indicies in a polygon
* loop according to proximity to a vertex.
*
* @tparam project_t A functor which converts vertices to a 2d
* projection.
* @tparam vert_t A vertex type.
*/
template<typename project_t, typename vert_t>
class heap_ordering {
const project_t &project;
const std::vector<vert_t> &loop;
const carve::geom2d::P2 p;
int axis;
public:
/**
*
* @param _project A functor which converts vertices to a 2d
* projection.
* @param _loop The polygon loop which indices address.
* @param _vert The vertex from which distance is measured.
*
*/
heap_ordering(const project_t &_project,
const std::vector<vert_t> &_loop,
vert_t _vert,
int _axis) :
project(_project),
loop(_loop),
p(_project(_vert)),
axis(_axis) {
}
bool operator()(size_t a, size_t b) const {
carve::geom2d::P2 pa = project(loop[a]);
carve::geom2d::P2 pb = project(loop[b]);
double da = carve::geom::distance2(p, pa);
double db = carve::geom::distance2(p, pb);
if (da > db) return true;
if (da < db) return false;
return axisOrdering(pa, pb, axis);
}
};
/**
* \brief Given a polygon loop and a hole loop, and attachment
* points, insert the hole loop vertices into the polygon loop.
*
* @param[in,out] f_loop The polygon loop to incorporate the
* hole into.
* @param f_loop_attach[in] The index of the vertex of the
* polygon loop that the hole is to be
* attached to.
* @param hole_attach[in] A pair consisting of a pointer to a
* hole container and an iterator into
* that container reflecting the point of
* attachment of the hole.
*/
template<typename vert_t>
void patchHoleIntoPolygon(std::vector<vert_t> &f_loop,
unsigned f_loop_attach,
const std::pair<const std::vector<vert_t> *,
typename std::vector<vert_t>::const_iterator> &hole_attach) {
// join the vertex curr of the polygon loop to the hole at
// h_loop_connect
f_loop.insert(f_loop.begin() + f_loop_attach + 1, hole_attach.first->size() + 2, NULL);
typename std::vector<vert_t>::iterator f = f_loop.begin() + f_loop_attach;
typename std::vector<vert_t>::const_iterator h = hole_attach.second;
while (h != hole_attach.first->end()) {
*++f = *h++;
}
h = hole_attach.first->begin();
typename std::vector<vert_t>::const_iterator he = hole_attach.second; ++he;
while (h != he) {
*++f = *h++;
}
*++f = f_loop[f_loop_attach];
}
struct vertex_info;
/**
* \brief Determine whether c is to the left of a->b.
*/
static inline bool isLeft(const vertex_info *a,
const vertex_info *b,
const vertex_info *c);
/**
* \brief Determine whether d is contained in the triangle abc.
*/
static inline bool pointInTriangle(const vertex_info *a,
const vertex_info *b,
const vertex_info *c,
const vertex_info *d);
/**
* \class vertex_info
* \brief Maintains a linked list of untriangulated vertices
* during a triangulation operation.
*/
struct vertex_info {
vertex_info *prev;
vertex_info *next;
carve::geom2d::P2 p;
size_t idx;
double score;
bool convex;
bool failed;
vertex_info(const carve::geom2d::P2 &_p, size_t _idx) :
prev(NULL), next(NULL),
p(_p), idx(_idx),
score(0.0), convex(false) {
}
static double triScore(const vertex_info *p, const vertex_info *v, const vertex_info *n);
double calcScore() const;
void recompute() {
score = calcScore();
convex = isLeft(prev, this, next);
failed = false;
}
bool isCandidate() const {
return convex && !failed;
}
void remove() {
next->prev = prev;
prev->next = next;
}
bool isClipable() const;
};
static inline bool isLeft(const vertex_info *a,
const vertex_info *b,
const vertex_info *c) {
if (a->idx < b->idx && b->idx < c->idx) {
return carve::geom2d::orient2d(a->p, b->p, c->p) > 0.0;
} else if (a->idx < c->idx && c->idx < b->idx) {
return carve::geom2d::orient2d(a->p, c->p, b->p) < 0.0;
} else if (b->idx < a->idx && a->idx < c->idx) {
return carve::geom2d::orient2d(b->p, a->p, c->p) < 0.0;
} else if (b->idx < c->idx && c->idx < a->idx) {
return carve::geom2d::orient2d(b->p, c->p, a->p) > 0.0;
} else if (c->idx < a->idx && a->idx < b->idx) {
return carve::geom2d::orient2d(c->p, a->p, b->p) > 0.0;
} else {
return carve::geom2d::orient2d(c->p, b->p, a->p) < 0.0;
}
}
static inline bool pointInTriangle(const vertex_info *a,
const vertex_info *b,
const vertex_info *c,
const vertex_info *d) {
return !isLeft(a, c, d) && !isLeft(b, a, d) && !isLeft(c, b, d);
}
size_t removeDegeneracies(vertex_info *&begin, std::vector<carve::triangulate::tri_idx> &result);
bool splitAndResume(vertex_info *begin, std::vector<carve::triangulate::tri_idx> &result);
bool doTriangulate(vertex_info *begin, std::vector<carve::triangulate::tri_idx> &result);
typedef std::pair<unsigned, unsigned> vert_edge_t;
struct hash_vert_edge_t {
size_t operator()(const vert_edge_t &e) const {
size_t r = (size_t)e.first;
size_t s = (size_t)e.second;
return r ^ ((s >> 16) | (s << 16));
}
};
static inline vert_edge_t ordered_vert_edge_t(unsigned a, unsigned b) {
return (a < b) ? vert_edge_t(a, b) : vert_edge_t(b, a);
}
struct tri_pair_t {
carve::triangulate::tri_idx *a, *b;
double score;
size_t idx;
tri_pair_t() : a(NULL), b(NULL), score(0.0) {
}
static inline unsigned N(unsigned i) { return (i+1)%3; }
static inline unsigned P(unsigned i) { return (i+2)%3; }
void findSharedEdge(unsigned &ai, unsigned &bi) const {
if (a->v[1] == b->v[0]) { if (a->v[0] == b->v[1]) { ai = 0; bi = 0; } else { ai = 1; bi = 2; } return; }
if (a->v[1] == b->v[1]) { if (a->v[0] == b->v[2]) { ai = 0; bi = 1; } else { ai = 1; bi = 0; } return; }
if (a->v[1] == b->v[2]) { if (a->v[0] == b->v[0]) { ai = 0; bi = 2; } else { ai = 1; bi = 1; } return; }
if (a->v[2] == b->v[0]) { ai = 2; bi = 2; return; }
if (a->v[2] == b->v[1]) { ai = 2; bi = 0; return; }
if (a->v[2] == b->v[2]) { ai = 2; bi = 1; return; }
CARVE_FAIL("should not be reached");
}
void flip(vert_edge_t &old_edge,
vert_edge_t &new_edge,
vert_edge_t perim[4]);
template<typename project_t, typename vert_t, typename distance_calc_t>
double calc(const project_t &project,
const std::vector<vert_t> &poly,
distance_calc_t dist) {
unsigned ai, bi;
unsigned cross_ai, cross_bi;
unsigned ea, eb;
findSharedEdge(ai, bi);
#if defined(CARVE_DEBUG)
if (carve::geom2d::signedArea(project(poly[a->v[0]]), project(poly[a->v[1]]), project(poly[a->v[2]])) > 0.0 ||
carve::geom2d::signedArea(project(poly[b->v[0]]), project(poly[b->v[1]]), project(poly[b->v[2]])) > 0.0) {
std::cerr << "warning: triangle pair " << this << " contains triangles with incorrect orientation" << std::endl;
}
#endif
cross_ai = P(ai);
cross_bi = P(bi);
ea = a->v[cross_ai];
eb = b->v[cross_bi];
double side_1 = carve::geom2d::orient2d(project(poly[ea]), project(poly[eb]), project(poly[a->v[ai]]));
double side_2 = carve::geom2d::orient2d(project(poly[ea]), project(poly[eb]), project(poly[a->v[N(ai)]]));
bool can_flip = (side_1 < 0.0 && side_2 > 0.0) || (side_1 > 0.0 && side_2 < 0.0);
if (!can_flip) {
score = -1;
} else {
score =
dist(poly[a->v[ai]], poly[b->v[bi]]) -
dist(poly[a->v[cross_ai]], poly[b->v[cross_bi]]);
}
return score;
}
template<typename project_t, typename vert_t, typename distance_calc_t>
double edgeLen(const project_t &project,
const std::vector<vert_t> &poly,
distance_calc_t dist) const {
unsigned ai, bi;
findSharedEdge(ai, bi);
return dist(poly[a->v[ai]], poly[b->v[bi]]);
}
};
struct max_score {
bool operator()(const tri_pair_t *a, const tri_pair_t *b) const { return a->score < b->score; }
};
struct tri_pairs_t {
typedef std::unordered_map<vert_edge_t, tri_pair_t *, hash_vert_edge_t> storage_t;
storage_t storage;
tri_pairs_t() : storage() {
};
~tri_pairs_t() {
for (storage_t::iterator i = storage.begin(); i != storage.end(); ++i) {
if ((*i).second) delete (*i).second;
}
}
void insert(unsigned a, unsigned b, carve::triangulate::tri_idx *t);
template<typename project_t, typename vert_t, typename distance_calc_t>
void updateEdge(tri_pair_t *tp,
const project_t &project,
const std::vector<vert_t> &poly,
distance_calc_t dist,
std::vector<tri_pair_t *> &edges,
size_t &n) {
double old_score = tp->score;
double new_score = tp->calc(project, poly, dist);
#if defined(CARVE_DEBUG)
std::cerr << "tp:" << tp << " old_score: " << old_score << " new_score: " << new_score << std::endl;
#endif
if (new_score > 0.0 && old_score <= 0.0) {
tp->idx = n;
edges[n++] = tp;
} else if (new_score <= 0.0 && old_score > 0.0) {
std::swap(edges[tp->idx], edges[--n]);
edges[tp->idx]->idx = tp->idx;
}
}
tri_pair_t *get(vert_edge_t &e) {
storage_t::iterator i;
i = storage.find(e);
if (i == storage.end()) return NULL;
return (*i).second;
}
template<typename project_t, typename vert_t, typename distance_calc_t>
void flip(const project_t &project,
const std::vector<vert_t> &poly,
distance_calc_t dist,
std::vector<tri_pair_t *> &edges,
size_t &n) {
vert_edge_t old_e, new_e;
vert_edge_t perim[4];
#if defined(CARVE_DEBUG)
std::cerr << "improvable edges: " << n << std::endl;
#endif
tri_pair_t *tp = *std::max_element(edges.begin(), edges.begin() + n, max_score());
#if defined(CARVE_DEBUG)
std::cerr << "improving tri-pair: " << tp << " with score: " << tp->score << std::endl;
#endif
tp->flip(old_e, new_e, perim);
#if defined(CARVE_DEBUG)
std::cerr << "old_e: " << old_e.first << "," << old_e.second << " -> new_e: " << new_e.first << "," << new_e.second << std::endl;
#endif
CARVE_ASSERT(storage.find(old_e) != storage.end());
storage.erase(old_e);
storage[new_e] = tp;
std::swap(edges[tp->idx], edges[--n]);
edges[tp->idx]->idx = tp->idx;
tri_pair_t *tp2;
tp2 = get(perim[0]);
if (tp2 != NULL) {
updateEdge(tp2, project, poly, dist, edges, n);
}
tp2 = get(perim[1]);
if (tp2 != NULL) {
CARVE_ASSERT(tp2->a == tp->b || tp2->b == tp->b);
if (tp2->a == tp->b) { tp2->a = tp->a; } else { tp2->b = tp->a; }
updateEdge(tp2, project, poly, dist, edges, n);
}
tp2 = get(perim[2]);
if (tp2 != NULL) {
updateEdge(tp2, project, poly, dist, edges, n);
}
tp2 = get(perim[3]);
if (tp2 != NULL) {
CARVE_ASSERT(tp2->a == tp->a || tp2->b == tp->a);
if (tp2->a == tp->a) { tp2->a = tp->b; } else { tp2->b = tp->b; }
updateEdge(tp2, project, poly, dist, edges, n);
}
}
template<typename project_t, typename vert_t, typename distance_calc_t>
size_t getInternalEdges(const project_t &project,
const std::vector<vert_t> &poly,
distance_calc_t dist,
std::vector<tri_pair_t *> &edges) {
size_t count = 0;
for (storage_t::iterator i = storage.begin(); i != storage.end();) {
tri_pair_t *tp = (*i).second;
if (tp->a && tp->b) {
tp->calc(project, poly, dist);
count++;
#if defined(CARVE_DEBUG)
std::cerr << "internal edge: " << (*i).first.first << "," << (*i).first.second << " -> " << tp << " " << tp->score << std::endl;
#endif
++i;
} else {
delete (*i).second;
storage.erase(i++);
}
}
edges.resize(count);
size_t fwd = 0;
size_t rev = count;
for (storage_t::iterator i = storage.begin(); i != storage.end(); ++i) {
tri_pair_t *tp = (*i).second;
if (tp && tp->a && tp->b) {
if (tp->score > 0.0) {
edges[fwd++] = tp;
} else {
edges[--rev] = tp;
}
}
}
CARVE_ASSERT(fwd == rev);
return fwd;
}
};
template<typename project_t, typename vert_t>
static bool
testCandidateAttachment(const project_t &project,
std::vector<vert_t> &current_f_loop,
size_t curr,
carve::geom2d::P2 hole_min) {
const size_t SZ = current_f_loop.size();
size_t prev, next;
if (curr == 0) {
prev = SZ - 1; next = 1;
} else if (curr == SZ - 1) {
prev = curr - 1; next = 0;
} else {
prev = curr - 1; next = curr + 1;
}
if (!carve::geom2d::internalToAngle(project(current_f_loop[next]),
project(current_f_loop[curr]),
project(current_f_loop[prev]),
hole_min)) {
return false;
}
if (hole_min == project(current_f_loop[curr])) {
return true;
}
carve::geom2d::LineSegment2 test(hole_min, project(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, project(current_f_loop[v1]));
double v2_side = 0;
while (v2 != current_f_loop.size()) {
v2_side = carve::geom2d::orient2d(test.v1, test.v2, project(current_f_loop[v2]));
if (v1_side != v2_side) {
// XXX: need to test vertices, not indices, because they may
// be duplicated.
if (project(current_f_loop[v1]) != project(current_f_loop[curr]) &&
project(current_f_loop[v2]) != project(current_f_loop[curr])) {
carve::geom2d::LineSegment2 test2(project(current_f_loop[v1]), project(current_f_loop[v2]));
if (carve::geom2d::lineSegmentIntersection_simple(test, test2)) {
// intersection; failed.
return false;
}
}
}
v1 = v2;
v1_side = v2_side;
++v2;
}
return true;
}
}
template<typename project_t, typename vert_t>
static std::vector<vert_t>
incorporateHolesIntoPolygon(const project_t &project,
const std::vector<vert_t> &f_loop,
const std::vector<std::vector<vert_t> > &h_loops) {
typedef std::vector<vert_t> hole_t;
typedef typename std::vector<vert_t>::const_iterator vert_iter;
typedef typename std::vector<std::vector<vert_t> >::const_iterator hole_iter;
size_t N = f_loop.size();
// work out how much space to reserve for the patched in holes.
for (hole_iter i = h_loops.begin(); i != h_loops.end(); ++i) {
N += 2 + (*i).size();
}
// this is the vector that we will build the result in.
std::vector<vert_t> current_f_loop;
current_f_loop.reserve(N);
std::vector<size_t> f_loop_heap;
f_loop_heap.reserve(N);
for (unsigned i = 0; i < f_loop.size(); ++i) {
current_f_loop.push_back(f_loop[i]);
}
std::vector<std::pair<const std::vector<vert_t> *, vert_iter> > h_loop_min_vertex;
h_loop_min_vertex.reserve(h_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.
bool first = true;
double min_x, min_y, max_x, max_y;
for (hole_iter i = h_loops.begin(); i != h_loops.end(); ++i) {
const hole_t &hole(*i);
for (vert_iter j = hole.begin(); j != hole.end(); ++j) {
carve::geom2d::P2 curr = project(*j);
if (first) {
min_x = max_x = curr.x;
min_y = max_y = curr.y;
first = false;
} else {
min_x = std::min(min_x, curr.x);
min_y = std::min(min_y, curr.y);
max_x = std::max(max_x, curr.x);
max_y = std::max(max_y, curr.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 (hole_iter i = h_loops.begin(); i != h_loops.end(); ++i) {
const hole_t &hole = *i;
vert_iter best_i = std::min_element(hole.begin(), hole.end(), detail::order_h_loops<project_t, vert_t>(project, axis));
h_loop_min_vertex.push_back(std::make_pair(&hole, best_i));
}
// sort the holes by the minimum vertex.
std::sort(h_loop_min_vertex.begin(), h_loop_min_vertex.end(), detail::order_h_loops<project_t, vert_t>(project, 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) {
const size_t N_f_loop = current_f_loop.size();
// the index of the vertex in the hole to connect.
vert_iter h_loop_connect = h_loop_min_vertex[i].second;
carve::geom2d::P2 hole_min = project(*h_loop_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
detail::heap_ordering<project_t, vert_t> _heap_ordering(project, current_f_loop, *h_loop_connect, axis);
for (size_t j = 0; j < N_f_loop; ++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 (project(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 (!detail::testCandidateAttachment(project, 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!");
}
detail::patchHoleIntoPolygon(current_f_loop, attachment_point, h_loop_min_vertex[i]);
}
return current_f_loop;
}
template<typename project_t, typename vert_t>
void triangulate(const project_t &project,
const std::vector<vert_t> &poly,
std::vector<tri_idx> &result) {
std::vector<detail::vertex_info *> vinfo;
const size_t N = poly.size();
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(project(poly[0]), 0);
for (size_t i = 1; i < N-1; ++i) {
vinfo[i] = new detail::vertex_info(project(poly[i]), i);
vinfo[i]->prev = vinfo[i-1];
vinfo[i-1]->next = vinfo[i];
}
vinfo[N-1] = new detail::vertex_info(project(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);
}
template<typename project_t, typename vert_t, typename distance_calc_t>
void improve(const project_t &project,
const std::vector<vert_t> &poly,
distance_calc_t dist,
std::vector<tri_idx> &result) {
detail::tri_pairs_t tri_pairs;
#if defined(CARVE_DEBUG)
bool warn = false;
for (size_t i = 0; i < result.size(); ++i) {
tri_idx &t = result[i];
if (carve::geom2d::signedArea(project(poly[t.a]), project(poly[t.b]), project(poly[t.c])) > 0) {
warn = true;
}
}
if (warn) {
std::cerr << "carve::triangulate::improve(): Some triangles are incorrectly oriented. Results may be incorrect." << std::endl;
}
#endif
for (size_t i = 0; i < result.size(); ++i) {
tri_idx &t = result[i];
tri_pairs.insert(t.a, t.b, &t);
tri_pairs.insert(t.b, t.c, &t);
tri_pairs.insert(t.c, t.a, &t);
}
std::vector<detail::tri_pair_t *> edges;
size_t n = tri_pairs.getInternalEdges(project, poly, dist, edges);
for (size_t i = 0; i < n; ++i) {
edges[i]->idx = i;
}
// procedure:
// while a tri pair with a positive score exists:
// p = pair with highest positive score
// flip p, rewriting its two referenced triangles.
// negate p's score
// for each q in the up-to-four adjoining tri pairs:
// update q's tri ptr, if changed, and its score.
#if defined(CARVE_DEBUG)
double initial_score = 0;
for (size_t i = 0; i < edges.size(); ++i) {
initial_score += edges[i]->edgeLen(project, poly, dist);
}
std::cerr << "initial score: " << initial_score << std::endl;
#endif
while (n) {
tri_pairs.flip(project, poly, dist, edges, n);
}
#if defined(CARVE_DEBUG)
double final_score = 0;
for (size_t i = 0; i < edges.size(); ++i) {
final_score += edges[i]->edgeLen(project, poly, dist);
}
std::cerr << "final score: " << final_score << std::endl;
#endif
#if defined(CARVE_DEBUG)
if (!warn) {
for (size_t i = 0; i < result.size(); ++i) {
tri_idx &t = result[i];
CARVE_ASSERT (carve::geom2d::signedArea(project(poly[t.a]), project(poly[t.b]), project(poly[t.c])) <= 0.0);
}
}
#endif
}
template<typename project_t, typename vert_t>
void improve(const project_t &project,
const std::vector<vert_t> &poly,
std::vector<tri_idx> &result) {
improve(project, poly, carve::geom::distance_functor(), result);
}
}
}