diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index 978a148a06..d516c51879 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -29,15 +29,154 @@ using namespace ifcopenshell::geometry; #include #include -// @todo use std::numbers::pi when upgrading to C++ 20 -static const double PI = boost::math::constants::pi(); - namespace { // @todo: rb is there a common math library these functions can be moved to? auto sign = [](double v) -> int { return v < 0 ? -1 : 1; }; // returns -1 or 1 auto binary_sign = [](double v) -> int { return v < 0 ? -1 : (0 < v ? 1 : 0); }; // returns -1, 0, or 1 + +// @todo change the calculation at end of this to std::lerp when upgrading to C++ 20 +auto interpolate = [](double u, double a, double b, double l) -> double { return l == 0.0 ? 0.0 : u * (b - a) / l; }; } // namespace + +// @todo use std::numbers::pi when upgrading to C++ 20 +static const double PI = boost::math::constants::pi(); + +// Current implementation uses the same segment_geometry_adjuster for all ParentCurve types. +// Comment/Uncomment to change the type of segment geometry adjuster +// Future implementations could use specialized adjusters based on ParentCurve type +//#define GEOMETRY_ADJUSTER segment_geometry_adjuster +#define GEOMETRY_ADJUSTER linear_segment_geometry_adjuster + +// Curve segments are evaluated using a parametric function over the curve length, u +// IfcCurveSegment.TransitionCode defines how the end of a segment connects to the next segment. +// When segments are continuously joined, the placement at u = length should be equal to the placement at u = 0 +// of the next segment. However, numerical errors can cause these two points to be slightly offset +// from one another (the tangents could be slightly different as well). +// +// The sources of these numerical errors include geometric approximations (series expansion versus integration +// for spiral curves), the IfcCurveSegment.SegmentStart or .SegmentLength parameters contain roundoff or +// truncation error, minor errors in placement at the start of a segment can magnify error at the end +// of the segment. There are probably others as well. +// +// The evaluation of the relative location of the end and start points of adjacent segments occurs +// after the IfcCurveSegment.Placement is applied to the ParentCurve. The ParentCurve can be defined in +// a convenient coordinate system, such as the center of a circle or the origin of a line at (0,0). The Placement +// them moves the computed geometry to its relative position. It is the geometry after applying the Placement +// that needs to be evaluated and any difference forms the bases for the adjustments made by segment_geometry_adjuster +// or one of its subclasses. +// +// This class applies the IfcCurveSegment.Placement to inst_. The placement at the start of next_inst_ can then be +// obtained from mapping->map and compared to the end placement of inst_ and the placement at u can be adjusted +// as needed. This default implementation doesn't make any adjustments. Subclass and override the transform_and_adjust +// function to specialize the refinement of the placement at u. +class segment_geometry_adjuster { + public: + segment_geometry_adjuster(mapping* mapping, const IfcSchema::IfcCurveSegment* inst, const IfcSchema::IfcCurveSegment* next_inst) : + end_of_inst_(Eigen::Matrix4d::Identity()), + start_of_next_inst_(Eigen::Matrix4d::Identity()), + transition_code_(inst->Transition()) + { + + transformation_matrix_ = taxonomy::cast(mapping->map(inst->Placement()))->ccomponents(); + length_ = fabs(*inst->SegmentLength()->as() * mapping->get_length_unit()); + + if (next_inst) { + // if there is a next segment, get the coordinates at the start. + // Note that mapping->map(next_inst) causes mapping to occur recursively + // through all of the curve segments until the end of curve is reached. + // Mapping of IfcCompositeCurve, IfcGradientCurve, and IfcSegmentedReferenceCurve may + // need to traverse the IfcCurveSegment objects in reverse order to avoid recursion. + auto next = taxonomy::cast(mapping->map(next_inst)); + start_of_next_inst_ = next->evaluate(0.0); + } + } + + // To determine the geometry adjustments the curve segment needs to be evaluated + // without adjustments. This function toggles the application of geometry adjustments + void enable_adjustments(bool adjustments) { adjustments_ = adjustments; } + + // This object doesn't have access to the eval_ property of the curve_segment_evaluator. + // The end point of the segment being adjusted, without adjustments, is computed externally + // and provided to the curve_segment_adjustor through this method + void set_segment_end_point(const Eigen::Matrix4d& end_of_inst) { + end_of_inst_ = end_of_inst; + } + + // Transforms the ParentCurve geometry with the IfcCurveSegment.Placement and + // applies geometric adjustments to the geometry, if enabled + Eigen::Matrix4d transform_and_adjust(double u, const Eigen::Matrix4d& parent_curve_point) const { + // transform the parent curve's value into the segment curve's coordinate system + Eigen::Matrix4d segment_curve_point = transformation_matrix_ * parent_curve_point; + if (adjustments_) { + apply_adjustments(u, segment_curve_point); + } + return segment_curve_point; + } + + protected: + // Applies geometric adjustment to the segment curve point evaluated at u + // This default implementation does nothing + virtual void apply_adjustments(double u, Eigen::Matrix4d& p) const { /* do nothing - override in subclass if needed */ } + + const Eigen::Matrix4d& get_end_of_segment() const { return end_of_inst_; } + const Eigen::Matrix4d& get_start_of_next_segment() const { return start_of_next_inst_; } + IfcSchema::IfcTransitionCode::Value get_transition_code() const { return transition_code_; } + double get_length() const { return length_; } + + private: + bool adjustments_ = true; + Eigen::Matrix4d transformation_matrix_; + Eigen::Matrix4d end_of_inst_; + Eigen::Matrix4d start_of_next_inst_; + double length_; + IfcSchema::IfcTransitionCode::Value transition_code_; +}; + +// This class refines the geometric adjustment along the segment by dividing the +// difference between the segment end point and the start point of the next segment +// into equal adjustments and applying the incremental adjustment to each position at u +class linear_segment_geometry_adjuster : public segment_geometry_adjuster { + public: + using segment_geometry_adjuster::segment_geometry_adjuster; + + virtual void apply_adjustments(double u, Eigen::Matrix4d& p) const override { + // make the adjustments based on the transition code + // all segments must connect end to end except for last segment IfcTransitionCode_DISCONTINUOUS for open curve + auto transition_code = get_transition_code(); + if (transition_code == IfcSchema::IfcTransitionCode::IfcTransitionCode_DISCONTINUOUS) + return; + + const auto& end_this = get_end_of_segment(); + const auto& start_next = get_start_of_next_segment(); + auto xe = end_this.col(3)(0); + auto ye = end_this.col(3)(1); + auto xs = start_next.col(3)(0); + auto ys = start_next.col(3)(1); + auto length = get_length(); + auto x = interpolate(u,xe,xs,length); + auto y = interpolate(u,ye,ys,length); + + p.col(3)(0) += x; + p.col(3)(1) += y; + + if (transition_code == IfcSchema::IfcTransitionCode::IfcTransitionCode_CONTSAMEGRADIENT or + transition_code == IfcSchema::IfcTransitionCode::IfcTransitionCode_CONTSAMEGRADIENTSAMECURVATURE) { + + for (int i = 0; i < 2; i++) { + auto dxe = end_this.col(i)(0); + auto dye = end_this.col(i)(1); + auto dxs = start_next.col(i)(0); + auto dys = start_next.col(i)(1); + auto dx = interpolate(u,dxe,dxs,length); + auto dy = interpolate(u,dye,dys,length); + p.col(i)(0) += dx; + p.col(i)(1) += dy; + } + } + } +}; + typedef boost::mpl::vector< IfcSchema::IfcLine #ifdef SCHEMA_HAS_IfcClothoid @@ -58,33 +197,66 @@ enum segment_type_t { class curve_segment_evaluator { private: mapping* mapping_; + const IfcSchema::IfcCurveSegment* inst_; + const IfcSchema::IfcCurveSegment* next_inst_; double length_unit_; double start_; double length_; segment_type_t segment_type_; - IfcSchema::IfcCurve* curve_; + const IfcSchema::IfcCurve* curve_; + + std::shared_ptr geometry_adjuster; std::optional> eval_; public: // First constructor, takes parameters from IfcCurveSegment - curve_segment_evaluator(mapping* mapping, double length_unit, segment_type_t segment_type, IfcSchema::IfcCurve* curve, IfcSchema::IfcCurveMeasureSelect* st, IfcSchema::IfcCurveMeasureSelect* le) + curve_segment_evaluator(mapping* mapping, const IfcSchema::IfcCurveSegment* inst, const IfcSchema::IfcCurveSegment* next_inst, double length_unit, segment_type_t segment_type) : mapping_(mapping), + inst_(inst), + next_inst_(next_inst), length_unit_(length_unit), segment_type_(segment_type), - curve_(curve) { + curve_(inst->ParentCurve()) { // @todo in IFC4X3_ADD2 this needs to be length measure - if (!st->as() || !le->as()) { + + if (!inst->SegmentStart()->as() || !inst->SegmentLength()->as()) { // @nb Parameter values are forbidden in the specification until parametrization is provided for all spirals throw std::runtime_error("Unsupported curve measure type"); } - start_ = *st->as() * length_unit; - length_ = *le->as() * length_unit; + start_ = *inst->SegmentStart()->as() * length_unit; + length_ = *inst->SegmentLength()->as() * length_unit; } - void set_spiral_functor(mapping* mapping_, IfcSchema::IfcSpiral* c, double s, std::function signX, std::function fnX, std::function signY, std::function fnY) { + void compute_segment_end_point() + { + // The segment_geometry_adjuster needs to have both the end point of this segment + // and the start point of the next segment. The start point of the next + // segment is easy to get and is handled by the segment_geometry_adjuster. + // The end point of this segment must be computed by calling the eval_ callback + // at u = length_. But things are a little more complicated than that. eval_ will + // use segment_geometry_adjuster to correct deviations between this segment's end point and + // the next segments start point. In order to compute those adjustments, the + // end point of this segment, without correction, must be known. The end point not known + // at this time because segment_geometry_adjuster doesn't have access to the eval_ callback. + // Additionally, the eval_ callback needs to know if it is evaluating the segment geometry + // with our without geometric adjustments. + // + // Solving that conundrum is the purpose of this function. The geometric adjustments + // of geometry_adjuster are disabled, eval_ is called to get the unadjusted end point + // of this segment, the geometry_adjuster is updated with the end point so it can + // compute and apply geometry adjustments. + if (eval_) { + geometry_adjuster->enable_adjustments(false); // disable adjustments + auto end_point = (*eval_)(fabs(length_)); // compute the end point without correction + geometry_adjuster->set_segment_end_point(end_point); // save the unadjusted end point it can be used to compute adjustments + geometry_adjuster->enable_adjustments(true); // enable adjustments + } + } + + void set_spiral_function(mapping* mapping_, const IfcSchema::IfcSpiral* c, double s, std::function signX, std::function fnX, std::function signY, std::function fnY) { // determine the length of the spiral from the local origin to the end point auto sign_s = binary_sign(start_); auto sign_l = binary_sign(length_); @@ -101,10 +273,11 @@ class curve_segment_evaluator { auto start = start_; + geometry_adjuster = std::make_shared(mapping_, inst_, next_inst_); + if (segment_type_ == ST_HORIZONTAL || segment_type_ == ST_VERTICAL) { auto segment_type = segment_type_; - eval_ = [L, start, s, signX, fnX, signY, fnY, transformation_matrix, segment_type](double u) { - using boost::math::quadrature::trapezoidal; + eval_ = [L, start, s, signX, fnX, signY, fnY, transformation_matrix, segment_type, geometry_adjuster = this->geometry_adjuster](double u) { u += start; @@ -112,6 +285,7 @@ class curve_segment_evaluator { auto a = 0.0; auto b = fabs(u / s); + using boost::math::quadrature::trapezoidal; auto x = signX(u) * trapezoidal(fnX, a, b); auto y = signY(u) * trapezoidal(fnY, a, b); @@ -142,19 +316,19 @@ class curve_segment_evaluator { m.col(3) = Eigen::Vector4d(0, 0, y, 1.0); // y is an elevation so store it as z } Eigen::Matrix4d result = transformation_matrix * m; - return result; + return geometry_adjuster->transform_and_adjust(u,result); }; } else if (segment_type_ == ST_CANT) { - eval_ = [](double u) { + eval_ = [geometry_adjuster = this->geometry_adjuster](double u) { Eigen::Matrix4d result; - return result; + return geometry_adjuster->transform_and_adjust(u, result); }; } else { Logger::Error(std::runtime_error("Unexpected segment type encountered")); } -} + } // Clothoid using Taylor Series approximation @@ -204,7 +378,7 @@ class curve_segment_evaluator { // Clothoid using numerical integration #ifdef SCHEMA_HAS_IfcClothoid // Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes - void operator()(IfcSchema::IfcClothoid* c) { + void operator()(const IfcSchema::IfcClothoid* c) { // see https://standards.buildingsmart.org/IFC/RELEASE/IFC4_3/HTML/lexical/IfcClothoid.htm // also see, https://standards.buildingsmart.org/IFC/RELEASE/IFC4_3/HTML/concepts/Partial_Templates/Geometry/Curve_Segment_Geometry/Clothoid_Transition_Segment/content.html, // which defines the clothoid constant as sqrt(L) and L is the length measured from the inflection point @@ -223,12 +397,12 @@ class curve_segment_evaluator { auto fn_x = [A,s](double t)->double {return s * cos(PI * fabs(A) * t * t / (2 * fabs(A))); }; auto fn_y = [A,s](double t)->double {return s * sin(PI * fabs(A) * t * t / (2 * fabs(A))); }; - set_spiral_functor(mapping_, c, s, sign_x, fn_x, sign_y, fn_y); + set_spiral_function(mapping_, c, s, sign_x, fn_x, sign_y, fn_y); } #endif #ifdef SCHEMA_HAS_IfcSecondOrderPolynomialSpiral - void operator()(IfcSchema::IfcSecondOrderPolynomialSpiral* c) + void operator()(const IfcSchema::IfcSecondOrderPolynomialSpiral* c) { // @todo: rb verify - this is an example implementation of a different kind of spiral - lots of clean up needed auto A0 = c->ConstantTerm(); @@ -250,23 +424,24 @@ class curve_segment_evaluator { auto fn_y = [theta](double t)->double {return sin(theta(t)); }; double s = 1.0; // @todo: rb - this is supposed to be the curve length when the parametric value u = 1.0 - set_spiral_functor(mapping_, c, s, sign_x, fn_x, sign_y, fn_y); + set_spiral_function(mapping_, c, s, sign_x, fn_x, sign_y, fn_y); } #endif - void operator()(IfcSchema::IfcCircle* c) + void operator()(const IfcSchema::IfcCircle* c) { auto R = c->Radius(); auto sign_l = sign(length_); auto start = start_; - //const auto& transformation_matrix = taxonomy::cast(mapping_->map(c->Position()))->ccomponents(); auto transformation_matrix = taxonomy::cast(mapping_->map(c->Position()))->ccomponents(); auto segment_type = segment_type_; - eval_ = [R, start, sign_l, transformation_matrix, segment_type](double u) + geometry_adjuster = std::make_shared(mapping_, inst_, next_inst_); + + eval_ = [R, start, sign_l, transformation_matrix, segment_type, geometry_adjuster = this->geometry_adjuster](double u) { auto angle = start + sign_l * u / R; @@ -297,11 +472,11 @@ class curve_segment_evaluator { Eigen::Matrix4d result = transformation_matrix * m; - return result; + return geometry_adjuster->transform_and_adjust(u, result); }; } - void operator()(IfcSchema::IfcPolyline* pl) + void operator()(const IfcSchema::IfcPolyline* pl) { struct Range { @@ -321,7 +496,7 @@ class curve_segment_evaluator { } auto std_compare = [](double u_start, double u, double u_end) {return u_start <= u && u < u_end; }; - auto end_compare = [](double u_start, double u, double u_end) {return u_start <= u && u <= (u_end + 0.001); }; + auto end_compare = [](double u_start, double u, double u_end) { return u_start <= u && u <= (u_end + 0.001); }; auto begin = p->begin(); auto iter = begin; @@ -391,7 +566,10 @@ class curve_segment_evaluator { u = u + l; } - eval_ = [fns](double u) { + + geometry_adjuster = std::make_shared(mapping_, inst_, next_inst_); + + eval_ = [fns, geometry_adjuster = this->geometry_adjuster](double u) { auto iter = std::find_if(fns.cbegin(), fns.cend(), [=](const auto& fn) { auto [u_start, u_end, compare] = fn.first; @@ -403,11 +581,11 @@ class curve_segment_evaluator { const auto& [u_start, u_end, compare] = iter->first; const auto& fn = iter->second; Eigen::Matrix4d m = fn(u - u_start); // (u - u_start) is distance from start of this segment of the polyline - return m; + return geometry_adjuster->transform_and_adjust(u, m); }; } - void operator()(IfcSchema::IfcLine* l) { + void operator()(const IfcSchema::IfcLine* l) { auto s = l->Pnt(); auto c = s->Coordinates(); auto v = l->Dir(); @@ -418,9 +596,10 @@ class curve_segment_evaluator { auto dx = dr[0] / m; auto dy = dr[1] / m; - if (segment_type_ == ST_HORIZONTAL) { + geometry_adjuster = std::make_shared(mapping_, inst_, next_inst_); + if (segment_type_ == ST_HORIZONTAL) { - eval_ = [px, py, dx, dy](double u) { + eval_ = [px, py, dx, dy, geometry_adjuster=this->geometry_adjuster](double u) { auto x = px + u * dx; auto y = py + u * dy; @@ -429,20 +608,18 @@ class curve_segment_evaluator { m.col(1) = Eigen::Vector4d(-dy, dx, 0, 0); // vector perpendicular to the curve, towards the left when looking from start to end along the curve (this is used for IfcAxis2PlacementLinear.RefDirection when it is not provided) m.col(2) = Eigen::Vector4d(0, 0, 1.0, 0); // cross product of x and y and will always be up (this is used for IfcAxis2PlacementLinear.Axis when it is not provided) m.col(3) = Eigen::Vector4d(x, y, 0.0, 1.0); - return m; - }; - + return geometry_adjuster->transform_and_adjust(u, m); + }; } else if (segment_type_ == ST_VERTICAL || segment_type_ == ST_CANT) { - - eval_ = [px, py, dx, dy](double u) { + eval_ = [py, dx, dy, geometry_adjuster = this->geometry_adjuster](double u) { // https://standards.buildingsmart.org/IFC/RELEASE/IFC4_3/HTML/lexical/IfcGradientCurve.htm // the parameter, u, is the parameter of the BaseCurve (u = plan view distance along base curve) // dx and dy are normalized so u needs to be scaled by dy/dx // Consider a 5% uphill grade defined by dr[0] = 1 and dr[1] = 0.05. // We would normally compute y = py + 0.05*u. - // However, m = sqrt(1*1 + 0.05*.0.05) = 1.0124922 we need to normalize the direction ratios as + // However, m = sqrt(1*1 + 0.05*0.05) = 1.0124922 we need to normalize the direction ratios as // dx = dr[0]/m and dy = dr[1]/m which makes dy = 0.05/1.0124922 = 0.0499376 // y = py + u * dy/dx = py + u * (dr[1]/m)*(m/dr[0]) = py + u * 0.05 auto y = py + u * dy/dx; @@ -452,31 +629,34 @@ class curve_segment_evaluator { m.col(1) = Eigen::Vector4d(0, 1, 0, 0); m.col(2) = Eigen::Vector4d(-dy, 0, dx, 0); m.col(3) = Eigen::Vector4d(0, 0, y, 1.0); // y is an elevation so store it as z - return m; - }; + return geometry_adjuster->transform_and_adjust(u, m); + }; } else { Logger::Error(std::runtime_error("Unexpected segment type encountered"), l); } - } + } - void operator()(IfcSchema::IfcPolynomialCurve* p) { + void operator()(const IfcSchema::IfcPolynomialCurve* p) { // see https://forums.buildingsmart.org/t/ifcpolynomialcurve-clarification/4716 for discussion on IfcPolynomialCurve auto coeffX = p->CoefficientsX().get_value_or(std::vector()); auto coeffY = p->CoefficientsY().get_value_or(std::vector()); auto coeffZ = p->CoefficientsZ().get_value_or(std::vector()); if (!coeffZ.empty()) - Logger::Warning("Expected IfcPolynomialCurve.CoefficientsZ to be undefined for alignment geometry", p); + Logger::Warning("Expected IfcPolynomialCurve.CoefficientsZ to be undefined for alignment geometry. Coefficients ignored.", p); auto transformation_matrix = taxonomy::cast(mapping_->map(p->Position()))->ccomponents(); auto segment_type = segment_type_; - eval_ = [coeffX, coeffY, coeffZ,transformation_matrix,segment_type](double u) { - std::array*, 3> coefficients{&coeffX, &coeffY, &coeffZ}; - std::array position{0.0, 0.0, 0.0}; // @todo: rb, use Eigen::VectorXd - I'm sure there is a way to do this with Eigen, but this is what I know - std::array slope{0.0, 0.0, 0.0}; // slope is derivative of the curve = SUM( coeff*pos*u^(pos-1) ) - for (int i = 0; i < 3; i++) { + geometry_adjuster = std::make_shared(mapping_, inst_, next_inst_); + + + eval_ = [coeffX, coeffY, transformation_matrix, segment_type, geometry_adjuster = this->geometry_adjuster](double u) { + std::array*, 2> coefficients{&coeffX, &coeffY}; + std::array position{0.0, 0.0}; + std::array slope{0.0, 0.0}; // slope is derivative of the curve = SUM( coeff*pos*u^(pos-1) ) + for (int i = 0; i < 2; i++) { auto begin = coefficients[i]->cbegin(); auto end = coefficients[i]->cend(); for (auto iter = begin; iter != end; iter++) { @@ -516,7 +696,7 @@ class curve_segment_evaluator { Logger::Error(std::runtime_error("Unexpected segment type encountered")); } - return m; + return geometry_adjuster->transform_and_adjust(u, m); }; } @@ -538,7 +718,27 @@ class curve_segment_evaluator { }; taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { - // @todo: rb figure out what to do with the zero length segments at the end of compound curves + // Find the next segment after inst + const IfcSchema::IfcCurveSegment* next_inst = nullptr; + auto composite_curves = inst->data().getInverse(&IfcSchema::IfcCompositeCurve::Class(), 0); + if (composite_curves) { + if (composite_curves->size() == 1) { + auto segments = (*composite_curves->begin())->as()->Segments(); + bool emit_next = false; + for (auto& s : *segments) { + if (emit_next) { + next_inst = s->as(); + break; + } + if (s == inst) { + emit_next = true; + } + } + } + else { + Logger::Warning("IfcCurveSegment belongs to multiple IfcCompositeCurve instances. Cannot determine the next segment. Geometry adjustments will not be made."); + } + } bool is_horizontal = false; bool is_vertical = false; @@ -569,25 +769,18 @@ 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,length_unit_, segment_type, inst->ParentCurve(), inst->SegmentStart(), inst->SegmentLength()); + curve_segment_evaluator cse(this, inst, next_inst, length_unit_, segment_type); boost::mpl::for_each>(std::ref(cse)); + cse.compute_segment_end_point(); auto& eval_fn = cse.evaluation_function(); if(!eval_fn) throw std::runtime_error(inst->ParentCurve()->declaration().name() + " not implemented"); auto fn = *eval_fn; auto length = fabs(cse.length()); - auto transformation_matrix = taxonomy::cast(map(inst->Placement()))->ccomponents(); - - auto fn_transformed = [fn, transformation_matrix](double u)->Eigen::Matrix4d { - Eigen::Matrix4d f = fn(u); - Eigen::Matrix4d result = transformation_matrix * f; - return result; - }; - // @todo it might be suboptimal that we no longer have the spans now auto pwf = taxonomy::make(); - pwf->spans.push_back({ length, fn_transformed }); + pwf->spans.push_back({ length, fn }); pwf->instance = inst; return pwf; }