diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index fe5a47cf2a..cec4f36beb 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -29,15 +29,168 @@ 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 +template +auto compute_adjustment = [](double u, const T& a, const T& 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; + init_adjustments(); + } + + // 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: + // precompute any values that are constant when applying geometry adjustments + //( subclasses to override. + virtual void init_adjustments() { /*do nothing*/ + } + // 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; + + protected: + virtual void init_adjustments() override { + // @todo: rb - implement to improve efficiency + // cache delta = (start_next - end_this)/length + // adjustment is then adj = u*delta + } + + 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 = compute_adjustment(u, xe, xs, length); + auto y = compute_adjustment(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 = compute_adjustment(u,dxe,dxs,length); + auto dy = compute_adjustment(u,dye,dys,length); + p.col(i)(0) += dx; + p.col(i)(1) += dy; + p.col(i).normalize(); + } + } + } +}; + typedef boost::mpl::vector< IfcSchema::IfcLine #ifdef SCHEMA_HAS_IfcClothoid @@ -56,98 +209,140 @@ enum segment_type_t { }; class curve_segment_evaluator { -private: - mapping* mapping_; - double length_unit_; - double start_; - double length_; - segment_type_t segment_type_; - IfcSchema::IfcCurve* curve_; + private: + mapping* mapping_; + const IfcSchema::IfcCurveSegment* inst_; + const IfcSchema::IfcCurveSegment* next_inst_; + double length_unit_; + double start_; + double length_; + segment_type_t segment_type_; + const IfcSchema::IfcCurve* curve_; - std::optional> eval_; + std::shared_ptr geometry_adjuster; -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) - : mapping_(mapping) - , length_unit_(length_unit) - , segment_type_(segment_type) - , curve_(curve) - { - // @todo in IFC4X3_ADD2 this needs to be length measure + std::optional> eval_; - if (!st->as() || !le->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; - } - - void set_spiral_functor(mapping* mapping_,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_); - double L = 0; - if (sign_s == 0) L = fabs(length_); // start_ is at zero so length_ is the L - else if (sign_s == sign_l) L = fabs(start_ + length_); // start_ and length_ are additive - else L = fabs(start_); // start_ and length_ are in opposite directions so start_ is furthest from the origin - - auto transformation_matrix = taxonomy::cast(mapping_->map(c->Position()))->ccomponents(); - - auto segment_type = segment_type_; - auto start = start_; - - eval_ = [L, start, s, signX, fnX, signY, fnY, transformation_matrix, segment_type](double u) { - using boost::math::quadrature::trapezoidal; - - u += start; - - // integration limits, integrate from a to b - auto a = 0.0; - auto b = fabs(u / s); - - auto x = signX(u) * trapezoidal(fnX, a, b); - auto y = signY(u) * trapezoidal(fnY, a, b); - - // From https://standards.buildingsmart.org/IFC/RELEASE/IFC4_3/HTML/lexical/IfcSpiral.htm, x = Integral(fnX du), y = Integral(fnY du) - // The tangent slope of a curve is the derivate of the curve, so the derivitive of an integral, is just the function - // Therefore, Dx/Du = fnX(u) and Dy/Du = fnY(u) which leads to du = Dx/fnX(u) and Dy = fnY(u)*Du = fnY(u)*Dx/fnX(u) so Dy/Dx = fnY(u)/fnX(u) - // However, Dx and Dy are not normalized. Recall that slope = rise/run - // If run = 1.0, then rise = Dy/Dx = fnY(u)/fnX(u) and l = sqrt((fnY(u)/fnX(u))^2 + 1.0^2) - // The direction ratios are dx = 1.0/l and dy = (fnY/fnX)/l; - auto rise = fnY(u) / fnX(u); - auto run = 1.0; - auto l = sqrt(run * run + rise * rise); - auto dx = run / l; - auto dy = rise / l; - - Eigen::Matrix4d m; - if (segment_type == ST_HORIZONTAL) { - // rotate about the Z-axis - m.col(0) = Eigen::Vector4d(dx, dy, 0, 0); // vector tangent to the curve, in the direction of the curve - 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); - } else if (segment_type == ST_VERTICAL) { - // rotate about the Y-axis (slope along u is dx, slope vertically is dy, vertical position is y) - m.col(0) = Eigen::Vector4d(dx, 0, dy, 0); - 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 - } else if (segment_type == ST_CANT) { - Logger::Warning(std::runtime_error("Use of IfcSpiral for cant is not supported")); - } else { - Logger::Error(std::runtime_error("Unexpected segment type encountered")); - } + public: + // First constructor, takes parameters from IfcCurveSegment + 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_(inst->ParentCurve()) { + // @todo in IFC4X3_ADD2 this needs to be length measure - Eigen::Matrix4d result = transformation_matrix * m; - return result; - }; - } + 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_ = *inst->SegmentStart()->as() * length_unit; + length_ = *inst->SegmentLength()->as() * length_unit; + } + + 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_); + double L = 0; + if (sign_s == 0) { + L = fabs(length_); // start_ is at zero so length_ is the L + } else if (sign_s == sign_l) { + L = fabs(start_ + length_); // start_ and length_ are additive + } else { + L = fabs(start_); // start_ and length_ are in opposite directions so start_ is furthest from the origin + } + + auto transformation_matrix = taxonomy::cast(mapping_->map(c->Position()))->ccomponents(); + + 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, geometry_adjuster = this->geometry_adjuster](double u) { + + u += start; + + // integration limits, integrate from a to b + 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); + + // From https://standards.buildingsmart.org/IFC/RELEASE/IFC4_3/HTML/lexical/IfcSpiral.htm, x = Integral(fnX du), y = Integral(fnY du) + // The tangent slope of a curve is the derivate of the curve, so the derivitive of an integral, is just the function + // Therefore, Dx/Du = fnX(u) and Dy/Du = fnY(u) which leads to du = Dx/fnX(u) and Dy = fnY(u)*Du = fnY(u)*Dx/fnX(u) so Dy/Dx = fnY(u)/fnX(u) + // However, Dx and Dy are not normalized. Recall that slope = rise/run + // If run = 1.0, then rise = Dy/Dx = fnY(u)/fnX(u) and l = sqrt((fnY(u)/fnX(u))^2 + 1.0^2) + // The direction ratios are dx = 1.0/l and dy = (fnY/fnX)/l; + auto rise = fnY(u) / fnX(u); + auto run = 1.0; + auto l = sqrt(run * run + rise * rise); + auto dx = run / l; + auto dy = rise / l; + + Eigen::Matrix4d m; + if (segment_type == ST_HORIZONTAL) { + // rotate about the Z-axis + m.col(0) = Eigen::Vector4d(dx, dy, 0, 0); // vector tangent to the curve, in the direction of the curve + 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); + } else if (segment_type == ST_VERTICAL) { + // rotate about the Y-axis (slope along u is dx, slope vertically is dy, vertical position is y) + m.col(0) = Eigen::Vector4d(dx, 0, dy, 0); + 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 + } + Eigen::Matrix4d result = transformation_matrix * m; + return geometry_adjuster->transform_and_adjust(u,result); + }; + } + else if (segment_type_ == ST_CANT) { + eval_ = [geometry_adjuster = this->geometry_adjuster](double u) { + Eigen::Matrix4d result; + return geometry_adjuster->transform_and_adjust(u, result); + }; + } + else { + Logger::Error(std::runtime_error("Unexpected segment type encountered")); + } + } // Clothoid using Taylor Series approximation @@ -197,7 +392,7 @@ public: // 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 @@ -216,12 +411,12 @@ public: 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(); @@ -243,23 +438,24 @@ public: 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; @@ -290,11 +486,11 @@ public: 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 { @@ -314,7 +510,7 @@ public: } 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; @@ -384,7 +580,10 @@ public: 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; @@ -396,11 +595,11 @@ public: 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(); @@ -411,9 +610,10 @@ public: 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; @@ -422,20 +622,19 @@ public: 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) { + 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; @@ -445,34 +644,34 @@ public: 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 if(segment_type_ == ST_CANT) { - Logger::Warning(std::runtime_error("Use of IfcLine for cant is not supported"), l); - } 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++) { @@ -487,11 +686,9 @@ public: auto x = position[0]; auto y = position[1]; - //auto z = position[2]; auto dx = slope[0]; auto dy = slope[1]; - //auto dz = slope[2]; Eigen::Matrix4d m; if (segment_type == ST_HORIZONTAL) { @@ -512,7 +709,7 @@ public: Logger::Error(std::runtime_error("Unexpected segment type encountered")); } - return m; + return geometry_adjuster->transform_and_adjust(u, m); }; } @@ -534,7 +731,27 @@ public: }; 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; @@ -565,25 +782,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; } diff --git a/src/ifcgeom/mapping/IfcGradientCurve.cpp b/src/ifcgeom/mapping/IfcGradientCurve.cpp index 85fa6f10af..a01d42f488 100644 --- a/src/ifcgeom/mapping/IfcGradientCurve.cpp +++ b/src/ifcgeom/mapping/IfcGradientCurve.cpp @@ -24,7 +24,10 @@ using namespace ifcopenshell::geometry; #ifdef SCHEMA_HAS_IfcGradientCurve taxonomy::ptr mapping::map_impl(const IfcSchema::IfcGradientCurve* inst) { - auto horizontal = taxonomy::cast(map(inst->BaseCurve())); + if (!inst->BaseCurve()->as()) + Logger::Warning("Expected IfcGradientCurve.BaseCurve to be IfcCompositeCurve", inst); // CT 4.1.7.1.1.2 + + auto horizontal = taxonomy::cast(map(inst->BaseCurve())); auto vertical = taxonomy::make(); auto segments = inst->Segments(); diff --git a/src/ifcgeom/taxonomy.cpp b/src/ifcgeom/taxonomy.cpp index bc04e9391e..45b0b1ad9d 100644 --- a/src/ifcgeom/taxonomy.cpp +++ b/src/ifcgeom/taxonomy.cpp @@ -472,6 +472,16 @@ ifcopenshell::geometry::taxonomy::item::ptr ifcopenshell::geometry::taxonomy::pi return polygon_from_points(polygon); } +Eigen::Matrix4d ifcopenshell::geometry::taxonomy::piecewise_function::evaluate(double u) const { + // @todo: rb optimize, assume monotonic evaluation and store last evaluated segment? + for (auto& [length, fn] : spans) { + if (u < length + 0.001) { // @todo: rb - need to use consistent tolerance + return fn(u); + } + u -= length; + } +} + ifcopenshell::geometry::taxonomy::collection::ptr ifcopenshell::geometry::flatten(taxonomy::collection::ptr deep) { auto flat = make(); ifcopenshell::geometry::visit(deep, [&flat](taxonomy::ptr i) { diff --git a/src/ifcgeom/taxonomy.h b/src/ifcgeom/taxonomy.h index 2fccbe760c..8d1358351c 100644 --- a/src/ifcgeom/taxonomy.h +++ b/src/ifcgeom/taxonomy.h @@ -151,15 +151,7 @@ typedef item const* ptr; virtual item::ptr evaluate() const; - Eigen::Matrix4d evaluate(double u) const { - // @todo: rb optimize, assume monotonic evaluation and store last evaluated segment? - for (auto& [length, fn] : spans) { - if (u < length+0.001) { // @todo: rb - need to use consistent tolerance - return fn(u); - } - u -= length; - } - } + Eigen::Matrix4d evaluate(double u) const; }; #ifdef TAXONOMY_USE_SHARED_PTR