Faster minkowski()

master
Oskar Linde 2014-06-03 13:18:40 +02:00
parent 8c977a45e5
commit cb2c5029c0
3 changed files with 215 additions and 0 deletions

View File

@ -110,6 +110,8 @@ GeometryEvaluator::ResultObject GeometryEvaluator::applyToChildren3D(const Abstr
// Only one child -> this is a noop
if (children.size() == 1) return ResultObject(children.front().second);
if (op == OPENSCAD_MINKOWSKI) return ResultObject(CGALUtils::applyMinkowski(children));
CGAL_Nef_polyhedron *N = new CGAL_Nef_polyhedron;
CGALUtils::applyOperator(children, *N, op);
return ResultObject(N);

View File

@ -65,6 +65,216 @@ namespace CGALUtils {
return false;
}
}
template<typename Polyhedron>
bool is_weakly_convex(Polyhedron const& p) {
for (typename Polyhedron::Edge_const_iterator i = p.edges_begin(); i != p.edges_end(); ++i) {
typename Polyhedron::Plane_3 p(i->opposite()->vertex()->point(), i->vertex()->point(), i->next()->vertex()->point());
if (p.has_on_positive_side(i->opposite()->next()->vertex()->point()) &&
CGAL::squared_distance(p, i->opposite()->next()->vertex()->point()) > 1e-8) {
return false;
}
}
return true;
}
Geometry const * applyMinkowski(const Geometry::ChildList &children)
{
CGAL::Timer t,t_tot;
assert(children.size() >= 2);
Geometry::ChildList::const_iterator it = children.begin();
t_tot.start();
Geometry const* operands[2] = {it->second.get(), nullptr};
try {
while (++it != children.end()) {
operands[1] = it->second.get();
typedef CGAL::Exact_predicates_inexact_constructions_kernel Hull_kernel;
std::list<CGAL_Polyhedron> P[2];
std::list<CGAL::Polyhedron_3<Hull_kernel> > result_parts;
for (int i = 0; i < 2; i++) {
CGAL_Polyhedron poly;
const PolySet * ps = dynamic_cast<const PolySet *>(operands[i]);
const CGAL_Nef_polyhedron * nef = dynamic_cast<const CGAL_Nef_polyhedron *>(operands[i]);
if (ps) createPolyhedronFromPolySet(*ps, poly);
else if (nef && nef->p3->is_simple()) nefworkaround::convert_to_Polyhedron<CGAL_Kernel3>(*nef->p3, poly);
else throw 0;
if (ps && ps->is_convex() || !ps && is_weakly_convex(poly)) {
PRINTDB("Minkowski: child %d is convex and %s",i % (ps?"PolySet":"Nef") );
P[i].push_back(poly);
} else {
CGAL_Nef_polyhedron3 decomposed_nef;
if (ps) {
PRINTDB("Minkowski: child %d is nonconvex PolySet, transforming to Nef and decomposing...", i);
CGAL_Nef_polyhedron *p = createNefPolyhedronFromGeometry(*ps);
decomposed_nef = *p->p3;
delete p;
} else {
PRINTDB("Minkowski: child %d is nonconvex Nef, decomposing...",i);
decomposed_nef = *nef->p3;
}
CGAL::convex_decomposition_3(decomposed_nef);
// the first volume is the outer volume, which ignored in the decomposition
CGAL_Nef_polyhedron3::Volume_const_iterator ci = ++decomposed_nef.volumes_begin();
for( ; ci != decomposed_nef.volumes_end(); ++ci) {
if(ci->mark()) {
CGAL_Polyhedron poly;
decomposed_nef.convert_inner_shell_to_polyhedron(ci->shells_begin(), poly);
P[i].push_back(poly);
}
}
PRINTDB("Minkowski: decomposed into %d convex parts", P[i].size());
}
}
std::vector<Hull_kernel::Point_3> points[2];
std::vector<Hull_kernel::Point_3> minkowski_points;
for (int i = 0; i < P[0].size(); i++) {
for (int j = 0; j < P[1].size(); j++) {
t.start();
points[0].clear();
points[1].clear();
for (int k = 0; k < 2; k++) {
std::list<CGAL_Polyhedron>::iterator it = P[k].begin();
std::advance(it, k==0?i:j);
CGAL_Polyhedron const& poly = *it;
points[k].reserve(poly.size_of_vertices());
for (CGAL_Polyhedron::Vertex_const_iterator pi = poly.vertices_begin(); pi != poly.vertices_end(); ++pi) {
CGAL_Polyhedron::Point_3 const& p = pi->point();
points[k].push_back(Hull_kernel::Point_3(to_double(p[0]),to_double(p[1]),to_double(p[2])));
}
}
minkowski_points.clear();
minkowski_points.reserve(points[0].size() * points[1].size());
for (int i = 0; i < points[0].size(); i++) {
for (int j = 0; j < points[1].size(); j++) {
minkowski_points.push_back(points[0][i]+(points[1][j]-CGAL::ORIGIN));
}
}
if (minkowski_points.size() <= 3) {
t.stop();
continue;
}
CGAL::Polyhedron_3<Hull_kernel> result;
t.stop();
PRINTDB("Minkowski: Point cloud creation (%d ⨉ %d -> %d) took %f ms", points[0].size() % points[1].size() % minkowski_points.size() % (t.time()*1000));
t.reset();
t.start();
CGAL::convex_hull_3(minkowski_points.begin(), minkowski_points.end(), result);
std::vector<Hull_kernel::Point_3> strict_points;
strict_points.reserve(minkowski_points.size());
for (CGAL::Polyhedron_3<Hull_kernel>::Vertex_iterator i = result.vertices_begin(); i != result.vertices_end(); ++i) {
Hull_kernel::Point_3 const& p = i->point();
CGAL::Polyhedron_3<Hull_kernel>::Vertex::Halfedge_handle h,e;
h = i->halfedge();
e = h;
bool collinear = false;
bool coplanar = true;
do {
Hull_kernel::Point_3 const& q = h->opposite()->vertex()->point();
if (coplanar && !CGAL::coplanar(p,q,
h->next_on_vertex()->opposite()->vertex()->point(),
h->next_on_vertex()->next_on_vertex()->opposite()->vertex()->point())) {
coplanar = false;
}
for (CGAL::Polyhedron_3<Hull_kernel>::Vertex::Halfedge_handle j = h->next_on_vertex();
j != h && !collinear && ! coplanar;
j = j->next_on_vertex()) {
Hull_kernel::Point_3 const& r = j->opposite()->vertex()->point();
if (CGAL::collinear(p,q,r)) {
collinear = true;
}
}
h = h->next_on_vertex();
} while (h != e && !collinear);
if (!collinear && !coplanar)
strict_points.push_back(p);
}
result.clear();
CGAL::convex_hull_3(strict_points.begin(), strict_points.end(), result);
t.stop();
PRINTDB("Minkowski: Computing convex hull took %f s", t.time());
t.reset();
result_parts.push_back(result);
}
}
if (it != std::next(children.begin()))
delete operands[0];
if (result_parts.size() == 1) {
PolySet *ps = new PolySet(3,true);
createPolySetFromPolyhedron(*result_parts.begin(), *ps);
operands[0] = ps;
} else if (!result_parts.empty()) {
t.start();
PRINTDB("Minkowski: Computing union of %d parts",result_parts.size());
Geometry::ChildList fake_children;
for (std::list<CGAL::Polyhedron_3<Hull_kernel> >::iterator i = result_parts.begin(); i != result_parts.end(); ++i) {
PolySet ps(3,true);
createPolySetFromPolyhedron(*i, ps);
fake_children.push_back(std::make_pair((const AbstractNode*)NULL,
shared_ptr<const Geometry>(createNefPolyhedronFromGeometry(ps))));
}
CGAL_Nef_polyhedron *N = new CGAL_Nef_polyhedron;
CGALUtils::applyOperator(fake_children, *N, OPENSCAD_UNION);
t.stop();
PRINTDB("Minkowski: Union done: %f s",t.time());
t.reset();
operands[0] = N;
} else {
return NULL;
}
}
t_tot.stop();
PRINTDB("Minkowski: Total execution time %f s", t_tot.time());
t_tot.reset();
return operands[0];
}
catch (...) {
// If anything throws we simply fall back to Nef Minkowski
PRINTD("Minkowski: Falling back to Nef Minkowski");
CGAL_Nef_polyhedron *N = new CGAL_Nef_polyhedron;
applyOperator(children, *N, OPENSCAD_MINKOWSKI);
return N;
}
}
/*!
Applies op to all children and stores the result in dest.
@ -451,6 +661,8 @@ bool createPolySetFromPolyhedron(const Polyhedron &p, PolySet &ps)
template bool createPolySetFromPolyhedron(const CGAL_Polyhedron &p, PolySet &ps);
template bool createPolySetFromPolyhedron(const CGAL::Polyhedron_3<CGAL::Epick> &p, PolySet &ps);
template bool createPolySetFromPolyhedron(const CGAL::Polyhedron_3<CGAL::Epeck> &p, PolySet &ps);
template bool createPolySetFromPolyhedron(const CGAL::Polyhedron_3<CGAL::Simple_cartesian<long> > &p, PolySet &ps);
/////// Tessellation begin

View File

@ -12,6 +12,7 @@ namespace CGALUtils {
Polygon2d *project(const CGAL_Nef_polyhedron &N, bool cut);
CGAL_Iso_cuboid_3 boundingBox(const CGAL_Nef_polyhedron3 &N);
bool is_approximately_convex(const PolySet &ps);
Geometry const* applyMinkowski(const Geometry::ChildList &children);
};
CGAL_Nef_polyhedron *createNefPolyhedronFromGeometry(const class Geometry &geom);