Glassvision/IfcOpenShell
1
0
Fork 0
mirror of https://github.com/IfcOpenShell/IfcOpenShell.git synced 2026-08-09 17:31:45 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
b301eb6520f9b9b68dbd0acea07e47adde0bfa74
IfcOpenShell/src/ifcgeom_schema_agnostic/boolean_utils.cpp
T
Thomas Krijnen b301eb6520 Update boolean_utils.cpp
2023-01-12 13:07:40 +01:00

1361 lines
38 KiB
C++
Raw Blame History

#include "boolean_utils.h"
#include "../ifcgeom_schema_agnostic/IfcGeomTree.h"
#include "../ifcgeom_schema_agnostic/base_utils.h"
#include <BRepBuilderAPI_Copy.hxx>
#include <TopExp_Explorer.hxx>
#include <GProp_GProps.hxx>
#include <BRepGProp.hxx>
#include <TopExp.hxx>
#include <TopoDS.hxx>
#include <Bnd_Box.hxx>
#include <Extrema_ExtPC.hxx>
#include <Geom_Plane.hxx>
#include <Geom_BSplineCurve.hxx>
#include <ShapeUpgrade_UnifySameDomain.hxx>
#include <GeomAPI_ExtremaCurveCurve.hxx>
#include <ShapeAnalysis_Surface.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <Standard_Version.hxx>
#include <BRepAlgoAPI_Fuse.hxx>
#include <BRepPrimAPI_MakePrism.hxx>
#include <BOPAlgo_PaveFiller.hxx>
#include <BOPAlgo_Alerts.hxx>
#include <ShapeFix_Shape.hxx>
#include <BRepCheck_Analyzer.hxx>
#include <BRepCheck_ListIteratorOfListOfStatus.hxx>
#include <BRepCheck.hxx>
#include <ShapeAnalysis_Edge.hxx>
#include <vector>
#include <thread>
void IfcGeom::util::copy_operand(const TopTools_ListOfShape & l, TopTools_ListOfShape & r) {
#if OCC_VERSION_HEX < 0x70000
r.Clear();
TopTools_ListIteratorOfListOfShape it(l);
for (; it.More(); it.Next()) {
r.Append(BRepBuilderAPI_Copy(it.Value()));
}
#else
// On OCCT 7.0 and higher BRepAlgoAPI_BuilderAlgo::SetNonDestructive(true) is
// called. Not entirely sure on the behaviour before 7.0, so overcautiously
// create copies.
r.Assign(l);
#endif
}
TopoDS_Shape IfcGeom::util::copy_operand(const TopoDS_Shape & s) {
#if OCC_VERSION_HEX < 0x70000
return BRepBuilderAPI_Copy(s);
#else
return s;
#endif
}
double IfcGeom::util::min_edge_length(const TopoDS_Shape & a) {
double min_edge_len = std::numeric_limits<double>::infinity();
TopExp_Explorer exp(a, TopAbs_EDGE);
for (; exp.More(); exp.Next()) {
const TopoDS_Edge& e = TopoDS::Edge(exp.Current());
TopoDS_Vertex v0, v1;
TopExp::Vertices(e, v0, v1);
if (!v0.IsNull() && !v1.IsNull() && v0.IsSame(v1)) {
// Don't consider a 3d-degenerate edge (for example cone apex)
// in calculating overall shape min edge length.
continue;
}
GProp_GProps prop;
BRepGProp::LinearProperties(e, prop);
double l = prop.Mass();
if (l < min_edge_len) {
min_edge_len = l;
}
}
return min_edge_len;
}
double IfcGeom::util::min_vertex_edge_distance(const TopoDS_Shape & a, double min_search, double max_search) {
double M = std::numeric_limits<double>::infinity();
TopTools_IndexedMapOfShape vertices, edges;
TopExp::MapShapes(a, TopAbs_VERTEX, vertices);
TopExp::MapShapes(a, TopAbs_EDGE, edges);
IfcGeom::impl::tree<int> tree;
// Add edges to tree
for (int i = 1; i <= edges.Extent(); ++i) {
tree.add(i, edges(i));
}
for (int j = 1; j <= vertices.Extent(); ++j) {
const TopoDS_Vertex& v = TopoDS::Vertex(vertices(j));
gp_Pnt p = BRep_Tool::Pnt(v);
Bnd_Box b;
b.Add(p);
b.Enlarge(max_search);
std::vector<int> edge_idxs = tree.select_box(b, false);
std::vector<int>::const_iterator it = edge_idxs.begin();
for (; it != edge_idxs.end(); ++it) {
const TopoDS_Edge& e = TopoDS::Edge(edges(*it));
TopoDS_Vertex v1, v2;
TopExp::Vertices(e, v1, v2);
if (v.IsSame(v1) || v.IsSame(v2)) {
continue;
}
BRepAdaptor_Curve crv(e);
Extrema_ExtPC ext(p, crv);
if (!ext.IsDone()) {
continue;
}
for (int i = 1; i <= ext.NbExt(); ++i) {
const double m = sqrt(ext.SquareDistance(i));
if (m < M && m > min_search) {
M = m;
}
}
}
}
return M;
}
bool IfcGeom::util::faces_overlap(const TopoDS_Face & f, const TopoDS_Face & g) {
points_on_planar_face_generator pgen(f);
BRep_Builder B;
gp_Pnt test;
double eps = BRep_Tool::Tolerance(f) + BRep_Tool::Tolerance(g);
BRepExtrema_DistShapeShape x;
x.LoadS1(g);
while (pgen(test)) {
TopoDS_Vertex V;
B.MakeVertex(V, test, Precision::Confusion());
x.LoadS2(V);
x.Perform();
if (x.IsDone() && x.NbSolution() == 1) {
if (x.Value() > eps) {
return false;
}
}
}
return true;
}
double IfcGeom::util::min_face_face_distance(const TopoDS_Shape & a, double max_search) {
/*
NB: This is currently only implemented for planar surfaces.
*/
double M = std::numeric_limits<double>::infinity();
TopTools_IndexedMapOfShape faces;
TopExp::MapShapes(a, TopAbs_FACE, faces);
IfcGeom::impl::tree<int> tree;
// Add faces to tree
for (int i = 1; i <= faces.Extent(); ++i) {
if (BRep_Tool::Surface(TopoDS::Face(faces(i)))->DynamicType() == STANDARD_TYPE(Geom_Plane)) {
tree.add(i, faces(i));
}
}
for (int j = 1; j <= faces.Extent(); ++j) {
const TopoDS_Face& f = TopoDS::Face(faces(j));
const Handle(Geom_Surface)& fs = BRep_Tool::Surface(f);
if (fs->DynamicType() != STANDARD_TYPE(Geom_Plane)) {
continue;
}
points_on_planar_face_generator pgen(f);
Bnd_Box b;
BRepBndLib::AddClose(f, b);
b.Enlarge(max_search);
std::vector<int> face_idxs = tree.select_box(b, false);
std::vector<int>::const_iterator it = face_idxs.begin();
for (; it != face_idxs.end(); ++it) {
if (*it == j) {
continue;
}
const TopoDS_Face& g = TopoDS::Face(faces(*it));
const Handle(Geom_Surface)& gs = BRep_Tool::Surface(g);
auto p0 = Handle(Geom_Plane)::DownCast(fs);
auto p1 = Handle(Geom_Plane)::DownCast(gs);
if (p0->Position().IsCoplanar(p1->Position(), max_search, asin(max_search))) {
pgen.reset();
BRepTopAdaptor_FClass2d cls(g, BRep_Tool::Tolerance(g));
gp_Pnt test;
while (pgen(test)) {
gp_Vec d = test.XYZ() - p1->Position().Location().XYZ();
double u = d.Dot(p1->Position().XDirection());
double v = d.Dot(p1->Position().YDirection());
// nb: TopAbs_ON is explicitly not considered to prevent matching adjacent faces
// with similar orientations.
if (cls.Perform(gp_Pnt2d(u, v)) == TopAbs_IN) {
gp_Pnt test2;
p1->D0(u, v, test2);
double w = std::abs(gp_Vec(p1->Position().Direction().XYZ()).Dot(test2.XYZ() - test.XYZ()));
if (w < M) {
M = w;
}
}
}
}
}
}
return M;
}
int IfcGeom::util::bounding_box_overlap(double p, const TopoDS_Shape & a, const TopTools_ListOfShape & b, TopTools_ListOfShape & c) {
int N = 0;
Bnd_Box A;
BRepBndLib::Add(a, A);
if (A.IsVoid()) {
return 0;
}
TopTools_ListIteratorOfListOfShape it(b);
for (; it.More(); it.Next()) {
Bnd_Box B;
BRepBndLib::Add(it.Value(), B);
if (B.IsVoid()) {
continue;
}
if (A.Distance(B) < p) {
c.Append(it.Value());
} else {
++N;
}
}
return N;
}
bool IfcGeom::util::get_edge_axis(const TopoDS_Edge & e, gp_Ax1 & ax) {
double _, __;
auto crv = BRep_Tool::Curve(e, _, __);
auto line = Handle_Geom_Line::DownCast(crv);
auto bsple = Handle_Geom_BSplineCurve::DownCast(crv);
if (line) {
ax = line->Position();
return true;
} else if (bsple) {
if (bsple->NbPoles() == 2 && bsple->Degree() == 1) {
gp_Dir V(bsple->Poles().Last().XYZ() - bsple->Poles().First().XYZ());
ax = gp_Ax1(bsple->Poles().First(), V);
return true;
}
}
return false;
}
bool IfcGeom::util::is_subset(const TopTools_IndexedMapOfShape & lhs, const TopTools_IndexedMapOfShape & rhs) {
if (rhs.Extent() < lhs.Extent()) {
return false;
}
for (int i = 1; i < lhs.Extent(); ++i) {
auto& s = lhs.FindKey(i);
if (!rhs.Contains(s)) {
return false;
}
}
return true;
}
bool IfcGeom::util::is_extrusion(const gp_Vec & v, const TopoDS_Shape & s, TopoDS_Face & base, std::pair<double, double>& interval) {
// This assumes UnifySameDomain has been processed on s, so that
// the extrusion top and bottom are a single face.
TopTools_IndexedDataMapOfShapeListOfShape mapping;
TopExp::MapShapesAndAncestors(s, TopAbs_EDGE, TopAbs_FACE, mapping);
TopExp::MapShapesAndAncestors(s, TopAbs_VERTEX, TopAbs_FACE, mapping);
TopTools_ListOfShape parallel;
TopTools_IndexedMapOfShape curved_orthogonal;
gp_Ax1 ax;
gp_Ax1 V(gp::Origin(), v);
// Segment edges in parallel to extrusion direction, and orthogonal or curved,
// where the latter two categories have to make the edges part of the base or
// top face. When neither of these categories the shape is not a extrusion
// or the extrusion direction is not orthogonal to its basis.
for (int i = 1; i < mapping.Extent(); ++i) {
auto& s = mapping.FindKey(i);
if (s.ShapeType() != TopAbs_EDGE) {
continue;
}
// @todo use a linear tolernace and the face extrimities, see #2218
const TopoDS_Edge& e = TopoDS::Edge(s);
if (!get_edge_axis(e, ax)) {
// curved
curved_orthogonal.Add(e);
} else if (ax.IsParallel(V, 1.e-7)) {
parallel.Append(e);
} else if (ax.IsNormal(V, 1.e-7)) {
// ortho
curved_orthogonal.Add(e);
} else {
return false;
}
}
// Select the two faces for which their edges are subsets
// of the ortho/curved edges
TopTools_IndexedMapOfShape ortho_faces;
for (TopExp_Explorer exp(s, TopAbs_FACE); exp.More(); exp.Next()) {
TopTools_IndexedMapOfShape face_edges;
TopExp::MapShapes(exp.Current(), TopAbs_EDGE, face_edges);
if (is_subset(face_edges, curved_orthogonal)) {
ortho_faces.Add(exp.Current());
}
}
// There should be a basis and top face
if (ortho_faces.Extent() != 2) {
return false;
}
// For the parallel edges assert that its two vertices are part
// of both the basis and the top face.
for (TopTools_ListIteratorOfListOfShape it(parallel);
it.More(); it.Next()) {
TopoDS_Vertex v01[2];
TopExp::Vertices(TopoDS::Edge(it.Value()), v01[0], v01[1]);
TopTools_IndexedMapOfShape v_ortho_faces;
int nb_ortho_faces[2] = { 0,0 };
for (int i = 0; i < 2; ++i) {
auto& faces = mapping.FindFromKey(v01[i]);
for (TopTools_ListIteratorOfListOfShape jt(faces);
jt.More(); jt.Next()) {
if (ortho_faces.Contains(jt.Value())) {
nb_ortho_faces[i] ++;
v_ortho_faces.Add(jt.Value());
}
}
}
bool sets_equal = v_ortho_faces.Size() == ortho_faces.Size() && is_subset(v_ortho_faces, ortho_faces);
if (!sets_equal) {
return false;
}
}
// Assert the base/top faces are planar and get the interval
// (dot products along axis) for which the extrusion is defined
// If necessary swap the two faces so that the basis face has
// the smallest dot product along the axis.
auto f0 = TopoDS::Face(ortho_faces.FindKey(1));
auto f1 = TopoDS::Face(ortho_faces.FindKey(2));
const Handle(Geom_Surface)& f0_s = BRep_Tool::Surface(f0);
const Handle(Geom_Surface)& f1_s = BRep_Tool::Surface(f1);
auto p0 = Handle(Geom_Plane)::DownCast(f0_s);
auto p1 = Handle(Geom_Plane)::DownCast(f1_s);
if (p0.IsNull() || p1.IsNull()) {
return false;
}
auto dot0 = p0->Location().XYZ().Dot(v.XYZ());
auto dot1 = p1->Location().XYZ().Dot(v.XYZ());
if (dot0 > dot1) {
std::swap(dot0, dot1);
std::swap(f0, f1);
}
base = f0;
interval = { dot0, dot1 };
return true;
}
int IfcGeom::util::eliminate_touching_operands(double prec, const TopoDS_Shape & a, const TopTools_ListOfShape & bs, TopTools_ListOfShape & c) {
TopTools_IndexedMapOfShape a_faces;
TopExp::MapShapes(a, TopAbs_FACE, a_faces);
// Check if any of the faces in a are non-planar, which is
// not supported by this quick check.
for (int i = 1; i <= a_faces.Extent(); ++i) {
auto surf = BRep_Tool::Surface(TopoDS::Face(a_faces(i)));
if (surf->DynamicType() != STANDARD_TYPE(Geom_Plane)) {
return 0;
}
}
TopTools_IndexedMapOfShape a_vertices;
TopExp::MapShapes(a, TopAbs_VERTEX, a_vertices);
IfcGeom::impl::tree<int> tree;
// Add faces to tree
for (int i = 1; i <= a_faces.Extent(); ++i) {
tree.add(i, a_faces(i));
}
int N = 0;
TopTools_ListIteratorOfListOfShape it(bs);
for (; it.More(); it.Next()) {
bool is_touching = false;
auto& b = it.Value();
TopTools_IndexedMapOfShape b_faces;
TopExp::MapShapes(b, TopAbs_FACE, b_faces);
// Check if any of the faces in b are non-planar, which is
// not supported by this quick check.
for (int i = 1; i <= b_faces.Extent(); ++i) {
auto surf = BRep_Tool::Surface(TopoDS::Face(b_faces(i)));
if (surf->DynamicType() != STANDARD_TYPE(Geom_Plane)) {
continue;
}
}
TopTools_IndexedMapOfShape b_vertices;
TopExp::MapShapes(b, TopAbs_VERTEX, b_vertices);
for (int k = 1; k <= b_faces.Extent(); ++k) {
const TopoDS_Face& f_b = TopoDS::Face(b_faces(k));
Bnd_Box B;
BRepBndLib::Add(f_b, B);
// Query tree using b_face bounding box
for (auto& i : tree.select_box(B, false)) {
const TopoDS_Face& f_a = TopoDS::Face(a_faces(i));
TopTools_IndexedMapOfShape f_a_vertices;
TopExp::MapShapes(f_a, TopAbs_VERTEX, f_a_vertices);
BRepGProp_Face prop_a(f_a);
BRepGProp_Face prop_b(f_b);
gp_Pnt p_a, p_b;
gp_Vec v_a, v_b;
double u0, u1, v0, v1;
prop_a.Bounds(u0, u1, v0, v1);
prop_a.Normal((u0 + u1) / 2., (u0 + u1) / 2., p_a, v_a);
prop_b.Bounds(u0, u1, v0, v1);
prop_b.Normal((u0 + u1) / 2., (u0 + u1) / 2., p_b, v_b);
bool all_vertices_behind_f_a = true;
// Check if all 'other' vertices in a are pointing
// away from the face in a, so that there is no geometry
// from a in front of the face that could participate
// in the boolean subtraction.
for (int j = 1; j <= a_vertices.Extent(); ++j) {
if (!f_a_vertices.Contains(a_vertices(j))) {
auto p = BRep_Tool::Pnt(TopoDS::Vertex(a_vertices(j)));
if ((p.XYZ() - p_a.XYZ()).Dot(v_a.XYZ()) > prec) {
all_vertices_behind_f_a = false;
break;
}
}
}
if (!all_vertices_behind_f_a) {
continue;
}
// Check if surface normals are opposite
if (v_a.IsOpposite(v_b, 1.e-5)) {
// Check if faces are co-planar
if ((p_b.XYZ() - p_a.XYZ()).Dot(v_a.XYZ()) <= prec) {
TopTools_IndexedMapOfShape f_b_vertices;
TopExp::MapShapes(f_b, TopAbs_VERTEX, f_b_vertices);
bool all_vertices_behind_f_b = true;
// Check if all 'other' vertices in b are pointing
// away from the face in a. So that a boolean subtraction
// would not alter a.
for (int j = 1; j <= b_vertices.Extent(); ++j) {
if (!f_b_vertices.Contains(b_vertices(j))) {
auto p = BRep_Tool::Pnt(TopoDS::Vertex(b_vertices(j)));
if ((p.XYZ() - p_a.XYZ()).Dot(v_a.XYZ()) < prec * 10.) {
all_vertices_behind_f_b = false;
break;
}
}
}
if (all_vertices_behind_f_b) {
is_touching = true;
break;
}
}
}
}
if (is_touching) {
break;
}
}
if (!is_touching) {
c.Append(it.Value());
} else {
++N;
}
}
return N;
}
TopoDS_Shape IfcGeom::util::unify(const TopoDS_Shape & s, double tolerance) {
tolerance = (std::min)(min_edge_length(s) / 2., tolerance);
ShapeUpgrade_UnifySameDomain usd(s);
#if OCC_VERSION_HEX >= 0x70200
usd.SetSafeInputMode(true);
#endif
#if OCC_VERSION_HEX >= 0x70100
usd.SetLinearTolerance(tolerance);
usd.SetAngularTolerance(1.e-3);
#endif
usd.Build();
return usd.Shape();
}
bool IfcGeom::util::boolean_subtraction_2d_using_builder(const TopoDS_Shape & a_input, const TopTools_ListOfShape & b_input, TopoDS_Shape & result, double eps) {
IfcGeom::impl::tree<int> edge_tree;
TopTools_ListOfShape ab_input = b_input;
ab_input.Prepend(a_input);
TopTools_ListIteratorOfListOfShape it(ab_input);
int shape_index = 0;
int edge_index = 0;
std::map<int, int> edge_index_to_shape_index;
std::vector<TopoDS_Shape> shapes;
std::vector<std::pair<size_t, TopoDS_Edge>> edges;
// First is the outer wire
std::vector<TopoDS_Wire> wires;
for (; it.More(); it.Next(), ++shape_index) {
if (it.Value().ShapeType() != TopAbs_FACE) {
return false;
}
const TopoDS_Face& f = TopoDS::Face(it.Value());
TopoDS_Wire outer_wire;
if (shape_index == 0) {
outer_wire = BRepTools::OuterWire(f);
wires.push_back(outer_wire);
}
size_t num_wires = 0;
TopoDS_Iterator it2(it.Value());
for (; it2.More(); it2.Next()) {
++num_wires;
if (outer_wire.IsNull() || !it2.Value().IsSame(outer_wire)) {
wires.push_back(TopoDS::Wire(it2.Value()));
if (shape_index == 0 && num_wires > 0) {
// An inner wire on the first operand face: reverse, because
// MakeFace expects inner boundaries to be added as bounded
// areas.
wires.back().Reverse();
}
}
}
if (num_wires > 1 && shape_index != 0) {
// The first operand can have inner wires, but the others
// can't because a inner wire would result in an additional
// outer wire for the result.
return false;
}
shapes.push_back(it.Value());
TopExp_Explorer exp(it.Value(), TopAbs_EDGE);
for (; exp.More(); exp.Next(), ++edge_index) {
edge_tree.add(edge_index, exp.Current());
edge_index_to_shape_index[edge_index] = shape_index;
edges.push_back({ shape_index, TopoDS::Edge(exp.Current()) });
}
}
shape_index = 0;
edge_index = 0;
it.Initialize(ab_input);
for (; it.More(); it.Next(), ++shape_index) {
TopExp_Explorer exp(it.Value(), TopAbs_EDGE);
for (; exp.More(); exp.Next(), ++edge_index) {
Bnd_Box b;
BRepBndLib::Add(exp.Current(), b);
b.Enlarge(eps);
for (auto& i : edge_tree.select_box(b)) {
if (i == edge_index) {
// Skip self-selection
continue;
}
if (edges[i].first == shape_index) {
// Skip edges of the same operand
continue;
}
const TopoDS_Edge& e0 = TopoDS::Edge(exp.Current());
const TopoDS_Edge& e1 = edges[i].second;
double u11, u12, u21, u22, U1, U2;
GeomAPI_ExtremaCurveCurve ecc(
BRep_Tool::Curve(e0, u11, u12),
BRep_Tool::Curve(e1, u21, u22)
);
// @todo: extend this to work in case of multiple extrema and curved segments.
const bool unbounded_intersects = (!ecc.Extrema().IsParallel() && ecc.NbExtrema() == 1 && ecc.Distance(1) < eps);
if (unbounded_intersects) {
ecc.Parameters(1, U1, U2);
if (u11 > u12) {
std::swap(u11, u12);
}
if (u21 > u22) {
std::swap(u21, u22);
}
/// @todo: tfk: probably need different thresholds on non-linear curves
u11 -= eps;
u12 += eps;
u21 -= eps;
u22 += eps;
if (u11 < U1 && U1 < u12 && u21 < U2 && U2 < u22) {
// Edge curves belonging to different operands intersect, don't process
// using builder.
Logger::Notice("Intersecting boundaries");
return false;
}
}
}
}
}
// Only inner wires are considered that are directly contained in the outer wire
// Redundant subtractions are eliminated.
std::vector<bool> redundant(wires.size(), false);
std::vector<TopoDS_Face> wire_faces;
wire_faces.reserve(wires.size());
std::vector<BRepTopAdaptor_FClass2d> wire_clss;
wire_clss.reserve(wires.size());
std::vector<std::unique_ptr<ShapeAnalysis_Surface>> sass;
sass.reserve(wires.size());
for (auto& w : wires) {
wire_faces.push_back(BRepBuilderAPI_MakeFace(w).Face());
wire_clss.emplace_back(wire_faces.back(), eps);
sass.push_back(std::make_unique<ShapeAnalysis_Surface>(BRep_Tool::Surface(wire_faces.back())));
}
// First check for containment in outer wire
for (auto it = ++wires.begin(); it != wires.end(); ++it) {
// Considering a single vertex is sufficient because we have already
// guaranteed that the edges of different operands do not cross.
TopoDS_Iterator it_ed(*it);
auto& ed = it_ed.Value();
TopoDS_Iterator it_v(ed);
auto& v = TopoDS::Vertex(it_v.Value());
auto pnt = BRep_Tool::Pnt(v);
auto p2d = sass[0]->ValueOfUV(pnt, eps);
if (wire_clss[0].Perform(p2d) != TopAbs_IN) {
// A wire is not contained in the outer wire, it's a subtraction without
// any effect and marked as redundant. Feeding it to the builder algo
// will likely cause problems.
redundant[std::distance(wires.begin(), it)] = true;
Logger::Notice("Subtraction operand outside of outer bound");
}
}
// Now build a tree to find inner wires contained in other inner wires
// NB first wire is *not* in this tree
IfcGeom::impl::tree<int> wire_tree;
for (size_t wire_index = 1; wire_index < wires.size(); ++wire_index) {
wire_tree.add(wire_index, wires[wire_index]);
}
for (size_t wire_index = 1; wire_index < wires.size(); ++wire_index) {
Bnd_Box b;
BRepBndLib::Add(wires[wire_index], b);
b.Enlarge(eps);
// We're only selecting operands completely within b because we
// have already guaranteed they do not intersect. So they are
// either fully in or out. Selecting with complete_within=true
// will filter out some unnecessary cases. It also means we need
// that due this asymmetry we need to process all pairs of wire
// indices and not just the pairs where the first element is less
// than the second element.
for (auto& other_index : wire_tree.select_box(b, true)) {
// other_index is fully contained in wire_index
if (wire_index == other_index) {
continue;
}
TopoDS_Iterator it_ed(wires[other_index]);
auto& ed = it_ed.Value();
TopoDS_Iterator it_v(ed);
auto& v = TopoDS::Vertex(it_v.Value());
auto pnt = BRep_Tool::Pnt(v);
auto p2d = sass[wire_index]->ValueOfUV(pnt, eps);
if (wire_clss[wire_index].Perform(p2d) == TopAbs_IN) {
// A wire is contained within another operand
redundant[other_index] = true;
Logger::Notice("Subtraction operand contained in other");
}
}
}
BRepBuilderAPI_MakeFace mf(wire_faces[0]);
for (size_t wire_index = 1; wire_index < wires.size(); ++wire_index) {
if (!redundant[wire_index]) {
mf.Add(TopoDS::Wire(wires[wire_index].Reversed()));
}
}
result = mf.Face();
return true;
}
void IfcGeom::util::points_on_planar_face_generator::reset() {
i = j = (int)inset_;
}
bool IfcGeom::util::points_on_planar_face_generator::operator()(gp_Pnt& p) {
while (j < N) {
double u = u0 + (u1 - u0) * i / N;
double v = v0 + (v1 - v0) * j / N;
i++;
if (i == N) {
i = 0;
j++;
}
// Specifically does not consider ON
if (cls_.Perform(gp_Pnt2d(u, v)) == TopAbs_IN) {
plane_->D0(u, v, p);
return true;
}
}
return false;
}
bool IfcGeom::util::boolean_operation(const boolean_settings& settings, const TopoDS_Shape& a_input, const TopTools_ListOfShape& b_input, BOPAlgo_Operation op, TopoDS_Shape& result, double fuzziness) {
using namespace std::string_literals;
const bool do_unify = true;
const bool do_subtraction_eliminate_disjoint_bbox = true;
const bool do_subtraction_eliminate_touching = true;
const bool do_attempt_2d_boolean = settings.attempt_2d;
const bool debug = settings.debug;
std::string debug_identifier;
if (debug) {
static my_thread_local size_t operation_counter_ = 0;
std::stringstream ss;
ss << "bool-" << std::this_thread::get_id() << "-" << (operation_counter_++);
debug_identifier = ss.str();
Logger::Notice("Boolean debug identifier: " + debug_identifier);
}
if (fuzziness < 0.) {
fuzziness = settings.precision / 10.;
}
// @todo, it does seem a bit odd, we first triangulate non-planar faces
// to later unify them again. Can we make this a bit more intelligent?
TopoDS_Shape a;
TopTools_ListOfShape b;
if (do_unify) {
PERF("boolean operation: unifying operands");
a = unify(a_input, fuzziness * 1000.);
Logger::Message(
Logger::LOG_DEBUG,
"Simplified operand A from "s +
std::to_string(count(a_input, TopAbs_FACE)) +
" to "s +
std::to_string(count(a, TopAbs_FACE))
);
{
TopTools_ListIteratorOfListOfShape it(b_input);
for (; it.More(); it.Next()) {
b.Append(unify(it.Value(), fuzziness));
Logger::Message(
Logger::LOG_DEBUG,
"Simplified operand B from "s +
std::to_string(count(it.Value(), TopAbs_FACE)) +
" to "s +
std::to_string(count(b.Last(), TopAbs_FACE))
);
}
}
} else {
a = a_input;
b = b_input;
}
bool success = false;
BRepAlgoAPI_BooleanOperation* builder;
TopTools_ListOfShape b_tmp;
if (op == BOPAlgo_CUT) {
builder = new BRepAlgoAPI_Cut();
if (do_subtraction_eliminate_disjoint_bbox) {
PERF("boolean subtraction: eliminate disjoint bbox");
auto N = bounding_box_overlap(fuzziness, a, b, b_tmp);
if (N) {
Logger::Notice("Eliminated " + std::to_string(N) + " disjoint operands");
std::swap(b, b_tmp);
}
}
if (do_subtraction_eliminate_touching) {
PERF("boolean subtraction: eliminate touching");
b_tmp.Clear();
auto N = eliminate_touching_operands(fuzziness, a, b, b_tmp);
if (N) {
Logger::Notice("Eliminated " + std::to_string(N) + " touching operands");
std::swap(b, b_tmp);
}
}
} else if (op == BOPAlgo_COMMON) {
builder = new BRepAlgoAPI_Common();
} else if (op == BOPAlgo_FUSE) {
builder = new BRepAlgoAPI_Fuse();
} else {
return false;
}
if (b.Extent() == 0) {
result = a;
return true;
}
if (Logger::LOG_NOTICE >= Logger::Verbosity()) {
PERF("preliminary manifoldness check");
Logger::Notice("Operand A is " + (is_manifold(a) ? ""s : "non-"s) + "manifold");
TopTools_ListIteratorOfListOfShape it(b);
for (int i = 0; it.More(); it.Next(), ++i) {
Logger::Notice("Operand B " + std::to_string(i) + " is " + (is_manifold(it.Value()) ? ""s : "non-"s) + "manifold");
}
}
// Find a sensible value for the fuzziness, based on precision
// and limited by edge lengths and vertex-edge distances.
double min_length_orig;
{
PERF("boolean operation: min edge length");
min_length_orig = min_edge_length(a);
TopTools_ListIteratorOfListOfShape it(b);
for (; it.More(); it.Next()) {
double d = min_edge_length(it.Value());
if (d < min_length_orig) {
min_length_orig = d;
}
}
}
{
PERF("boolean operation: min vertex-edge dist");
double d = min_vertex_edge_distance(a, settings.precision, min_length_orig);
if (d < min_length_orig) {
min_length_orig = d;
}
TopTools_ListIteratorOfListOfShape it(b);
for (; it.More(); it.Next()) {
d = min_vertex_edge_distance(it.Value(), settings.precision, min_length_orig);
if (d < min_length_orig) {
min_length_orig = d;
}
}
}
const double fuzz = (std::min)(min_length_orig / 3., fuzziness);
Logger::Notice("Used fuzziness: " + std::to_string(fuzz));
TopTools_ListOfShape s1s;
s1s.Append(copy_operand(a));
if (debug) {
TopTools_ListOfShape* lists[2] = { &s1s, &b };
static std::string operand_names[2] = { "a", "b" };
for (int i = 0; i < 2; ++i) {
TopTools_ListIteratorOfListOfShape it(*lists[i]);
for (int j = 0; it.More(); it.Next(), ++j) {
std::string fn = debug_identifier + "-" + operand_names[i] + "-" + std::to_string(j) + ".brep";
BRepTools::Write(it.Value(), fn.c_str());
}
}
}
if (op == BOPAlgo_CUT) {
TopoDS_Face a_face;
std::pair<double, double> a_interval;
TopTools_ListOfShape b_faces, b_remainder_3d;
bool is_extrusion_a = false;
if (do_attempt_2d_boolean) {
PERF("boolean subtraction: extrusion check");
is_extrusion_a = is_extrusion(gp::DY(), a, a_face, a_interval);
}
if (is_extrusion_a) {
Logger::Notice("Operand A 1/1 is an extrusion");
TopTools_ListIteratorOfListOfShape it(b);
for (int nb = 1; it.More(); it.Next(), ++nb) {
bool process_2d = false;
TopoDS_Face b_face;
std::pair<double, double> b_interval;
bool is_extrusion_b;
{
PERF("boolean subtraction: extrusion check");
is_extrusion_b = is_extrusion(gp::DY(), it.Value(), b_face, b_interval);
}
if (is_extrusion_b) {
Logger::Notice("Operand B " + std::to_string(nb) + "/" + std::to_string(b.Extent()) + " is an extrusion");
if (b_interval.first < a_interval.first + fuzz && b_interval.second > a_interval.second - fuzz) {
Logger::Notice("Operand B creates a through hole");
// Align b with a operand
gp_Trsf trsf;
trsf.SetTranslation(gp_Vec(gp::DY()) * (a_interval.first - b_interval.first));
b_faces.Append(b_face.Moved(trsf));
process_2d = true;
}
}
if (!process_2d) {
b_remainder_3d.Append(it.Value());
}
}
if (b_faces.Extent()) {
TopoDS_Shape face_result;
bool boolean_op_2d_success;
{
PERF("boolean operation: 2d builder");
// First try using face builder
boolean_op_2d_success = boolean_subtraction_2d_using_builder(a_face, b_faces, face_result, fuzziness);
}
if (!boolean_op_2d_success) {
PERF("boolean operation: 2d");
// Retry using generic 2d using boolean algo on faces
boolean_op_2d_success = boolean_operation(settings, a_face, b_faces, op, face_result, fuzziness);
}
if (boolean_op_2d_success) {
PERF("boolean operation: 2d to 3d");
BRepPrimAPI_MakePrism mp(face_result, gp_Vec(gp::DY()) * (a_interval.second - a_interval.first));
if (mp.IsDone()) {
if (b_remainder_3d.Extent()) {
Logger::Notice(std::to_string(b_remainder_3d.Extent()) + " operands remaining to process in 3D");
b = b_remainder_3d;
s1s.Clear();
s1s.Append(mp.Shape());
} else {
Logger::Notice("Processed fully in 2D");
result = mp.Shape();
return true;
}
} else {
Logger::Notice("Failed to extrude 2D boolean result. Retrying in 3D.");
}
} else {
Logger::Notice("Failed to perform 2D boolean operation. Retrying in 3D.");
}
} else {
Logger::Notice("No second operands can be processed as 2D inner bounds. Retrying in 3D.");
}
}
}
#if OCC_VERSION_HEX >= 0x70000
builder->SetNonDestructive(true);
#endif
builder->SetFuzzyValue(fuzz);
builder->SetArguments(s1s);
copy_operand(b, b_tmp);
std::swap(b, b_tmp);
builder->SetTools(b);
{
PERF("boolean operation: build");
builder->Build();
}
if (builder->IsDone()) {
if (builder->DSFiller()->HasWarning(STANDARD_TYPE(BOPAlgo_AlertAcquiredSelfIntersection))) {
Logger::Notice("Builder reports self-intersection in output");
success = false;
} else {
TopoDS_Shape r = *builder;
{
PERF("boolean operation: shape healing");
ShapeFix_Shape fix(r);
try {
fix.SetMaxTolerance(fuzz);
fix.Perform();
r = fix.Shape();
} catch (...) {
Logger::Error("Shape healing failed on boolean result");
}
}
{
PERF("boolean operation: shape analysis");
BRepCheck_Analyzer ana(r);
success = ana.IsValid() != 0;
if (!success) {
Logger::Notice("Boolean operation yields invalid result");
std::stringstream str;
bool any_emitted = false;
std::function<void(const TopoDS_Shape&)> dump;
dump = [&ana, &str, &dump, &any_emitted](const TopoDS_Shape& s) {
if (!ana.Result(s).IsNull()) {
BRepCheck_ListIteratorOfListOfStatus itl;
itl.Initialize(ana.Result(s)->Status());
for (; itl.More(); itl.Next()) {
if (itl.Value() != BRepCheck_NoError) {
if (any_emitted) {
str << ", ";
}
BRepCheck::Print(itl.Value(), str);
str.seekp(str.tellp() - (std::streamoff)1);
str << " on ";
TopAbs::Print(s.ShapeType(), str);
any_emitted = true;
}
}
}
for (TopoDS_Iterator it(s); it.More(); it.Next()) {
dump(it.Value());
}
};
dump(r);
Logger::Notice(str.str());
}
}
if (success) {
{
PERF("boolean operation: manifoldness check");
success = !is_manifold(a) || is_manifold(r);
}
if (!success) {
PERF("boolean operation: manifoldness check excemption");
// An excemption for the requirement to be manifold: When the cut operands have overlapping edge belonging to faces that do not overlap.
bool operands_nonmanifold = false;
if (op == BOPAlgo_CUT) {
TopTools_IndexedMapOfShape edges;
TopTools_IndexedDataMapOfShapeListOfShape map;
for (TopTools_ListIteratorOfListOfShape it2(b); it2.More(); it2.Next()) {
auto& bb = it2.Value();
TopExp::MapShapes(bb, TopAbs_EDGE, edges);
TopExp::MapShapesAndAncestors(bb, TopAbs_EDGE, TopAbs_FACE, map);
}
IfcGeom::impl::tree<int> tree;
for (int i = 1; i <= edges.Extent(); ++i) {
tree.add(i, edges.FindKey(i));
}
for (int i = 1; i <= edges.Extent(); ++i) {
const TopoDS_Edge& ei = TopoDS::Edge(edges.FindKey(i));
Bnd_Box bb;
BRepBndLib::Add(ei, bb);
bb.Enlarge(fuzziness);
auto ii = tree.select_box(bb, false);
for (int j : ii) {
if (j != i) {
const TopoDS_Edge& ej = TopoDS::Edge(edges.FindKey(j));
ShapeAnalysis_Edge sae;
double f = fuzziness;
bool edges_overlapping = sae.CheckOverlapping(ei, ej, f, 0.) ||
sae.CheckOverlapping(ej, ei, f, 0.);
if (edges_overlapping) {
auto faces_i = map.FindFromKey(edges.FindKey(i));
auto faces_j = map.FindFromKey(edges.FindKey(j));
bool overlap = false;
for (TopTools_ListIteratorOfListOfShape it4(faces_i); it4.More(); it4.Next()) {
auto& fi = it4.Value();
for (TopTools_ListIteratorOfListOfShape it2(faces_j); it2.More(); it2.Next()) {
auto& fj = it2.Value();
if (faces_overlap(TopoDS::Face(fi), TopoDS::Face(fj))) {
overlap = true;
}
}
if (overlap) {
break;
}
}
operands_nonmanifold = !overlap;
break;
}
}
}
if (operands_nonmanifold) {
break;
}
}
}
success = operands_nonmanifold;
}
if (success) {
bool all_faces_included_in_result = true;
bool has_open_shells = false;
if (op == BOPAlgo_CUT) {
PERF("boolean operation: open shell face adition check");
for (TopExp_Explorer exp(a, TopAbs_SHELL); exp.More(); exp.Next()) {
if (!exp.Current().Closed()) {
// This 'face addition check' is only done when the first operand
// contains open shells (which was initially the aim of this check
// see #1472).
// Later in #1914 we found that the logic to apply openings in groups
// of similar edge lengths can create a situation of inner voids, which
// trigger a false positive in this check. This could have also been
// solved below by checking whether the opening(s) are included as a
// unmodified (interior) shell within a solid of multiple shells.
// Checking for open shells in first operand was quicker and more
// straightforward. The question still is whether in cases like #1472
// we need to first try the boolean union as solid/solid interference
// to trigger this case or whether we can immediately proceed to a face/
// solid operation.
has_open_shells = true;
break;
}
}
if (has_open_shells) {
TopTools_IndexedMapOfShape faces;
TopExp::MapShapes(r, TopAbs_FACE, faces);
for (TopExp_Explorer exp(a, TopAbs_FACE); exp.More(); exp.Next()) {
auto& f = TopoDS::Face(exp.Current());
if (!faces.Contains(f)) {
all_faces_included_in_result = false;
break;
}
}
} else {
all_faces_included_in_result = false;
}
}
int result_n_faces = count(r, TopAbs_FACE);
int first_op_n_faces = count(a, TopAbs_FACE);
if (op == BOPAlgo_CUT && has_open_shells && all_faces_included_in_result && result_n_faces > first_op_n_faces) {
success = false;
Logger::Notice("Boolean result discarded because subtractions results in only the addition of faces");
} else {
// when there are edges or vertex-edge distances close to the used fuzziness, the
// output is not trusted and the operation is attempted with a higher fuzziness.
int reason = 0;
double v;
{
PERF("boolean operation: result min edge length check");
if ((v = min_edge_length(r)) < fuzziness * 3.) {
reason = 0;
success = false;
goto skip_further_checks;
}
}
{
PERF("boolean operation: result min vertex-edge dist check");
if ((v = min_vertex_edge_distance(r, settings.precision, fuzziness * 3.)) < fuzziness * 3.) {
reason = 1;
success = false;
goto skip_further_checks;
}
}
{
PERF("boolean operation: result min face-face dist check");
if ((v = min_face_face_distance(r, 1.e-4)) < 1.e-4) {
// #2095 Check if this distance wasn't already realized in the input first operand.
if (v < min_face_face_distance(a, 1.e-4)) {
reason = 2;
success = false;
}
}
}
skip_further_checks:
if (!success) {
static const char* const reason_strings[] = { "edge length", "vertex-edge", "face-face" };
std::stringstream str;
str << "Boolean operation result failing " << reason_strings[reason] << " interference check, with fuzziness " << fuzziness << " with length " << v;
Logger::Notice(str.str());
}
}
if (success) {
result = r;
}
} else {
Logger::Notice("Boolean operation yields non-manifold result");
}
}
}
} else {
std::stringstream str;
#if OCC_VERSION_HEX >= 0x70200
if (builder->HasError(STANDARD_TYPE(BOPAlgo_AlertBOPNotAllowed))) {
Logger::Error("Invalid operands. Using first operand");
result = a;
success = true;
}
#endif
#if OCC_VERSION_HEX >= 0x70000
builder->DumpErrors(str);
#else
str << "Error code: " << builder->ErrorStatus();
#endif
std::string str_str = str.str();
if (str_str.size()) {
Logger::Notice(str_str);
}
}
delete builder;
if (!success) {
const double new_fuzziness = fuzziness * 10.;
if (new_fuzziness - 1e-15 <= settings.precision * 10000. && new_fuzziness < min_length_orig) {
return boolean_operation(settings, a, b, op, result, new_fuzziness);
} else {
Logger::Notice("No longer attempting boolean operation with higher fuzziness");
}
}
return success && !result.IsNull();
}
bool IfcGeom::util::boolean_operation(const boolean_settings& settings, const TopoDS_Shape& a, const TopoDS_Shape& b, BOPAlgo_Operation op, TopoDS_Shape& result, double fuzziness) {
TopTools_ListOfShape bs;
bs.Append(b);
return boolean_operation(settings, a, bs, op, result, fuzziness);
}
const TopoDS_Shape& IfcGeom::util::ensure_fit_for_subtraction(const TopoDS_Shape& shape, TopoDS_Shape& solid, double tol) {
const bool is_comp = is_compound(shape);
if (!is_comp) {
return solid = shape;
}
if (!create_solid_from_compound(shape, solid, tol)) {
return solid = shape;
}
return solid;
}
Reference in New Issue View Git Blame Copy Permalink