Spatial Querying: Implementation of BVH tree for IFC files using NCollection_UBTree (#130)

* Implementation of BVH tree for IFC files using NCollection_UBTree

* Fixes for gcc and OCCT < v6.8

* Enable BVH selection by TopoDS_Shape

* Small fixes to tree

* Merge conflicted files and other changes

* Eliminate one warning (some remain) + add tree test

src/ifcgeom/IfcGeom.h

@@ -205,7 +205,7 @@ public:
 	bool closest(const gp_Pnt&, const std::vector<gp_Pnt>&, gp_Pnt&);
 	bool project(const Handle_Geom_Curve&, const gp_Pnt&, gp_Pnt& p, double& u, double& d);
 	bool project(const Handle_Geom_Surface&, const TopoDS_Shape&, double& u1, double& v1, double& u2, double& v2, double widen=0.1);
-	int count(const TopoDS_Shape&, TopAbs_ShapeEnum);
+	static int count(const TopoDS_Shape&, TopAbs_ShapeEnum);
 
 	bool find_wall_end_points(const IfcSchema::IfcWall*, gp_Pnt& start, gp_Pnt& end);
 

src/ifcgeom/IfcGeomIterator.h

@@ -433,7 +433,7 @@ namespace IfcGeom {
 											IfcSchema::IfcProductRepresentation::list::ptr prodreps = representation_mapped_to->OfProductRepresentation();
 
 											bool all_product_without_openings = true;
-											IfcSchema::IfcProduct::list::ptr products;
+											IfcSchema::IfcProduct::list::ptr products(new IfcSchema::IfcProduct::list);
 
 											for (IfcSchema::IfcProductRepresentation::list::it it = prodreps->begin(); it != prodreps->end(); ++it) {
 												IfcSchema::IfcProduct::list::ptr products_of_prodrep = (*it)->entity->getInverse(IfcSchema::Type::IfcProduct, -1)->as<IfcSchema::IfcProduct>();

src/ifcgeom/IfcGeomRepresentation.cpp

@@ -31,13 +31,8 @@ IfcGeom::Representation::Serialization::Serialization(const BRep& brep)
 	: Representation(brep.settings())
 	, id_(brep.id())
 {
-	TopoDS_Compound compound;
-	BRep_Builder builder;
-	builder.MakeCompound(compound);
+	TopoDS_Compound compound = brep.as_compound();
 	for (IfcGeom::IfcRepresentationShapeItems::const_iterator it = brep.begin(); it != brep.end(); ++ it) {
-		const TopoDS_Shape& s = it->Shape();
-		gp_GTrsf trsf = it->Placement();
-		
 		if (it->hasStyle() && it->Style().Diffuse()) {
 			const IfcGeom::SurfaceStyle::ColorComponent& clr = *it->Style().Diffuse();
 			surface_styles_.push_back(clr.R());
@@ -53,18 +48,28 @@ IfcGeom::Representation::Serialization::Serialization(const BRep& brep)
 		} else {
 			surface_styles_.push_back(1.);
 		}
-		
+	}
+	std::stringstream sstream;
+	BRepTools::Write(compound,sstream);
+	brep_data_ = sstream.str();
+}
+
+TopoDS_Compound IfcGeom::Representation::BRep::as_compound() const {
+	TopoDS_Compound compound;
+	BRep_Builder builder;
+	builder.MakeCompound(compound);
+	for (IfcGeom::IfcRepresentationShapeItems::const_iterator it = begin(); it != end(); ++it) {
+		const TopoDS_Shape& s = it->Shape();
+		gp_GTrsf trsf = it->Placement();
+
 		if (settings().get(IteratorSettings::CONVERT_BACK_UNITS)) {
 			gp_Trsf scale;
 			scale.SetScaleFactor(1.0 / settings().unit_magnitude());
 			trsf.PreMultiply(scale);
 		}
-		
-		const TopoDS_Shape moved_shape = IfcGeom::Kernel::apply_transformation(s, trsf);
 
+		const TopoDS_Shape moved_shape = IfcGeom::Kernel::apply_transformation(s, trsf);
 		builder.Add(compound, moved_shape);
 	}
-	std::stringstream sstream;
-	BRepTools::Write(compound,sstream);
-	brep_data_ = sstream.str();
+	return compound;
 }

src/ifcgeom/IfcGeomRepresentation.h

@@ -38,6 +38,8 @@
 #include "../ifcgeom/IfcGeomMaterial.h"
 #include "../ifcgeom/IfcRepresentationShapeItem.h"
 
+#include <TopoDS_Compound.hxx>
+
 namespace IfcGeom {
 
 	namespace Representation {
@@ -72,6 +74,7 @@ namespace IfcGeom {
 			IfcGeom::IfcRepresentationShapeItems::const_iterator end() const { return shapes_.end(); }
 			const IfcGeom::IfcRepresentationShapeItems& shapes() const { return shapes_; }
 			const std::string& id() const { return id_; }
+			TopoDS_Compound as_compound() const;
 		};
 
 		class IFC_GEOM_API Serialization : public Representation  {

src/ifcgeom/IfcGeomTree.h

@@ -0,0 +1,278 @@
+/********************************************************************************
+*                                                                              *
+* This file is part of IfcOpenShell.                                           *
+*                                                                              *
+* IfcOpenShell is free software: you can redistribute it and/or modify         *
+* it under the terms of the Lesser GNU General Public License as published by  *
+* the Free Software Foundation, either version 3.0 of the License, or          *
+* (at your option) any later version.                                          *
+*                                                                              *
+* IfcOpenShell is distributed in the hope that it will be useful,              *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of               *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the                 *
+* Lesser GNU General Public License for more details.                          *
+*                                                                              *
+* You should have received a copy of the Lesser GNU General Public License     *
+* along with this program. If not, see <http://www.gnu.org/licenses/>.         *
+*                                                                              *
+********************************************************************************/
+
+#ifndef IFCGEOMTREE_H
+#define IFCGEOMTREE_H
+
+#include "../ifcparse/IfcFile.h"
+#include "../ifcgeom/IfcGeomIterator.h"
+
+#include <NCollection_UBTree.hxx>
+#include <BRepBndLib.hxx>
+#include <Bnd_Box.hxx>
+#include <BRepAlgoAPI_Common.hxx>
+#include <BRepAlgoAPI_Cut.hxx>
+#include <BRepClass3d_SolidClassifier.hxx>
+
+namespace IfcGeom {
+
+	namespace impl {
+		template <typename T>
+		class tree {
+
+		public:
+
+			void add(const T& t, const Bnd_Box& b) {
+				tree_.Add(t, b);
+			}
+
+			void add(const T& t, const TopoDS_Shape& s) {
+				Bnd_Box b;
+				BRepBndLib::AddClose(s, b);
+				add(t, b);
+				shapes_[t] = s;
+			}
+
+			std::vector<T> select_box(const T& t, bool completely_within = false, double extend=-1.e-5) const {
+				typename map_t::const_iterator it = shapes_.find(t);
+				if (it == shapes_.end()) {
+					return std::vector<T>();
+				}
+
+				Bnd_Box b;
+				BRepBndLib::AddClose(it->second, b);
+
+				// Gap is assumed to be positive throughout the codebase,
+				// but at least for IsOut() in the selector a negative
+				// Gap should work as well.
+				b.SetGap(b.GetGap() + extend); 
+				
+				return select_box(b, completely_within);
+			}
+
+			std::vector<T> select_box(const gp_Pnt& p) const {
+				Bnd_Box b;
+				b.Add(p);
+				return select_box(b);
+			}
+
+			std::vector<T> select_box(const Bnd_Box& b, bool completely_within = false) const {
+				selector s(b);
+				tree_.Select(s);
+				if (completely_within) {
+					std::vector<T> ts = s.results();
+					std::vector<T> ts_filtered;
+					ts_filtered.reserve(ts.size());
+					typename std::vector<T>::const_iterator it = ts.begin();
+					for (; it != ts.end(); ++it) {
+						const TopoDS_Shape& shp = shapes_.find(*it)->second;
+						Bnd_Box B;
+						BRepBndLib::AddClose(shp, B);
+
+						// BndBox::CornerMin() /-Max() introduced in OCCT 6.8
+						double x1, y1, z1, x2, y2, z2;
+						b.Get(x1, y1, z1, x2, y2, z2);
+						double gap = B.GetGap();
+						gp_Pnt p1(x1 - gap, y1 - gap, z1 - gap);
+						gp_Pnt p2(x2 + gap, y2 + gap, z2 + gap);
+						
+						if (!b.IsOut(p1) && !b.IsOut(p2)) {
+							ts_filtered.push_back(*it);
+						}
+					}
+					return ts_filtered;
+				} else {
+					return s.results();
+				}
+			}
+
+			std::vector<T> select(const T& t, bool completely_within = false) const {
+				std::vector<T> ts = select_box(t);
+				if (ts.empty()) {
+					return ts;
+				}
+
+				std::vector<T> ts_filtered;
+
+				const TopoDS_Shape& A = shapes_.find(t)->second;
+				if (IfcGeom::Kernel::count(A, TopAbs_SHELL) == 0) {
+					return ts_filtered;
+				}
+
+				ts_filtered.reserve(ts.size());
+
+				typename std::vector<T>::const_iterator it = ts.begin();
+				for (it = ts.begin(); it != ts.end(); ++it) {
+					const TopoDS_Shape& B = shapes_.find(*it)->second;
+					if (IfcGeom::Kernel::count(B, TopAbs_SHELL) == 0) {
+						continue;
+					}
+
+					if (completely_within) {
+						BRepAlgoAPI_Cut cut(B, A);
+						if (cut.IsDone()) {
+							if (IfcGeom::Kernel::count(cut.Shape(), TopAbs_SHELL) == 0) {
+								ts_filtered.push_back(*it);
+							}
+						}
+					} else {
+						BRepAlgoAPI_Common common(A, B);
+						if (common.IsDone()) {
+							if (IfcGeom::Kernel::count(common.Shape(), TopAbs_SHELL) > 0) {
+								ts_filtered.push_back(*it);
+							}
+						}
+					}
+				}
+
+				return ts_filtered;
+			}
+
+			std::vector<T> select(const TopoDS_Shape& s) const {
+				Bnd_Box bb;
+				BRepBndLib::AddClose(s, bb);
+
+				std::vector<T> ts;
+
+				if (IfcGeom::Kernel::count(s, TopAbs_SHELL) == 0) {
+					return ts;
+				}
+
+				ts = select_box(bb);
+
+				if (ts.empty()) {
+					return ts;
+				}
+
+				std::vector<T> ts_filtered;
+				ts_filtered.reserve(ts.size());
+
+				typename std::vector<T>::const_iterator it = ts.begin();
+				for (it = ts.begin(); it != ts.end(); ++it) {
+					const TopoDS_Shape& B = shapes_.find(*it)->second;
+					
+					if (IfcGeom::Kernel::count(B, TopAbs_SHELL) == 0) {
+						continue;
+					}
+
+					BRepAlgoAPI_Common common(s, B);
+					if (common.IsDone()) {
+						if (IfcGeom::Kernel::count(common.Shape(), TopAbs_SHELL) > 0) {
+							ts_filtered.push_back(*it);
+						}
+					}
+				}
+
+				return ts_filtered;
+			}
+
+			std::vector<T> select(const gp_Pnt& p) const {
+				std::vector<T> ts = select_box(p);
+				if (ts.empty()) {
+					return ts;
+				}
+
+				std::vector<T> ts_filtered;
+				ts_filtered.reserve(ts.size());
+
+				typename std::vector<T>::const_iterator it = ts.begin();
+				for (it = ts.begin(); it != ts.end(); ++it) {
+					const TopoDS_Shape& B = shapes_.find(*it)->second;
+					TopExp_Explorer exp(B, TopAbs_SOLID);
+					for (; exp.More(); exp.Next()) {
+						BRepClass3d_SolidClassifier cls(exp.Current(), p, 1e-5);
+						if (cls.State() != TopAbs_OUT) {
+							ts_filtered.push_back(*it);
+							break;
+						}
+					}					
+				}
+
+				return ts_filtered;
+			}
+
+		protected:
+
+			typedef NCollection_UBTree<T, Bnd_Box> tree_t;
+			typedef std::map<T, TopoDS_Shape> map_t;
+			tree_t tree_;
+			map_t shapes_;
+
+			class selector : public tree_t::Selector
+			{
+			public:
+				selector(const Bnd_Box& b)
+					: tree_t::Selector()
+					, bounds_(b)
+				{}
+
+				Standard_Boolean Reject(const Bnd_Box& b) const {
+					return bounds_.IsOut(b);
+				}
+
+				Standard_Boolean Accept(const T& o) {
+					results_.push_back(o);
+					return Standard_True;
+				}
+
+				const std::vector<T>& results() const {
+					return results_;
+				}
+
+			private:
+				std::vector<T> results_;
+				const Bnd_Box& bounds_;
+			};
+
+		};
+	}
+
+	class tree : public impl::tree<IfcSchema::IfcProduct*> {
+	public:
+
+		tree() {};
+
+		tree(IfcParse::IfcFile& f) {
+			add_file(f, IfcGeom::IteratorSettings());
+		}
+
+		tree(IfcParse::IfcFile& f, const IfcGeom::IteratorSettings& settings) {
+			add_file(f, settings);
+		}
+
+		void add_file(IfcParse::IfcFile& f, const IfcGeom::IteratorSettings& settings) {
+			IfcGeom::IteratorSettings settings_ = settings;
+			settings_.set(IfcGeom::IteratorSettings::DISABLE_TRIANGULATION, true);
+			settings_.set(IfcGeom::IteratorSettings::USE_WORLD_COORDS, true);
+			settings_.set(IfcGeom::IteratorSettings::SEW_SHELLS, true);
+
+			IfcGeom::Iterator<double> it(settings_, &f);
+
+			if (it.initialize()) {
+				do {
+					IfcGeom::BRepElement<double>* elem = (IfcGeom::BRepElement<double>*)it.get();
+					add((IfcSchema::IfcProduct*)f.entityById(elem->id()), elem->geometry().as_compound());
+				} while (it.next());
+			}
+		}
+	};
+
+}
+
+#endif

src/ifcopenshell-python/ifcopenshell/geom/main.py

@@ -72,6 +72,50 @@ class iterator(_iterator):
             return wrap_shape_creation(self.settings, _iterator.get(self))
 
 
+class tree(ifcopenshell_wrapper.tree):
+    
+    def __init__(self, file=None, settings=None):
+        args = [self]
+        if file is not None: 
+            args.append(file.wrapped_data)
+            if settings is not None:
+                args.append(settings)
+        ifcopenshell_wrapper.tree.__init__(*args)
+        
+    def add_file(self, file, settings):
+        ifcopenshell_wrapper.tree.add_file(self, file.wrapped_data, settings)
+        
+    def select(self, value, **kwargs):
+        def unwrap(value):
+            if isinstance(value, entity_instance):
+                return value.wrapped_data
+            elif all(map(lambda v: hasattr(value, v), "XYZ")):
+                return value.X(), value.Y(), value.Z()
+            return value
+        args = [self, unwrap(value)]
+        if isinstance(value, entity_instance):
+            args.append(kwargs.get("completely_within", False))
+        elif has_occ:
+            import OCC.TopoDS
+            if isinstance(value, OCC.TopoDS.TopoDS_Shape):
+                args[1] = utils.serialize_shape(value)
+        return [entity_instance(e) for e in ifcopenshell_wrapper.tree.select(*args)]
+        
+    def select_box(self, value, **kwargs):
+        def unwrap(value):
+            if isinstance(value, entity_instance):
+                return value.wrapped_data
+            elif hasattr(value, "Get"):
+                return value.Get()[:3], value.Get()[3:]
+            return value
+        args = [self, unwrap(value)]
+        if "extend" in kwargs or "completely_within" in kwargs:
+            args.append(kwargs.get("completely_within", False))
+        if "extend" in kwargs:
+            args.append(kwargs.get("extend", -1.e-5))
+        return [entity_instance(e) for e in ifcopenshell_wrapper.tree.select_box(*args)]
+
+
 def create_shape(settings, inst, repr=None): 
     return wrap_shape_creation(
         settings,

src/ifcwrap/IfcGeomWrapper.i

@@ -36,6 +36,8 @@
 	$result = PyBool_FromLong(static_cast<long>(*$1));
 }
 
+%ignore IfcGeom::impl::tree::selector;
+
 %include "../ifcgeom/ifc_geom_api.h"
 %include "../ifcgeom/IfcGeomIteratorSettings.h"
 %include "../ifcgeom/IfcGeomElement.h"
@@ -43,6 +45,63 @@
 %include "../ifcgeom/IfcGeomRepresentation.h"
 %include "../ifcgeom/IfcGeomIterator.h"
 
+// A Template instantantation should be defined before it is used as a base class. 
+// But frankly I don't care as most methods are subtlely different anyway.
+%include "../ifcgeom/IfcGeomTree.h"
+
+%extend IfcGeom::tree {
+
+	static IfcEntityList::ptr vector_to_list(const std::vector<IfcSchema::IfcProduct*>& ps) const {
+		IfcEntityList::ptr r(new IfcEntityList);
+		for (std::vector<IfcSchema::IfcProduct*>::const_iterator it = ps.begin(); it != ps.end(); ++it) {
+			r->push(*it);
+		}
+		return r;
+	}
+
+	IfcEntityList::ptr select_box(IfcUtil::IfcBaseClass* e, bool completely_within = false, double extend=-1.e-5) const {
+		if (!e->is(IfcSchema::Type::IfcProduct)) {
+			throw IfcParse::IfcException("Instance should be an IfcProduct");
+		}
+		std::vector<IfcSchema::IfcProduct*> ps = $self->select_box((IfcSchema::IfcProduct*)e, completely_within, extend);
+		return IfcGeom_tree_vector_to_list(ps);
+	}
+
+	IfcEntityList::ptr select_box(const gp_Pnt& p) const {
+		std::vector<IfcSchema::IfcProduct*> ps = $self->select_box(p);
+		return IfcGeom_tree_vector_to_list(ps);
+	}
+
+	IfcEntityList::ptr select_box(const Bnd_Box& b, bool completely_within = false) const {
+		std::vector<IfcSchema::IfcProduct*> ps = $self->select_box(b, completely_within);
+		return IfcGeom_tree_vector_to_list(ps);
+	}
+
+	IfcEntityList::ptr select(IfcUtil::IfcBaseClass* e, bool completely_within = false) const {
+		if (!e->is(IfcSchema::Type::IfcProduct)) {
+			throw IfcParse::IfcException("Instance should be an IfcProduct");
+		}
+		std::vector<IfcSchema::IfcProduct*> ps = $self->select((IfcSchema::IfcProduct*)e, completely_within);
+		return IfcGeom_tree_vector_to_list(ps);
+	}
+
+	IfcEntityList::ptr select(const gp_Pnt& p) const {
+		std::vector<IfcSchema::IfcProduct*> ps = $self->select(p);
+		return IfcGeom_tree_vector_to_list(ps);
+	}
+
+	IfcEntityList::ptr select(const std::string& shape_serialization) const {
+		std::stringstream stream(shape_serialization);
+		BRepTools_ShapeSet shapes;
+		shapes.Read(stream);
+		const TopoDS_Shape& shp = shapes.Shape(shapes.NbShapes());
+
+		std::vector<IfcSchema::IfcProduct*> ps = $self->select(shp);
+		return IfcGeom_tree_vector_to_list(ps);
+	}
+
+}
+
 // Using RTTI return a more specialized type of Element
 // Note that these elements are not to be owned by SWIG/Python as they will be freed automatically upon the next iteration
 // except for the IfcGeom::Element instances which are returned by Iterator::getObject() calls

src/ifcwrap/IfcPython.i

@@ -72,6 +72,7 @@
 %module ifcopenshell_wrapper %{
 	#include "../ifcgeom/IfcGeom.h"
 	#include "../ifcgeom/IfcGeomIterator.h"
+	#include "../ifcgeom/IfcGeomTree.h"
 
 	#include "../ifcparse/IfcBaseClass.h"
 	#include "../ifcparse/IfcFile.h"

src/ifcwrap/utils/typemaps_in.i

@@ -108,4 +108,39 @@ CREATE_VECTOR_TYPEMAP_IN(std::string, STRING, str)
 	} else {
 		SWIG_exception(SWIG_TypeError, "Attribute of type AGGREGATE OF AGGREGATE OF ENTITY INSTANCE needs a python sequence of sequence of entity instances");
 	}
+}
+
+%typemap(in) const gp_Pnt& {
+	if (!check_aggregate_of_type($input, get_python_type<double>())) {
+		SWIG_exception(SWIG_TypeError, "<Point> type needs a python sequence of 3 floats");
+	}
+	std::vector<double> ds = python_sequence_as_vector<double>($input);
+	if (ds.size() != 3) {
+		SWIG_exception(SWIG_TypeError, "<Point> type needs a python sequence of 3 floats");
+	}
+	$1 = new gp_Pnt(ds[0], ds[1], ds[2]);
+}
+
+%typemap(in) const Bnd_Box& {
+	if (!check_aggregate_of_aggregate_of_type($input, get_python_type<double>())) {
+		SWIG_exception(SWIG_TypeError, "<AABB> type needs a python sequence of 2 x 3 floats");
+	}
+	std::vector< std::vector<double> > ds = python_sequence_as_vector_of_vector<double>($input);
+	if (ds.size() != 2) {
+		SWIG_exception(SWIG_TypeError, "<AABB> type needs a python sequence of 2 x 3 floats");
+	}
+	if (ds[0].size() != 3 || ds[1].size() != 3) {
+		SWIG_exception(SWIG_TypeError, "<AABB> type needs a python sequence of 2 x 3 floats");
+	}
+	$1 = new Bnd_Box();
+	$1->Add(gp_Pnt(ds[0][0], ds[0][1], ds[0][2]));
+	$1->Add(gp_Pnt(ds[1][0], ds[1][1], ds[1][2]));
+}
+
+%typemap(typecheck, precedence=SWIG_TYPECHECK_INTEGER) const gp_Pnt& {
+   $1 = check_aggregate_of_type($input, get_python_type<double>());
+}
+
+%typemap(typecheck, precedence=SWIG_TYPECHECK_INTEGER) const Bnd_Box& {
+   $1 = check_aggregate_of_aggregate_of_type($input, get_python_type<double>());
 }

test/tests.py

@@ -23,6 +23,7 @@ import os
 import uuid
 
 import ifcopenshell
+import ifcopenshell.geom
 import ifcopenshell.guid
 
 f = ifcopenshell.open("input/acad2010_walls.ifc")
@@ -91,6 +92,13 @@ assert f[288].ConnectedTo == (rel,)
 # Some operations on ifcopenshell.guid
 assert len(ifcopenshell.guid.compress(uuid.uuid1().hex)) == 22
 
+# Test the BVH tree
+tree_settings = ifcopenshell.geom.settings()
+tree_settings.set(tree_settings.DISABLE_OPENING_SUBTRACTIONS, True)
+t = ifcopenshell.geom.tree(f, tree_settings)
+# This wall is connected to two other walls
+assert len(t.select_box(f[48], extend=0.1)) == 3
+
 # Test serialization
 f.write("output.ifc")
 with open("output.ifc") as txt: