Updates taxonomy::piecewise_function so it's immutable and implements internal caching. Adds new evaluate2 method that returns loop and vector of distance along for each loop point.

src/ifcgeom/mapping/IfcCompositeCurve.cpp

@@ -23,7 +23,7 @@ using namespace ifcopenshell::geometry;
 
 taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) {
 	auto loop = taxonomy::make<taxonomy::loop>();
-	auto pwf = taxonomy::make<taxonomy::piecewise_function>(&settings_);
+	std::vector<taxonomy::piecewise_function::ptr> pwfs;
 
 #ifdef SCHEMA_HAS_IfcSegment
 	// 4x3
@@ -31,8 +31,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) {
 #else
 	IfcSchema::IfcCompositeCurveSegment::list::ptr segments = inst->Segments();
 #endif
-	current_segment_count_ = segments->size();
+	
-
 	for (auto& segment : *segments) {
 		if (segment->as<IfcSchema::IfcCompositeCurveSegment>() && segment->as<IfcSchema::IfcCompositeCurveSegment>()->ParentCurve()->as<IfcSchema::IfcLine>()) {
 			Logger::Notice("Infinite IfcLine used as ParentCurve of segment, treating as a segment", segment);
@@ -78,8 +77,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) {
 				}
 			}
 			else if (crv && crv->kind() == taxonomy::PIECEWISE_FUNCTION) {
-				auto seg = taxonomy::cast<taxonomy::piecewise_function>(crv);
+            pwfs.push_back(taxonomy::cast<taxonomy::piecewise_function>(crv));
-				pwf->spans.insert(pwf->spans.end(), seg->spans.begin(), seg->spans.end());
 			} else if (!crv) {
 				return nullptr;
 			}
@@ -87,7 +85,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) {
 #endif
 	}
 
-	if (pwf->spans.empty()) {
+	if (pwfs.empty()) {
 		aggregate_of_instance::ptr profile = inst->data().getInverse(&IfcSchema::IfcProfileDef::Class(), -1);
 		const bool force_close = profile && profile->size() > 0;
 		loop->closed = force_close;
@@ -95,7 +93,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) {
 		return loop;
 	}
 	else {
-		pwf->instance = inst;
+      auto pwf = taxonomy::make<taxonomy::piecewise_function>(pwfs,&settings_,inst);
 		return pwf;
 	}
 }

src/ifcgeom/mapping/IfcCurveSegment.cpp

