SectionedSolidHorizontal: sort sections, bounds check

src/ifcgeom/mapping/IfcSectionedSolidHorizontal.cpp

@@ -33,13 +33,32 @@ namespace {
 	T lerp(const T& a, const T& b, double t) {
 		return a + t * (b - a);
 	}
+
+	struct cross_section {
+		double dist_along;
+		taxonomy::face::ptr section_geometry;
+		Eigen::Vector3d offset;
+
+		bool operator <(const cross_section& other) const {
+			return dist_along < other.dist_along;
+		}
+	};
 }
 
 taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSolidHorizontal* inst) {
+	std::vector<cross_section> cross_sections;
+
 	auto dir = map(inst->Directrix());
+	auto pwf = taxonomy::dcast<taxonomy::piecewise_function>(dir);
+	if (!pwf) {
+		// Only implement on alignment curves
+		return nullptr;
+	}
+
+	{	
 	auto css = inst->CrossSections();
 	auto csps = inst->CrossSectionPositions();
-	std::vector<taxonomy::face::ptr> cross_sections;
+	std::vector<taxonomy::face::ptr> faces;
 
 	// The PointByDistanceExpressesions are factored out into (a) a cartesian offset relative to the
 	// reference frame along a certain curve location (b) the longitude.
@@ -49,19 +68,13 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSolidHorizontal* in
 	std::vector<Eigen::Vector3d> profile_offsets;
 	std::vector<double> longitudes;
 
-	auto pwf = taxonomy::dcast<taxonomy::piecewise_function>(dir);
-	if (!pwf) {
-		// Only implement on alignment curves
-		return nullptr;
-	}
-
 	for (auto& cs : *css) {
-		cross_sections.push_back(std::move(taxonomy::cast<taxonomy::face>(map(cs))));
+		faces.push_back(std::move(taxonomy::cast<taxonomy::face>(map(cs))));
 	}
 #ifdef SCHEMA_HAS_IfcPointByDistanceExpression
 	for (auto& csp : *csps) {
 		auto pbde = csp->Location()->as<IfcSchema::IfcPointByDistanceExpression>(true);
-		
+
 		longitudes.push_back(*pbde->DistanceAlong()->as<IfcSchema::IfcLengthMeasure>(true) * length_unit_);
 
 		// Corresponds to the profile X, Y directions (hopefully).
@@ -77,23 +90,35 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSolidHorizontal* in
 #else
 	return nullptr;
 #endif
-	if (cross_sections.size() != profile_offsets.size()) {
-		Logger::Error("Expected CrossSections and CrossSectionPositions to be equal length, but got " + std::to_string(cross_sections.size()) + " and " + std::to_string(profile_offsets.size())  + " respectively", inst);
+	if (faces.size() != profile_offsets.size()) {
+		Logger::Warning("Expected CrossSections and CrossSectionPositions to be equal length, but got " + std::to_string(faces.size()) + " and " + std::to_string(profile_offsets.size()) + " respectively", inst);
 		return nullptr;
 	}
-	if (cross_sections.size() < 2) {
-		Logger::Error("Expected at least two cross sections, but got " + std::to_string(cross_sections.size()), inst);
+	if (faces.size() < 2) {
+		Logger::Warning("Expected at least two cross sections, but got " + std::to_string(faces.size()), inst);
 		return nullptr;
 	}
 
+	for (size_t i = 0; i < faces.size(); ++i) {
+		cross_sections.push_back({ longitudes[i], faces[i], profile_offsets[i] });
+	}
+	}
+
+	std::sort(cross_sections.begin(), cross_sections.end());
+
 	auto loft = taxonomy::make<taxonomy::loft>();
 	// @todo intialize as default
 	loft->axis = nullptr;
 
 	// @todo currently only the case is handled where directrix returns a piecewise_function
 	if (pwf) {
-		double start = std::max(0., longitudes.front());
-		double end = std::min(pwf->length(), longitudes.back());
+		double start = std::max(0., cross_sections.front().dist_along);
+		double end = std::min(pwf->length(), cross_sections.back().dist_along);
+
+		if (end - start < 1.e-9) {
+			Logger::Warning("Empty sweep domain with start at " + std::to_string(cross_sections.front().dist_along) + " end at " + std::to_string(cross_sections.back().dist_along) + " and curve domain length " + std::to_string(pwf->length()), inst);
+			return nullptr;
+		}
 
 		auto curve_length = end - start;
 		auto param_type = settings_.get<ifcopenshell::geometry::settings::PiecewiseStepType>().get();
@@ -106,6 +131,10 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSolidHorizontal* in
 			// parameter is minimum number of steps
 			num_steps = (size_t) std::ceil(param);
 		}
+		std::vector<double> longitudes;
+		for (auto& x : cross_sections) {
+			longitudes.push_back(x.dist_along);
+		}
 		longitudes.push_back(std::numeric_limits<double>::infinity());
 		auto profile_index = longitudes.begin();
 		for (size_t i = 0; i <= num_steps; ++i) {
@@ -118,8 +147,8 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSolidHorizontal* in
 			}
 
 			auto relative_dist_along = (dist_along - *profile_index) / (*(profile_index+1) - *profile_index);
-			const auto& profile_a = cross_sections[std::distance(longitudes.begin(), profile_index)];
-			const auto& offset_a = profile_offsets[std::distance(longitudes.begin(), profile_index)];
+			const auto& profile_a = cross_sections[std::distance(longitudes.begin(), profile_index)].section_geometry;
+			const auto& offset_a = cross_sections[std::distance(longitudes.begin(), profile_index)].offset;
 
 			taxonomy::face::ptr interpolated = nullptr;
 
@@ -134,8 +163,8 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSolidHorizontal* in
 				taxonomy::face::ptr profile_b;
 				Eigen::Vector3d offset_b;
 				if ((profile_index + 1 < longitudes.end())) {
-					profile_b = cross_sections[std::distance(longitudes.begin(), profile_index) + 1];
-					offset_b = profile_offsets[std::distance(longitudes.begin(), profile_index) + 1];
+					profile_b = cross_sections[std::distance(longitudes.begin(), profile_index) + 1].section_geometry;
+					offset_b = cross_sections[std::distance(longitudes.begin(), profile_index) + 1].offset;
 				} else {
 					profile_b = profile_a;
 					offset_b = offset_a;
@@ -179,6 +208,9 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcSectionedSolidHorizontal* in
 			}
 
 			auto m4 = pwf->evaluate(dist_along);
+			/* {
+				std::wcout << "#" << pwf->instance->data().id() << " " << dist_along << ": " << m4.col(3).row(2).value() << std::endl;
+			}*/
 
 			Eigen::Matrix4d m4b = Eigen::Matrix4d::Identity();
 			m4b.col(0).head<3>() = m4.col(1).head<3>().normalized();