vitalif
/
openscad
mirror of https://github.com/vitalif/openscad
1
0
Fork 0
Code Issues Releases Wiki Activity
openscad/cgal/decompose.cpp

668 lines
21 KiB
C++
Raw Permalink Blame History

#include <boost/foreach.hpp>
#include <boost/regex.hpp>
#include <sstream>
#include <iostream>
#include <locale.h>
#include "cgalutils.h"
#include "export.h"
#include "polyset.h"
#include "CGAL_Nef_polyhedron.h"
using namespace CGALUtils;
// Nef polyhedron are using CGAL_Kernel3 (Cartesian<Gmpq>)
// Triangulation will use Epick
typedef CGAL::Epick K;
typedef CGAL::Polyhedron_3<K> PolyhedronK;
#include <boost/assign/std/vector.hpp>
#include <boost/assign/list_of.hpp>
using namespace boost::assign; // bring 'operator+=()' into scope
std::vector<Color4f> colors = boost::assign::list_of
(Color4f(240, 248, 255))
(Color4f(250, 235, 215))
(Color4f(0, 255, 255))
(Color4f(127, 255, 212))
(Color4f(240, 255, 255))
(Color4f(245, 245, 220))
(Color4f(255, 228, 196))
(Color4f(0, 0, 0))
(Color4f(255, 235, 205))
(Color4f(0, 0, 255))
(Color4f(138, 43, 226))
(Color4f(165, 42, 42))
(Color4f(222, 184, 135))
(Color4f(95, 158, 160))
(Color4f(127, 255, 0))
(Color4f(210, 105, 30))
(Color4f(255, 127, 80))
(Color4f(100, 149, 237))
(Color4f(255, 248, 220))
(Color4f(220, 20, 60))
(Color4f(0, 255, 255))
(Color4f(0, 0, 139))
(Color4f(0, 139, 139))
(Color4f(184, 134, 11))
(Color4f(169, 169, 169))
(Color4f(0, 100, 0))
(Color4f(169, 169, 169))
(Color4f(189, 183, 107))
(Color4f(139, 0, 139))
(Color4f(85, 107, 47))
(Color4f(255, 140, 0))
(Color4f(153, 50, 204))
(Color4f(139, 0, 0))
(Color4f(233, 150, 122))
(Color4f(143, 188, 143))
(Color4f(72, 61, 139))
(Color4f(47, 79, 79))
(Color4f(47, 79, 79))
(Color4f(0, 206, 209))
(Color4f(148, 0, 211))
(Color4f(255, 20, 147))
(Color4f(0, 191, 255))
(Color4f(105, 105, 105))
(Color4f(105, 105, 105))
(Color4f(30, 144, 255))
(Color4f(178, 34, 34))
(Color4f(255, 250, 240))
(Color4f(34, 139, 34))
(Color4f(255, 0, 255))
(Color4f(220, 220, 220))
(Color4f(248, 248, 255))
(Color4f(255, 215, 0))
(Color4f(218, 165, 32))
(Color4f(128, 128, 128))
(Color4f(0, 128, 0))
(Color4f(173, 255, 47))
(Color4f(128, 128, 128))
(Color4f(240, 255, 240))
(Color4f(255, 105, 180))
(Color4f(205, 92, 92))
(Color4f(75, 0, 130))
(Color4f(255, 255, 240))
(Color4f(240, 230, 140))
(Color4f(230, 230, 250))
(Color4f(255, 240, 245))
(Color4f(124, 252, 0))
(Color4f(255, 250, 205))
(Color4f(173, 216, 230))
(Color4f(240, 128, 128))
(Color4f(224, 255, 255))
(Color4f(250, 250, 210))
(Color4f(211, 211, 211))
(Color4f(144, 238, 144))
(Color4f(211, 211, 211))
(Color4f(255, 182, 193))
(Color4f(255, 160, 122))
(Color4f(32, 178, 170))
(Color4f(135, 206, 250))
(Color4f(119, 136, 153))
(Color4f(119, 136, 153))
(Color4f(176, 196, 222))
(Color4f(255, 255, 224))
(Color4f(0, 255, 0))
(Color4f(50, 205, 50))
(Color4f(250, 240, 230))
(Color4f(255, 0, 255))
(Color4f(128, 0, 0))
(Color4f(102, 205, 170))
(Color4f(0, 0, 205))
(Color4f(186, 85, 211))
(Color4f(147, 112, 219))
(Color4f(60, 179, 113))
(Color4f(123, 104, 238))
(Color4f(0, 250, 154))
(Color4f(72, 209, 204))
(Color4f(199, 21, 133))
(Color4f(25, 25, 112))
(Color4f(245, 255, 250))
(Color4f(255, 228, 225))
(Color4f(255, 228, 181))
(Color4f(255, 222, 173))
(Color4f(0, 0, 128))
(Color4f(253, 245, 230))
(Color4f(128, 128, 0))
(Color4f(107, 142, 35))
(Color4f(255, 165, 0))
(Color4f(255, 69, 0))
(Color4f(218, 112, 214))
(Color4f(238, 232, 170))
(Color4f(152, 251, 152))
(Color4f(175, 238, 238))
(Color4f(219, 112, 147))
(Color4f(255, 239, 213))
(Color4f(255, 218, 185))
(Color4f(205, 133, 63))
(Color4f(255, 192, 203))
(Color4f(221, 160, 221))
(Color4f(176, 224, 230))
(Color4f(128, 0, 128))
(Color4f(255, 0, 0))
(Color4f(188, 143, 143))
(Color4f(65, 105, 225))
(Color4f(139, 69, 19))
(Color4f(250, 128, 114))
(Color4f(244, 164, 96))
(Color4f(46, 139, 87))
(Color4f(255, 245, 238))
(Color4f(160, 82, 45))
(Color4f(192, 192, 192))
(Color4f(135, 206, 235))
(Color4f(106, 90, 205))
(Color4f(112, 128, 144))
(Color4f(112, 128, 144))
(Color4f(255, 250, 250))
(Color4f(0, 255, 127))
(Color4f(70, 130, 180))
(Color4f(210, 180, 140))
(Color4f(0, 128, 128))
(Color4f(216, 191, 216))
(Color4f(255, 99, 71))
(Color4f(0, 0, 0, 0))
(Color4f(64, 224, 208))
(Color4f(238, 130, 238))
(Color4f(245, 222, 179))
(Color4f(255, 255, 255))
(Color4f(245, 245, 245))
(Color4f(255, 255, 0))
(Color4f(154, 205, 50));
#include <boost/unordered_set.hpp>
#include <CGAL/convex_hull_3.h>
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;
}
}
// Also make sure that there is only one shell:
boost::unordered_set<typename Polyhedron::Facet_const_handle, typename CGAL::Handle_hash_function> visited;
// c++11
// visited.reserve(p.size_of_facets());
std::queue<typename Polyhedron::Facet_const_handle> to_explore;
to_explore.push(p.facets_begin()); // One arbitrary facet
visited.insert(to_explore.front());
while (!to_explore.empty()) {
typename Polyhedron::Facet_const_handle f = to_explore.front();
to_explore.pop();
typename Polyhedron::Facet::Halfedge_around_facet_const_circulator he, end;
end = he = f->facet_begin();
CGAL_For_all(he,end) {
typename Polyhedron::Facet_const_handle o = he->opposite()->facet();
if (!visited.count(o)) {
visited.insert(o);
to_explore.push(o);
}
}
}
return visited.size() == p.size_of_facets();
}
class Shell_explorer
{
public:
std::vector<K::Point_3> vertices;
Shell_explorer() {}
void visit(CGAL_Nef_polyhedron3::Vertex_const_handle v) {
vertices.push_back(K::Point_3(to_double(v->point()[0]),
to_double(v->point()[1]),
to_double(v->point()[2])));
}
void visit(CGAL_Nef_polyhedron3::Halfedge_const_handle ) {}
void visit(CGAL_Nef_polyhedron3::Halffacet_const_handle ) {}
void visit(CGAL_Nef_polyhedron3::SHalfedge_const_handle ) {}
void visit(CGAL_Nef_polyhedron3::SHalfloop_const_handle ) {}
void visit(CGAL_Nef_polyhedron3::SFace_const_handle ) {}
};
template<class Output>
void decompose(const CGAL_Nef_polyhedron3 *N, Output out_iter)
{
int parts = 0;
assert(N);
CGAL_Polyhedron poly;
if (N->is_simple()) {
nefworkaround::convert_to_Polyhedron<CGAL_Kernel3>(*N, poly);
}
if (is_weakly_convex(poly)) {
PRINTD("Minkowski: Object is convex and Nef");
PolyhedronK poly2;
CGALUtils::copyPolyhedron(poly, poly2);
*out_iter++ = poly2;
return;
}
else {
PRINTD("Minkowski: Object is nonconvex Nef, decomposing...");
CGAL_Nef_polyhedron3 decomposed_nef = *N;
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();
// Convert each convex volume to a Polyhedron
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.push_back(poly);
auto s = CGAL_Nef_polyhedron3::SFace_const_handle(ci->shells_begin());
CGAL_Nef_polyhedron3::SFace_const_iterator sf = ci->shells_begin();
Shell_explorer SE;
decomposed_nef.visit_shell_objects(CGAL_Nef_polyhedron3::SFace_const_handle(sf),SE);
PolyhedronK poly;
CGAL::convex_hull_3(SE.vertices.begin(), SE.vertices.end(), poly);
*out_iter++ = poly;
parts++;
}
}
PRINTDB("Minkowski: decomposed into %d convex parts", parts);
}
}
Geometry const * minkowskitest(const Geometry::ChildList &children)
{
CGAL::Timer t,t_tot;
assert(children.size() >= 2);
// Iterate over children, perform pairwise minkowski on children:
// operands = [ch, ch+1]
Geometry::ChildList::const_iterator minkowski_ch_it = children.begin();
t_tot.start();
Geometry const *operands[2] = {minkowski_ch_it->second.get(), NULL};
try {
while (++minkowski_ch_it != children.end()) {
operands[1] = minkowski_ch_it->second.get();
std::vector<PolyhedronK> convexP[2]; // Stores decomposed operands
std::list<PolyhedronK> result_parts;
for (int i = 0; i < 2; i++) {
shared_ptr<const CGAL_Nef_polyhedron> N;
if (const PolySet *ps = dynamic_cast<const PolySet *>(operands[i])) {
if (ps->is_convex()) {
PRINTDB("Minkowski: child %d is convex and PolySet", i);
PolyhedronK poly;
CGALUtils::createPolyhedronFromPolySet(*ps, poly);
convexP[i].push_back(poly);
}
else {
PRINTDB("Minkowski: child %d is nonconvex PolySet, transforming to Nef", i);
N.reset(createNefPolyhedronFromGeometry(*ps));
}
}
else if (const CGAL_Nef_polyhedron *n = dynamic_cast<const CGAL_Nef_polyhedron *>(operands[i])) {
CGAL_Polyhedron poly;
if (n->p3->is_simple()) {
nefworkaround::convert_to_Polyhedron<CGAL_Kernel3>(*n->p3, poly);
// FIXME: Can we calculate weakly_convex on a PolyhedronK instead?
if (is_weakly_convex(poly)) {
PRINTDB("Minkowski: child %d is convex and Nef", i);
PolyhedronK poly2;
CGALUtils::copyPolyhedron(poly, poly2);
convexP[i].push_back(poly2);
}
else {
PRINTDB("Minkowski: child %d is nonconvex Nef",i);
N.reset(n);
}
}
else throw 0; // We cannot handle this, fall back to CGAL's minkowski
}
// If not convex...
if (N && N->p3) {
PRINTD("Decomposing...");
decompose(N->p3.get(), std::back_inserter(convexP[i]));
}
PRINTD("Hulling convex parts...");
std::vector<K::Point_3> points[2];
std::vector<K::Point_3> minkowski_points;
// For each permutation of convex operands..
BOOST_FOREACH(const PolyhedronK &p0, convexP[0]) {
BOOST_FOREACH(const PolyhedronK &p1, convexP[1]) {
t.start();
// Create minkowski pointcloud
minkowski_points.clear();
minkowski_points.reserve(p0.size_of_vertices() * p0.size_of_vertices());
BOOST_FOREACH(const K::Point_3 &p0p, std::make_pair(p0.points_begin(), p0.points_end())) {
BOOST_FOREACH(const K::Point_3 &p1p, std::make_pair(p1.points_begin(), p1.points_end())) {
minkowski_points.push_back(p0p+(p1p-CGAL::ORIGIN));
}
}
t.stop();
// Ignore empty volumes
if (minkowski_points.size() <= 3) continue;
// Hull point cloud
PolyhedronK result;
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<K::Point_3> strict_points;
strict_points.reserve(minkowski_points.size());
for (PolyhedronK::Vertex_iterator i = result.vertices_begin(); i != result.vertices_end(); ++i) {
K::Point_3 const &p = i->point();
PolyhedronK::Vertex::Halfedge_handle h,e;
h = i->halfedge();
e = h;
bool collinear = false;
bool coplanar = true;
do {
K::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 (PolyhedronK::Vertex::Halfedge_handle j = h->next_on_vertex();
j != h && !collinear && ! coplanar;
j = j->next_on_vertex()) {
K::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 (minkowski_ch_it != boost::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<PolyhedronK>::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 = CGALUtils::applyOperator(fake_children, OPENSCAD_UNION);
t.stop();
if (N) PRINTDB("Minkowski: Union done: %f s",t.time());
else PRINTDB("Minkowski: Union failed: %f s",t.time());
t.reset();
operands[0] = N;
} else {
operands[0] = new CGAL_Nef_polyhedron();
}
}
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 = applyOperator(children, OPENSCAD_MINKOWSKI);
return N;
}
}
#define STL_FACET_NUMBYTES 4*3*4+2
// as there is no 'float32_t' standard, we assume the systems 'float'
// is a 'binary32' aka 'single' standard IEEE 32-bit floating point type
union stl_facet {
uint8_t data8[ STL_FACET_NUMBYTES ];
uint32_t data32[4*3];
struct facet_data {
float i, j, k;
float x1, y1, z1;
float x2, y2, z2;
float x3, y3, z3;
uint16_t attribute_byte_count;
} data;
};
void uint32_byte_swap( uint32_t &x )
{
#if __GNUC__ >= 4 && __GNUC_MINOR__ >= 3
x = __builtin_bswap32( x );
#elif defined(__clang__)
x = __builtin_bswap32( x );
#elif defined(_MSC_VER)
x = _byteswap_ulong( x );
#else
uint32_t b1 = ( 0x000000FF & x ) << 24;
uint32_t b2 = ( 0x0000FF00 & x ) << 8;
uint32_t b3 = ( 0x00FF0000 & x ) >> 8;
uint32_t b4 = ( 0xFF000000 & x ) >> 24;
x = b1 | b2 | b3 | b4;
#endif
}
void read_stl_facet( std::ifstream &f, stl_facet &facet )
{
f.read( (char*)facet.data8, STL_FACET_NUMBYTES );
#ifdef BOOST_BIG_ENDIAN
for ( int i = 0; i < 12; i++ ) {
uint32_byte_swap( facet.data32[ i ] );
}
// we ignore attribute byte count
#endif
}
PolySet *import_stl(const std::string &filename)
{
PolySet *p = new PolySet(3);
// Open file and position at the end
std::ifstream f(filename.c_str(), std::ios::in | std::ios::binary | std::ios::ate);
if (!f.good()) {
PRINTB("WARNING: Can't open import file '%s'.", filename);
return NULL;
}
boost::regex ex_sfe("solid|facet|endloop");
boost::regex ex_outer("outer loop");
boost::regex ex_vertex("vertex");
boost::regex ex_vertices("\\s*vertex\\s+([^\\s]+)\\s+([^\\s]+)\\s+([^\\s]+)");
bool binary = false;
std::streampos file_size = f.tellg();
f.seekg(80);
if (f.good() && !f.eof()) {
uint32_t facenum = 0;
f.read((char *)&facenum, sizeof(uint32_t));
#ifdef BOOST_BIG_ENDIAN
uint32_byte_swap( facenum );
#endif
if (file_size == static_cast<std::streamoff>(80 + 4 + 50*facenum)) {
binary = true;
}
}
f.seekg(0);
char data[5];
f.read(data, 5);
if (!binary && !f.eof() && f.good() && !memcmp(data, "solid", 5)) {
int i = 0;
double vdata[3][3];
std::string line;
std::getline(f, line);
while (!f.eof()) {
std::getline(f, line);
boost::trim(line);
if (boost::regex_search(line, ex_sfe)) {
continue;
}
if (boost::regex_search(line, ex_outer)) {
i = 0;
continue;
}
boost::smatch results;
if (boost::regex_search(line, results, ex_vertices)) {
try {
for (int v=0;v<3;v++) {
vdata[i][v] = boost::lexical_cast<double>(results[v+1]);
}
}
catch (const boost::bad_lexical_cast &blc) {
PRINTB("WARNING: Can't parse vertex line '%s'.", line);
i = 10;
continue;
}
if (++i == 3) {
p->append_poly();
p->append_vertex(vdata[0][0], vdata[0][1], vdata[0][2]);
p->append_vertex(vdata[1][0], vdata[1][1], vdata[1][2]);
p->append_vertex(vdata[2][0], vdata[2][1], vdata[2][2]);
}
}
}
}
else if (binary && !f.eof() && f.good())
{
f.ignore(80-5+4);
while (1) {
stl_facet facet;
read_stl_facet( f, facet );
if (f.eof()) break;
p->append_poly();
p->append_vertex(facet.data.x1, facet.data.y1, facet.data.z1);
p->append_vertex(facet.data.x2, facet.data.y2, facet.data.z2);
p->append_vertex(facet.data.x3, facet.data.y3, facet.data.z3);
}
}
return p;
}
/*!
file format:
1. polygon coordinates (x,y,z) are comma separated (+/- spaces) and
each coordinate is on a separate line
2. each polygon is separated by one or more blank lines
*/
bool import_polygon(PolyholeK &polyhole, const std::string &filename)
{
std::ifstream ifs(filename.c_str());
if (!ifs) return false;
std::string line;
PolygonK polygon;
while (std::getline(ifs, line)) {
std::stringstream ss(line);
double X = 0.0, Y = 0.0, Z = 0.0;
if (!(ss >> X)) {
//ie blank lines => flag start of next polygon
if (polygon.size() > 0) polyhole.push_back(polygon);
polygon.clear();
continue;
}
char c = ss.peek();
while (c == ' ') {ss.read(&c, 1); c = ss.peek();} //gobble spaces before comma
if (c == ',') {ss.read(&c, 1); c = ss.peek();} //gobble comma
while (c == ' ') {ss.read(&c, 1); c = ss.peek();} //gobble spaces after comma
if (!(ss >> Y)) {
std::cerr << "Y error\n";
return false;
}
c = ss.peek();
while (c == ' ') {ss.read(&c, 1); c = ss.peek();} //gobble spaces before comma
if (c == ',') {ss.read(&c, 1); c = ss.peek();} //gobble comma
while (c == ' ') {ss.read(&c, 1); c = ss.peek();} //gobble spaces after comma
if (!(ss >> Z)) {
std::cerr << "Z error\n";
return false;
}
polygon.push_back(Vertex3K(X, Y, Z));
}
if (polygon.size() > 0) polyhole.push_back(polygon);
ifs.close();
return true;
}
//------------------------------------------------------------------------------
int main(int argc, char *argv[])
{
OpenSCAD::debug = "decompose";
PolySet *ps = NULL;
if (argc == 2) {
if (!(ps = import_stl(argv[1]))) {
std::cerr << "Error importing STL " << argv[1] << std::endl;
exit(1);
}
std::cerr << "Imported " << ps->numPolygons() << " polygons" << std::endl;
}
else {
std::cerr << "Usage: " << argv[0] << " <file.stl> <file.stl>" << std::endl;
exit(1);
}
Geometry::ChildList children;
CGAL_Nef_polyhedron *N = createNefPolyhedronFromGeometry(*ps);
std::vector<PolyhedronK> result;
decompose(N->p3.get(), std::back_inserter(result));
std::cerr << "Decomposed into " << result.size() << " convex parts" << std::endl;
int idx = 0;
BOOST_FOREACH(const PolyhedronK &P, result) {
PolySet result_ps(3);
if (CGALUtils::createPolySetFromPolyhedron(P, result_ps)) {
std::cerr << "Error converting to PolySet\n";
}
else {
std::stringstream ss;
ss << "out" << idx++ << ".stl";
exportFileByName(&result_ps, OPENSCAD_STL, ss.str().c_str(), ss.str().c_str());
std::cout << "color([" << colors[idx%147][0] << "," << colors[idx%147][1] << "," << colors[idx%147][2] << "]) " << "import(\"" << ss.str() << "\");\n";
}
}
std::cerr << "Done." << std::endl;
}
Reference in New Issue View Git Blame Copy Permalink