@@ -103,7 +103,6 @@ class curve_segment_evaluator {
     double length_; // length along the curve, as provided from the IfcCurveSegment
     segment_type_t segment_type_;
     const IfcSchema::IfcCurve* parent_curve_ = nullptr;
-    size_t current_segment_count_;
 
     double projected_length_; // for vertical segments, this is the length of curve projected onto the "Distance Along" axis
 
@@ -111,13 +110,12 @@ class curve_segment_evaluator {
     std::optional<Eigen::Matrix4d> parent_curve_placement_;                 // placement matrix for the parent curve
 
   public:
-    curve_segment_evaluator(mapping* mapping, const IfcSchema::IfcCurveSegment* inst, double length_unit, segment_type_t segment_type, size_t current_segment_count)
+    curve_segment_evaluator(mapping* mapping, const IfcSchema::IfcCurveSegment* inst, double length_unit, segment_type_t segment_type)
         : mapping_(mapping),
           inst_(inst),
           length_unit_(length_unit),
           segment_type_(segment_type),
-          parent_curve_(inst->ParentCurve()),
+          parent_curve_(inst->ParentCurve()) {
-          current_segment_count_(current_segment_count) {
         start_ = translate_if_param_value(inst->ParentCurve(), inst->SegmentStart()) * length_unit;
         length_ = translate_if_param_value(inst->ParentCurve(), inst->SegmentLength()) * length_unit;
     }
@@ -734,7 +732,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) {
 
     auto segment_type = is_horizontal ? ST_HORIZONTAL : is_vertical ? ST_VERTICAL : ST_CANT;
 
-    curve_segment_evaluator cse(this, inst, length_unit_, segment_type, current_segment_count_);
+    curve_segment_evaluator cse(this, inst, length_unit_, segment_type);
     boost::mpl::for_each<curve_seg_types, boost::type<boost::mpl::_>>(std::ref(cse));
     const auto& parent_curve_fn = cse.parent_curve_function();
     const auto& parent_curve_placement = cse.parent_curve_placement();
@@ -772,11 +770,11 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) {
         return curve_segment_placement * rotation * translation * p;
     };
 
-    // @todo it might be suboptimal that we no longer have the spans now
+    auto length = cse.length();
-    auto pwf = taxonomy::make<taxonomy::piecewise_function>(&settings_);
+
-    auto length = fabs(cse.length());
+    taxonomy::piecewise_function::spans spans;
-    pwf->spans.push_back({length, fn});
+    spans.emplace_back(fabs(length), fn);
-    pwf->instance = inst;
+    auto pwf = taxonomy::make<taxonomy::piecewise_function>(spans,&settings_,inst);
     return pwf;
 }
 

src/ifcgeom/mapping/IfcGradientCurve.cpp

@@ -28,18 +28,17 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcGradientCurve* inst) {
        Logger::Warning("Expected IfcGradientCurve.BaseCurve to be IfcCompositeCurve", inst); // CT 4.1.7.1.1.2
 
    auto horizontal = taxonomy::cast<taxonomy::piecewise_function>(map(inst->BaseCurve()));
-	auto vertical = taxonomy::make<taxonomy::piecewise_function>(&settings_);
 
 	auto segments = inst->Segments();
-	current_segment_count_ = segments->size();
+
+    std::vector<taxonomy::piecewise_function::ptr> pwfs;
 
 	for (auto& segment : *segments) {
 		if (segment->as<IfcSchema::IfcCurveSegment>()) {
 			// @todo check that we don't get a mixture of implicit and explicit definitions
 			auto crv = map(segment->as<IfcSchema::IfcCurveSegment>());
 			if (crv && crv->kind() == taxonomy::PIECEWISE_FUNCTION) {
-				auto seg = taxonomy::cast<taxonomy::piecewise_function>(crv);
+				pwfs.push_back(taxonomy::cast<taxonomy::piecewise_function>(crv));
-				vertical->spans.insert(vertical->spans.end(), seg->spans.begin(), seg->spans.end());
 			} else {
 				Logger::Error("Unsupported");
 				return nullptr;
@@ -49,6 +48,8 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcGradientCurve* inst) {
 			return nullptr;
 		}
 	}
+  
+	auto vertical = taxonomy::make<taxonomy::piecewise_function>(pwfs,&settings_);
 
 	auto composition = [horizontal, vertical](double u)->Eigen::Matrix4d {
 		auto xy = horizontal->evaluate(u);
@@ -63,26 +64,12 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcGradientCurve* inst) {
       return m;
 	};
 
-	std::array<taxonomy::piecewise_function::ptr, 2> both = { horizontal , vertical };
+   double min_length = std::min(horizontal->length(), vertical->length());
-	// @todo: rb - this constrains the range of u to the minimum of horizontal and vertical
-	// we discussed using the maximum for the range of us and then using std::numeric_limits<double>::NAN
-	// for values of u that horizontal or vertical cannot be computed.
-	// Review and decide what to do.
-	double min_length = std::numeric_limits<double>::infinity();
-	for (auto i = 0; i < 2; ++i) {
-		double l = 0;
-		for (auto& s : both[i]->spans) {
-			l += s.first;
-		}
-		if (l < min_length) {
-			min_length = l;
-		}
-	}
 
-	auto pwf = taxonomy::make<taxonomy::piecewise_function>(&settings_);
+	taxonomy::piecewise_function::spans spans;
-	pwf->spans.emplace_back( min_length, composition );
+   spans.emplace_back(min_length, composition);
-	pwf->instance = inst;
+   auto pwf = taxonomy::make<taxonomy::piecewise_function>(spans, &settings_, inst);
-	return pwf;
+   return pwf;
 }
 
 #endif

src/ifcgeom/mapping/IfcOffsetCurveByDistance.cpp

@@ -42,7 +42,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcOffsetCurveByDistances* inst
     auto pw_curve = taxonomy::dcast<taxonomy::piecewise_function>(map(basis_curve));
     double basis_curve_length = pw_curve->length();
 
-    auto offsets = taxonomy::make<taxonomy::piecewise_function>(&settings_);
+    taxonomy::piecewise_function::spans offset_spans;
 
 #if defined SCHEMA_HAS_IfcDistanceExpression
    double first_distance = first_offset_value->DistanceAlong();
@@ -63,12 +63,12 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcOffsetCurveByDistances* inst
         py *= length_unit_;
         pz *= length_unit_;
         
-        auto fn = [py, pz](double u) -> Eigen::Matrix4d { 
+        auto fn = [py, pz](double /*u*/) -> Eigen::Matrix4d { 
            Eigen::Matrix4d m = Eigen::Matrix4d::Identity(); 
            m.col(3)(1) = py; 
            m.col(3)(2) = pz; 
            return m; };
-        offsets->spans.push_back({first_distance, fn});
+        offset_spans.emplace_back(first_distance, fn);
 	}
 
     auto iter = offset_values->begin();
@@ -112,7 +112,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcOffsetCurveByDistances* inst
             m.col(3)(2) = (l == 0.0 ? zp : (zp + (zn - zp) * u / l));
             return m;
         };
