Keep moving schema agnostic code out of kernel

2026-09-27 02:31:09 +00:00 · 2022-09-12 10:34:04 +02:00
parent 4cd0932d7d
commit a63d6475b2
8 changed files with 938 additions and 1026 deletions
@@ -38,6 +38,7 @@
 #include "../ifcgeom/IfcGeom.h"
 #include "../ifcgeom_schema_agnostic/face_definition.h"
 #include "../ifcgeom_schema_agnostic/wire_utils.h"
 #define Kernel MAKE_TYPE_NAME(Kernel)
@@ -156,7 +157,7 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcFace* l, TopoDS_Shape& result)
 				}
 			} else {
 				gp_Pln pln;
-				if (approximate_plane_through_wire(wire, pln)) {
+				if (util::approximate_plane_through_wire(wire, pln, getValue(GV_PRECISION))) {
 					fd.surface() = new Geom_Plane(pln);
 				}
 			}
@@ -185,11 +186,17 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcFace* l, TopoDS_Shape& result)
 		if (fd.all_outer()) {
 			for (const auto& w : fd.wires()) {
 				TopTools_ListOfShape fl;
-				triangulate_wire({ w }, fl);
+				auto r = util::triangulate_wire({ w }, fl);
 				face_list.Append(fl);
 				if (faceset_helper_ && r == util::TRIANGULATE_WIRE_NON_MANIFOLD) {
 					faceset_helper_->non_manifold() = true;
 				}
 			}
 		} else {
-			triangulate_wire(fd.wires(), face_list);			
+			auto r = util::triangulate_wire(fd.wires(), face_list);
 			if (faceset_helper_ && r == util::TRIANGULATE_WIRE_NON_MANIFOLD) {
 				faceset_helper_->non_manifold() = true;
 			}
 		}
 	} else if (!fd.all_outer()) {
 		BRepBuilderAPI_MakeFace mf(fd.surface(), fd.outer_wire());
@@ -40,10 +40,6 @@
 #include <gp_Ax3.hxx>
 #include <gp_Pln.hxx>
 #include <boost/range/irange.hpp>
 #include <boost/range/algorithm_ext/push_back.hpp>
 #include <Geom_Line.hxx>
 #include <Geom_Circle.hxx>
@@ -143,6 +139,7 @@
 #include "../ifcgeom_schema_agnostic/IfcGeomTree.h"
 #include "../ifcgeom_schema_agnostic/boolean_utils.h"
 #include "../ifcgeom_schema_agnostic/wire_utils.h"
 #include <memory>
 #include <thread>
@@ -217,10 +214,6 @@ void MAKE_INIT_FN(KernelImplementation_)(IfcGeom::impl::KernelFactoryImplementat
 #define Kernel MAKE_TYPE_NAME(Kernel)
 namespace {
 }
 void IfcGeom::Kernel::set_offset(const std::array<double, 3> &p_offset) {
    offset = gp_Vec(p_offset[0], p_offset[1], p_offset[2]);
@@ -897,13 +890,24 @@ bool IfcGeom::Kernel::convert_openings_fast(const IfcSchema::IfcProduct* entity,
 }
 #endif
 double IfcGeom::Kernel::get_wire_intersection_tolerance(const TopoDS_Wire& wire) const {
 	return getValue(GV_NO_WIRE_INTERSECTION_TOLERANCE) > 0.
 		? 0
 		: faceset_helper_
 			// eps is added to both ends of the parametric domain, so 3. is chosen to be on the safe side here.
 			? (faceset_helper_->epsilon() / 3.)
 			// @todo re-evaluate 2. here for the reasons above:
 			: (std::min)(util::min_edge_length(wire) / 2., getValue(GV_PRECISION) * 10.);
 }
 bool IfcGeom::Kernel::convert_wire_to_face(const TopoDS_Wire& w, TopoDS_Face& face) {
 	TopoDS_Wire wire = w;
 	TopTools_ListOfShape results;
-	if (wire_intersections(wire, results)) {
+
 	if (getValue(GV_NO_WIRE_INTERSECTION_CHECK) == 0. && util::wire_intersections(wire, results, get_wire_intersection_tolerance(wire), getValue(GV_PRECISION))) {
 		Logger::Warning("Self-intersections with " + boost::lexical_cast<std::string>(results.Extent()) + " cycles detected");
-		select_largest(results, wire);
+		util::select_largest(results, wire);
 	}
 	bool is_2d = true;
@@ -958,7 +962,7 @@ bool IfcGeom::Kernel::convert_wire_to_faces(const TopoDS_Wire& w, TopoDS_Compoun
 	}
 	TopTools_ListOfShape results;
-	if (wire_intersections(w, results)) {
+	if (getValue(GV_NO_WIRE_INTERSECTION_CHECK) == 0. && util::wire_intersections(w, results, get_wire_intersection_tolerance(w), getValue(GV_PRECISION))) {
 		Logger::Warning("Self-intersections with " + boost::lexical_cast<std::string>(results.Extent()) + " cycles detected");
 	} else {
 		results.Clear();
@@ -3286,295 +3290,6 @@ bool IfcGeom::Kernel::is_identity_transform(IfcUtil::IfcBaseInterface* l) {
 	}
 }
 bool IfcGeom::Kernel::approximate_plane_through_wire(const TopoDS_Wire& wire, gp_Pln& plane, double eps) {
 	// Newell's Method is used for the normal calculation
 	// as a simple edge cross product can give opposite results
 	// for a concave face boundary.
 	// Reference: Graphics Gems III p. 231
 	const double eps_ = eps < 1. ? getValue(GV_PRECISION) : eps;
 	const double eps2 = eps_ * eps_;
 	double x = 0, y = 0, z = 0;
 	gp_Pnt current, previous, first;
 	gp_XYZ center;
 	int n = 0;
 	BRepTools_WireExplorer exp(wire);
 	for (;; exp.Next()) {
 		const bool has_more = exp.More() != 0;
 		if (has_more) {
 			const TopoDS_Vertex& v = exp.CurrentVertex();
 			current = BRep_Tool::Pnt(v);
 			center += current.XYZ();
 		} else {
 			current = first;
 		}
 		if (n) {
 			const double& xn = previous.X();
 			const double& yn = previous.Y();
 			const double& zn = previous.Z();
 			const double& xn1 = current.X();
 			const double& yn1 = current.Y();
 			const double& zn1 = current.Z();
 			x += (yn - yn1)*(zn + zn1);
 			y += (xn + xn1)*(zn - zn1);
 			z += (xn - xn1)*(yn + yn1);
 		} else {
 			first = current;
 		}
 		if (!has_more) {
 			break;
 		}
 		previous = current;
 		++n;
 	}
 	if (n < 3) {
 		return false;
 	}
 	plane = gp_Pln(center / n, gp_Dir(x, y, z));
 	exp.Init(wire);
 	for (; exp.More(); exp.Next()) {
 		const TopoDS_Vertex& v = exp.CurrentVertex();
 		current = BRep_Tool::Pnt(v);
 		if (plane.SquareDistance(current) > eps2) {
 			return false;
 		}
 	}
 	return true;
 }
 bool IfcGeom::Kernel::flatten_wire(TopoDS_Wire& wire) {
 	gp_Pln pln;
 	if (!approximate_plane_through_wire(wire, pln)) {
 		return false;
 	}
 	TopoDS_Face face = BRepBuilderAPI_MakeFace(pln).Face();
 	BRepAlgo_NormalProjection proj(face);
 	proj.Add(wire);
 	proj.Build();
 	if (!proj.IsDone()) {
 		return false;
 	}
 	TopTools_ListOfShape list;
 	proj.BuildWire(list);
 	if (list.Extent() != 1) {
 		return false;
 	}
 	wire = TopoDS::Wire(list.First());
 	return true;
 }
 bool IfcGeom::Kernel::triangulate_wire(const std::vector<TopoDS_Wire>& wires, TopTools_ListOfShape& faces) {
 	// This is a bit of a precarious approach, but seems to work for the
 	// versions of OCCT tested for. OCCT has a Delaunay triangulation function
 	// BRepMesh_Delaun, but it is notoriously hard to interpret the results
 	// (due to the Bowyer-Watson super triangle perhaps?). Therefore
 	// alternatively we use the regular OCCT incremental mesher on a new face
 	// created from the UV coordinates of the original wire. Pray to our gods
 	// that the vertex coordinates are unaffected by the meshing algorithm and
 	// map them back to 3d coordinates when iterating over the mesh triangles.
 	// In addition, to maintain a manifold shell, we need to make sure that
 	// every edge from the input wire is used exactly once in the list of
 	// resulting faces. And that other internal edges are used twice.
 	typedef std::pair<double, double> uv_node;
 	gp_Pln pln;
 	if (!approximate_plane_through_wire(wires.front(), pln, std::numeric_limits<double>::infinity())) {
 		return false;
 	}
 	const gp_XYZ& udir = pln.Position().XDirection().XYZ();
 	const gp_XYZ& vdir = pln.Position().YDirection().XYZ();
 	const gp_XYZ& pnt = pln.Position().Location().XYZ();
 	std::map<uv_node, TopoDS_Vertex> mapping;
 	std::map<std::pair<uv_node, uv_node>, TopoDS_Edge> existing_edges, new_edges;
 	std::unique_ptr<BRepBuilderAPI_MakeFace> mf;
 	for (auto it = wires.begin(); it != wires.end(); ++it) {
 		const TopoDS_Wire& wire = *it;
 		BRepTools_WireExplorer exp(wire);
 		BRepBuilderAPI_MakePolygon mp;
 		// Add UV coordinates to a newly created polygon
 		for (; exp.More(); exp.Next()) {
 			// Project onto plane
 			const TopoDS_Vertex& V = exp.CurrentVertex();
 			gp_Pnt p = BRep_Tool::Pnt(V);
 			double u = (p.XYZ() - pnt).Dot(udir);
 			double v = (p.XYZ() - pnt).Dot(vdir);
 			mp.Add(gp_Pnt(u, v, 0.));
 			mapping.insert(std::make_pair(std::make_pair(u, v), V));
 			// Store existing edges in a map so that triangles can
 			// actually reference the preexisting edges.
 			const TopoDS_Edge& e = exp.Current();
 			TopoDS_Vertex V0, V1;
 			TopExp::Vertices(e, V0, V1, true);
 			gp_Pnt p0 = BRep_Tool::Pnt(V0);
 			gp_Pnt p1 = BRep_Tool::Pnt(V1);
 			double u0 = (p0.XYZ() - pnt).Dot(udir);
 			double v0 = (p0.XYZ() - pnt).Dot(vdir);
 			double u1 = (p1.XYZ() - pnt).Dot(udir);
 			double v1 = (p1.XYZ() - pnt).Dot(vdir);
 			uv_node uv0 = std::make_pair(u0, v0);
 			uv_node uv1 = std::make_pair(u1, v1);
 			existing_edges.insert(std::make_pair(std::make_pair(uv0, uv1), e));
 			existing_edges.insert(std::make_pair(std::make_pair(uv1, uv0), TopoDS::Edge(e.Reversed())));
 		}
 		// Not closed by default
 		mp.Close();
 		if (mf) {
 			if (it - 1 == wires.begin()) {
 				// @todo is this necessary?
 				TopoDS_Face f = mf->Face();
 				mf->Init(f);
 			}
 			mf->Add(mp.Wire());
 		} else {
 			mf.reset(new BRepBuilderAPI_MakeFace(mp.Wire()));
 		}
 	}
 	const TopoDS_Face& face = mf->Face();
 	// Create a triangular mesh from the face
 	BRepMesh_IncrementalMesh(face, Precision::Confusion());
 	int n123[3];
 	TopLoc_Location loc;
 	Handle_Poly_Triangulation tri = BRep_Tool::Triangulation(face, loc);
 	if (!tri.IsNull()) {
 		const Poly_Array1OfTriangle& triangles = tri->Triangles();
 		for (int i = 1; i <= triangles.Length(); ++i) {
 			if (face.Orientation() == TopAbs_REVERSED)
 				triangles(i).Get(n123[2], n123[1], n123[0]);
 			else triangles(i).Get(n123[0], n123[1], n123[2]);
 			// Create polygons from the mesh vertices
 			BRepBuilderAPI_MakeWire mp2;
 			for (int j = 0; j < 3; ++j) {
 				uv_node uvnodes[2];
 				TopoDS_Vertex vs[2];
 				for (int k = 0; k < 2; ++k) {
 					const gp_Pnt& uv = tri->Node(n123[(j + k) % 3]);
 					uvnodes[k] = std::make_pair(uv.X(), uv.Y());
 					auto it = mapping.find(uvnodes[k]);
 					if (it == mapping.end()) {
 						Logger::Error("Internal error: unable to unproject uv-mesh");
 						return false;
 					}
 					vs[k] = it->second;
 				}
 				auto it = existing_edges.find(std::make_pair(uvnodes[0], uvnodes[1]));
 				if (it != existing_edges.end()) {
 					// This is a boundary edge, reuse existing edge from wire
 					mp2.Add(it->second);
 				} else {
 					auto jt = new_edges.find(std::make_pair(uvnodes[0], uvnodes[1]));
 					if (jt != new_edges.end()) {
 						// We have already added the reverse as part of another
 						// triangle, reuse this edge.
 						mp2.Add(TopoDS::Edge(jt->second));
 					} else {
 						// This is a new internal edge. Register the reverse
 						// for reuse later. We need to be sure to reuse vertices
 						// for the edge construction because otherwise the wire
 						// builder will use geometrical proximity for vertex
 						// connections in which case the edge will be copied
 						// and no longer partner with other edges from the shell.
 						TopoDS_Edge ne = BRepBuilderAPI_MakeEdge(vs[0], vs[1]);
 						mp2.Add(ne);
 						// Store the reverse to be picked up later.
 						new_edges.insert(std::make_pair(std::make_pair(uvnodes[1], uvnodes[0]), TopoDS::Edge(ne.Reversed())));
 					}
 				}
 			}
 			BRepBuilderAPI_MakeFace mft(mp2.Wire());
 			if (mft.IsDone()) {
 				TopoDS_Face triangle_face = mft.Face();
 				TopoDS_Iterator jt(triangle_face, false);
 				for (; jt.More(); jt.Next()) {
 					const TopoDS_Wire& w = TopoDS::Wire(jt.Value());
 					if (w.Orientation() != wires.front().Orientation()) {
 						triangle_face.Reverse();
 					}
 				}
 				faces.Append(triangle_face);
 			} else {
 				Logger::Error("Internal error: missing face");
 				return false;
 			}
 		}
 	}
 	TopTools_IndexedDataMapOfShapeListOfShape mape, mapn;
 	for (auto& wire : wires) {
 		TopExp::MapShapesAndAncestors(wire, TopAbs_EDGE, TopAbs_WIRE, mape);
 	}
 	TopTools_ListIteratorOfListOfShape it(faces);
 	for (; it.More(); it.Next()) {
 		TopExp::MapShapesAndAncestors(it.Value(), TopAbs_EDGE, TopAbs_WIRE, mapn);
 	}
 	// Validation
 	for (int i = 1; i <= mape.Extent(); ++i) {
 #if OCC_VERSION_HEX >= 0x70000
 		TopTools_ListOfShape val;
 		if (!mapn.FindFromKey(mape.FindKey(i), val)) {
 #else
 		bool contains = false;
 		try {
 			TopTools_ListOfShape val = mapn.FindFromKey(mape.FindKey(i));
 			contains = true;
 		} catch (Standard_NoSuchObject&) {}
 		if (!contains) {
 #endif
 			// All existing edges need to exist in the new faces
 			Logger::Error("Internal error, missing edge from triangulation");
 			if (faceset_helper_ != nullptr) {
 				faceset_helper_->non_manifold() = true;
 			}
 		}
 	}
 	for (int i = 1; i <= mapn.Extent(); ++i) {
 		const TopoDS_Shape& v = mapn.FindKey(i);
 		int n = mapn.FindFromIndex(i).Extent();
 		// Existing edges are boundaries with use 1
 		// New edges are internal with use 2
 		if (n != (mape.Contains(v) ? 1 : 2)) {
 			Logger::Error("Internal error, non-manifold result from triangulation");
 			if (faceset_helper_ != nullptr) {
 				faceset_helper_->non_manifold() = true;
 			}
 		}
 	}
 	return true;
 }
 TopoDS_Shape IfcGeom::Kernel::apply_transformation(const TopoDS_Shape& s, const gp_Trsf& t) {
 	if (t.Form() == gp_Identity) {
 		return s;
@@ -3597,260 +3312,6 @@ TopoDS_Shape IfcGeom::Kernel::apply_transformation(const TopoDS_Shape& s, const
 	}
 }
 namespace {
 	/*
 	 * A small helper utility to wrap around a numeric range
 	*/
 	class bounded_int {
 	private:
 		int i;
 		size_t n;
 	public:
 		bounded_int(int i, size_t n) : i(i), n(n) {}
 		bounded_int& operator--() {
 			--i;
 			if (i == -1) {
 				i = (int) n - 1;
 			}
 			return *this;
 		}
 		bounded_int& operator++() {
 			++i;
 			if (i == (int) n) {
 				i = 0;
 			}
 			return *this;
 		}
 		operator int() { return i; }
 	};
 }
 bool IfcGeom::Kernel::wire_intersections(const TopoDS_Wire& wire, TopTools_ListOfShape& wires) {
 	if (getValue(GV_NO_WIRE_INTERSECTION_CHECK) > 0.) {
 		return false;
 	}
 	if (!wire.Closed()) {
 		wires.Append(wire);
 		return false;
 	}
 	int n = count(wire, TopAbs_EDGE);
 	if (n < 3) {
 		wires.Append(wire);
 		return false;
 	}
 	// Note: initialize empty
 	Handle(ShapeExtend_WireData) wd = new ShapeExtend_WireData();
 	// ... to be sure to get consecutive edges
 	BRepTools_WireExplorer exp(wire);
 	IfcGeom::impl::tree<int> tree;
 	int edge_idx = 0;
 	for (; exp.More(); exp.Next()) {
 		wd->Add(exp.Current());
 		if (n > 64) {
 			// tfk: indices in tree are 0-based vd 1-based in wiredata
 			tree.add(edge_idx++, exp.Current());
 		}
 	}
 	if (wd->NbEdges() != n) {
 		// If the number of edges differs, BRepTools_WireExplorer did not
 		// reach every edge, probably due to loops exactly at vertex locations.
 		// This is not supported by this algorithm which only elimates loops
 		// due to edge crossings.
 		throw geometry_exception("Invalid loop");
 	}
 	bool intersected = false;
 	// tfk: Extrema on infinite curves proved to be more robust.
 	// TopoDS_Face face = BRepBuilderAPI_MakeFace(wire, true).Face();
 	// ShapeAnalysis_Wire saw(wd, face, getValue(GV_PRECISION));
 	double eps = 0;
 	if (getValue(GV_NO_WIRE_INTERSECTION_TOLERANCE) < 0.) {
 		eps = faceset_helper_
 			  // eps is added to both ends of the parametric domain, so 3. is chosen to be on the safe side here.
 			  ? (faceset_helper_->epsilon() / 3.)
 			  // @todo re-evaluate 2. here for the reasons above:
 			  : (std::min)(util::min_edge_length(wire) / 2., getValue(GV_PRECISION) * 10.);
 	}
 	// @todo: should this start from 0 in case of n > 64?
 	for (int i = 2; i < n; ++i) {
 		std::vector<int> js;
 		if (n > 64) {
 			Bnd_Box b;
 			BRepBndLib::Add(wd->Edge(i + 1), b);
 			b.Enlarge(eps);
 			js = tree.select_box(b, false);
 		} else {
 			boost::push_back(js, boost::irange(0, i - 1));
 		}
 		for(std::vector<int>::const_iterator it = js.begin(); it != js.end(); ++it) {
            int j = *it;
 			if (n > 64) {
 				if (j > i) {
 					continue;
 				}
 				if ((std::max)(i, j) - (std::min)(i, j) <= 1) {
 					continue;
 				}
 			}
 			// Only check non-consecutive edges
 			if (i == n - 1 && j == 0) continue;
 			double u11, u12, u21, u22, U1, U2;
 			GeomAPI_ExtremaCurveCurve ecc(
 				BRep_Tool::Curve(wd->Edge(i + 1), u11, u12),
 				BRep_Tool::Curve(wd->Edge(j + 1), 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;
 				// tfk: code below is for ShapeAnalysis_Wire::CheckIntersectingEdges()
 				// IntRes2d_SequenceOfIntersectionPoint points2d;
 				// TColgp_SequenceOfPnt points3d;
 				// TColStd_SequenceOfReal errors;
 				// if (saw.CheckIntersectingEdges(i + 1, j + 1, points2d, points3d, errors)) {
 				if (u11 < U1 && U1 < u12 && u21 < U2 && U2 < u22) {
 					intersected = true;
 					// Explore a forward and backward cycle from the intersection point
 					for (int fb = 0; fb <= 1; ++fb) {
 						const bool forward = fb == 0;
 						BRepBuilderAPI_MakeWire mw;
 						bool first = true;
 						for (bounded_int k(j, n);;) {
 							bool intersecting = k == j || k == i;
 							if (intersecting) {
 								TopoDS_Edge e = wd->Edge(k + 1);
 								TopoDS_Vertex v1, v2;
 								TopExp::Vertices(e, v1, v2, true);
 								const TopoDS_Vertex* v = first == forward ? &v2 : &v1;
 								// gp_Pnt p2 = points3d.Value(1);
 								gp_Pnt p1 = BRep_Tool::Pnt(*v);
 								gp_Pnt pp1, pp2;
 								ecc.Points(1, pp1, pp2);
 								const gp_Pnt& p2 = k == i ? pp1 : pp2;
 								// Substitute with a new edge from/to the intersection point
 								if (p1.Distance(p2) > getValue(GV_PRECISION) * 2) {
 									double _, __;
 									Handle_Geom_Curve crv = BRep_Tool::Curve(e, _, __);
 									BRepBuilderAPI_MakeEdge me(crv, p1, p2);
 									TopoDS_Edge ed = me.Edge();
 									mw.Add(ed);
 								}
 								first = false;
 							} else {
 								// Re-use original edge
 								mw.Add(wd->Edge(k + 1));
 							}
 							if (k == i) {
 								break;
 							}
 							if (forward) {
 								++k;
 							} else {
 								--k;
 							}
 						}
 						// Recursively process both cuts
 						wire_intersections(mw.Wire(), wires);
 					}
 					return true;
 				}
 			}
 		}
 	}
 	// No intersections found, append original wire
 	if (!intersected) {
 		wires.Append(wire);
 	}
 	return intersected;
 }
 void IfcGeom::Kernel::select_largest(const TopTools_ListOfShape& shapes, TopoDS_Shape& largest) {
 	double mass = 0.;
 	TopTools_ListIteratorOfListOfShape it(shapes);
 	for (; it.More(); it.Next()) {
 		/*
 		// tfk: bounding box is more efficient probably
 		const TopoDS_Wire& w = TopoDS::Wire(it.Value());
 		TopoDS_Face face = BRepBuilderAPI_MakeFace(w).Face();
 		const double m = face_area(face);
 		*/
 		Bnd_Box bb;
 		BRepBndLib::AddClose(it.Value(), bb);
 		double xyz_min[3], xyz_max[3];
 		bb.Get(xyz_min[0], xyz_min[1], xyz_min[2], xyz_max[0], xyz_max[1], xyz_max[2]);
 		const double eps = getValue(GV_PRECISION);
 		double m = 1.;
 		for (int i = 0; i < 3; ++i) {
 			if (Precision::IsNegativeInfinite(xyz_min[i])) {
 				xyz_min[i] = 0.;
 			}
 			if (Precision::IsInfinite(xyz_max[i])) {
 				xyz_max[i] = 0.;
 			}
 			m *= (xyz_max[i] + eps) - (xyz_min[i] - eps);
 		}
 		if (m > mass) {
 			mass = m;
 			largest = it.Value();
 		}
 	}
 }
 bool IfcGeom::Kernel::fit_halfspace(const TopoDS_Shape& a, const TopoDS_Shape& b, TopoDS_Shape& box, double& height) {
 	TopExp_Explorer exp(b, TopAbs_FACE);
 	if (!exp.More()) {
@@ -4105,6 +3566,16 @@ bool IfcGeom::Kernel::boolean_operation(const TopoDS_Shape& a_input, const TopTo
 		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.
@@ -4543,361 +4014,6 @@ bool IfcGeom::Kernel::boolean_operation(const TopoDS_Shape& a, const TopoDS_Shap
 }
 #endif
 namespace {
 	void find_neighbours(IfcGeom::impl::tree<int>& tree, std::vector<std::unique_ptr<gp_Pnt>>& pnts, std::set<int>& visited, int p, double eps) {
 		visited.insert(p);
 		Bnd_Box b;
 		b.Set(*pnts[p].get());
 		b.Enlarge(eps);
 		std::vector<int> js = tree.select_box(b, false);
 		for (int j : js) {
 			visited.insert(j);
 #ifdef FACESET_HELPER_RECURSIVE
 			if (visited.find(j) == visited.end()) {
 				// @todo, making this recursive removes the dependence on the initial ordering, but will
 				// likely result in empty results when all vertices are within 1 eps from another point.
 				find_neighbours(tree, pnts, visited, j, eps);
 			}
 #endif
 		}
 	}
 }
 template <typename CP, typename LP>
 IfcGeom::Kernel::faceset_helper<CP, LP>::~faceset_helper() {
 	// @todo this is super ugly, but how else can we be notified that the unique_ptr goes out of scope?
 	// Perhaps just supply a custom std::deleter?
 	kernel_->faceset_helper_ = nullptr;
 }
 template <typename CP, typename LP>
 bool IfcGeom::Kernel::faceset_helper<CP, LP>::construct(const IfcSchema::IfcCartesianPoint* cp, gp_Pnt* l) {
 	return kernel_->convert(cp, *l);
 }
 template <typename CP, typename LP>
 bool IfcGeom::Kernel::faceset_helper<CP, LP>::construct(const std::vector<double>& cp, gp_Pnt* l) {
 	if (cp.size() != 3) {
 		return false;
 	}
 	auto LU = kernel_->getValue(GV_LENGTH_UNIT);
 	l->SetCoord(cp[0] * LU, cp[1] * LU, cp[2] * LU);
 	return true;
 }
 /*
 template <typename CP, typename LP>
 IfcGeom::Kernel::faceset_helper<CP, LP>::faceset_helper(Kernel* kernel, const IfcSchema::IfcConnectedFaceSet* l)
 	: kernel_(kernel)
 	, non_manifold_(false)
 {
 	kernel->faceset_helper_ = this;
 	IfcSchema::IfcCartesianPoint::list::ptr points = IfcParse::traverse((IfcUtil::IfcBaseClass*) l)->as<IfcSchema::IfcCartesianPoint>();
 	std::vector<std::unique_ptr<gp_Pnt>> pnts(std::distance(points->begin(), points->end()));
 	std::vector<TopoDS_Vertex> vertices(pnts.size());
 	IfcGeom::impl::tree<int> tree;
 	BRep_Builder B;
 	Bnd_Box box;
 	for (size_t i = 0; i < points->size(); ++i) {
 		gp_Pnt* p = new gp_Pnt();
 		if (kernel->convert(*(points->begin() + i), *p)) {
 			pnts[i].reset(p);
 			B.MakeVertex(vertices[i], *p, Precision::Confusion());
 			tree.add(i, vertices[i]);
 			box.Add(*p);
 		} else {
 			delete p;
 		}
 	}
 	// Use the bbox diagonal to influence local epsilon
 	// double bdiff = std::sqrt(box.SquareExtent());
 	// @todo the bounding box diagonal is not used (see above)
 	// because we're explicitly interested in the miminal
 	// dimension of the element to limit the tolerance (for sheet-
 	// like elements for example). But the way below is very
 	// dependent on orientation due to the usage of the
 	// axis-aligned bounding box. Use PCA to find three non-aligned
 	// set of dimensions and use the one with the smallest eigenvalue.
 	// Find the minimal bounding box edge
 	double bmin[3], bmax[3];
 	box.Get(bmin[0], bmin[1], bmin[2], bmax[0], bmax[1], bmax[2]);
 	double bdiff = std::numeric_limits<double>::infinity();
 	for (size_t i = 0; i < 3; ++i) {
 		const double d = bmax[i] - bmin[i];
 		if (d > kernel->getValue(GV_PRECISION) * 10. && d < bdiff) {
 			bdiff = d;
 		}
 	}
 	eps_ = kernel->getValue(GV_PRECISION) * 10. * (std::min)(1.0, bdiff);
 	// @todo, there a tiny possibility that the duplicate faces are triggered
 	// for an internal boundary, that is also present as an external boundary.
 	// This will result in non-manifold configuration then, but this is deemed
 	// such as corner-case that it is not considered.
 	IfcSchema::IfcPolyLoop::list::ptr loops = IfcParse::traverse((IfcUtil::IfcBaseClass*)l)->as<IfcSchema::IfcPolyLoop>();
 	size_t loops_removed, non_manifold, duplicate_faces;
 	std::map<std::pair<int, int>, int> edge_use;
 	for (int i = 0; i < 3; ++i) {
 		// Some times files, have large tolerance values specified collapsing too many vertices.
 		// This case we detect below and re-run the loop with smaller epsilon. Normally
 		// the body of this loop would only be executed once.
 		loops_removed = 0;
 		non_manifold = 0;
 		duplicate_faces = 0;
 		vertex_mapping_.clear();
 		duplicates_.clear();
 		edge_use.clear();
 		if (eps_ < Precision::Confusion()) {
 			// occt uses some hard coded precision values, don't go smaller than that.
 			// @todo, can be reset though with BRepLib::Precision(double)
 			eps_ = Precision::Confusion();
 		}
 		for (int pnt_i = 0; pnt_i < (int)pnts.size(); ++pnt_i) {
 			if (pnts[pnt_i]) {
 				std::set<int> vs;
 				find_neighbours(tree, pnts, vs, pnt_i, eps_);
 				for (int v : vs) {
 					auto pt = *(points->begin() + v);
 					// NB: insert() ignores duplicate keys
 					vertex_mapping_.insert({ get_idx(pt), pnt_i });
 				}
 			}
 		}
 		typedef std::array<int, 2> edge_t;
 		typedef std::set<edge_t> edge_set_t;
 		std::set<edge_set_t> edge_sets;
 		for (auto& loop : *loops) {
 			auto ps = loop->Polygon();
 			std::vector<std::pair<int, int> > segments;
 			edge_set_t segment_set;
 			loop_(ps, [&segments, &segment_set](int C, int D, bool) {
 				segment_set.insert(edge_t{C,D});
 				segments.push_back(std::make_pair(C, D));
 			});
 			if (edge_sets.find(segment_set) != edge_sets.end()) {
 				duplicate_faces++;
 				duplicates_.insert(loop);
 				continue;
 			}
 			edge_sets.insert(segment_set);
 			if (segments.size() >= 3) {
 				for (auto& p : segments) {
 					edge_use[p] ++;
 				}
 			} else {
 				loops_removed += 1;
 			}
 		}
 		if (edge_use.size() != 0) {
 			break;
 		} else {
 			eps_ /= 10.;
 		}
 	}
 	for (auto& p : edge_use) {
 		int a, b;
 		std::tie(a, b) = p.first;
 		edges_[p.first] = BRepBuilderAPI_MakeEdge(vertices[a], vertices[b]);
 		if (p.second != 2) {
 			non_manifold += 1;
 		}
 	}
 	if (loops_removed || (non_manifold && l->declaration().is(IfcSchema::IfcClosedShell::Class()))) {
 		Logger::Warning(boost::lexical_cast<std::string>(duplicate_faces) + " duplicate faces removed, " + boost::lexical_cast<std::string>(loops_removed) + " loops removed and " + boost::lexical_cast<std::string>(non_manifold) + " non-manifold edges for:", l);
 	}
 }
 */
 namespace {
 	const std::vector<std::vector<double>>* store_cache(const std::vector<std::vector<double>>& p) {
 		return &p;
 	}
 	const std::vector<std::vector<double>>* store_cache(const std::vector<const IfcSchema::IfcCartesianPoint*>& /*p*/) {
 		return nullptr;
 	}
 }
 template <typename CP, typename LP>
 IfcGeom::Kernel::faceset_helper<CP, LP>::faceset_helper(
 	Kernel* kernel,
 	const std::vector<CP>& points,
 	const std::vector<LP>& indices,
 	bool should_be_closed
 )
 	: kernel_(kernel)
 	, non_manifold_(false)
 	, points_(store_cache(points))
 {
 	std::vector<std::unique_ptr<gp_Pnt>> pnts(std::distance(points.begin(), points.end()));
 	std::vector<TopoDS_Vertex> vertices(pnts.size());
 	IfcGeom::impl::tree<int> tree;
 	BRep_Builder B;
 	Bnd_Box box;
 	for (size_t i = 0; i < points.size(); ++i) {
 		gp_Pnt* p = new gp_Pnt;
 		if (construct(points[i], p)) {
 			pnts[i].reset(p);
 			B.MakeVertex(vertices[i], *p, Precision::Confusion());
 			tree.add((int) i, vertices[i]);
 			box.Add(*p);
 		} else {
 			delete p;
 		}
 	}
 	// Use the bbox diagonal to influence local epsilon
 	// double bdiff = std::sqrt(box.SquareExtent());
 	// @todo the bounding box diagonal is not used (see above)
 	// because we're explicitly interested in the miminal
 	// dimension of the element to limit the tolerance (for sheet-
 	// like elements for example). But the way below is very
 	// dependent on orientation due to the usage of the
 	// axis-aligned bounding box. Use PCA to find three non-aligned
 	// set of dimensions and use the one with the smallest eigenvalue.
 	// Find the minimal bounding box edge
 	double bmin[3], bmax[3];
 	box.Get(bmin[0], bmin[1], bmin[2], bmax[0], bmax[1], bmax[2]);
 	double bdiff = std::numeric_limits<double>::infinity();
 	for (size_t i = 0; i < 3; ++i) {
 		const double d = bmax[i] - bmin[i];
 		if (d > kernel->getValue(GV_PRECISION) * 10. && d < bdiff) {
 			bdiff = d;
 		}
 	}
 	eps_ = kernel->getValue(GV_PRECISION) * 10. * (std::min)(1.0, bdiff);
 	size_t loops_removed, non_manifold, duplicate_faces;
 	std::map<std::pair<int, int>, int> edge_use;
 	for (int i = 0; i < 3; ++i) {
 		// Some times files, have large tolerance values specified collapsing too many vertices.
 		// This case we detect below and re-run the loop with smaller epsilon. Normally
 		// the body of this loop would only be executed once.
 		loops_removed = 0;
 		non_manifold = 0;
 		duplicate_faces = 0;
 		vertex_mapping_.clear();
 		duplicates_.clear();
 		edge_use.clear();
 		if (eps_ < Precision::Confusion()) {
 			// occt uses some hard coded precision values, don't go smaller than that.
 			// @todo, can be reset though with BRepLib::Precision(double)
 			eps_ = Precision::Confusion();
 		}
 		for (int pnt_i = 0; pnt_i < (int)pnts.size(); ++pnt_i) {
 			if (pnts[pnt_i]) {
 				std::set<int> vs;
 				find_neighbours(tree, pnts, vs, pnt_i, eps_);
 				for (int v : vs) {
 					// NB: insert() ignores duplicate keys
 					// v-1?
 					vertex_mapping_.insert({ get_idx(points[v]), pnt_i });
 				}
 			}
 		}
 		typedef std::array<int, 2> edge_t;
 		typedef std::set<edge_t> edge_set_t;
 		std::set<edge_set_t> edge_sets;
 		for (auto ps = indices.begin(); ps != indices.end(); ++ps) {
 			std::vector<std::pair<int, int> > segments;
 			edge_set_t segment_set;
 			loop_(*ps, [&segments, &segment_set](int C, int D, bool) {
 				segment_set.insert(edge_t{ C,D });
 				segments.push_back(std::make_pair(C, D));
 			});
 			if (edge_sets.find(segment_set) != edge_sets.end()) {
 				duplicate_faces++;
 				duplicates_.insert(util::conditional_address_of(*ps));
 				continue;
 			}
 			edge_sets.insert(segment_set);
 			if (segments.size() >= 3) {
 				for (auto& p : segments) {
 					edge_use[p] ++;
 				}
 			}
 			else {
 				loops_removed += 1;
 			}
 		}
 		if (edge_use.size() != 0) {
 			break;
 		}
 		else {
 			eps_ /= 10.;
 		}
 	}
 	for (auto& p : edge_use) {
 		int a, b;
 		std::tie(a, b) = p.first;
 		edges_[p.first] = BRepBuilderAPI_MakeEdge(vertices[a], vertices[b]);
 		if (p.second != 2) {
 			non_manifold += 1;
 		}
 	}
 	if (duplicates_.size() || loops_removed || (non_manifold && should_be_closed)) {
 		Logger::Warning(boost::lexical_cast<std::string>(duplicate_faces) + " duplicate faces removed, " + boost::lexical_cast<std::string>(loops_removed) + " degenerate loops eliminated and " + boost::lexical_cast<std::string>(non_manifold) + " non-manifold edges");
 	}
 }
 template class IfcGeom::Kernel::faceset_helper<const IfcSchema::IfcCartesianPoint*, const IfcSchema::IfcPolyLoop*>;
 template class IfcGeom::Kernel::faceset_helper<std::vector<double>, std::vector<int>>;
 void IfcGeom::Kernel::set_conversion_placement_rel_to_type(const IfcParse::declaration* type) {
 	placement_rel_to_type_ = type;
 }
@@ -119,28 +119,8 @@ private:
 		const std::vector<std::vector<double>>* points_ = nullptr;
 		double eps_;
 		bool non_manifold_;
-
+		
-		template <typename Fn>
+		void loop_(const LP& lp, const std::function<void(int, int, bool)>& callback);
 		void loop_(const LP& lp, const Fn& callback) {
 			auto ps = get_idxs(lp);
 			if (ps.size() < 3) {
 				return;
 			}
 			auto A = ps.back();
 			for (auto& B : ps) {
 				auto C = vertex_mapping_[A], D = vertex_mapping_[B];
 				bool fwd = C < D;
 				if (!fwd) {
 					std::swap(C, D);
 				}
 				if (C != D) {
 					callback(C, D, fwd);
 					A = B;
 				}
 			}
 		}
 		bool construct(const IfcSchema::IfcCartesianPoint* cp, gp_Pnt* l);
 		bool construct(const std::vector<double>& cp, gp_Pnt* l);
@@ -153,32 +133,9 @@ private:
 			return &cp;
 		}
-		std::vector<const void*> get_idxs(const IfcSchema::IfcPolyLoop* lp) {
+		std::vector<const void*> get_idxs(const IfcSchema::IfcPolyLoop* lp);
-			auto poly = lp->Polygon();
+		std::vector<const void*> get_idxs(const std::vector<int>& it);
 			std::vector<const void*> idxs;
 			std::transform(poly->begin(), poly->end(), std::back_inserter(idxs), [this](const IfcSchema::IfcCartesianPoint* p) {return get_idx(p); });
 			return idxs;
 		}
 		std::vector<const void*> get_idxs(const std::vector<int>& it) {
 			std::vector<const void*> idxs;
 			std::transform(it.begin(), it.end(), std::back_inserter(idxs), [this](int i) { return get_idx((*points_)[i - 1]); });
 			return idxs;
 		}
 		/*
 		std::vector<uintptr_t> get_idxs(std::vector<std::vector<int>>::const_iterator it) {
 			std::vector<uintptr_t> idxs;
 			std::transform(it->begin(), it->end(), std::back_inserter(idxs), [this](int i) {return get_idx(i); });
 			return idxs;
 		}
 		*/
 	public:
 		/*
 		faceset_helper(MAKE_TYPE_NAME(Kernel)* kernel, const IfcSchema::IfcConnectedFaceSet* l);
 		*/
 		faceset_helper(
 			MAKE_TYPE_NAME(Kernel)* kernel, 
 			const std::vector<CP>& points,
@@ -189,64 +146,11 @@ private:
 		bool non_manifold() const { return non_manifold_; }
 		bool& non_manifold() { return non_manifold_; }
 		double epsilon() const { return eps_; }
 		bool edge(int A, int B, TopoDS_Edge& e);
-		/*
+		bool wire(const LP& loop, TopoDS_Wire& wire);		
 		bool edge(CP a, CP b, TopoDS_Edge& e) {
 			int A = vertex_mapping_[get_idx(a)];
 			int B = vertex_mapping_[get_idx(b)];
 			if (A == B) {
 				return false;
 			}
 			return edge(A, B, e);
 		}
 		*/
 		bool edge(int A, int B, TopoDS_Edge& e) {
 			auto it = edges_.find({A, B});
 			if (it == edges_.end()) {
 				return false;
 			}
 			e = it->second;
 			return true;
 		}
 		bool wire(const LP& loop, TopoDS_Wire& wire) {
 			if (duplicates_.find(util::conditional_address_of(loop)) != duplicates_.end()) {
 				return false;
 			}
 			BRep_Builder builder;
 			builder.MakeWire(wire);
 			int count = 0;
 			loop_(loop, [this, &builder, &wire, &count](int A, int B, bool fwd) {
 				TopoDS_Edge e;
 				if (edge(A, B, e)) {
 					if (!fwd) {
 						e.Reverse();
 					}
 					builder.Add(wire, e);
 					count += 1;
 				}
 			});
 			if (count >= 3) {
 				wire.Closed(true);
 				TopTools_ListOfShape results;
 				if (kernel_->wire_intersections(wire, results)) {
 					Logger::Warning("Self-intersections with " + boost::lexical_cast<std::string>(results.Extent()) + " cycles detected");
 					kernel_->select_largest(results, wire);
 					non_manifold_ = true;
 				}
 				return true;
 			} else {
 				return false;
 			}
 		}
 		double epsilon() const {
 			return eps_;
 		}
 	};
 	double deflection_tolerance;
@@ -355,6 +259,7 @@ public:
 	void set_offset(const std::array<double, 3>& offset);
 	void set_rotation(const std::array<double, 4>& rotation);
 	double get_wire_intersection_tolerance(const TopoDS_Wire&) const;
 	bool convert_wire_to_face(const TopoDS_Wire& wire, TopoDS_Face& face);
 	bool convert_wire_to_faces(const TopoDS_Wire& wire, TopoDS_Compound& face);
@@ -419,12 +324,6 @@ public:
 	void remove_collinear_points_from_loop(TColgp_SequenceOfPnt& polygon, bool closed, double tol=-1.);
 	bool wire_to_sequence_of_point(const TopoDS_Wire&, TColgp_SequenceOfPnt&);
 	void sequence_of_point_to_wire(const TColgp_SequenceOfPnt&, TopoDS_Wire&, bool closed);
 	bool approximate_plane_through_wire(const TopoDS_Wire&, gp_Pln&, double eps=-1.);
 	bool flatten_wire(TopoDS_Wire&);
 	/// Triangulate the set of wires. The firstmost wire is assumed to be the outer wire.
 	bool triangulate_wire(const std::vector<TopoDS_Wire>&, TopTools_ListOfShape&);
 	bool wire_intersections(const TopoDS_Wire & wire, TopTools_ListOfShape & wires);
 	void select_largest(const TopTools_ListOfShape& shapes, TopoDS_Shape& largest);
 	static double shape_volume(const TopoDS_Shape& s);
 	static double face_area(const TopoDS_Face& f);
@@ -17,11 +17,13 @@
 *                                                                              *
 ********************************************************************************/
 #include "../ifcgeom/IfcGeom.h"
 #include "../ifcgeom_schema_agnostic/wire_utils.h"
 #include <gp_Pnt.hxx>
 #include <BRepBuilderAPI_MakePolygon.hxx>
 #include <TopoDS_Wire.hxx>
 #include <TopTools_ListOfShape.hxx>
 #include "../ifcgeom/IfcGeom.h"
 #define _USE_MATH_DEFINES
 #define Kernel MAKE_TYPE_NAME(Kernel)
@@ -68,9 +70,9 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcPolyLoop* l, TopoDS_Wire& resu
 	result = w.Wire();
 	TopTools_ListOfShape results;
-	if (wire_intersections(result, results)) {
+	if (getValue(GV_NO_WIRE_INTERSECTION_CHECK) == 0. && util::wire_intersections(result, results, get_wire_intersection_tolerance(result), getValue(GV_PRECISION))) {
 		Logger::Error("Self-intersections with " + boost::lexical_cast<std::string>(results.Extent()) + " cycles detected", l);
-		select_largest(results, result);
+		util::select_largest(results, result);
 	}
 	return true;
@@ -22,6 +22,7 @@
 #include <TopoDS_Wire.hxx>
 #include <TopoDS_Face.hxx>
 #include "../ifcgeom/IfcGeom.h"
 #include "../ifcgeom_schema_agnostic/wire_utils.h"
 #define Kernel MAKE_TYPE_NAME(Kernel)
@@ -111,7 +112,7 @@ bool IfcGeom::Kernel::convert(const IfcSchema::IfcPolygonalFaceSet* pfs, TopoDS_
 		if (not_planar) {
 			TopTools_ListOfShape fs;
-			if (triangulate_wire(ws, fs)) {
+			if (util::triangulate_wire(ws, fs)) {
 				Logger::Warning("Triangulated face boundary:", pfs);
 				TopTools_ListIteratorOfListOfShape it(fs);
 				for (; it.More(); it.Next()) {
@@ -0,0 +1,301 @@
 #include "IfcGeom.h"
 #include "../ifcgeom_schema_agnostic/IfcGeomTree.h"
 #include "../ifcgeom_schema_agnostic/wire_utils.h"
 #define Kernel MAKE_TYPE_NAME(Kernel)
 namespace {
 	void find_neighbours(IfcGeom::impl::tree<int>& tree, std::vector<std::unique_ptr<gp_Pnt>>& pnts, std::set<int>& visited, int p, double eps) {
 		visited.insert(p);
 		Bnd_Box b;
 		b.Set(*pnts[p].get());
 		b.Enlarge(eps);
 		std::vector<int> js = tree.select_box(b, false);
 		for (int j : js) {
 			visited.insert(j);
 #ifdef FACESET_HELPER_RECURSIVE
 			if (visited.find(j) == visited.end()) {
 				// @todo, making this recursive removes the dependence on the initial ordering, but will
 				// likely result in empty results when all vertices are within 1 eps from another point.
 				find_neighbours(tree, pnts, visited, j, eps);
 			}
 #endif
 		}
 	}
 }
 namespace {
 	const std::vector<std::vector<double>>* store_cache(const std::vector<std::vector<double>>& p) {
 		return &p;
 	}
 	const std::vector<std::vector<double>>* store_cache(const std::vector<const IfcSchema::IfcCartesianPoint*>& /*p*/) {
 		return nullptr;
 	}
 }
 template <typename CP, typename LP>
 IfcGeom::Kernel::faceset_helper<CP, LP>::faceset_helper(
 	Kernel* kernel,
 	const std::vector<CP>& points,
 	const std::vector<LP>& indices,
 	bool should_be_closed
 )
 	: kernel_(kernel)
 	, non_manifold_(false)
 	, points_(store_cache(points))
 {
 	std::vector<std::unique_ptr<gp_Pnt>> pnts(std::distance(points.begin(), points.end()));
 	std::vector<TopoDS_Vertex> vertices(pnts.size());
 	IfcGeom::impl::tree<int> tree;
 	BRep_Builder B;
 	Bnd_Box box;
 	for (size_t i = 0; i < points.size(); ++i) {
 		gp_Pnt* p = new gp_Pnt;
 		if (construct(points[i], p)) {
 			pnts[i].reset(p);
 			B.MakeVertex(vertices[i], *p, Precision::Confusion());
 			tree.add((int)i, vertices[i]);
 			box.Add(*p);
 		} else {
 			delete p;
 		}
 	}
 	// Use the bbox diagonal to influence local epsilon
 	// double bdiff = std::sqrt(box.SquareExtent());
 	// @todo the bounding box diagonal is not used (see above)
 	// because we're explicitly interested in the miminal
 	// dimension of the element to limit the tolerance (for sheet-
 	// like elements for example). But the way below is very
 	// dependent on orientation due to the usage of the
 	// axis-aligned bounding box. Use PCA to find three non-aligned
 	// set of dimensions and use the one with the smallest eigenvalue.
 	// Find the minimal bounding box edge
 	double bmin[3], bmax[3];
 	box.Get(bmin[0], bmin[1], bmin[2], bmax[0], bmax[1], bmax[2]);
 	double bdiff = std::numeric_limits<double>::infinity();
 	for (size_t i = 0; i < 3; ++i) {
 		const double d = bmax[i] - bmin[i];
 		if (d > kernel->getValue(GV_PRECISION) * 10. && d < bdiff) {
 			bdiff = d;
 		}
 	}
 	eps_ = kernel->getValue(GV_PRECISION) * 10. * (std::min)(1.0, bdiff);
 	size_t loops_removed, non_manifold, duplicate_faces;
 	std::map<std::pair<int, int>, int> edge_use;
 	for (int i = 0; i < 3; ++i) {
 		// Some times files, have large tolerance values specified collapsing too many vertices.
 		// This case we detect below and re-run the loop with smaller epsilon. Normally
 		// the body of this loop would only be executed once.
 		loops_removed = 0;
 		non_manifold = 0;
 		duplicate_faces = 0;
 		vertex_mapping_.clear();
 		duplicates_.clear();
 		edge_use.clear();
 		if (eps_ < Precision::Confusion()) {
 			// occt uses some hard coded precision values, don't go smaller than that.
 			// @todo, can be reset though with BRepLib::Precision(double)
 			eps_ = Precision::Confusion();
 		}
 		for (int pnt_i = 0; pnt_i < (int)pnts.size(); ++pnt_i) {
 			if (pnts[pnt_i]) {
 				std::set<int> vs;
 				find_neighbours(tree, pnts, vs, pnt_i, eps_);
 				for (int v : vs) {
 					// NB: insert() ignores duplicate keys
 					// v-1?
 					vertex_mapping_.insert({ get_idx(points[v]), pnt_i });
 				}
 			}
 		}
 		std::set<std::tuple<double, double, double>> unique;
 		for (int pnt_i = 0; pnt_i < (int)pnts.size(); ++pnt_i) {
 			if (pnts[pnt_i]) {
 				unique.insert(std::make_tuple(
 					(*pnts[pnt_i]).X(),
 					(*pnts[pnt_i]).Y(),
 					(*pnts[pnt_i]).Z()
 				));
 			}
 		}
 		Logger::Notice("Collapsed vertices from " + std::to_string(pnts.size()) + " (" + std::to_string(unique.size()) + " unique) to " + std::to_string(vertex_mapping_.size()));
 		typedef std::array<int, 2> edge_t;
 		typedef std::set<edge_t> edge_set_t;
 		std::set<edge_set_t> edge_sets;
 		for (auto ps = indices.begin(); ps != indices.end(); ++ps) {
 			std::vector<std::pair<int, int> > segments;
 			edge_set_t segment_set;
 			loop_(*ps, [&segments, &segment_set](int C, int D, bool) {
 				segment_set.insert(edge_t{ C,D });
 				segments.push_back(std::make_pair(C, D));
 			});
 			if (edge_sets.find(segment_set) != edge_sets.end()) {
 				duplicate_faces++;
 				duplicates_.insert(util::conditional_address_of(*ps));
 				continue;
 			}
 			edge_sets.insert(segment_set);
 			if (segments.size() >= 3) {
 				for (auto& p : segments) {
 					edge_use[p] ++;
 				}
 			} else {
 				loops_removed += 1;
 			}
 		}
 		if (edge_use.size() != 0) {
 			break;
 		} else {
 			eps_ /= 10.;
 		}
 	}
 	for (auto& p : edge_use) {
 		int a, b;
 		std::tie(a, b) = p.first;
 		edges_[p.first] = BRepBuilderAPI_MakeEdge(vertices[a], vertices[b]);
 		if (p.second != 2) {
 			non_manifold += 1;
 		}
 	}
 	if (duplicates_.size() || loops_removed || (non_manifold && should_be_closed)) {
 		Logger::Warning(boost::lexical_cast<std::string>(duplicate_faces) + " duplicate faces removed, " + boost::lexical_cast<std::string>(loops_removed) + " degenerate loops eliminated and " + boost::lexical_cast<std::string>(non_manifold) + " non-manifold edges");
 	}
 }
 template <typename CP, typename LP>
 void IfcGeom::Kernel::faceset_helper<CP, LP>::loop_(const LP& lp, const std::function<void(int, int, bool)>& callback) {
 	auto ps = get_idxs(lp);
 	if (ps.size() < 3) {
 		return;
 	}
 	auto A = ps.back();
 	for (auto& B : ps) {
 		auto C = vertex_mapping_[A], D = vertex_mapping_[B];
 		bool fwd = C < D;
 		if (!fwd) {
 			std::swap(C, D);
 		}
 		if (C != D) {
 			callback(C, D, fwd);
 			A = B;
 		}
 	}
 }
 template <typename CP, typename LP>
 std::vector<const void*> IfcGeom::Kernel::faceset_helper<CP, LP>::get_idxs(const IfcSchema::IfcPolyLoop* lp) {
 	auto poly = lp->Polygon();
 	std::vector<const void*> idxs;
 	std::transform(poly->begin(), poly->end(), std::back_inserter(idxs), [this](const IfcSchema::IfcCartesianPoint* p) {return get_idx(p); });
 	return idxs;
 }
 template <typename CP, typename LP>
 std::vector<const void*> IfcGeom::Kernel::faceset_helper<CP, LP>::get_idxs(const std::vector<int>& it) {
 	std::vector<const void*> idxs;
 	std::transform(it.begin(), it.end(), std::back_inserter(idxs), [this](int i) { return get_idx((*points_)[i - 1]); });
 	return idxs;
 }
 template <typename CP, typename LP>
 bool IfcGeom::Kernel::faceset_helper<CP, LP>::edge(int A, int B, TopoDS_Edge& e) {
 	auto it = edges_.find({ A, B });
 	if (it == edges_.end()) {
 		return false;
 	}
 	e = it->second;
 	return true;
 }
 template <typename CP, typename LP>
 bool IfcGeom::Kernel::faceset_helper<CP, LP>::wire(const LP& loop, TopoDS_Wire& wire) {
 	if (duplicates_.find(util::conditional_address_of(loop)) != duplicates_.end()) {
 		return false;
 	}
 	BRep_Builder builder;
 	builder.MakeWire(wire);
 	int count = 0;
 	loop_(loop, [this, &builder, &wire, &count](int A, int B, bool fwd) {
 		TopoDS_Edge e;
 		if (edge(A, B, e)) {
 			if (!fwd) {
 				e.Reverse();
 			}
 			builder.Add(wire, e);
 			count += 1;
 		}
 	});
 	if (count >= 3) {
 		wire.Closed(true);
 		TopTools_ListOfShape results;
 		if (kernel_->getValue(GV_NO_WIRE_INTERSECTION_CHECK) == 0. && util::wire_intersections(wire, results, kernel_->get_wire_intersection_tolerance(wire), kernel_->getValue(IfcGeom::Kernel::GV_PRECISION))) {
 			Logger::Warning("Self-intersections with " + boost::lexical_cast<std::string>(results.Extent()) + " cycles detected");
 			util::select_largest(results, wire);
 			non_manifold_ = true;
 		}
 		return true;
 	} else {
 		return false;
 	}
 }
 template <typename CP, typename LP>
 IfcGeom::Kernel::faceset_helper<CP, LP>::~faceset_helper() {
 	// @todo this is super ugly, but how else can we be notified that the unique_ptr goes out of scope?
 	// Perhaps just supply a custom std::deleter?
 	kernel_->faceset_helper_ = nullptr;
 }
 template <typename CP, typename LP>
 bool IfcGeom::Kernel::faceset_helper<CP, LP>::construct(const IfcSchema::IfcCartesianPoint* cp, gp_Pnt* l) {
 	return kernel_->convert(cp, *l);
 }
 template <typename CP, typename LP>
 bool IfcGeom::Kernel::faceset_helper<CP, LP>::construct(const std::vector<double>& cp, gp_Pnt* l) {
 	if (cp.size() != 3) {
 		return false;
 	}
 	auto LU = kernel_->getValue(GV_LENGTH_UNIT);
 	l->SetCoord(cp[0] * LU, cp[1] * LU, cp[2] * LU);
 	return true;
 }
 template class IfcGeom::Kernel::faceset_helper<const IfcSchema::IfcCartesianPoint*, const IfcSchema::IfcPolyLoop*>;
 template class IfcGeom::Kernel::faceset_helper<std::vector<double>, std::vector<int>>;
@@ -0,0 +1,558 @@
 #include "wire_utils.h"
 #include "../ifcparse/IfcLogger.h"
 #include "../ifcgeom_schema_agnostic/Kernel.h"
 #include "../ifcgeom_schema_agnostic/IfcGeomTree.h"
 #include <TopExp.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Iterator.hxx>
 #include <BRep_Tool.hxx>
 #include <BRepTools_WireExplorer.hxx>
 #include <BRepBuilderAPI_MakeFace.hxx>
 #include <BRepAlgo_NormalProjection.hxx>
 #include <BRepMesh_IncrementalMesh.hxx>
 #include <BRepBuilderAPI_MakePolygon.hxx>
 #include <BRepBuilderAPI_MakeWire.hxx>
 #include <BRepBuilderAPI_MakeEdge.hxx>
 #include <TopTools_ListOfShape.hxx>
 #include <ShapeExtend_WireData.hxx>
 #include <Standard_Version.hxx>
 #include <GeomAPI_ExtremaCurveCurve.hxx>
 #include <boost/range/irange.hpp>
 #include <boost/range/algorithm_ext/push_back.hpp>
 #include <map>
 bool IfcGeom::util::approximate_plane_through_wire(const TopoDS_Wire& wire, gp_Pln& plane, double eps_) {
 	// Newell's Method is used for the normal calculation
 	// as a simple edge cross product can give opposite results
 	// for a concave face boundary.
 	// Reference: Graphics Gems III p. 231
 	const double eps2 = eps_ * eps_;
 	double x = 0, y = 0, z = 0;
 	gp_Pnt current, previous, first;
 	gp_XYZ center;
 	int n = 0;
 	BRepTools_WireExplorer exp(wire);
 	for (;; exp.Next()) {
 		const bool has_more = exp.More() != 0;
 		if (has_more) {
 			const TopoDS_Vertex& v = exp.CurrentVertex();
 			current = BRep_Tool::Pnt(v);
 			center += current.XYZ();
 		} else {
 			current = first;
 		}
 		if (n) {
 			const double& xn = previous.X();
 			const double& yn = previous.Y();
 			const double& zn = previous.Z();
 			const double& xn1 = current.X();
 			const double& yn1 = current.Y();
 			const double& zn1 = current.Z();
 			x += (yn - yn1)*(zn + zn1);
 			y += (xn + xn1)*(zn - zn1);
 			z += (xn - xn1)*(yn + yn1);
 		} else {
 			first = current;
 		}
 		if (!has_more) {
 			break;
 		}
 		previous = current;
 		++n;
 	}
 	if (n < 3) {
 		return false;
 	}
 	plane = gp_Pln(center / n, gp_Dir(x, y, z));
 	exp.Init(wire);
 	for (; exp.More(); exp.Next()) {
 		const TopoDS_Vertex& v = exp.CurrentVertex();
 		current = BRep_Tool::Pnt(v);
 		if (plane.SquareDistance(current) > eps2) {
 			return false;
 		}
 	}
 	return true;
 }
 bool IfcGeom::util::flatten_wire(TopoDS_Wire& wire, double eps) {
 	gp_Pln pln;
 	if (!approximate_plane_through_wire(wire, pln, eps)) {
 		return false;
 	}
 	TopoDS_Face face = BRepBuilderAPI_MakeFace(pln).Face();
 	BRepAlgo_NormalProjection proj(face);
 	proj.Add(wire);
 	proj.Build();
 	if (!proj.IsDone()) {
 		return false;
 	}
 	TopTools_ListOfShape list;
 	proj.BuildWire(list);
 	if (list.Extent() != 1) {
 		return false;
 	}
 	wire = TopoDS::Wire(list.First());
 	return true;
 }
 IfcGeom::util::triangulate_wire_result IfcGeom::util::triangulate_wire(const std::vector<TopoDS_Wire>& wires, TopTools_ListOfShape& faces) {
 	// This is a bit of a precarious approach, but seems to work for the
 	// versions of OCCT tested for. OCCT has a Delaunay triangulation function
 	// BRepMesh_Delaun, but it is notoriously hard to interpret the results
 	// (due to the Bowyer-Watson super triangle perhaps?). Therefore
 	// alternatively we use the regular OCCT incremental mesher on a new face
 	// created from the UV coordinates of the original wire. Pray to our gods
 	// that the vertex coordinates are unaffected by the meshing algorithm and
 	// map them back to 3d coordinates when iterating over the mesh triangles.
 	// In addition, to maintain a manifold shell, we need to make sure that
 	// every edge from the input wire is used exactly once in the list of
 	// resulting faces. And that other internal edges are used twice.
 	typedef std::pair<double, double> uv_node;
 	gp_Pln pln;
 	if (!approximate_plane_through_wire(wires.front(), pln, std::numeric_limits<double>::infinity())) {
 		return TRIANGULATE_WIRE_FAIL;
 	}
 	const gp_XYZ& udir = pln.Position().XDirection().XYZ();
 	const gp_XYZ& vdir = pln.Position().YDirection().XYZ();
 	const gp_XYZ& pnt = pln.Position().Location().XYZ();
 	std::map<uv_node, TopoDS_Vertex> mapping;
 	std::map<std::pair<uv_node, uv_node>, TopoDS_Edge> existing_edges, new_edges;
 	std::unique_ptr<BRepBuilderAPI_MakeFace> mf;
 	for (auto it = wires.begin(); it != wires.end(); ++it) {
 		const TopoDS_Wire& wire = *it;
 		BRepTools_WireExplorer exp(wire);
 		BRepBuilderAPI_MakePolygon mp;
 		// Add UV coordinates to a newly created polygon
 		for (; exp.More(); exp.Next()) {
 			// Project onto plane
 			const TopoDS_Vertex& V = exp.CurrentVertex();
 			gp_Pnt p = BRep_Tool::Pnt(V);
 			double u = (p.XYZ() - pnt).Dot(udir);
 			double v = (p.XYZ() - pnt).Dot(vdir);
 			mp.Add(gp_Pnt(u, v, 0.));
 			mapping.insert(std::make_pair(std::make_pair(u, v), V));
 			// Store existing edges in a map so that triangles can
 			// actually reference the preexisting edges.
 			const TopoDS_Edge& e = exp.Current();
 			TopoDS_Vertex V0, V1;
 			TopExp::Vertices(e, V0, V1, true);
 			gp_Pnt p0 = BRep_Tool::Pnt(V0);
 			gp_Pnt p1 = BRep_Tool::Pnt(V1);
 			double u0 = (p0.XYZ() - pnt).Dot(udir);
 			double v0 = (p0.XYZ() - pnt).Dot(vdir);
 			double u1 = (p1.XYZ() - pnt).Dot(udir);
 			double v1 = (p1.XYZ() - pnt).Dot(vdir);
 			uv_node uv0 = std::make_pair(u0, v0);
 			uv_node uv1 = std::make_pair(u1, v1);
 			existing_edges.insert(std::make_pair(std::make_pair(uv0, uv1), e));
 			existing_edges.insert(std::make_pair(std::make_pair(uv1, uv0), TopoDS::Edge(e.Reversed())));
 		}
 		// Not closed by default
 		mp.Close();
 		if (mf) {
 			if (it - 1 == wires.begin()) {
 				// @todo is this necessary?
 				TopoDS_Face f = mf->Face();
 				mf->Init(f);
 			}
 			mf->Add(mp.Wire());
 		} else {
 			mf.reset(new BRepBuilderAPI_MakeFace(mp.Wire()));
 		}
 	}
 	const TopoDS_Face& face = mf->Face();
 	// Create a triangular mesh from the face
 	BRepMesh_IncrementalMesh(face, Precision::Confusion());
 	int n123[3];
 	TopLoc_Location loc;
 	Handle_Poly_Triangulation tri = BRep_Tool::Triangulation(face, loc);
 	if (!tri.IsNull()) {
 		const Poly_Array1OfTriangle& triangles = tri->Triangles();
 		for (int i = 1; i <= triangles.Length(); ++i) {
 			if (face.Orientation() == TopAbs_REVERSED)
 				triangles(i).Get(n123[2], n123[1], n123[0]);
 			else triangles(i).Get(n123[0], n123[1], n123[2]);
 			// Create polygons from the mesh vertices
 			BRepBuilderAPI_MakeWire mp2;
 			for (int j = 0; j < 3; ++j) {
 				uv_node uvnodes[2];
 				TopoDS_Vertex vs[2];
 				for (int k = 0; k < 2; ++k) {
 					const gp_Pnt& uv = tri->Node(n123[(j + k) % 3]);
 					uvnodes[k] = std::make_pair(uv.X(), uv.Y());
 					auto it = mapping.find(uvnodes[k]);
 					if (it == mapping.end()) {
 						Logger::Error("Internal error: unable to unproject uv-mesh");
 						return TRIANGULATE_WIRE_FAIL;
 					}
 					vs[k] = it->second;
 				}
 				auto it = existing_edges.find(std::make_pair(uvnodes[0], uvnodes[1]));
 				if (it != existing_edges.end()) {
 					// This is a boundary edge, reuse existing edge from wire
 					mp2.Add(it->second);
 				} else {
 					auto jt = new_edges.find(std::make_pair(uvnodes[0], uvnodes[1]));
 					if (jt != new_edges.end()) {
 						// We have already added the reverse as part of another
 						// triangle, reuse this edge.
 						mp2.Add(TopoDS::Edge(jt->second));
 					} else {
 						// This is a new internal edge. Register the reverse
 						// for reuse later. We need to be sure to reuse vertices
 						// for the edge construction because otherwise the wire
 						// builder will use geometrical proximity for vertex
 						// connections in which case the edge will be copied
 						// and no longer partner with other edges from the shell.
 						TopoDS_Edge ne = BRepBuilderAPI_MakeEdge(vs[0], vs[1]);
 						mp2.Add(ne);
 						// Store the reverse to be picked up later.
 						new_edges.insert(std::make_pair(std::make_pair(uvnodes[1], uvnodes[0]), TopoDS::Edge(ne.Reversed())));
 					}
 				}
 			}
 			BRepBuilderAPI_MakeFace mft(mp2.Wire());
 			if (mft.IsDone()) {
 				TopoDS_Face triangle_face = mft.Face();
 				TopoDS_Iterator jt(triangle_face, false);
 				for (; jt.More(); jt.Next()) {
 					const TopoDS_Wire& w = TopoDS::Wire(jt.Value());
 					if (w.Orientation() != wires.front().Orientation()) {
 						triangle_face.Reverse();
 					}
 				}
 				faces.Append(triangle_face);
 			} else {
 				Logger::Error("Internal error: missing face");
 				return TRIANGULATE_WIRE_FAIL;
 			}
 		}
 	}
 	TopTools_IndexedDataMapOfShapeListOfShape mape, mapn;
 	for (auto& wire : wires) {
 		TopExp::MapShapesAndAncestors(wire, TopAbs_EDGE, TopAbs_WIRE, mape);
 	}
 	TopTools_ListIteratorOfListOfShape it(faces);
 	for (; it.More(); it.Next()) {
 		TopExp::MapShapesAndAncestors(it.Value(), TopAbs_EDGE, TopAbs_WIRE, mapn);
 	}
 	// Validation
 	bool non_manifold = false;
 	for (int i = 1; i <= mape.Extent(); ++i) {
 #if OCC_VERSION_HEX >= 0x70000
 		TopTools_ListOfShape val;
 		if (!mapn.FindFromKey(mape.FindKey(i), val)) {
 #else
 		bool contains = false;
 		try {
 			TopTools_ListOfShape val = mapn.FindFromKey(mape.FindKey(i));
 			contains = true;
 		} catch (Standard_NoSuchObject&) {}
 		if (!contains) {
 #endif
 			// All existing edges need to exist in the new faces
 			Logger::Error("Internal error, missing edge from triangulation");
 			non_manifold = true;
 		}
 		}
 	for (int i = 1; i <= mapn.Extent(); ++i) {
 		const TopoDS_Shape& v = mapn.FindKey(i);
 		int n = mapn.FindFromIndex(i).Extent();
 		// Existing edges are boundaries with use 1
 		// New edges are internal with use 2
 		if (n != (mape.Contains(v) ? 1 : 2)) {
 			Logger::Error("Internal error, non-manifold result from triangulation");
 			non_manifold = true;
 		}
 	}
 	return non_manifold ? TRIANGULATE_WIRE_NON_MANIFOLD : TRIANGULATE_WIRE_OK;
 }
 namespace {
 	/*
 	 * A small helper utility to wrap around a numeric range
 	*/
 	class bounded_int {
 	private:
 		int i;
 		size_t n;
 	public:
 		bounded_int(int i, size_t n) : i(i), n(n) {}
 		bounded_int& operator--() {
 			--i;
 			if (i == -1) {
 				i = (int)n - 1;
 			}
 			return *this;
 		}
 		bounded_int& operator++() {
 			++i;
 			if (i == (int)n) {
 				i = 0;
 			}
 			return *this;
 		}
 		operator int() { return i; }
 	};
 }
 bool IfcGeom::util::wire_intersections(const TopoDS_Wire& wire, TopTools_ListOfShape& wires, double eps, double eps_real) {
 	if (!wire.Closed()) {
 		wires.Append(wire);
 		return false;
 	}
 	int n = IfcGeom::Kernel::count(wire, TopAbs_EDGE);
 	if (n < 3) {
 		wires.Append(wire);
 		return false;
 	}
 	// Note: initialize empty
 	Handle(ShapeExtend_WireData) wd = new ShapeExtend_WireData();
 	// ... to be sure to get consecutive edges
 	BRepTools_WireExplorer exp(wire);
 	IfcGeom::impl::tree<int> tree;
 	int edge_idx = 0;
 	for (; exp.More(); exp.Next()) {
 		wd->Add(exp.Current());
 		if (n > 64) {
 			// tfk: indices in tree are 0-based vd 1-based in wiredata
 			tree.add(edge_idx++, exp.Current());
 		}
 	}
 	if (wd->NbEdges() != n) {
 		// If the number of edges differs, BRepTools_WireExplorer did not
 		// reach every edge, probably due to loops exactly at vertex locations.
 		// This is not supported by this algorithm which only elimates loops
 		// due to edge crossings.
 		throw geometry_exception("Invalid loop");
 	}
 	bool intersected = false;
 	// tfk: Extrema on infinite curves proved to be more robust.
 	// TopoDS_Face face = BRepBuilderAPI_MakeFace(wire, true).Face();
 	// ShapeAnalysis_Wire saw(wd, face, getValue(GV_PRECISION));
 	// @todo: should this start from 0 in case of n > 64?
 	for (int i = 2; i < n; ++i) {
 		std::vector<int> js;
 		if (n > 64) {
 			Bnd_Box b;
 			BRepBndLib::Add(wd->Edge(i + 1), b);
 			b.Enlarge(eps);
 			js = tree.select_box(b, false);
 		} else {
 			boost::push_back(js, boost::irange(0, i - 1));
 		}
 		for (std::vector<int>::const_iterator it = js.begin(); it != js.end(); ++it) {
 			int j = *it;
 			if (n > 64) {
 				if (j > i) {
 					continue;
 				}
 				if ((std::max)(i, j) - (std::min)(i, j) <= 1) {
 					continue;
 				}
 			}
 			// Only check non-consecutive edges
 			if (i == n - 1 && j == 0) continue;
 			double u11, u12, u21, u22, U1, U2;
 			GeomAPI_ExtremaCurveCurve ecc(
 				BRep_Tool::Curve(wd->Edge(i + 1), u11, u12),
 				BRep_Tool::Curve(wd->Edge(j + 1), 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;
 				// tfk: code below is for ShapeAnalysis_Wire::CheckIntersectingEdges()
 				// IntRes2d_SequenceOfIntersectionPoint points2d;
 				// TColgp_SequenceOfPnt points3d;
 				// TColStd_SequenceOfReal errors;
 				// if (saw.CheckIntersectingEdges(i + 1, j + 1, points2d, points3d, errors)) {
 				if (u11 < U1 && U1 < u12 && u21 < U2 && U2 < u22) {
 					intersected = true;
 					// Explore a forward and backward cycle from the intersection point
 					for (int fb = 0; fb <= 1; ++fb) {
 						const bool forward = fb == 0;
 						BRepBuilderAPI_MakeWire mw;
 						bool first = true;
 						for (bounded_int k(j, n);;) {
 							bool intersecting = k == j || k == i;
 							if (intersecting) {
 								TopoDS_Edge e = wd->Edge(k + 1);
 								TopoDS_Vertex v1, v2;
 								TopExp::Vertices(e, v1, v2, true);
 								const TopoDS_Vertex* v = first == forward ? &v2 : &v1;
 								// gp_Pnt p2 = points3d.Value(1);
 								gp_Pnt p1 = BRep_Tool::Pnt(*v);
 								gp_Pnt pp1, pp2;
 								ecc.Points(1, pp1, pp2);
 								const gp_Pnt& p2 = k == i ? pp1 : pp2;
 								// Substitute with a new edge from/to the intersection point
 								if (p1.Distance(p2) > eps_real * 2) {
 									double _, __;
 									Handle_Geom_Curve crv = BRep_Tool::Curve(e, _, __);
 									BRepBuilderAPI_MakeEdge me(crv, p1, p2);
 									TopoDS_Edge ed = me.Edge();
 									mw.Add(ed);
 								}
 								first = false;
 							} else {
 								// Re-use original edge
 								mw.Add(wd->Edge(k + 1));
 							}
 							if (k == i) {
 								break;
 							}
 							if (forward) {
 								++k;
 							} else {
 								--k;
 							}
 						}
 						// Recursively process both cuts
 						// @todo this is a change in behaviour with eps precomputed from the kernel
 						// instead of adaptively calculated for the successive iterations.
 						wire_intersections(mw.Wire(), wires, eps, eps_real);
 					}
 					return true;
 				}
 			}
 		}
 	}
 	// No intersections found, append original wire
 	if (!intersected) {
 		wires.Append(wire);
 	}
 	return intersected;
 }
 void IfcGeom::util::select_largest(const TopTools_ListOfShape& shapes, TopoDS_Shape& largest) {
 	double mass = 0.;
 	TopTools_ListIteratorOfListOfShape it(shapes);
 	for (; it.More(); it.Next()) {
 		/*
 		// tfk: bounding box is more efficient probably
 		const TopoDS_Wire& w = TopoDS::Wire(it.Value());
 		TopoDS_Face face = BRepBuilderAPI_MakeFace(w).Face();
 		const double m = face_area(face);
 		*/
 		Bnd_Box bb;
 		BRepBndLib::AddClose(it.Value(), bb);
 		double xyz_min[3], xyz_max[3];
 		bb.Get(xyz_min[0], xyz_min[1], xyz_min[2], xyz_max[0], xyz_max[1], xyz_max[2]);
 		// @todo hard coded precision.
 		// @todo this is a really strange measure for wire size. Why not use newell's
 		// method to project to plane and then calculate size of the 2d bbox?
 		const double eps = 1.e-5;
 		double m = 1.;
 		for (int i = 0; i < 3; ++i) {
 			if (Precision::IsNegativeInfinite(xyz_min[i])) {
 				xyz_min[i] = 0.;
 			}
 			if (Precision::IsInfinite(xyz_max[i])) {
 				xyz_max[i] = 0.;
 			}
 			m *= (xyz_max[i] + eps) - (xyz_min[i] - eps);
 		}
 		if (m > mass) {
 			mass = m;
 			largest = it.Value();
 		}
 	}
 }
@@ -0,0 +1,28 @@
 #include <gp_Pln.hxx>
 #include <TopoDS_Wire.hxx>
 #include <TopTools_ListOfShape.hxx>
 #include <vector>
 namespace IfcGeom {
 	namespace util {
 		bool approximate_plane_through_wire(const TopoDS_Wire& wire, gp_Pln& plane, double eps);
 		bool flatten_wire(TopoDS_Wire& wire, double eps);
 		enum triangulate_wire_result {
 			TRIANGULATE_WIRE_FAIL,
 			TRIANGULATE_WIRE_OK,
 			TRIANGULATE_WIRE_NON_MANIFOLD,
 		};
 		/// Triangulate the set of wires. The firstmost wire is assumed to be the outer wire.
 		triangulate_wire_result triangulate_wire(const std::vector<TopoDS_Wire>& wires, TopTools_ListOfShape& faces);
 		// eps: tolerance added to wire intersection checks, can be zero
 		// eps_real: tolerance used to construct new edge geometry around intersection points, cannot be zero
 		bool wire_intersections(const TopoDS_Wire& wire, TopTools_ListOfShape& wires, double eps, double eps_real);
 		void select_largest(const TopTools_ListOfShape& shapes, TopoDS_Shape& largest);
 	}
 }