From 852311279d8826e7c4569f7332f7285b033008a8 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Mon, 25 May 2026 10:29:38 -0700 Subject: [PATCH] Simplifies line and circle parent curves and parent curve normalization --- src/ifcgeom/mapping/IfcCurveSegment.cpp | 157 +++++++++++++++--------- 1 file changed, 97 insertions(+), 60 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index 6dfd831e6d..8e6fac1259 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -133,23 +133,21 @@ struct spiral_parent_curve : public parent_curve_function { // this is the piecewise curve segment function for horizontal and vertical struct curve_segment_function { - curve_segment_function(const Eigen::Matrix4d& curve_segment_placement, const Eigen::Matrix4d& remove_parent_curve_rotation, const Eigen::Matrix4d& remove_parent_curve_translation, std::shared_ptr parent_curve_fn) : + curve_segment_function(const Eigen::Matrix4d& curve_segment_placement, const Eigen::Matrix4d& parent_curve_normalization, std::shared_ptr parent_curve_fn) : curve_segment_placement_(curve_segment_placement), - remove_parent_curve_rotation_(remove_parent_curve_rotation), - remove_parent_curve_translation_(remove_parent_curve_translation), + parent_curve_normalization_(parent_curve_normalization), parent_curve_fn_(parent_curve_fn) { } Eigen::Matrix4d operator()(double u) const { Eigen::Matrix4d parent_curve_point = (*parent_curve_fn_)(u); - Eigen::Matrix4d curve_segment_point = curve_segment_placement_ * remove_parent_curve_rotation_ * remove_parent_curve_translation_ * parent_curve_point; + Eigen::Matrix4d curve_segment_point = curve_segment_placement_ * parent_curve_normalization_ * parent_curve_point; return curve_segment_point + parent_curve_fn_->curvature(u); } private: Eigen::Matrix4d curve_segment_placement_; - Eigen::Matrix4d remove_parent_curve_rotation_; - Eigen::Matrix4d remove_parent_curve_translation_; + Eigen::Matrix4d parent_curve_normalization_; std::shared_ptr parent_curve_fn_; }; @@ -166,9 +164,46 @@ struct cant_curve_segment_function { // Subtract the parent_curve_start_point to get the incremental cant rotation and superelevation // Add the incremental cant rotation and superelevation to curve_segment_placement to get the curve_segment_point Eigen::Matrix4d parent_curve_point = (*parent_curve_fn_)(u); - Eigen::Matrix4d cant_increment = parent_curve_point - parent_curve_start_point_; - Eigen::Matrix4d curve_segment_point = curve_segment_placement_ + cant_increment; + + Eigen::Matrix3d parent_curve_start_rotation_ = parent_curve_start_point_.block<3, 3>(0, 0); + Eigen::Matrix3d parent_curve_point_rotation = parent_curve_point.block<3, 3>(0, 0); + Eigen::Matrix3d incremental_rotation = parent_curve_point_rotation * parent_curve_start_rotation_.transpose(); + Eigen::Matrix3d placement_rotation_ = curve_segment_placement_.block<3, 3>(0, 0); + Eigen::Matrix3d curve_segment_rotation = incremental_rotation * placement_rotation_; + + Eigen::Vector3d parent_curve_start_translation_ = parent_curve_start_point_.block<3, 1>(0, 3); + Eigen::Vector3d parent_curve_point_translation = parent_curve_point.block<3, 1>(0, 3); + Eigen::Vector3d incremental_translation = parent_curve_point_translation - parent_curve_start_translation_; + Eigen::Vector3d placement_translation_ = curve_segment_placement_.block<3, 1>(0, 3); + Eigen::Vector3d curve_segment_translation = incremental_translation + placement_translation_; + + Eigen::Matrix4d curve_segment_point = Eigen::Matrix4d::Identity(); + curve_segment_point.block<3, 3>(0, 0) = curve_segment_rotation; + curve_segment_point.block<3, 1>(0, 3) = curve_segment_translation; + + //if (0.0 < u) { + // Eigen::IOFormat latexFormat( + // Eigen::FullPrecision, // full precision + // 0, // no alignment flags + // " & ", // coeff separator + // " \\\\ \n", // row separator + // "", // row prefix + // "", // row suffix + // "", // matrix prefix + // "" // matrix suffix + // ); + // std::cout << "Placement (M_CSP)" << std::endl; + // std::cout << curve_segment_placement_.format(latexFormat) << std::endl; + // std::cout << "Parent curve start point (M_PCS)" << std::endl; + // std::cout << parent_curve_start_point_.format(latexFormat) << std::endl; + // std::cout << "Parent curve point (M_PCl)" << std::endl; + // std::cout << parent_curve_point.format(latexFormat) << std::endl; + // std::cout << "Curve segment point (M_c)" << std::endl; + // std::cout << curve_segment_point.format(latexFormat) << std::endl; + //} + return curve_segment_point + parent_curve_fn_->curvature(u); + } private: @@ -314,33 +349,26 @@ class curve_segment_evaluator { return taxonomy::make(length, fn); } else { // The parent curve function returns the 4x4 matrix for the parent curve. - // Subtract the parent curve start point (remove the translation and rotation) - // to get the incremental translation and rotation. Apply the incremental + // Normalize the parent curve so that the trim start point and tangent direction at the start point + // are aligned with the origin. This is accomplished with a normalization matrix that subtracts the + // incremental parent curve start point and applies a rotation. Apply the incremental // translation and rotation to the curve_segment_placement to get the curve_segment_point - // Do a negative translation of the parent curve point relative to the start of the parent curve. - // This moves parent_curve_fn(u=0.0) to coordinate (0,0). - // This is done so the curve_segment_placement is applied relative to (0,0) - Eigen::Matrix4d remove_parent_curve_translation = Eigen::Matrix4d::Identity(); - remove_parent_curve_translation.col(3) = -1.0 * (*parent_curve_start_point_).col(3); - remove_parent_curve_translation(3, 3) = 1.0; + auto rotation = (*parent_curve_start_point_).block<3, 3>(0, 0); + auto dxo = rotation(0, 0); + auto dyo = rotation(1, 0); + rotation(0, 1) *= -1.0; + rotation(1, 0) *= -1.0; + auto xo = (*parent_curve_start_point_)(0, 3); + auto yo = (*parent_curve_start_point_)(1, 3); + auto xn = -xo*dxo - yo*dyo; + auto yn = xo*dyo - yo*dxo; + Eigen::Matrix4d parent_curve_normalization = Eigen::Matrix4d::Identity(); + parent_curve_normalization.block<3, 3>(0, 0) = rotation; + parent_curve_normalization(0, 3) = xn; + parent_curve_normalization(1, 3) = yn; - // Do a rotation so that the tangent of the parent curve is in the direction (1,0) - // Example: if the parent curve IfcLine is at a 30 degree clockwise angle, this does - // a 30 degree counter-clockwise rotation - // Clockwise rotation matrix = [cos(angle) -sin(angle)] - // [sin(angle) cos(angle)] - // - // Counter-clockwise rotation = [ cos(angle) sin(angle)] - // [-sin(angle) cos(angle)] - // - // That's just a sign flip in positions (0,1) and (1,0) - Eigen::Matrix4d remove_parent_curve_rotation = (*parent_curve_start_point_); - remove_parent_curve_rotation(0, 1) *= -1.0; - remove_parent_curve_rotation(1, 0) *= -1.0; - remove_parent_curve_rotation.col(3) = Eigen::Vector4d(0, 0, 0, 1); // remove the parent curve placement point - - auto fn = curve_segment_function(*curve_segment_placement_, remove_parent_curve_rotation, remove_parent_curve_translation, parent_curve_fn_); + auto fn = curve_segment_function(*curve_segment_placement_, parent_curve_normalization, parent_curve_fn_); return taxonomy::make(length, fn); } } @@ -826,6 +854,8 @@ class curve_segment_evaluator { auto R = c.Radius() * length_unit_; auto parent_curve_position = taxonomy::cast(mapping_->map(c.Position()))->ccomponents(); + auto sign_l = sign(length_); + // center point of the parent curve auto pcCenterX = parent_curve_position(0, 3); auto pcCenterY = parent_curve_position(1, 3); @@ -839,7 +869,8 @@ class curve_segment_evaluator { // angle from X = 0 to the first point on the trimmed curve auto start_angle = pc_axis_angle + sweep_start_angle; - auto sign_l = sign(length_); + auto pcStartX = pcCenterX + R * cos(start_angle); + auto pcStartY = pcCenterY + R * sin(start_angle); projected_length_ = length_; @@ -851,34 +882,27 @@ class curve_segment_evaluator { #ifdef SCHEMA_IfcCurveSegment_HAS_Placement curve_segment_placement = taxonomy::cast(mapping_->map(inst_.Placement()))->ccomponents(); #endif - auto csStartX = curve_segment_placement(0, 3); - auto csStartY = curve_segment_placement(1, 3); - auto csStartDx = curve_segment_placement(0, 0); - auto csStartDy = curve_segment_placement(1, 0); - auto csCenterX = csStartX - sign_l * csStartDy * R; - auto csCenterY = csStartY + sign_l * csStartDx * R; - // determine projected length along the x-axis auto subtended_angle = R ? length_ / R : 0.0; auto end_angle = start_angle + subtended_angle; - auto csEndX = csCenterX + R * cos(end_angle); - projected_length_ = csEndX - csStartX; + auto pcEndX = pcCenterX + R * cos(end_angle); + projected_length_ = pcEndX - pcStartX; - convert_u = [csStartX, csStartY, csCenterX, csCenterY, R, sign_l](double u) { + convert_u = [pcStartX, pcStartY, pcCenterX, pcCenterY, R, sign_l](double u) { // for vertical, u is measured along the horizonal but we need it to be an arc length // x and y are coordinates on the curve segment for horizontal distance u from the start point // u is a horizontal distance so x = csStartX + u // Recognizing the triangle - // R^2 = (u + csStartX - csCenterX)^2 + (y - csCenterY)^2 + // R^2 = (u + pcStartX - pcCenterX)^2 + (y - pcCenterY)^2 // solve for y - // (y - csCenterY) = sqrt( R^2 - (u + csStartX - csCenterX)^2 ) - // y = csCenterY + sqrt( R^2 - (u + csStartX - csCenterX)^2 ) - auto x = csStartX + u; - auto y = csCenterY - sign_l * sqrt(pow(R, 2) - pow(u + csStartX - csCenterX, 2)); + // (y - pcCenterY) = sqrt( R^2 - (u + pcStartX - pcCenterX)^2 ) + // y = pcCenterY + sqrt( R^2 - (u + pcStartX - pcCenterX)^2 ) + auto x = pcStartX + u; + auto y = pcCenterY - sign_l * sqrt(pow(R, 2) - pow(u + pcStartX - pcCenterX, 2)); // compute the chord distance between the start point and (x,y) - auto c = sqrt(pow(x - csStartX, 2.0) + pow(y - csStartY, 2.0)); + auto c = sqrt(pow(x - pcStartX, 2.0) + pow(y - pcStartY, 2.0)); // compute the subtended angle // c = 2R*sin(delta/2) @@ -981,18 +1005,18 @@ class curve_segment_evaluator { double m_squared = std::inner_product(dr.begin(), dr.end(), dr.begin(), 0.0); double m = sqrt(m_squared); std::transform(dr.begin(), dr.end(), dr.begin(), [m](auto& d) { return d / m; }); - auto pcDx = dr[0]; - auto pcDy = dr[1]; + auto pcDXx = dr[0]; + auto pcDXy = dr[1]; if (segment_type_ == ST_VERTICAL && curve_segment_placement_) { // the general algorithm for mapping parent curve onto curve segment doesn't // exactly work for IfcLine. This is easily overcome by using the curve segment // placement for the IfcLine direction - pcDx = (*curve_segment_placement_)(0, 0); - pcDy = (*curve_segment_placement_)(1, 0); + pcDXx = (*curve_segment_placement_)(0, 0); + pcDXy = (*curve_segment_placement_)(1, 0); // projected length along the x-axis is the 'i' component of the total length - projected_length_ = length_ * pcDx; + projected_length_ = length_ * pcDXx; } if (segment_type_ == ST_HORIZONTAL || segment_type_ == ST_VERTICAL || segment_type_ == ST_CANT) { @@ -1001,19 +1025,32 @@ class curve_segment_evaluator { convert_u = [](double u) { return u; }; // u is along curve } else { // u is along horizontal, convert to along curve - convert_u = [pcDx](double u) { return u/pcDx; }; + convert_u = [pcDXx](double u) { return u/pcDXx; }; } + auto pcDZy = curve_segment_placement_ ? (*curve_segment_placement_)(1, 2) : 0.; + auto pcDZz = curve_segment_placement_ ? (*curve_segment_placement_)(2, 2) : 1.; + parent_curve_fn_ = std::make_shared( - [pcX, pcY, pcDx, pcDy, convert_u](double u)->Eigen::Matrix4d { + [segment_type = segment_type_,pcX, pcY, pcDXx, pcDXy, pcDZy, pcDZz, convert_u](double u)->Eigen::Matrix4d { u = convert_u(u); - auto x = pcX + pcDx * u; - auto y = pcY + pcDy * u; + auto x = pcX + pcDXx * u; + auto y = pcY + pcDXy * u; Eigen::Matrix4d m = Eigen::Matrix4d::Identity(); - m.col(0) = Eigen::Vector4d(pcDx, pcDy, 0, 0); - m.col(1) = Eigen::Vector4d(-pcDy, pcDx, 0, 0); + Eigen::Vector3d X(pcDXx, pcDXy, 0); + Eigen::Vector3d Z(0, 0, 1); + + if (segment_type == ST_CANT) { + Z = Eigen::Vector3d(0, pcDZy, pcDZz); + } + + Eigen::Vector3d Y = Z.cross(X).normalized(); + + m.col(0) = Eigen::Vector4d(X[0], X[1], X[2], 0); + m.col(1) = Eigen::Vector4d(Y[0], Y[1], Y[2], 0); + m.col(2) = Eigen::Vector4d(Z[0], Z[1], Z[2], 0); m.col(3) = Eigen::Vector4d(x, y, 0.0, 1.0); return m; },