-		  offsets->spans.push_back({l, fn});
+        offset_spans.emplace_back(l, fn);
     }
 
 	 // at this point, next == end and prev == end-1
@@ -134,15 +134,17 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcOffsetCurveByDistances* inst
          py *= length_unit_;
          pz *= length_unit_;
          double l = basis_curve_length - last_distance;
-         auto fn = [py, pz](double u) -> Eigen::Matrix4d { 
+         auto fn = [py, pz](double /*u*/) -> Eigen::Matrix4d { 
             Eigen::Matrix4d m = Eigen::Matrix4d::Identity(); 
             m.col(3)(1) = py; 
             m.col(3)(2) = pz; 
             return m; };
 
-         offsets->spans.push_back({l, fn});
+         offset_spans.emplace_back(l, fn);
     }
 
+   auto offsets = taxonomy::make<taxonomy::piecewise_function>(offset_spans,&settings_);
+
 	auto composition = [pw_curve, offsets](double u) -> Eigen::Matrix4d {
       auto p = pw_curve->evaluate(u);
 		auto offset = offsets->evaluate(u);
@@ -152,9 +154,9 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcOffsetCurveByDistances* inst
 
    // current implementation assumes that composition is equal to the full length of basis curve
    // this may change depending on decisions in the bSI-IF
-	auto pwf = taxonomy::make<taxonomy::piecewise_function>(&settings_);
+   taxonomy::piecewise_function::spans spans;
-	pwf->spans.emplace_back( basis_curve_length, composition );
+   spans.emplace_back(basis_curve_length, composition);
-	pwf->instance = inst;
+	auto pwf = taxonomy::make<taxonomy::piecewise_function>(spans,&settings_,inst);
 	return pwf;
 }
 

src/ifcgeom/mapping/IfcSegmentedReferenceCurve.cpp

@@ -28,18 +28,16 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSegmentedReferenceCurve* ins
        Logger::Warning("Expected IfcSegmentedReferenceCurve.BaseCurve to be IfcGradient", inst); // CT 4.1.7.1.1.3
 		  
 	auto gradient = taxonomy::cast<taxonomy::piecewise_function>(map(inst->BaseCurve()));
