diff --git a/src/examples/IfcAdvancedHouse.cpp b/src/examples/IfcAdvancedHouse.cpp index 466bb2df89..93371cf7ae 100644 --- a/src/examples/IfcAdvancedHouse.cpp +++ b/src/examples/IfcAdvancedHouse.cpp @@ -38,18 +38,13 @@ #include -#ifdef USE_IFC4 -#include "../ifcparse/Ifc4.h" -#define IfcSchema Ifc4 -#else -#include "../ifcparse/Ifc2x3.h" #define IfcSchema Ifc2x3 -#endif - +#include "../ifcparse/macros.h" +#include "../ifcparse/Ifc2x3.h" #include "../ifcparse/IfcBaseClass.h" #include "../ifcparse/IfcHierarchyHelper.h" -#include "../ifcgeom/IfcGeom.h" -#include "../ifcgeom_schema_agnostic/Serialization.h" + +#include "../ifcgeom/schema_agnostic/Serialization.h" #if USE_VLD #include diff --git a/src/examples/IfcOpenHouse.cpp b/src/examples/IfcOpenHouse.cpp index 189dd4320e..6ea1eb3dff 100644 --- a/src/examples/IfcOpenHouse.cpp +++ b/src/examples/IfcOpenHouse.cpp @@ -33,18 +33,15 @@ #include #include -#ifdef USE_IFC4 -#include "../ifcparse/Ifc4.h" -#define IfcSchema Ifc4 -#else -#include "../ifcparse/Ifc2x3.h" -#define IfcSchema Ifc2x3 -#endif +#include +#define IfcSchema Ifc2x3 +#include "../ifcparse/macros.h" +#include "../ifcparse/Ifc2x3.h" #include "../ifcparse/IfcBaseClass.h" #include "../ifcparse/IfcHierarchyHelper.h" -#include "../ifcgeom/IfcGeom.h" -#include "../ifcgeom_schema_agnostic/Serialization.h" + +#include "../ifcgeom/schema_agnostic/Serialization.h" #if USE_VLD #include diff --git a/src/ifcconvert/validate_space_boundaries.cpp b/src/ifcconvert/validate_space_boundaries.cpp new file mode 100644 index 0000000000..0f660e5901 --- /dev/null +++ b/src/ifcconvert/validate_space_boundaries.cpp @@ -0,0 +1,185 @@ +#include "validation_utils.h" + +using namespace ifcopenshell::geometry; + +#include +#include +#include +#include + +typedef Kernel_::FT FT; +typedef Kernel_::Point_3 Point; +typedef Kernel_::Segment_3 Segment; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef CGAL::AABB_face_graph_triangle_primitive Primitive; +typedef CGAL::AABB_traits Traits; +typedef CGAL::AABB_tree Tree; +typedef Tree::Point_and_primitive_id Point_and_primitive_id; + +void fix_spaceboundaries(IfcParse::IfcFile& f, bool no_progress, bool quiet, bool stderr_progress) { + intersection_validator v(f, { "IfcWall", "IfcSpace", "IfcSlab", "IfcCovering" }, 1.e-5, no_progress, quiet, stderr_progress); + + auto rels = f.instances_by_type("IfcRelSpaceBoundary"); + + std::map, const IfcUtil::IfcBaseClass*> rel_by_space_elem; + + + if (rels) { + std::for_each(rels->begin(), rels->end(), [&rel_by_space_elem](const IfcUtil::IfcBaseClass* rel) { + auto x = ((IfcUtil::IfcBaseEntity*)rel)->get_value("RelatingSpace"); + try { + auto y = ((IfcUtil::IfcBaseEntity*)rel)->get_value("RelatedBuildingElement"); + rel_by_space_elem.insert({ { x,y }, rel }); + } catch (IfcParse::IfcException&) { + // RelatedBuildingElement can be NULL + } + }); + } + + std::set rels_encounted; + + IfcParse::IfcFile f2("boundaries-triangulated.ifc"); + if (!f2.good()) { + return; + } + + ifcopenshell::geometry::settings settings; + + settings.set(ifcopenshell::geometry::settings::USE_WORLD_COORDS, false); + settings.set(ifcopenshell::geometry::settings::WELD_VERTICES, false); + settings.set(ifcopenshell::geometry::settings::SEW_SHELLS, true); + settings.set(ifcopenshell::geometry::settings::CONVERT_BACK_UNITS, true); + settings.set(ifcopenshell::geometry::settings::DISABLE_TRIANGULATION, true); + settings.set(ifcopenshell::geometry::settings::DISABLE_OPENING_SUBTRACTIONS, true); + + ifcopenshell::geometry::Converter c("cgal", &f2, settings); + + std::map, std::vector> elem_to_space_boundary_coords; + + for (auto& i : *f2.instances_by_type("IfcProduct")) { + auto n = ((IfcUtil::IfcBaseEntity*)i)->get_value("Name"); + auto g1 = n.substr(0, 22); + auto g2 = n.substr(23); + auto item = c.mapping()->map(i); + if (((ifcopenshell::geometry::taxonomy::collection*) item)->children[0] == nullptr) { + continue; + } + auto shell = (taxonomy::shell*) ((taxonomy::collection*)((taxonomy::collection*) item)->children[0])->children[0]; + for (auto& f : shell->children) { + auto face = (taxonomy::face*) f; + for (auto& w : face->children) { + auto wire = (taxonomy::loop*) w; + for (auto& e : wire->children) { + auto edge = (taxonomy::edge*) e; + auto p3 = boost::get(edge->start); + auto p4 = *((taxonomy::geom_item*)item)->matrix.components * p3.components->homogeneous(); + Kernel_::Point_3 P(p4(0), p4(1), p4(2)); + elem_to_space_boundary_coords[{g1, g2}].emplace_back(P); + } + } + } + } + + std::set< std::set > guid_pairs_visited; + + v([&rel_by_space_elem, &elem_to_space_boundary_coords, &guid_pairs_visited](const intersection_validator::Box& a, const intersection_validator::Box& b) { + std::ostringstream ss; + // ss << id_map[a.id()]->first->data().toString() << "x" << id_map[b.id()]->first->data().toString() << std::endl; + // auto x = id_map[a.id()]->second * id_map[b.id()]->second; + + auto A = a.handle()->first; + auto B = b.handle()->first; + + auto Aguid = A->get_value("GlobalId"); + auto Bguid = B->get_value("GlobalId"); + + int space_count = 0; + if (A->declaration().name() == "IfcSpace") { + space_count += 1; + } + if (B->declaration().name() == "IfcSpace") { + space_count += 1; + } + if (space_count != 1) { + return; + } + + ss << a.handle()->first->data().toString() << "x" << a.handle()->first->data().toString() << std::endl; + auto x = a.handle()->second * b.handle()->second; + + if (x.is_empty()) { + return; + } + + guid_pairs_visited.insert({ Aguid, Bguid }); + + cgal_shape_t x_poly; + x.convert_to_polyhedron(x_poly); + + { + std::string fn = "computed_boundaries_" + Aguid + "_" + Bguid + ".off"; + std::ofstream computed_boundaries(fn.c_str()); + computed_boundaries.precision(17); + computed_boundaries << x_poly; + } + + Tree tree(faces(x_poly).first, faces(x_poly).second, x_poly); + tree.accelerate_distance_queries(); + + auto itelem = elem_to_space_boundary_coords.find({ Aguid, Bguid }); + + if (itelem == elem_to_space_boundary_coords.end()) { + Logger::Error("Missing space boundary relationship " + Aguid + " " + Bguid); + return; + } + + const auto& coords = itelem->second; + std::vector distances; + std::transform(coords.begin(), coords.end(), std::back_inserter(distances), [&tree](const Kernel_::Point_3& p) { + return std::sqrt(CGAL::to_double(tree.squared_distance(p))); + }); + + bool valid = *std::max_element(distances.begin(), distances.end()) < 0.4; + + if (!valid) { + Logger::Error("Wrong connection geometry " + Aguid + " " + Bguid); + } + + /*{ + remove_thickness r(x_poly); + std::string fn = "thin_computed_boundaries_" + Aguid + "_" + Bguid + ".off"; + std::ofstream computed_boundaries(fn.c_str()); + computed_boundaries.precision(17); + computed_boundaries << r.flattened; + }*/ + + /* + { + auto FN = s0 + "-" + s1 + "-" + std::to_string(i0) + "-" + std::to_string(i1) + "-sides-sb.off"; + std::ofstream os(FN.c_str()); + os.precision(17); + os << r.polyhedron2; + } + { + auto FN = s0 + "-" + s1 + "-" + std::to_string(i0) + "-" + std::to_string(i1) + "-flat-sb.off"; + std::ofstream os(FN.c_str()); + os.precision(17); + os << r.flattened; + } + */ + }); + + auto is_wall_space_or_slab = [&f](const std::string& g) { + auto decl = f.instance_by_guid(g)->declaration(); + return decl.is("IfcWall") || decl.is("IfcSpace") || decl.is("IfcSlab"); + }; + + for (auto& i : *f2.instances_by_type("IfcProduct")) { + auto n = ((IfcUtil::IfcBaseEntity*)i)->get_value("Name"); + auto g1 = n.substr(0, 22); + auto g2 = n.substr(23); + if (is_wall_space_or_slab(g1) && is_wall_space_or_slab(g2) && guid_pairs_visited.find({ g1, g2 }) == guid_pairs_visited.end()) { + Logger::Error("Space boundary for non-bounding geometry " + g1 + " " + g2); + } + } +} diff --git a/src/ifcconvert/validate_storey_containment.cpp b/src/ifcconvert/validate_storey_containment.cpp new file mode 100644 index 0000000000..c003291230 --- /dev/null +++ b/src/ifcconvert/validate_storey_containment.cpp @@ -0,0 +1,233 @@ +#include "../ifcgeom/kernels/cgal/CgalKernel.h" +#include "../ifcgeom/schema_agnostic/IfcGeomFilter.h" +#include "../ifcgeom/schema_agnostic/IfcGeomIterator.h" + +#include +#include + +#include + +void fix_storeycontainment(IfcParse::IfcFile& f, bool no_progress, bool quiet, bool stderr_progress) { + ifcopenshell::geometry::settings settings; + settings.set(ifcopenshell::geometry::settings::USE_WORLD_COORDS, false); + settings.set(ifcopenshell::geometry::settings::WELD_VERTICES, false); + settings.set(ifcopenshell::geometry::settings::SEW_SHELLS, true); + settings.set(ifcopenshell::geometry::settings::CONVERT_BACK_UNITS, true); + settings.set(ifcopenshell::geometry::settings::DISABLE_TRIANGULATION, true); + settings.set(ifcopenshell::geometry::settings::DISABLE_OPENING_SUBTRACTIONS, true); + + std::vector no_openings_and_spaces = { + IfcGeom::entity_filter(false, false, {"IfcOpeningElement", "IfcSpace"}) + }; + + ifcopenshell::geometry::Iterator context_iterator("cgal", settings, &f, no_openings_and_spaces); + + auto get_elevation = [](IfcUtil::IfcBaseClass* a) { + return ((IfcUtil::IfcBaseEntity*)a)->get_value_or("Elevation", 0.); + }; + + // latebound inverse attribute lookup not working + auto rels = f.instances_by_type("IfcRelContainedInSpatialStructure"); + std::map elem_to_storey; + std::for_each(rels->begin(), rels->end(), [&elem_to_storey](IfcUtil::IfcBaseClass* r) { + auto elems = ((IfcUtil::IfcBaseEntity*)r)->get_value("RelatedElements"); + auto storey = ((IfcUtil::IfcBaseEntity*)r)->get_value("RelatingStructure"); + + if (storey->declaration().name() == "IfcBuildingStorey") { + for (auto it = elems->begin(); it != elems->end(); ++it) { + elem_to_storey[*it] = storey; + } + } + }); + + auto storeys = f.instances_by_type("IfcBuildingStorey"); + std::vector storeys_sorted(storeys->begin(), storeys->end()); + std::sort(storeys_sorted.begin(), storeys_sorted.end(), [&get_elevation](IfcUtil::IfcBaseClass* a, IfcUtil::IfcBaseClass* b) { + return get_elevation(a) < get_elevation(b); + }); + + /* + std::wcout << "Storeys "; + for (auto& s : storeys_sorted) { + auto n = ((IfcUtil::IfcBaseEntity*)s)->get_value("Name"); + std::wcout << n.c_str() << " "; + } + std::wcout << std::endl; + */ + + std::vector elevations; + std::transform(storeys_sorted.begin(), storeys_sorted.end(), std::back_inserter(elevations), get_elevation); + + double LARGE = 1e4; + + std::vector> elevation_slices; + for (size_t i = 0; i < elevations.size(); ++i) { + elevation_slices.push_back({ + i == 0 ? -LARGE : elevations[i], + i + 1 == elevations.size() ? LARGE : elevations[i + 1] + }); + } + + std::for_each(elevation_slices.begin(), elevation_slices.end(), [](std::pair& p) { + p.first -= 0.3; + p.second += 0.3; + }); + + std::vector> nefs; + std::transform(elevation_slices.begin(), elevation_slices.end(), std::back_inserter(nefs), [&LARGE](const std::pair& p) { + // std::wcout << p.first << " - " << p.second << std::endl; + Kernel_::Point_3 p1(-LARGE, -LARGE, p.first); + Kernel_::Point_3 p2(+LARGE, +LARGE, p.second); + auto poly = ifcopenshell::geometry::utils::create_cube(p1, p2); + return ifcopenshell::geometry::utils::create_nef_polyhedron(poly); + }); + + /* + for (auto& n : nefs) { + auto poly = ifcopenshell::geometry::utils::create_polyhedron(n); + auto bounds = CGAL::Polygon_mesh_processing::bbox_3(poly); + for (int i = 0; i < 3; ++i) { + std::wcout << bounds.min(i) << std::endl; + } + for (int i = 0; i < 3; ++i) { + std::wcout << bounds.max(i) << std::endl; + } + std::wcout << "---" << std::endl; + } + */ + + if (!context_iterator.initialize()) { + return; + } + + size_t num_created = 0; + int old_progress = quiet ? 0 : -1; + + for (;; ++num_created) { + bool has_more = true; + if (num_created) { + has_more = context_iterator.next(); + } + ifcopenshell::geometry::NativeElement* geom_object = nullptr; + if (has_more) { + geom_object = context_iterator.get_native(); + } + if (!geom_object) { + break; + } + + /* + std::stringstream ss; + ss << geom_object->product()->data().toString(); + auto sss = ss.str(); + std::wcout << sss.c_str() << std::endl; + */ + + if (elem_to_storey.find(geom_object->product()) == elem_to_storey.end()) { + // std::wcout << "not associated to storey" << std::endl; + continue; + } + + std::vector intersection_volumes(nefs.size()); + + for (auto& g : geom_object->geometry()) { + auto s = ((ifcopenshell::geometry::CgalShape*) g.Shape())->shape(); + const auto& m = *g.Placement().components; + const auto& n = *geom_object->transformation().data().components; + + const cgal_placement_t trsf( + m(0, 0), m(0, 1), m(0, 2), m(0, 3), + m(1, 0), m(1, 1), m(1, 2), m(1, 3), + m(2, 0), m(2, 1), m(2, 2), m(2, 3)); + + const cgal_placement_t trsf2( + n(0, 0), n(0, 1), n(0, 2), n(0, 3), + n(1, 0), n(1, 1), n(1, 2), n(1, 3), + n(2, 0), n(2, 1), n(2, 2), n(2, 3)); + + // Apply transformation + for (auto &vertex : vertices(s)) { + vertex->point() = vertex->point().transform(trsf).transform(trsf2); + } + + /* + { + auto bounds = CGAL::Polygon_mesh_processing::bbox_3(s); + for (int i = 0; i < 3; ++i) { + std::wcout << bounds.min(i) << std::endl; + } + for (int i = 0; i < 3; ++i) { + std::wcout << bounds.max(i) << std::endl; + } + std::wcout << "---" << std::endl; + } + */ + + CGAL::Nef_polyhedron_3 part_nef = ifcopenshell::geometry::utils::create_nef_polyhedron(s); + + if (!part_nef.is_simple()) { + // std::wcout << "not simple" << std::endl; + continue; + } + + std::vector::iterator accumulator = intersection_volumes.begin(); + std::for_each(nefs.begin(), nefs.end(), [&accumulator, &part_nef](const CGAL::Nef_polyhedron_3& storey_nef) { + auto poly = ifcopenshell::geometry::utils::create_polyhedron(part_nef * storey_nef); + CGAL::Polygon_mesh_processing::triangulate_faces(poly); + *accumulator += CGAL::to_double(CGAL::Polygon_mesh_processing::volume(poly)); + accumulator++; + }); + } + + /* + std::wcout << "volumes: "; + for (auto& v : intersection_volumes) { + std::wcout << v << " "; + } + std::wcout << std::endl; + */ + + auto calc_idx = std::max_element(intersection_volumes.begin(), intersection_volumes.end()) - intersection_volumes.begin(); + auto calc_overlap = intersection_volumes[calc_idx]; + auto assigned_idx = std::distance(storeys_sorted.begin(), std::find(storeys_sorted.begin(), storeys_sorted.end(), elem_to_storey[geom_object->product()])); + auto assigned_overlap = intersection_volumes[assigned_idx]; + if (calc_overlap > 0 && assigned_overlap < calc_overlap * 0.9) { + auto s = geom_object->product()->get_value("GlobalId"); + auto s1 = ((IfcUtil::IfcBaseEntity*)storeys_sorted[calc_idx])->get_value("GlobalId"); + auto s2 = ((IfcUtil::IfcBaseEntity*)elem_to_storey[geom_object->product()])->get_value("GlobalId"); + Logger::Error("Element " + s + " contained in " + s2 + " located on " + s1); + } + + if (!no_progress) { + if (quiet) { + const int progress = context_iterator.progress(); + for (; old_progress < progress; ++old_progress) { + std::cout << "."; + if (stderr_progress) + std::cerr << "."; + } + std::cout << std::flush; + if (stderr_progress) + std::cerr << std::flush; + } else { + const int progress = context_iterator.progress() / 2; + if (old_progress != progress) Logger::ProgressBar(progress); + old_progress = progress; + } + } + } + + if (!no_progress && quiet) { + for (; old_progress < 100; ++old_progress) { + std::cout << "."; + if (stderr_progress) + std::cerr << "."; + } + std::cout << std::flush; + if (stderr_progress) + std::cerr << std::flush; + } else { + Logger::Status("\rDone fixing space boundaries for " + boost::lexical_cast(num_created) + + " objects "); + } +} diff --git a/src/ifcconvert/validate_wall_connectivity.cpp b/src/ifcconvert/validate_wall_connectivity.cpp new file mode 100644 index 0000000000..28de1f22e3 --- /dev/null +++ b/src/ifcconvert/validate_wall_connectivity.cpp @@ -0,0 +1,198 @@ +#include "validation_utils.h" + +#include +#include + +#include + +using namespace ifcopenshell::geometry; + +void fix_wallconnectivity(IfcParse::IfcFile& f, bool no_progress, bool quiet, bool stderr_progress) { + intersection_validator v(f, { "IfcWall" }, 1.e-3, no_progress, quiet, stderr_progress); + + + ifcopenshell::geometry::settings settings; + settings.set(ifcopenshell::geometry::settings::USE_WORLD_COORDS, false); + settings.set(ifcopenshell::geometry::settings::WELD_VERTICES, false); + settings.set(ifcopenshell::geometry::settings::SEW_SHELLS, true); + settings.set(ifcopenshell::geometry::settings::CONVERT_BACK_UNITS, true); + settings.set(ifcopenshell::geometry::settings::DISABLE_TRIANGULATION, true); + settings.set(ifcopenshell::geometry::settings::DISABLE_OPENING_SUBTRACTIONS, true); + + settings.set(ifcopenshell::geometry::settings::INCLUDE_CURVES, true); + settings.set(ifcopenshell::geometry::settings::EXCLUDE_SOLIDS_AND_SURFACES, true); + + ifcopenshell::geometry::Converter c("cgal", &f, settings); + + auto rels = f.instances_by_type("IfcRelConnectsPathElements"); + std::map, const IfcUtil::IfcBaseClass*> rel_by_elem; + std::for_each(rels->begin(), rels->end(), [&rel_by_elem](const IfcUtil::IfcBaseClass* rel) { + auto x = ((IfcUtil::IfcBaseEntity*)rel)->get_value("RelatingElement"); + auto y = ((IfcUtil::IfcBaseEntity*)rel)->get_value("RelatedElement"); + rel_by_elem.insert({{ x,y }, rel}); + }); + + std::set rels_encounted; + + double total_nef_intersection_time = 0.; + double conversion_to_poly = 0.; + + v([&c, &rel_by_elem, &rels_encounted, &total_nef_intersection_time, &conversion_to_poly](const intersection_validator::Box& a, const intersection_validator::Box& b) { + auto A = a.handle()->first; + auto B = b.handle()->first; + + const IfcUtil::IfcBaseClass* rel = nullptr; + std::string a_type, b_type; + + auto rit = rel_by_elem.find({ A, B }); + if (rit != rel_by_elem.end()) { + rel = rit->second; + const bool a_is_relating = A == ((IfcUtil::IfcBaseEntity*)rel)->get_value("RelatingElement"); + a_type = ((IfcUtil::IfcBaseEntity*)rel)->get_value("RelatingConnectionType"); + b_type = ((IfcUtil::IfcBaseEntity*)rel)->get_value("RelatedConnectionType"); + if (!a_is_relating) { + std::swap(a_type, b_type); + } + } + +#if 0 + auto a_poly = ifcopenshell::geometry::utils::create_polyhedron(a.handle()->second); + auto b_poly = ifcopenshell::geometry::utils::create_polyhedron(b.handle()->second); + + std::wcout << "a" << std::endl; + for (auto& v : vertices(a_poly)) { + for (int i = 0; i < 3; ++i) { + std::wcout << CGAL::to_double(v->point().cartesian(i)) << " "; + } + std::wcout << std::endl; + } + + std::wcout << "b" << std::endl; + for (auto& v : vertices(b_poly)) { + for (int i = 0; i < 3; ++i) { + std::wcout << CGAL::to_double(v->point().cartesian(i)) << " "; + } + std::wcout << std::endl; + } +#endif + + std::ostringstream ss; + ss << A->data().toString() << "x" << B->data().toString() << std::endl; + std::clock_t intersection_begin = std::clock(); + auto x = a.handle()->second * b.handle()->second; + std::clock_t intersection_end = std::clock(); + + total_nef_intersection_time += (intersection_end - intersection_begin) / (double) CLOCKS_PER_SEC; + + if (x.is_empty()) { + return; + } + + std::clock_t poly_begin = std::clock(); + cgal_shape_t x_poly; + x.convert_to_polyhedron(x_poly); + std::clock_t poly_end = std::clock(); + conversion_to_poly += (poly_end - poly_begin) / (double)CLOCKS_PER_SEC; + + auto dza = a.bbox().zmax() - a.bbox().zmin(); + auto dzb = b.bbox().zmax() - b.bbox().zmin(); + auto bb = CGAL::Polygon_mesh_processing::bbox_3(x_poly); + if (bb.zmax() - bb.zmin() < std::min(dza, dzb) / 3.) { + return; + } + + CGAL::Polygon_mesh_processing::triangulate_faces(x_poly); + if (CGAL::Polygon_mesh_processing::area(x_poly) > 4.0) { + return; + } + + auto get_axis_parameter_min_max = [&c, &x_poly](IfcUtil::IfcBaseEntity* inst) { + auto item = c.mapping()->map(inst); + auto shaperep = ((taxonomy::collection*) item)->children[0]; + auto loop = ((taxonomy::collection*) shaperep)->children[0]; + + if (loop->kind() != taxonomy::LOOP) { + // std::wcout << "no suitable axis" << std::endl; + } else { + auto first_vertex = ((taxonomy::edge*) ((taxonomy::loop*) loop)->children.front())->start; + auto last_vertex = ((taxonomy::edge*) ((taxonomy::loop*) loop)->children.back())->end; + + if (first_vertex.which() != 0 || last_vertex.which() != 0) { + // std::wcout << "trims not supported" << std::endl; + } else { + auto p0 = boost::get(first_vertex); + auto p1 = boost::get(last_vertex); + + auto v0 = *((taxonomy::geom_item*)item)->matrix.components * p0.components->homogeneous(); + auto v1 = *((taxonomy::geom_item*)item)->matrix.components * p1.components->homogeneous(); + + auto P0 = Kernel_::Point_3(v0(0), v0(1), v0(2)); + auto P1 = Kernel_::Point_3(v1(0), v1(1), v1(2)); + + auto D = P1 - P0; + auto len = std::sqrt(CGAL::to_double(D.squared_length())); + D /= len; + + std::vector parameters; + + std::transform(vertices(x_poly).begin(), vertices(x_poly).end(), std::back_inserter(parameters), [&P0, D](cgal_vertex_descriptor_t& v) { + return (v->point() - P0) * D; + }); + + auto pit = std::minmax_element(parameters.begin(), parameters.end()); + return std::make_pair(len, std::make_pair(CGAL::to_double(*pit.first), CGAL::to_double(*pit.second))); + } + } + const auto& nan = std::numeric_limits::quiet_NaN(); + return std::make_pair(nan, std::make_pair(nan, nan)); + }; + + auto qualify_connection_type = [](double l, const std::pair& p) { + if (p.first < 1.e-3) { + return "ATSTART"; + } else if (p.second > l - 1.e-3) { + return "ATEND"; + } else { + return "ATPATH"; + } + }; + + auto alu0u1 = get_axis_parameter_min_max(A); + auto blu0u1 = get_axis_parameter_min_max(B); + + auto atype_computed = qualify_connection_type(alu0u1.first, alu0u1.second); + auto btype_computed = qualify_connection_type(blu0u1.first, blu0u1.second); + + rels_encounted.insert(rel); + + if (a_type != atype_computed || b_type != btype_computed) { + if (rel) { + Logger::Error(std::string("Connection type ") + atype_computed + " " + btype_computed + " for:", rel); + } else { + auto A_str = A->get_value("GlobalId"); + auto B_str = B->get_value("GlobalId"); + Logger::Error("No connection for adjacent " + A_str + " " + B_str); + } + } + }); + + std::for_each(rels->begin(), rels->end(), [&rels_encounted, &v](const IfcUtil::IfcBaseClass* rel) { + if (rels_encounted.find(rel) == rels_encounted.end()) { + auto x = (IfcUtil::IfcBaseEntity*)((IfcUtil::IfcBaseEntity*)rel)->get_value("RelatingElement"); + auto y = (IfcUtil::IfcBaseEntity*)((IfcUtil::IfcBaseEntity*)rel)->get_value("RelatedElement"); + if (v.succesfully_processed.find(x) != v.succesfully_processed.end() && v.succesfully_processed.find(y) != v.succesfully_processed.end()) { + Logger::Error("Connection for non-adjacent walls", rel); + } + } + }); + + std::wcout << std::setprecision(14); + std::wcout << "total_map_time " << v.total_map_time << std::endl; + std::wcout << "total_geom_time " << v.total_geom_time << std::endl; + std::wcout << "total_nef_time " << v.total_nef_time << std::endl; + std::wcout << "total_minkowsky_time " << v.total_minkowsky_time << std::endl; + std::wcout << "total_box_time " << v.total_box_time << std::endl; + std::wcout << "total_nef_intersection_time " << total_nef_intersection_time << std::endl; + std::wcout << "total_conversion_to_poly_time " << conversion_to_poly << std::endl; +} + diff --git a/src/ifcconvert/validation_utils.cpp b/src/ifcconvert/validation_utils.cpp new file mode 100644 index 0000000000..512f6163d6 --- /dev/null +++ b/src/ifcconvert/validation_utils.cpp @@ -0,0 +1,28 @@ +#include "validation_utils.h" + +double facet_area(const cgal_shape_t::Facet_handle& f) { + auto p0 = f->facet_begin()->vertex()->point(); + auto p1 = f->facet_begin()->next()->vertex()->point(); + auto p2 = f->facet_begin()->next()->next()->vertex()->point(); + return std::sqrt(CGAL::to_double(CGAL::cross_product(p0 - p1, p2 - p1).squared_length())); +} + +void dump_facet(const cgal_shape_t::Facet_handle& f) { + auto p0 = f->facet_begin()->vertex()->point(); + auto p1 = f->facet_begin()->next()->vertex()->point(); + auto p2 = f->facet_begin()->next()->next()->vertex()->point(); + auto V = CGAL::cross_product(p0 - p1, p2 - p1); + auto d = std::sqrt(CGAL::to_double(V.squared_length())); + if (d > 1.e-20) { + V /= d; + } + + std::ostringstream oss; + oss.precision(8); + oss << "Facet with area " << facet_area(f) << " and normal (" + << CGAL::to_double(V.cartesian(0)) << " " << CGAL::to_double(V.cartesian(1)) << " " + << CGAL::to_double(V.cartesian(2)) << ")"; + + auto osss = oss.str(); + std::wcout << osss.c_str() << std::endl; +} \ No newline at end of file diff --git a/src/ifcconvert/validation_utils.h b/src/ifcconvert/validation_utils.h new file mode 100644 index 0000000000..ecf4a876a4 --- /dev/null +++ b/src/ifcconvert/validation_utils.h @@ -0,0 +1,581 @@ +#include "../ifcgeom/kernels/cgal/CgalKernel.h" +#include "../ifcgeom/schema_agnostic/IfcGeomFilter.h" +#include "../ifcgeom/schema_agnostic/IfcGeomIterator.h" + +#include +#include + +#include +#include +#include +#include + +#include +#include + +template +T enlarge(const T& t, double d = 1.e-5) { + typename T::NT min[3]; + typename T::NT max[3]; + for (int i = 0; i < t.dimension(); ++i) { + min[i] = t.min_coord(i) - d; + max[i] = t.max_coord(i) + d; + } + return T(min, max, t.handle()); +} + +template +struct Build_Offset : public CGAL::Modifier_base { + std::list input; + + void operator()(HDS& hds) { + // Postcondition: hds is a valid polyhedral surface. + CGAL::Polyhedron_incremental_builder_3 B(hds); + + int Nv = 0, Nf = 0; + for (auto& f : input) { + Nv += 3; + Nf += 1; + } + + B.begin_surface(Nv, Nf); + + for (auto& f : input) { + auto p0 = f->facet_begin()->vertex()->point(); + auto p1 = f->facet_begin()->next()->vertex()->point(); + auto p2 = f->facet_begin()->next()->next()->vertex()->point(); + + auto O = CGAL::centroid(p0, p1, p2); + + Kernel_::Point_3* p012[3] = { &p0, &p1, &p2 }; + for (int i = 0; i < 3; ++i) { + *p012[i] = CGAL::ORIGIN + (((*(p012[i])) - CGAL::ORIGIN) + ((*(p012[i])) - O)); + B.add_vertex(*p012[i]); + } + } + + Nv = 0; + for (int i = 0; i < Nf; ++i) { + B.begin_facet(); + B.add_vertex_to_facet(Nv++); + B.add_vertex_to_facet(Nv++); + B.add_vertex_to_facet(Nv++); + B.end_facet(); + } + + B.end_surface(); + } +}; + +template +std::list connected_faces(cgal_shape_t::Facet_handle f, const Ts& excluded) { + std::set fs = { f }; + + std::function process; + process = [&fs, &process, &excluded](cgal_shape_t::Facet_handle& f) { + cgal_shape_t::Halfedge_around_facet_circulator circ = f->facet_begin(), end(circ); + do { + auto ff = circ->opposite()->facet(); + if (excluded.find(ff) == excluded.end()) { + auto p = fs.insert(ff); + if (p.second) { + process(ff); + } + } + } while (++circ != end); + }; + + process(f); + return std::list(fs.begin(), fs.end()); +} + +template +struct Builder_With_Map : public CGAL::Modifier_base { + std::list input; + std::map mapping; + + void operator()(HDS& hds) { + // Postcondition: hds is a valid polyhedral surface. + CGAL::Polyhedron_incremental_builder_3 B(hds); + + std::set used_points; + + for (auto& f : input) { + cgal_shape_t::Halfedge_around_facet_circulator circ = f->facet_begin(), end(circ); + do { + auto P = circ->vertex()->point(); + auto it = mapping.find(P); + if (it == mapping.end()) { + std::wcout << "WARNING unprojected point :(" << std::endl; + } else { + P = it->second; + } + used_points.insert(P); + } while (++circ != end); + } + + B.begin_surface(used_points.size(), input.size()); + + for (auto& p : used_points) { + B.add_vertex(p); + } + + for (auto& f : input) { + B.begin_facet(); + cgal_shape_t::Halfedge_around_facet_circulator circ = f->facet_begin(), end(circ); + do { + auto P = circ->vertex()->point(); + auto it = mapping.find(P); + if (it == mapping.end()) { + std::wcout << "WARNING unprojected point :(" << std::endl; + } else { + P = it->second; + } + + auto jt = used_points.find(P); + if (jt == used_points.end()) { + throw std::runtime_error("Unable to map point"); + } + size_t idx = std::distance(used_points.begin(), jt); + std::wcout << "idx " << idx << std::endl; + B.add_vertex_to_facet(idx); + } while (++circ != end); + + B.end_facet(); + } + + B.end_surface(); + } +}; + +double facet_area(const cgal_shape_t::Facet_handle& f); + +void dump_facet(const cgal_shape_t::Facet_handle& f); + +struct remove_thickness { + typedef Kernel_::Point_3 Point; + typedef Kernel_::Plane_3 Plane; + typedef Kernel_::Vector_3 Vector; + typedef Kernel_::Segment_3 Segment; + typedef Kernel_::Ray_3 Ray; + + typedef CGAL::Polyhedron_3 Polyhedron; + typedef CGAL::AABB_face_graph_triangle_primitive Primitive; + typedef CGAL::AABB_traits Traits; + typedef CGAL::AABB_tree Tree; + typedef boost::optional::Type> Ray_intersection; + + cgal_shape_t polyhedron, polyhedron2, flattened; + + remove_thickness(const cgal_shape_t& p) + // edge_collapse(p) still does not work :( + : polyhedron(p) + , polyhedron2(p) { + CGAL::Polygon_mesh_processing::triangulate_faces(polyhedron); + CGAL::Polygon_mesh_processing::triangulate_faces(polyhedron2); + + std::list non_degenerate, degenerate, longitudonal; + std::set thin_sides; + + std::wcout << "ALL FACES:" << std::endl; + + for (auto& f : faces(polyhedron)) { + dump_facet(f); + if (facet_area(f) > 1.e-20) { + non_degenerate.push_back(f); + } else { + degenerate.push_front(f); + std::wcout << "Degenerate, area: " << facet_area(f) << std::endl; + } + } + + std::wcout << "NON DEGENERATE:" << std::endl; + for (auto& f : non_degenerate) { + dump_facet(f); + } + + cgal_shape_t enlarged_non_degenerate_triangles; + Build_Offset bo; + bo.input = non_degenerate; + enlarged_non_degenerate_triangles.delegate(bo); + + // @todo, first on non-enlarged faces, then on enlarged; to fix projection on concave surfaces where the enlarging operation shortens projection distances. + + Tree tree(faces(enlarged_non_degenerate_triangles).first, faces(enlarged_non_degenerate_triangles).second, enlarged_non_degenerate_triangles); + + std::map face_normals; + boost::associative_property_map> face_normals_map(face_normals); + CGAL::Polygon_mesh_processing::compute_face_normals(polyhedron, face_normals_map); + + for (auto& f : non_degenerate) { + auto O = CGAL::centroid( + f->facet_begin()->vertex()->point(), + f->facet_begin()->next()->vertex()->point(), + f->facet_begin()->next()->next()->vertex()->point() + ); + + Ray ray(O, -face_normals_map[f]); + + std::list intersections; + tree.all_intersections(ray, std::back_inserter(intersections)); + double N = std::numeric_limits::infinity(); + Point P; + for (auto& intersection : intersections) { + if (boost::get(&(intersection->first))) { + const Point* p = boost::get(&(intersection->first)); + const double d = std::sqrt(CGAL::to_double((*p - O).squared_length())); + if (d > 1.e-20 && d < N) { + N = d; + } + } + } + if (N != std::numeric_limits::infinity() && N > 1.e-4) { + thin_sides.insert(f); + } + } + + std::wcout << "THIN SIDES:" << std::endl; + for (auto& f : thin_sides) { + dump_facet(f); + } + + for (auto& f : non_degenerate) { + if (thin_sides.find(f) == thin_sides.end()) { + longitudonal.push_back(f); + } + } + + std::wcout << "LONGITUDONAL:" << std::endl; + for (auto& f : longitudonal) { + dump_facet(f); + } + + std::wcout << "faces " << faces(polyhedron).size() << "long " << longitudonal.size() << "thin " << thin_sides.size() << "non-degen " << non_degenerate.size() << std::endl; + + cgal_shape_t enlarged_indiv_triangles; + Build_Offset bo2; + bo2.input = longitudonal; + enlarged_indiv_triangles.delegate(bo2); + + { + std::ofstream ofs("enlarged.off"); + ofs.precision(17); + ofs << enlarged_indiv_triangles; + } + + Tree tree2(faces(enlarged_indiv_triangles).begin(), faces(enlarged_indiv_triangles).end(), enlarged_indiv_triangles); + + std::map new_points; + + for (Polyhedron::Facet_iterator fit = polyhedron.facets_begin(); + fit != polyhedron.facets_end(); + ++fit) { + if (CGAL::collinear( + fit->halfedge()->vertex()->point(), + fit->halfedge()->next()->vertex()->point(), + fit->halfedge()->opposite()->vertex()->point())) { + std::wcout << "degenerate triangle" << std::endl; + } + } + + for (auto& v : vertices(polyhedron)) { + auto O = v->point(); + + Kernel_::Vector_3 norm; + Kernel_::Vector_3 accum; + int count = 0; + CGAL::Face_around_target_circulator it(v->halfedge(), polyhedron), end(it); + do { + cgal_shape_t::Facet_handle fh = (*it)->halfedge()->facet(); + + auto jt = std::find(non_degenerate.begin(), non_degenerate.end(), fh); + std::wcout << "non degen: " << (jt != non_degenerate.end()) << std::endl; + auto kt = std::find(thin_sides.begin(), thin_sides.end(), fh); + std::wcout << "thin side: " << (kt != thin_sides.end()) << std::endl; + if (jt != non_degenerate.end() && kt == thin_sides.end()) { + // else degenerate, prevent div by zero, do not incorporate in vnorm. + // or else part of thin side + + auto p0 = (*it)->facet_begin()->vertex()->point(); + auto p1 = (*it)->facet_begin()->next()->vertex()->point(); + auto p2 = (*it)->facet_begin()->next()->next()->vertex()->point(); + + { + std::ostringstream oss; + oss.precision(8); + oss << "p0 " << p0.cartesian(0) << " " << p0.cartesian(1) << " " << p0.cartesian(2) << "\n"; + oss << "p1 " << p1.cartesian(0) << " " << p1.cartesian(1) << " " << p1.cartesian(2) << "\n"; + oss << "p2 " << p2.cartesian(0) << " " << p2.cartesian(1) << " " << p2.cartesian(2) << "\n"; + auto osss = oss.str(); + std::wcout << osss.c_str() << std::endl; + } + + auto fnorm = CGAL::cross_product(p0 - p1, p2 - p1); + fnorm /= std::sqrt(CGAL::to_double(fnorm.squared_length())); + + // const auto& fnorm = face_normals_map_2[*it]; + std::ostringstream oss; + oss.precision(8); + oss << fnorm.cartesian(0) << " " << fnorm.cartesian(1) << " " << fnorm.cartesian(2); + auto osss = oss.str(); + std::wcout << osss.c_str() << std::endl; + accum += fnorm; + + ++count; + } + + ++it; + } while (it != end); + + norm = accum / count; + std::wcout << "count " << count << std::endl; + + if (count == 0) { + // part of only degenerate or only thin sides + continue; + } + + // v->vertex_begin(); + Ray ray(O, norm); + std::ostringstream oss; + oss.precision(8); + oss << O << " -> " << norm; + auto osss = oss.str(); + std::wcout << osss.c_str() << std::endl; + + std::list intersections; + tree2.all_intersections(ray, std::back_inserter(intersections)); + double N = std::numeric_limits::infinity(); + Point P; + + bool used_intersection = false; + + if (intersections.size()) { + for (auto& intersection : intersections) { + if (boost::get(&(intersection->first))) { + const Point* p = boost::get(&(intersection->first)); + const double d = std::sqrt(CGAL::to_double((*p - O).squared_length())); + if (d < N && d > 1.e-20) { + N = d; + P = *p; + std::wcout << "intersection @ " << d << std::endl; + } + } + } + std::wcout << "-----------" << std::endl; + + // average the new point + new_points[O] = CGAL::ORIGIN + (((O - CGAL::ORIGIN) + (P - CGAL::ORIGIN))) / 2; + used_intersection = true; + } + + if (!used_intersection) { + std::wcout << "no intersection :(" << std::endl; + } + } + + auto thin_sides_degenerate = thin_sides; + thin_sides_degenerate.insert(degenerate.begin(), degenerate.end()); + + // @todo choose connected / connected_opposing based on largest combined area of facets? + + if (longitudonal.size() == 0) { + std::wcout << "no longitudonal faces detected :(" << std::endl; + return; + } + + auto connected = connected_faces(*longitudonal.begin(), thin_sides_degenerate); + decltype(connected) connected_opposing; + + for (auto& f : longitudonal) { + if (std::find(connected.begin(), connected.end(), f) == connected.end()) { + connected_opposing = connected_faces(f, thin_sides_degenerate); + + std::set longi(longitudonal.begin(), longitudonal.end()); + std::set both_sides(connected.begin(), connected.end()); + both_sides.insert(connected_opposing.begin(), connected_opposing.end()); + + if (longi == both_sides) { + std::wcout << "Facet connection functioning properly" << std::endl; + } else { + std::wcout << "Facet connection functioning incorrectly" << std::endl; + } + + break; + } + } + + Builder_With_Map b2; + b2.input = connected; + b2.mapping = new_points; + + flattened.delegate(b2); + } +}; + + +struct intersection_validator { + typedef std::list> > nefs_t; + typedef CGAL::Box_intersection_d::Box_with_handle_d Box; + + std::vector boxes; + nefs_t nefs; + + double total_map_time = 0.; + double total_geom_time = 0.; + double total_nef_time = 0.; + double total_minkowsky_time = 0.; + double total_box_time = 0.; + + std::set succesfully_processed; + + intersection_validator(IfcParse::IfcFile& f, std::initializer_list entities, double eps, bool no_progress, bool quiet, bool stderr_progress) { + + ifcopenshell::geometry::settings settings; + settings.set(ifcopenshell::geometry::settings::USE_WORLD_COORDS, false); + settings.set(ifcopenshell::geometry::settings::WELD_VERTICES, false); + settings.set(ifcopenshell::geometry::settings::SEW_SHELLS, true); + settings.set(ifcopenshell::geometry::settings::CONVERT_BACK_UNITS, true); + settings.set(ifcopenshell::geometry::settings::DISABLE_TRIANGULATION, true); + settings.set(ifcopenshell::geometry::settings::DISABLE_OPENING_SUBTRACTIONS, true); + + std::vector spaces_and_walls = { + IfcGeom::entity_filter(true, false, entities) + }; + + ifcopenshell::geometry::Iterator context_iterator("cgal", settings, &f, spaces_and_walls); + + if (!context_iterator.initialize()) { + return; + } + + auto polycube = ifcopenshell::geometry::utils::create_cube(eps); + auto cube = ifcopenshell::geometry::utils::create_nef_polyhedron(polycube); + + size_t num_created = 0; + int old_progress = quiet ? 0 : -1; + + for (;; ++num_created) { + bool has_more = true; + if (num_created) { + has_more = context_iterator.next(); + } + ifcopenshell::geometry::NativeElement* geom_object = nullptr; + if (has_more) { + geom_object = context_iterator.get_native(); + } + if (!geom_object) { + break; + } + + std::stringstream ss; + ss << geom_object->product()->data().toString(); + auto sss = ss.str(); + std::wcout << sss.c_str() << std::endl; + + for (auto& g : geom_object->geometry()) { + auto s = ((ifcopenshell::geometry::CgalShape*) g.Shape())->shape(); + const auto& m = *g.Placement().components; + const auto& n = *geom_object->transformation().data().components; + + const cgal_placement_t trsf( + m(0, 0), m(0, 1), m(0, 2), m(0, 3), + m(1, 0), m(1, 1), m(1, 2), m(1, 3), + m(2, 0), m(2, 1), m(2, 2), m(2, 3)); + + const cgal_placement_t trsf2( + n(0, 0), n(0, 1), n(0, 2), n(0, 3), + n(1, 0), n(1, 1), n(1, 2), n(1, 3), + n(2, 0), n(2, 1), n(2, 2), n(2, 3)); + + // Apply transformation + for (auto &vertex : vertices(s)) { + vertex->point() = vertex->point().transform(trsf).transform(trsf2); + } + + std::clock_t nef_begin = std::clock(); + CGAL::Nef_polyhedron_3 nef = ifcopenshell::geometry::utils::create_nef_polyhedron(s); + std::clock_t nef_end = std::clock(); + total_nef_time += (nef_end - nef_begin) / (double) CLOCKS_PER_SEC; + if (nef.is_empty()) { + std::wcout << "Failed to create nef" << std::endl; + continue; + } + + succesfully_processed.insert(geom_object->product()); + + nef = CGAL::minkowski_sum_3(nef, cube); + std::clock_t minkowski_end = std::clock(); + total_minkowsky_time += (minkowski_end - nef_end) / (double) CLOCKS_PER_SEC; + + std::wcout << "product: " << geom_object->product() << std::endl; + nefs.push_back({ geom_object->product(), nef }); + + Box b(&*(nefs.rbegin())); + // id_map[b.id()] = ; + + for (auto &vertex : vertices(s)) { + double p[3] = { + CGAL::to_double(vertex->point().cartesian(0)), + CGAL::to_double(vertex->point().cartesian(1)), + CGAL::to_double(vertex->point().cartesian(2)) + }; + b.extend(p); + } + + boxes.push_back(enlarge(b)); + + /* + std::ostringstream ss; + ss << geom_object->product()->data().toString() << std::endl << b.min_coord(0) << " - " << b.max_coord(0) << std::endl; + auto sss = ss.str(); + std::wcout << sss.c_str(); + */ + } + + if (!no_progress) { + if (quiet) { + const int progress = context_iterator.progress(); + for (; old_progress < progress; ++old_progress) { + std::cout << "."; + if (stderr_progress) + std::cerr << "."; + } + std::cout << std::flush; + if (stderr_progress) + std::cerr << std::flush; + } else { + const int progress = context_iterator.progress() / 2; + if (old_progress != progress) Logger::ProgressBar(progress); + old_progress = progress; + } + } + } + + if (!no_progress && quiet) { + for (; old_progress < 100; ++old_progress) { + std::cout << "."; + if (stderr_progress) + std::cerr << "."; + } + std::cout << std::flush; + if (stderr_progress) + std::cerr << std::flush; + } else { + Logger::Status("\rDone fixing space boundaries for " + boost::lexical_cast(num_created) + + " objects "); + } + + total_geom_time = context_iterator.converter().total_geom_time; + total_map_time = context_iterator.converter().total_map_time; + } + + template + void operator()(Fn fn) { + std::clock_t box_overlap_begin = std::clock(); + CGAL::box_self_intersection_d(boxes.begin(), boxes.end(), [](Box& x, Box& y) {}); + std::clock_t box_overlap_end = std::clock(); + total_box_time += (box_overlap_end - box_overlap_begin) / (double) CLOCKS_PER_SEC; + CGAL::box_self_intersection_d(boxes.begin(), boxes.end(), fn); + } +}; diff --git a/src/ifcgeom/abstract_mapping.cpp b/src/ifcgeom/abstract_mapping.cpp new file mode 100644 index 0000000000..215b9d5caf --- /dev/null +++ b/src/ifcgeom/abstract_mapping.cpp @@ -0,0 +1,35 @@ +#include "abstract_mapping.h" + +#include "../ifcparse/IfcFile.h" + +ifcopenshell::geometry::impl::MappingFactoryImplementation& ifcopenshell::geometry::impl::mapping_implementations() { + static MappingFactoryImplementation impl; + return impl; +} + +extern void init_MappingImplementation_Ifc2x3(ifcopenshell::geometry::impl::MappingFactoryImplementation*); +extern void init_MappingImplementation_Ifc4(ifcopenshell::geometry::impl::MappingFactoryImplementation*); +extern void init_MappingImplementation_Ifc4x1(ifcopenshell::geometry::impl::MappingFactoryImplementation*); +extern void init_MappingImplementation_Ifc4x2(ifcopenshell::geometry::impl::MappingFactoryImplementation*); + +ifcopenshell::geometry::impl::MappingFactoryImplementation::MappingFactoryImplementation() { + init_MappingImplementation_Ifc2x3(this); + init_MappingImplementation_Ifc4(this); + init_MappingImplementation_Ifc4x1(this); + init_MappingImplementation_Ifc4x2(this); +} + +void ifcopenshell::geometry::impl::MappingFactoryImplementation::bind(const std::string& schema_name, ifcopenshell::geometry::impl::mapping_fn fn) { + const std::string schema_name_lower = boost::to_lower_copy(schema_name); + this->insert(std::make_pair(schema_name_lower, fn)); +} + +ifcopenshell::geometry::abstract_mapping* ifcopenshell::geometry::impl::MappingFactoryImplementation::construct(IfcParse::IfcFile* file, settings& s) { + const std::string schema_name_lower = boost::to_lower_copy(file->schema()->name()); + std::map::const_iterator it; + it = this->find(schema_name_lower); + if (it == end()) { + throw IfcParse::IfcException("No geometry mapping registered for " + schema_name_lower); + } + return it->second(file, s); +} diff --git a/src/ifcgeom/abstract_mapping.h b/src/ifcgeom/abstract_mapping.h new file mode 100644 index 0000000000..9a39fb55d0 --- /dev/null +++ b/src/ifcgeom/abstract_mapping.h @@ -0,0 +1,55 @@ +#ifndef ABSTRACT_MAPPING_H +#define ABSTRACT_MAPPING_H + +#include "../ifcparse/IfcBaseClass.h" +#include "../ifcparse/IfcEntityList.h" +#include "../ifcgeom/taxonomy.h" +#include "../ifcgeom/settings.h" + +#include + +#include +#include + +namespace ifcopenshell { + +namespace geometry { + + struct geometry_conversion_task { + int index; + IfcUtil::IfcBaseEntity* representation; + IfcEntityList::ptr products; + }; + + typedef boost::function filter_t; + + class abstract_mapping { + protected: + settings settings_; + public: + abstract_mapping(settings& s) : settings_(s) {} + + virtual ifcopenshell::geometry::taxonomy::item* map(const IfcUtil::IfcBaseClass*) = 0; + virtual void get_representations(std::vector& tasks, std::vector& filters, settings& s) = 0; + virtual IfcUtil::IfcBaseEntity* get_decomposing_entity(IfcUtil::IfcBaseEntity* product, bool include_openings = true) = 0; + virtual std::map get_layers(IfcUtil::IfcBaseEntity*) = 0; + }; + + namespace impl { + typedef boost::function2 mapping_fn; + + class MappingFactoryImplementation : public std::map { + public: + MappingFactoryImplementation(); + void bind(const std::string& schema_name, mapping_fn); + abstract_mapping* construct(IfcParse::IfcFile*, settings&); + }; + + MappingFactoryImplementation& mapping_implementations(); + } + +} + +} + +#endif diff --git a/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp b/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp new file mode 100644 index 0000000000..f0217c9c0e --- /dev/null +++ b/src/ifcgeom/kernels/cgal/CgalConversionResult.cpp @@ -0,0 +1,102 @@ +#include "CgalConversionResult.h" + +#include "../../../ifcparse/IfcLogger.h" +#include "../../../ifcgeom/schema_agnostic/IfcGeomRepresentation.h" + +void ifcopenshell::geometry::CgalShape::Triangulate(const settings& settings, const ifcopenshell::geometry::taxonomy::matrix4& place, Representation::Triangulation* t, int surface_style_id) const { + // Copy is made because triangulate_faces() does not accept a const argument + cgal_shape_t s = shape_; + + if (!place.components->isIdentity()) { + const auto& m = *place.components; + + // @todo check + const cgal_placement_t trsf( + m(0, 0), m(0, 1), m(0, 2), m(0, 3), + m(1, 0), m(1, 1), m(1, 2), m(1, 3), + m(2, 0), m(2, 1), m(2, 2), m(2, 3)); + + // Apply transformation + for (auto &vertex : vertices(s)) { + vertex->point() = vertex->point().transform(trsf); + } + } + + if (!s.is_valid()) { + Logger::Message(Logger::LOG_ERROR, "Invalid Polyhedron_3 in object (before triangulation)"); + return; + } + + // Triangulate the shape and compute the normals +// std::map vertex_normals; +// boost::associative_property_map> vertex_normals_map(vertex_normals); + std::map face_normals; + boost::associative_property_map> face_normals_map(face_normals); + + bool success = false; + try { + success = CGAL::Polygon_mesh_processing::triangulate_faces(s); + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "Triangulation crashed"); + return; + } + + if (!success) { + Logger::Message(Logger::LOG_ERROR, "Triangulation failed"); + return; + } + // std::cout << "Triangulated model: " << s.size_of_facets() << " facets and " << s.size_of_vertices() << " vertices" << std::endl; + + if (!s.is_valid()) { + Logger::Message(Logger::LOG_ERROR, "Invalid Polyhedron_3 in object (after triangulation)"); + return; + } + +// CGAL::Polygon_mesh_processing::compute_normals(s, vertex_normals_map, face_normals_map); + try { + CGAL::Polygon_mesh_processing::compute_face_normals(s, face_normals_map); + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "Face normal calculation failed"); + return; + } + + int num_faces = 0, num_vertices = 0; + for (auto &face: faces(s)) { + if (!face->is_triangle()) { + std::cout << "Warning: non-triangular face!" << std::endl; + continue; + } + CGAL::Polyhedron_3::Halfedge_around_facet_const_circulator current_halfedge = face->facet_begin(); + int vertexidx[3]; + int i = 0; + do { + vertexidx[i++] = t->addVertex(surface_style_id, + CGAL::to_double(current_halfedge->vertex()->point().cartesian(0)), + CGAL::to_double(current_halfedge->vertex()->point().cartesian(1)), + CGAL::to_double(current_halfedge->vertex()->point().cartesian(2))); + + double nx = 0.; + double ny = 0.; + double nz = 1.; + // @todo normal calculation throws divide by zero? + // try { + if (false) { + nx = CGAL::to_double(face_normals_map[face].cartesian(0)); + ny = CGAL::to_double(face_normals_map[face].cartesian(1)); + nz = CGAL::to_double(face_normals_map[face].cartesian(2)); + } + // catch (...) { + // Logger::Error("Error during normal calculation"); + // } + t->addNormal(nx, ny, nz); + + ++num_vertices; + ++current_halfedge; + } while (current_halfedge != face->facet_begin()); + + t->addFace(surface_style_id, vertexidx[0], vertexidx[1], vertexidx[2]); + + ++num_faces; + } + +} diff --git a/src/ifcgeom/kernels/cgal/CgalConversionResult.h b/src/ifcgeom/kernels/cgal/CgalConversionResult.h new file mode 100644 index 0000000000..35aaa5d61d --- /dev/null +++ b/src/ifcgeom/kernels/cgal/CgalConversionResult.h @@ -0,0 +1,91 @@ +/******************************************************************************** +* * +* This file is part of IfcOpenShell. * +* * +* IfcOpenShell is free software: you can redistribute it and/or modify * +* it under the terms of the Lesser GNU General Public License as published by * +* the Free Software Foundation, either version 3.0 of the License, or * +* (at your option) any later version. * +* * +* IfcOpenShell is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* Lesser GNU General Public License for more details. * +* * +* You should have received a copy of the Lesser GNU General Public License * +* along with this program. If not, see . * +* * +********************************************************************************/ + +#ifndef CGALCONVERSIONRESULT_H +#define CGALCONVERSIONRESULT_H + +#include "../../../ifcgeom/schema_agnostic/IfcGeomElement.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel_; + +typedef Kernel_::Aff_transformation_3 cgal_placement_t; +typedef Kernel_::Point_3 cgal_point_t; +typedef Kernel_::Vector_3 cgal_direction_t; +typedef Kernel_::Vector_3 cgal_vector_t; +typedef Kernel_::Plane_3 cgal_plane_t; +typedef std::vector cgal_curve_t; +typedef std::vector cgal_wire_t; + +struct cgal_face_t { + cgal_wire_t outer; + std::vector inner; +}; + +typedef CGAL::Polyhedron_3 cgal_shape_t; +typedef boost::graph_traits>::vertex_descriptor cgal_vertex_descriptor_t; +typedef boost::graph_traits>::face_descriptor cgal_face_descriptor_t; + +#include "../../../ifcgeom/schema_agnostic/ConversionResult.h" + +namespace ifcopenshell { namespace geometry { + + class CgalShape : public ConversionResultShape { + public: + CgalShape(const cgal_shape_t& shape) + : shape_(shape) + {} + + const cgal_shape_t& shape() const { return shape_; } + operator const cgal_shape_t& () { return shape_; } + + virtual void Triangulate(const settings& settings, const ifcopenshell::geometry::taxonomy::matrix4& place, Representation::Triangulation* t, int surface_style_id) const; + + virtual void Serialize(std::string&) const { + throw std::runtime_error("Not implemented"); + } + + virtual ConversionResultShape* clone() const { + return new CgalShape(shape_); + } + + virtual bool is_manifold() const { + throw std::runtime_error("Not implemented"); + } + + virtual int surface_genus() const { + throw std::runtime_error("Not implemented"); + } + private: + cgal_shape_t shape_; + }; + +}} + +#endif diff --git a/src/ifcgeom/kernels/cgal/CgalKernel.cpp b/src/ifcgeom/kernels/cgal/CgalKernel.cpp new file mode 100644 index 0000000000..a99c09dccb --- /dev/null +++ b/src/ifcgeom/kernels/cgal/CgalKernel.cpp @@ -0,0 +1,863 @@ +/******************************************************************************** + * * + * This file is part of IfcOpenShell. * + * * + * IfcOpenShell is free software: you can redistribute it and/or modify * + * it under the terms of the Lesser GNU General Public License as published by * + * the Free Software Foundation, either version 3.0 of the License, or * + * (at your option) any later version. * + * * + * IfcOpenShell is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * Lesser GNU General Public License for more details. * + * * + * You should have received a copy of the Lesser GNU General Public License * + * along with this program. If not, see . * + * * + ********************************************************************************/ +#define _USE_MATH_DEFINES +#include + +#include "CgalKernel.h" + +#include "../../../ifcparse/IfcLogger.h" +#include "../../../ifcgeom/kernels/cgal/CgalConversionResult.h" + +#include + +using namespace ifcopenshell::geometry; +using namespace ifcopenshell::geometry::kernels; + +void CgalKernel::remove_duplicate_points_from_loop(cgal_wire_t& polygon) { + std::set points; + for (int i = 0; i < polygon.size(); ++i) { + if (points.count(polygon[i])) { + polygon.erase(polygon.begin() + i); + --i; + } else points.insert(polygon[i]); + } +} + +CGAL::Polyhedron_3 ifcopenshell::geometry::utils::create_polyhedron(std::list &face_list) { + + // Naive creation + CGAL::Polyhedron_3 polyhedron; + PolyhedronBuilder builder(&face_list); + polyhedron.delegate(builder); + + // Stitch edges + // std::cout << "Before: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl; + CGAL::Polygon_mesh_processing::stitch_borders(polyhedron); + if (!polyhedron.is_valid()) { + Logger::Message(Logger::LOG_ERROR, "create_polyhedron: Polyhedron not valid!"); + // std::ofstream fresult; + // fresult.open("/Users/ken/Desktop/invalid.off"); + // fresult << polyhedron << std::endl; + // fresult.close(); + return CGAL::Polyhedron_3(); + } if (polyhedron.is_closed()) { + if (!CGAL::Polygon_mesh_processing::is_outward_oriented(polyhedron)) { + CGAL::Polygon_mesh_processing::reverse_face_orientations(polyhedron); + } + } + + // std::cout << "After: " << polyhedron.size_of_vertices() << " vertices and " << polyhedron.size_of_facets() << " facets" << std::endl; + + return polyhedron; +} + +CGAL::Polyhedron_3 ifcopenshell::geometry::utils::create_polyhedron(const CGAL::Nef_polyhedron_3& nef_polyhedron) { + if (nef_polyhedron.is_simple()) { + try { + CGAL::Polyhedron_3 polyhedron; + nef_polyhedron.convert_to_polyhedron(polyhedron); + return polyhedron; + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "Conversion from Nef to polyhedron failed!"); + return CGAL::Polyhedron_3(); + } + } else { + Logger::Message(Logger::LOG_ERROR, "Nef polyhedron not simple: cannot create polyhedron!"); + return CGAL::Polyhedron_3(); + } +} + +CGAL::Nef_polyhedron_3 ifcopenshell::geometry::utils::create_nef_polyhedron(std::list &face_list) { + CGAL::Polyhedron_3 polyhedron = create_polyhedron(face_list); + CGAL::Polygon_mesh_processing::triangulate_faces(polyhedron); + CGAL::Nef_polyhedron_3 nef_polyhedron; + try { + nef_polyhedron = CGAL::Nef_polyhedron_3(polyhedron); + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "Conversion to Nef polyhedron failed!"); + } + return nef_polyhedron; +} + +CGAL::Nef_polyhedron_3 ifcopenshell::geometry::utils::create_nef_polyhedron(CGAL::Polyhedron_3 &polyhedron) { + if (polyhedron.is_valid() && polyhedron.is_closed()) { + // @todo is it necessary to triangulat? + CGAL::Polygon_mesh_processing::triangulate_faces(polyhedron); + CGAL::Nef_polyhedron_3 nef_polyhedron; + try { + nef_polyhedron = CGAL::Nef_polyhedron_3(polyhedron); + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "Conversion to Nef polyhedron failed!"); + } + return nef_polyhedron; + } else { + Logger::Message(Logger::LOG_ERROR, "Polyhedron not valid: cannot create Nef polyhedron!"); + return CGAL::Nef_polyhedron_3(); + } +} + +bool CgalKernel::convert(const taxonomy::shell* l, cgal_shape_t& shape) { + auto faces = l->children_as(); + + std::list face_list; + for (auto& f : faces) { + bool success = false; + cgal_face_t face; + + try { + success = convert(f, face); + } catch (...) {} + + if (!success) { + Logger::Message(Logger::LOG_WARNING, "Failed to convert face:", f->instance); + continue; + } + + // std::cout << "Face in ConnectedFaceSet: " << std::endl; + // for (auto &point: face.outer) { + // std::cout << "\tPoint(" << point << ")" << std::endl; + // } + + face_list.push_back(face); + } + + shape = utils::create_polyhedron(face_list); + return true; +} + +bool CgalKernel::convert(const taxonomy::face* face, cgal_face_t& result) { + auto bounds = face->children_as(); + + int num_outer_bounds = 0; + + for (auto& bound : bounds) { + if (bound->external.get_value_or(false)) num_outer_bounds++; + } + + if (num_outer_bounds != 1) { + Logger::Message(Logger::LOG_ERROR, "Invalid configuration of boundaries for:", face->instance); + return false; + } + + cgal_face_t mf; + + for (auto& bound : bounds) { + + const bool is_interior = !bound->external.get_value_or(false); + + cgal_wire_t wire; + if (!convert(bound, wire)) { + Logger::Message(Logger::LOG_ERROR, "Failed to process face boundary loop", bound->instance); + return false; + } + + if (!is_interior) { + mf.outer = wire; + } else { + mf.inner.push_back(wire); + } + } + + result = mf; + + // std::cout << "Face: " << std::endl; + // for (auto &point: face.outer) { + // std::cout << "\tPoint(" << point << ")" << std::endl; + // } + + return true; +} + +namespace { + // @todo obsolete? + bool convert_curve(CgalKernel* kernel, const taxonomy::item* curve, cgal_wire_t& builder) { + if (curve->kind() == taxonomy::EDGE) { + auto e = (taxonomy::edge*) curve; + if (true || e->basis == nullptr) { + if (builder.empty()) { + const auto& p = boost::get(e->start); + cgal_point_t pnt((*p.components)(0), (*p.components)(1), (*p.components)(2)); + builder.push_back(pnt); + } + const auto& p = boost::get(e->end); + cgal_point_t pnt((*p.components)(0), (*p.components)(1), (*p.components)(2)); + builder.push_back(pnt); + } else if (e->basis->kind() == taxonomy::CIRCLE) { + // @todo + } else if (e->basis->kind() == taxonomy::ELLIPSE) { + + } else { + throw std::runtime_error("Not implemented basis kind"); + } + } else if (curve->kind() == taxonomy::LOOP) { + const auto& edges = ((taxonomy::loop*) curve)->children; + for (auto& c : edges) { + convert_curve(kernel, c, builder); + } + } else { + throw std::runtime_error("Not implemented curve"); + } + } +} + +namespace { + typedef std::pair parameter_range; + + static const parameter_range unbounded = { + -std::numeric_limits::infinity(), + +std::numeric_limits::infinity() + }; + + void evaluate_curve(const taxonomy::line& c, double u, taxonomy::point3& p) { + Eigen::Vector4d xy{ u, 0, 0, 1. }; + *p.components = (*c.matrix.components * xy).head<3>(); + } + + void evaluate_curve(const taxonomy::circle& c, double u, taxonomy::point3& p) { + Eigen::Vector4d xy{ c.radius * std::cos(u), c.radius * std::sin(u), 0, 1. }; + *p.components = (*c.matrix.components * xy).head<3>(); + } + + void evaluate_curve(const taxonomy::ellipse& c, double u, taxonomy::point3& p) { + Eigen::Vector4d xy{ c.radius * std::cos(u), c.radius2 * std::sin(u), 0, 1. }; + *p.components = (*c.matrix.components * xy).head<3>(); + } + + // ---- + + void project_onto_curve(const taxonomy::line& c, const taxonomy::point3& p, double& u) { + u = (c.matrix.components->inverse() * p.components->homogeneous())(0); + } + + void project_onto_curve(const taxonomy::circle& c, const taxonomy::point3& p, double& u) { + Eigen::Vector2d xy = (c.matrix.components->inverse() * p.components->homogeneous()).head<2>(); + u = std::atan2(xy(1), xy(0)); + } + + void project_onto_curve(const taxonomy::ellipse& c, const taxonomy::point3& p, double& u) { + Eigen::Vector2d xy = (c.matrix.components->inverse() * p.components->homogeneous()).head<2>(); + u = std::atan2(xy(1), xy(0)); + } + + struct point_projection_visitor_ { + taxonomy::point3 p; + double u; + + void operator()(const taxonomy::line& c) { + project_onto_curve(c, p, u); + } + + void operator()(const taxonomy::circle& c) { + project_onto_curve(c, p, u); + } + + void operator()(const taxonomy::ellipse& c) { + project_onto_curve(c, p, u); + } + + void operator()(const taxonomy::item& c) { + throw std::runtime_error("Point projection not implemented on this geometry type"); + } + }; + + struct point_projection_visitor { + taxonomy::item* curve; + double u; + + void operator()(const taxonomy::point3& p) { + point_projection_visitor_ v{ p }; + dispatch_curve_creation::dispatch(curve, v); + u = v.u; + } + + void operator()(const double& u) { + this->u = u; + } + }; + + struct cgal_curve_creation_visitor { + static const int FULL_CIRCLE_NUM_SEGMENTS = 32; + parameter_range param; + + std::vector points; + + cgal_curve_creation_visitor() : param(unbounded) {} + cgal_curve_creation_visitor(const parameter_range& p) : param(p) {} + + void operator()(const taxonomy::line& l) { + if (param == unbounded) { + throw std::runtime_error("Cannot represent infinite line segment"); + } + taxonomy::point3 start, end; + evaluate_curve(l, param.first, start); + evaluate_curve(l, param.second, end); + points.push_back(start); + points.push_back(end); + } + + template + void evaluate_conic(const T& t) { + double a, b; + if (param == unbounded) { + a = 0.; + b = 2 * M_PI; + } else { + std::tie(a, b) = param; + } + int num_segments = (int)std::ceil(std::fabs(a - b) / (2 * M_PI) * FULL_CIRCLE_NUM_SEGMENTS); + double du = (b - a) / num_segments; + taxonomy::point3 P; + // @nb for loop is not inclusive of the both end points + evaluate_curve(t, a, P); + points.push_back(P); + for (int i = 1; i < num_segments; ++i) { + double u = a + du * i; + evaluate_curve(t, u, P); + points.push_back(P); + } + evaluate_curve(t, b, P); + points.push_back(P); + } + + void operator()(const taxonomy::circle& c) { + evaluate_conic(c); + } + + void operator()(const taxonomy::ellipse& e) { + evaluate_conic(e); + } + + void operator()(const taxonomy::trimmed_curve& e) { + point_projection_visitor v1, v2; + boost::apply_visitor(v1, e.start); + boost::apply_visitor(v2, e.end); + + cgal_curve_creation_visitor v({ v1.u, v2.u }); + + dispatch_curve_creation::dispatch(e.basis, v); + this->points = v.points; + } + + void operator()(const taxonomy::item& e) { + throw std::runtime_error("Not supported"); + } + }; + + void convert_curve(taxonomy::item* i, std::vector& points) { + cgal_curve_creation_visitor v; + dispatch_curve_creation::dispatch(i, v); + points = v.points; + } + + // @nb mutates a + void extend_wire(std::vector& a, const std::vector& b) { + if (a.empty()) { + a = b; + } + if (b.empty()) { + return; + } + double d = (*a.back().components - *b.front().components).norm(); + size_t offset = d < 1.e-5 ? 1 : 0; + a.insert(a.end(), b.begin() + offset, b.end()); + } +} + +bool CgalKernel::convert(const taxonomy::loop* loop, cgal_wire_t& result) { + // @todo only implement polygonal loops + + auto edges = loop->children_as(); + std::vector points; + + for (auto& e : edges) { + if (e->basis) { + std::vector edge; + convert_curve(e->basis, points); + extend_wire(points, edge); + } else { + extend_wire(points, { + boost::get(e->start), + boost::get(e->end) + }); + } + } + + if (points.size() >= 2) { + // the edges -> conversion left us with a duplicate global begin,end point. + double d = (*points.back().components - *points.front().components).norm(); + points.erase(points.end() - 1); + } + + // Parse and store the points in a sequence + cgal_wire_t polygon = std::vector(); + for (auto& p : points) { + cgal_point_t pnt((*p.components)(0), (*p.components)(1), (*p.components)(2)); + polygon.push_back(pnt); + } + + // A loop should consist of at least three vertices + std::size_t original_count = polygon.size(); + if (original_count < 3) { + Logger::Message(Logger::LOG_ERROR, "Not enough edges for:", loop->instance); + return false; + } + + // Remove points that are too close to one another + remove_duplicate_points_from_loop(polygon); + + std::size_t count = polygon.size(); + if (original_count - count != 0) { + std::stringstream ss; ss << (original_count - count) << " edges removed for:"; + Logger::Message(Logger::LOG_WARNING, ss.str(), loop->instance); + } + + if (count < 3) { + Logger::Message(Logger::LOG_ERROR, "Not enough edges for:", loop->instance); + return false; + } + + result = polygon; + + // std::cout << "PolyLoop: " << std::endl; + // for (auto &point: polygon) { + // std::cout << "\tPoint(" << point << ")" << std::endl; + // } + + return true; +} + + +bool CgalKernel::convert_impl(const taxonomy::shell *shell, ifcopenshell::geometry::ConversionResults& results) { + cgal_shape_t shape; + if (!convert(shell, shape)) { + return false; + } + results.emplace_back(ConversionResult( + shell->instance->data().id(), + shell->matrix, + new CgalShape(shape), + shell->surface_style + )); + return true; +} + +bool CgalKernel::convert_impl(const taxonomy::extrusion* extrusion, ifcopenshell::geometry::ConversionResults& results) { + cgal_shape_t shape; + if (!convert(extrusion, shape)) { + return false; + } + results.emplace_back(ConversionResult( + extrusion->instance->data().id(), + extrusion->matrix, + new CgalShape(shape), + extrusion->surface_style + )); + return true; +} + +bool CgalKernel::convert(const taxonomy::extrusion* extrusion, cgal_shape_t &shape) { + const double& height = extrusion->depth; + if (height < precision_) { + Logger::Message(Logger::LOG_ERROR, "Non-positive extrusion height encountered for:", extrusion->instance); + return false; + } + + // Outer + cgal_face_t bottom_face; + if (!convert(&extrusion->basis, bottom_face)) { + return false; + } + // std::cout << "Face vertices: " << face.outer.size() << std::endl; + + auto fs = *extrusion->direction.components; + cgal_direction_t dir(fs(0), fs(1), fs(2)); + // std::cout << "Direction: " << dir << std::endl; + + std::list face_list; + face_list.push_back(bottom_face); + + for (std::vector::const_iterator current_vertex = bottom_face.outer.begin(); + current_vertex != bottom_face.outer.end(); + ++current_vertex) { + std::vector::const_iterator next_vertex = current_vertex; + ++next_vertex; + if (next_vertex == bottom_face.outer.end()) { + next_vertex = bottom_face.outer.begin(); + } cgal_face_t side_face; + side_face.outer.push_back(*next_vertex); + side_face.outer.push_back(*current_vertex); + side_face.outer.push_back(*current_vertex + height * dir); + side_face.outer.push_back(*next_vertex + height * dir); + face_list.push_back(side_face); + } + + cgal_face_t top_face; + for (std::vector::const_reverse_iterator vertex = bottom_face.outer.rbegin(); + vertex != bottom_face.outer.rend(); + ++vertex) { + top_face.outer.push_back(*vertex + height * dir); + } face_list.push_back(top_face); + + if (bottom_face.inner.empty()) { + shape = utils::create_polyhedron(face_list); + // if (has_position) for (auto &vertex : vertices(shape)) vertex->point() = vertex->point().transform(trsf); + return true; + } + + CGAL::Nef_polyhedron_3 nef_shape = utils::create_nef_polyhedron(face_list); + + // Inner + // TODO: Would be faster to triangulate top/bottom face template rather than use Nef polyhedra for subtraction + for (auto &inner : bottom_face.inner) { + // std::cout << "Inner wire" << std::endl; + face_list.clear(); + + cgal_face_t hole_bottom_face; + hole_bottom_face.outer = inner; + remove_duplicate_points_from_loop(hole_bottom_face.outer); + face_list.push_back(hole_bottom_face); + + for (std::vector::const_iterator current_vertex = inner.begin(); + current_vertex != inner.end(); + ++current_vertex) { + std::vector::const_iterator next_vertex = current_vertex; + ++next_vertex; + if (next_vertex == inner.end()) { + next_vertex = inner.begin(); + } cgal_face_t hole_side_face; + hole_side_face.outer.push_back(*next_vertex); + hole_side_face.outer.push_back(*current_vertex); + hole_side_face.outer.push_back(*current_vertex + height * dir); + hole_side_face.outer.push_back(*next_vertex + height * dir); + face_list.push_back(hole_side_face); + } + + cgal_face_t hole_top_face; + for (std::vector::const_reverse_iterator vertex = inner.rbegin(); + vertex != inner.rend(); + ++vertex) { + hole_top_face.outer.push_back(*vertex + height * dir); + } face_list.push_back(hole_top_face); + + try { + nef_shape -= utils::create_nef_polyhedron(face_list); + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "IfcExtrudedAreaSolid: cannot subtract opening for:", extrusion->instance); + return false; + } + } + + /*if (has_position) { + // IfcSweptAreaSolid.Position (trsf) is an IfcAxis2Placement3D + // and therefore has a unit scale factor + nef_shape.transform(trsf); + }*/ + + try { + nef_shape.convert_to_polyhedron(shape); + return true; + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "IfcExtrudedAreaSolid: cannot convert Nef to polyhedron for:", extrusion->instance); + return false; + } + +} + +CGAL::Polyhedron_3 ifcopenshell::geometry::utils::create_cube(double d) { + cgal_face_t bottom_face; + bottom_face.outer.push_back(Kernel_::Point_3(-d, -d, -d)); + bottom_face.outer.push_back(Kernel_::Point_3(+d, -d, -d)); + bottom_face.outer.push_back(Kernel_::Point_3(+d, +d, -d)); + bottom_face.outer.push_back(Kernel_::Point_3(-d, +d, -d)); + + cgal_direction_t dir(0, 0, 2 * d); + + std::list face_list = { bottom_face }; + + for (std::vector::const_iterator current_vertex = bottom_face.outer.begin(); + current_vertex != bottom_face.outer.end(); + ++current_vertex) + { + std::vector::const_iterator next_vertex = current_vertex; + ++next_vertex; + + if (next_vertex == bottom_face.outer.end()) { + next_vertex = bottom_face.outer.begin(); + } + + cgal_face_t side_face; + + side_face.outer.push_back(*next_vertex); + side_face.outer.push_back(*current_vertex); + side_face.outer.push_back(*current_vertex + dir); + side_face.outer.push_back(*next_vertex + dir); + + face_list.push_back(side_face); + } + + cgal_face_t top_face; + + for (std::vector::const_reverse_iterator vertex = bottom_face.outer.rbegin(); + vertex != bottom_face.outer.rend(); + ++vertex) + { + top_face.outer.push_back(*vertex + dir); + } + + face_list.push_back(top_face); + + return create_polyhedron(face_list); +} + + +CGAL::Polyhedron_3 ifcopenshell::geometry::utils::create_cube(const Kernel_::Point_3& lower, const Kernel_::Point_3& upper) { + cgal_face_t bottom_face; + + auto a0 = lower.cartesian(0); + auto a1 = lower.cartesian(1); + auto a2 = lower.cartesian(2); + + auto b0 = upper.cartesian(0); + auto b1 = upper.cartesian(1); + auto b2 = upper.cartesian(2); + + bottom_face.outer.push_back(Kernel_::Point_3(a0, a1, a2)); + bottom_face.outer.push_back(Kernel_::Point_3(b0, a1, a2)); + bottom_face.outer.push_back(Kernel_::Point_3(b0, b1, a2)); + bottom_face.outer.push_back(Kernel_::Point_3(a0, b1, a2)); + + cgal_direction_t dir(0, 0, b2 - a2); + + std::list face_list = { bottom_face }; + + for (std::vector::const_iterator current_vertex = bottom_face.outer.begin(); + current_vertex != bottom_face.outer.end(); + ++current_vertex) + { + std::vector::const_iterator next_vertex = current_vertex; + ++next_vertex; + + if (next_vertex == bottom_face.outer.end()) { + next_vertex = bottom_face.outer.begin(); + } + + cgal_face_t side_face; + + side_face.outer.push_back(*next_vertex); + side_face.outer.push_back(*current_vertex); + side_face.outer.push_back(*current_vertex + dir); + side_face.outer.push_back(*next_vertex + dir); + + face_list.push_back(side_face); + } + + cgal_face_t top_face; + + for (std::vector::const_reverse_iterator vertex = bottom_face.outer.rbegin(); + vertex != bottom_face.outer.rend(); + ++vertex) + { + top_face.outer.push_back(*vertex + dir); + } + + face_list.push_back(top_face); + + return create_polyhedron(face_list); +} + +bool CgalKernel::thin_solid(const CGAL::Nef_polyhedron_3& a, CGAL::Nef_polyhedron_3& result) { + // @todo this should be possible as a minkowski sum of facet & cube. rather than a set of boolean ops. + + auto a_nonconst = a; + auto ax = CGAL::minkowski_sum_3(a_nonconst, precision_cube_); + auto x = ax - a; + + result = x; + return true; + + auto yxy = CGAL::minkowski_sum_3(x, precision_cube_); + auto y = yxy * a; + auto zyz = CGAL::minkowski_sum_3(y, precision_cube_); + result = yxy * zyz; + + return true; +} + +bool CgalKernel::preprocess_boolean_operand(const IfcUtil::IfcBaseClass* log_reference, const cgal_shape_t& shape_const, CGAL::Nef_polyhedron_3& result, bool dilate) { + cgal_shape_t shape = shape_const; + + if (!shape.is_valid()) { + Logger::Message(Logger::LOG_ERROR, "Conversion to Nef will fail. Invalid geometry:", log_reference); + return false; + } + + if (!shape.is_closed()) { + // TODO: There can be substractions to remove parts of non-volumetric objects. Maybe iterate over all faces of an entity and put them in a Nef_polyhedron_3 through Boolean union? Highly inefficient but maybe desirable... + Logger::Message(Logger::LOG_ERROR, "Subtraction of openings not supported for non-closed geometry:", log_reference); + return false; + } + + bool success = false; + + try { + success = CGAL::Polygon_mesh_processing::triangulate_faces(shape); + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "Triangulation of geometry crashed:", log_reference); + return false; + } + + if (!success) { + Logger::Message(Logger::LOG_ERROR, "Triangulation of geometry failed:", log_reference); + return false; + } + + if (CGAL::Polygon_mesh_processing::does_self_intersect(shape)) { + Logger::Message(Logger::LOG_ERROR, "Conversion to Nef will fail. Self-intersecting geometry:", log_reference); + return false; + } + + try { + result = CGAL::Nef_polyhedron_3(shape); + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "Could not convert geometry to Nef:", log_reference); + return false; + } + + if (dilate) { + try { + // @todo don't dilate in 3 dimensions but only in the XY plane, orthogonal to wall axis. + result = CGAL::minkowski_sum_3(result, precision_cube_); + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "Could not dilate boolean operand", log_reference); + return false; + } + } + + try { + cgal_shape_t convert_back; + result.convert_to_polyhedron(convert_back); + } catch (...) { + Logger::Message(Logger::LOG_WARNING, "Final conversion will likely fail. Could not convert geometry from Nef:", log_reference); + } + + return true; +} + +namespace { + bool convert_placement(const ifcopenshell::geometry::taxonomy::matrix4& place, cgal_placement_t& trsf) { + const auto& m = *place.components; + + // @todo check + trsf = cgal_placement_t( + m(0, 0), m(0, 1), m(0, 2), m(0, 3), + m(1, 0), m(1, 1), m(1, 2), m(1, 3), + m(2, 0), m(2, 1), m(2, 2), m(2, 3)); + + return true; + } +} + + + +bool CgalKernel::convert_impl(const taxonomy::boolean_result* br, ifcopenshell::geometry::ConversionResults& results) { + bool first = true; + + CGAL::Nef_polyhedron_3 a; + + taxonomy::style first_item_style; + + for (auto& c : br->children) { + // AbstractKernel::convert(c, results); + // continue; + + ifcopenshell::geometry::ConversionResults cr; + // @todo half-space detection + AbstractKernel::convert(c, cr); + + if (first && br->operation == taxonomy::boolean_result::SUBTRACTION) { + first_item_style = ((taxonomy::geom_item*)c)->surface_style; + if (!first_item_style.diffuse && c->kind() == taxonomy::COLLECTION) { + first_item_style = ((taxonomy::geom_item*) ((taxonomy::collection*)c)->children[0])->surface_style; + } + } + + for (auto it = cr.begin(); it != cr.end(); ++it) { + const cgal_shape_t& entity_shape_unlocated(((CgalShape*)it->Shape())->shape()); + cgal_shape_t entity_shape(entity_shape_unlocated); + if (!it->Placement().components->isIdentity()) { + cgal_placement_t trsf; + convert_placement(it->Placement(), trsf); + for (auto &vertex : vertices(entity_shape)) { + if (false) { + auto x = CGAL::to_double(vertex->point().x()); + auto y = CGAL::to_double(vertex->point().y()); + auto z = CGAL::to_double(vertex->point().z()); + std::wcout << x << " " << y << " " << z << std::endl; + } + vertex->point() = vertex->point().transform(trsf); + if (false) { + auto x = CGAL::to_double(vertex->point().x()); + auto y = CGAL::to_double(vertex->point().y()); + auto z = CGAL::to_double(vertex->point().z()); + std::wcout << x << " " << y << " " << z << std::endl; + } + } + } + + CGAL::Nef_polyhedron_3 nef; + preprocess_boolean_operand(c->instance, entity_shape, nef, + // Dilate boolean subtraction operands + (!first && br->operation == taxonomy::boolean_result::SUBTRACTION)); + + if (first) { + a = nef; + } else { + if (br->operation == taxonomy::boolean_result::SUBTRACTION) { + a -= nef; + } else if (br->operation == taxonomy::boolean_result::INTERSECTION) { + a *= nef; + } else if (br->operation == taxonomy::boolean_result::UNION) { + a += nef; + } + } + } + + first = false; + } + + cgal_shape_t a_poly, b_poly; + + // CGAL::Nef_polyhedron_3 b; + // thin_solid(a, b); + + try { + a.convert_to_polyhedron(a_poly); + } catch (...) { + Logger::Message(Logger::LOG_ERROR, "Could not convert geometry with openings from Nef:", br->instance); + return false; + } + + results.emplace_back(ConversionResult( + br->instance->data().id(), + br->matrix, + new CgalShape(a_poly), + br->surface_style.diffuse ? br->surface_style : first_item_style + )); + return true; +} \ No newline at end of file diff --git a/src/ifcgeom/kernels/cgal/CgalKernel.h b/src/ifcgeom/kernels/cgal/CgalKernel.h new file mode 100644 index 0000000000..588fd5b3a8 --- /dev/null +++ b/src/ifcgeom/kernels/cgal/CgalKernel.h @@ -0,0 +1,142 @@ +/******************************************************************************** +* * +* This file is part of IfcOpenShell. * +* * +* IfcOpenShell is free software: you can redistribute it and/or modify * +* it under the terms of the Lesser GNU General Public License as published by * +* the Free Software Foundation, either version 3.0 of the License, or * +* (at your option) any later version. * +* * +* IfcOpenShell is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* Lesser GNU General Public License for more details. * +* * +* You should have received a copy of the Lesser GNU General Public License * +* along with this program. If not, see . * +* * +********************************************************************************/ + +#ifndef CGAL_KERNEL_H +#define CGAL_KERNEL_H + +/* +#ifdef NO_CACHE + +#define IN_CACHE(T,E,t,e) +#define CACHE(T,E,e) + +#else + +#define IN_CACHE(T,E,t,e) std::map::const_iterator it = cache.T.find(E->entity->id());\ +if ( it != cache.T.end() ) { e = it->second; return true; } +#define CACHE(T,E,e) cache.T[E->entity->id()] = e; + +#endif +*/ + +#include + +#include "../../../ifcparse/macros.h" + +#include "../../../ifcgeom/kernel_agnostic/AbstractKernel.h" + +#include "../../../ifcgeom/schema_agnostic/IfcGeomElement.h" +#include "../../../ifcgeom/kernels/cgal/CgalConversionResult.h" + +struct PolyhedronBuilder : public CGAL::Modifier_base::HalfedgeDS> { +private: + std::list *face_list; +public: + PolyhedronBuilder(std::list *face_list) { + this->face_list = face_list; + } + + void operator()(CGAL::Polyhedron_3::HalfedgeDS &hds) { + std::list points; + std::list> facet_vertices; + CGAL::Polyhedron_incremental_builder_3::HalfedgeDS> builder(hds, true); + + for (auto &face: *face_list) { + facet_vertices.push_back(std::list()); + for (auto &point: face.outer) { + facet_vertices.back().push_back(points.size()); + points.push_back(point); + } + } + + builder.begin_surface(points.size(), facet_vertices.size()); + + for (auto &point: points) { +// std::cout << "Adding point " << point << std::endl; + builder.add_vertex(point); + } + + for (auto &facet: facet_vertices) { + builder.begin_facet(); +// std::cout << "Adding facet "; + for (auto &vertex: facet) { +// std::cout << vertex << " "; + builder.add_vertex_to_facet(vertex); + } +// std::cout << std::endl; + builder.end_facet(); + } + + builder.end_surface(); + } +}; + +namespace ifcopenshell { +namespace geometry { +namespace utils { + IFC_GEOM_API CGAL::Polyhedron_3 create_cube(double d); + IFC_GEOM_API CGAL::Polyhedron_3 create_cube(const Kernel_::Point_3& lower, const Kernel_::Point_3& upper); + IFC_GEOM_API CGAL::Polyhedron_3 create_polyhedron(std::list &face_list); + IFC_GEOM_API CGAL::Polyhedron_3 create_polyhedron(const CGAL::Nef_polyhedron_3 &nef_polyhedron); + IFC_GEOM_API CGAL::Nef_polyhedron_3 create_nef_polyhedron(std::list &face_list); + IFC_GEOM_API CGAL::Nef_polyhedron_3 create_nef_polyhedron(CGAL::Polyhedron_3 &polyhedron); +} + +namespace kernels { + + class IFC_GEOM_API CgalKernel : public AbstractKernel { + private: + double precision_; + size_t circle_segments_; + CGAL::Nef_polyhedron_3 precision_cube_; + + bool preprocess_boolean_operand(const IfcUtil::IfcBaseClass* log_reference, const cgal_shape_t& shape_const, CGAL::Nef_polyhedron_3& result, bool dilate); + bool thin_solid(const CGAL::Nef_polyhedron_3& a, CGAL::Nef_polyhedron_3& result); + public: + + CgalKernel() + : AbstractKernel("cgal") + // @todo + , precision_(1.e-5) + , circle_segments_(16) + { + auto cc = utils::create_cube(precision_); + precision_cube_ = CGAL::Nef_polyhedron_3(cc); + } + + void remove_duplicate_points_from_loop(cgal_wire_t& polygon); + + bool convert(const taxonomy::extrusion*, cgal_shape_t&); + bool convert(const taxonomy::face*, cgal_face_t&); + bool convert(const taxonomy::loop*, cgal_wire_t&); + // bool convert(const taxonomy::matrix4*, cgal_placement_t&); + bool convert(const taxonomy::shell*, cgal_shape_t&); + + // virtual bool convert_impl(const taxonomy::face*, ifcopenshell::geometry::ConversionResults&); + virtual bool convert_impl(const taxonomy::shell*, ifcopenshell::geometry::ConversionResults&); + virtual bool convert_impl(const taxonomy::extrusion*, ifcopenshell::geometry::ConversionResults&); + virtual bool convert_impl(const taxonomy::boolean_result*, ifcopenshell::geometry::ConversionResults&); + + const CGAL::Nef_polyhedron_3& precision_cube() const { return precision_cube_; } + }; + +} +} +} +#endif \ No newline at end of file