-	auto cant = taxonomy::make<taxonomy::piecewise_function>(&settings_);
+	
-
 	auto segments = inst->Segments();
-	current_segment_count_ = segments->size();
 
-	for (auto& segment : *segments) {
+    std::vector<taxonomy::piecewise_function::ptr> pwfs;
+    for (auto& segment : *segments) {
 		if (segment->as<IfcSchema::IfcCurveSegment>()) {
 			// @todo check that we don't get a mixture of implicit and explicit definitions
 			auto crv = map(segment->as<IfcSchema::IfcCurveSegment>());
 			if (crv && crv->kind() == taxonomy::PIECEWISE_FUNCTION) {
-				auto seg = taxonomy::cast<taxonomy::piecewise_function>(crv);
+				pwfs.push_back(taxonomy::cast<taxonomy::piecewise_function>(crv));
-				cant->spans.insert(cant->spans.end(), seg->spans.begin(), seg->spans.end());
 			} else {
 				Logger::Error("Unsupported");
 				return nullptr;
@@ -49,6 +47,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSegmentedReferenceCurve* ins
 			return nullptr;
 		}
 	}
+   auto cant = taxonomy::make<taxonomy::piecewise_function>(pwfs,&settings_);
 
 	auto composition = [gradient, cant](double u)->Eigen::Matrix4d {
 		auto g = gradient->evaluate(u);
@@ -63,26 +62,12 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSegmentedReferenceCurve* ins
       return m;
 	};
 
-	std::array<taxonomy::piecewise_function::ptr, 2> both = { gradient , cant };
+   double min_length = std::min(gradient->length(), cant->length());
-	// @todo: rb - this constrains the range of u to the minimum of gradient and cant
-	// we discussed using the maximum for the range of us and then using std::numeric_limits<double>::NAN
-	// for values of u that gradient or cant cannot be computed.
-	// Review and decide what to do.
-	double min_length = std::numeric_limits<double>::infinity();
-	for (auto i = 0; i < 2; ++i) {
-		double l = 0;
-		for (auto& s : both[i]->spans) {
-			l += s.first;
-		}
-		if (l < min_length) {
-			min_length = l;
-		}
-	}
 
-	auto pwf = taxonomy::make<taxonomy::piecewise_function>(&settings_);
+	taxonomy::piecewise_function::spans spans;
-	pwf->spans.emplace_back( min_length, composition );
+   spans.emplace_back(min_length, composition);
-	pwf->instance = inst;
+   auto pwf = taxonomy::make<taxonomy::piecewise_function>(spans, &settings_, inst);
-	return pwf;
+   return pwf;
 }
 
 #endif

src/ifcgeom/mapping/mapping.h

@@ -22,7 +22,6 @@ namespace geometry {
 		IfcParse::IfcFile* file_;
 		double length_unit_, angle_unit_;
 		std::string length_unit_name_;
-		size_t current_segment_count_;
 
 		std::map<uint32_t, ifcopenshell::geometry::taxonomy::ptr> cache_;
 

src/ifcgeom/taxonomy.cpp

@@ -471,37 +471,67 @@ ifcopenshell::geometry::taxonomy::item::ptr ifcopenshell::geometry::taxonomy::pi
    } else {
 	   // parameter is minimum number of steps
        num_steps = (unsigned)std::ceil(param);
-    }
+   }
 
 	return evaluate(0.0, curve_length, num_steps);
 }
 
 item::ptr ifcopenshell::geometry::taxonomy::piecewise_function::evaluate(double ustart, double uend,unsigned nsteps) const {
+    return evaluate2(ustart, uend, nsteps).first;
+}
+
+std::pair<item::ptr, std::vector<double>> ifcopenshell::geometry::taxonomy::piecewise_function::evaluate2(double ustart, double uend, unsigned nsteps) const {
     double curve_length = length();
     ustart = std::max(0.0, ustart);
     uend = std::min(uend, curve_length);
 
- 	 auto resolution = (uend - ustart)/ nsteps;
+    auto resolution = (uend - ustart) / nsteps;
 
     std::vector<taxonomy::point3::ptr> polygon;
+    std::vector<double> u_values;
+    polygon.reserve(nsteps);
+    u_values.reserve(nsteps);
 
-    for (int i = 0; i <= nsteps; ++i) {
+    for (unsigned i = 0; i <= nsteps; ++i) {
         auto u = resolution * i + ustart;
+        u_values.push_back(u);
         Eigen::Matrix4d m = evaluate(u);
         polygon.push_back(taxonomy::make<taxonomy::point3>(m.col(3)(0), m.col(3)(1), m.col(3)(2)));
     }
 
-    return polygon_from_points(polygon);
+    return {polygon_from_points(polygon), u_values};
 }
 
 Eigen::Matrix4d ifcopenshell::geometry::taxonomy::piecewise_function::evaluate(double u) const {
-   // @todo: rb optimize, assume monotonic evaluation and store last evaluated segment?
+	// assume monotonic evaluation and store last evaluated segment
-   for (auto& [length, fn] : spans) {
+    if (current_span_fn_ == nullptr || (u < current_span_start_ || current_span_end_ < u)) {
-   if (u < length + 0.001) { // @todo: rb - need to use consistent tolerance
+       // there isn't a current span or u is outside the range of the current span
-         return fn(u);
+		 // get a new "current span"
-      }
+        std::tie(current_span_start_,current_span_end_, current_span_fn_) = get_span(u);
-      u -= length;
+    }
-   }
+
+	 u -= current_span_start_; // make u relative to start of span
+	 return (*current_span_fn_)(u);
+}
+
+std::tuple<double, double, const std::function<Eigen::Matrix4d(double u)>*> ifcopenshell::geometry::taxonomy::piecewise_function::get_span(double u) const {
+	 // force u to be within bounds of the curve
+    double curve_length = length();
+    u = std::max(0.0, u);
+    u = std::min(u, curve_length);
+
+	double start = 0;
+    for (auto& [length, fn] : spans_) {
+		 auto tolerance = settings_ ? settings_->get<ifcopenshell::geometry::settings::Precision>().get() : 0.001;
+       if (u < length + tolerance) {
+           return {start, start+length, &fn} ;
+       }
+       start += length;
+       u -= length;
+    }
+
+    Logger::Error("taxonomy::piecewise_function::get_span span not found.");
+    return {0, 0, nullptr};
 }
 
 ifcopenshell::geometry::taxonomy::collection::ptr ifcopenshell::geometry::flatten(const taxonomy::collection::ptr& deep) {

src/ifcgeom/taxonomy.h

@@ -137,21 +137,32 @@ typedef item const* ptr;
 			};
 
 			struct piecewise_function : public implicit_item {
-				DECLARE_PTR(piecewise_function)
+            DECLARE_PTR(piecewise_function)
 
-				piecewise_function(const IfcUtil::IfcBaseInterface* instance = nullptr) : implicit_item(instance){};
+				using spans = std::vector<std::pair<double, std::function<Eigen::Matrix4d(double u)>>>;
-				piecewise_function(ifcopenshell::geometry::Settings* settings) : settings_(settings){};
+
+            piecewise_function(const spans& s, ifcopenshell::geometry::Settings* settings = nullptr, const IfcUtil::IfcBaseInterface* instance = nullptr) : 
+					implicit_item(instance), settings_(settings), spans_(s){};
+            piecewise_function(const std::vector<piecewise_function::ptr>& pwfs, ifcopenshell::geometry::Settings* settings = nullptr, const IfcUtil::IfcBaseInterface* instance = nullptr) : 
+					implicit_item(instance), settings_(settings)
+				{
+                for (auto& pwf : pwfs) {
+                    spans_.insert(spans_.end(), pwf->spans_.begin(), pwf->spans_.end());
+                }
+				};
             piecewise_function(piecewise_function&&) = default;
             piecewise_function(const piecewise_function&) = default;
 
-				ifcopenshell::geometry::Settings* settings_ = nullptr;
+				const ifcopenshell::geometry::Settings* settings_ = nullptr;
 
-				// length of span, function to evaluate span
+				bool is_empty() const { return spans_.empty(); }
-				std::vector<std::pair<double, std::function<Eigen::Matrix4d(double u)>>> spans;
 
 				double length() const {
-                    return std::accumulate(spans.begin(), spans.end(), 0.0, [](const auto& v,const auto& s) { return v + s.first; });
+					if (!length_.has_value()) {
-                }
+						length_ = std::accumulate(spans_.begin(), spans_.end(), 0.0, [](const auto& v, const auto& s) { return v + s.first; });
+					}
+               return *length_;
+            }
 
 				void print(std::ostream& o, int = 0) const {
 					o << "piecewise_function" << std::endl;
@@ -165,10 +176,34 @@ typedef item const* ptr;
 					return boost::hash<decltype(v)>{}(v);
 				}
 
-				virtual item::ptr evaluate() const;
+				item::ptr evaluate() const override;
-            item::ptr evaluate(double ustart, double uend,unsigned nsteps) const;
 
+            /// @brief evaluates the piecewise function between ustart and uend
+            /// @param ustart starting location - taken as 0.0 if before start
+            /// @param uend ending location - taken as length if beyond end
+            /// @param nsteps number of steps to evaluate
+            /// @return taxonomy::loop::ptr
+            item::ptr evaluate(double ustart, double uend, unsigned nsteps) const;
+
+            /// @brief evaluates the piecewise function between ustart and uend
+            /// @param ustart starting location - taken as 0.0 if before start
+            /// @param uend ending location - taken as length if beyond end
+            /// @param nsteps number of steps to evaluate
+            /// @return taxonomy::loop::ptr and vector of u values
+            std::pair<item::ptr, std::vector<double>> evaluate2(double ustart, double uend, unsigned nsteps) const;
+
+				/// @brief evaluates the piecewise function at u
+				/// @param u u is constrained to be between 0 and length
+				/// @return 4x4 placement matrix
 				Eigen::Matrix4d evaluate(double u) const;
+
+            private:
+               std::tuple<double, double, const std::function<Eigen::Matrix4d(double u)>*> get_span(double u) const;
+               spans spans_;
+               mutable double current_span_start_ = 0;
+				   mutable double current_span_end_ = 0;
+				   mutable const std::function<Eigen::Matrix4d(double u)>* current_span_fn_ = nullptr;
+               mutable boost::optional<double> length_;
 			};
 
 #ifdef TAXONOMY_USE_SHARED_PTR
@@ -1156,7 +1191,8 @@ typedef item const* ptr;
                         if (loop->pwf.is_initialized()) {
                             pwf_ = loop->pwf;
                         } else {
-                            pwf_ = taxonomy::make<taxonomy::piecewise_function>();
+                            taxonomy::piecewise_function::spans spans;
+                            spans.reserve(loop->children.size());
                             for (auto& edge : loop->children) {
                                 // the edge could be an arc or trimmed circle in the case of IfcIndexPolyCurve - support for this isn't implemented yet
                                 if (edge->basis) {
@@ -1174,8 +1210,9 @@ typedef item const* ptr;
                                     axis = refDirection.cross(Y).normalized();
                                     return taxonomy::make<taxonomy::matrix4>(o, axis, refDirection)->components();
                                 };
-                                (*pwf_)->spans.emplace_back(l, fn);
+                                spans.emplace_back(l, fn);
                             }
+                            pwf_ = taxonomy::make<taxonomy::piecewise_function>(spans);
                             loop->pwf = pwf_;
                         }
                     }