From 53028d9c97d5e9d3728e7bde9975b9971ab8c9b9 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Thu, 21 Sep 2023 15:48:12 -0700 Subject: [PATCH 01/15] Updated build scripts to use boost_1.79.0 --- win/build-deps.cmd | 2 +- win/run-cmake.bat | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/win/build-deps.cmd b/win/build-deps.cmd index 07acffd6ec..9d5ac6f07a 100644 --- a/win/build-deps.cmd +++ b/win/build-deps.cmd @@ -114,7 +114,7 @@ if "%CMAKE_VERSION%" LSS "cmake version 3.11.4" ( ) :: NOTE Boost < 1.64 doesn't work without tricks if the user has only VS 2017 installed and no earlier versions. -set BOOST_VERSION=1.74.0 +set BOOST_VERSION=1.79.0 :: Version string with underscores instead of dots. set BOOST_VER=%BOOST_VERSION:.=_% diff --git a/win/run-cmake.bat b/win/run-cmake.bat index 92741419bc..b8c81f03f9 100755 --- a/win/run-cmake.bat +++ b/win/run-cmake.bat @@ -80,7 +80,7 @@ IF NOT EXIST ..\%BUILD_DIR%. mkdir ..\%BUILD_DIR% pushd ..\%BUILD_DIR% :: tfk: todo remove duplication -set BOOST_VERSION=1.74.0 +set BOOST_VERSION=1.79.0 set BOOST_VER=%BOOST_VERSION:.=_% set BOOST_ROOT=%DEPS_DIR%\boost_%BOOST_VER% From 04b70efd7fb6706c8e05d0247ade83dac131a0e7 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Thu, 21 Sep 2023 15:47:40 -0700 Subject: [PATCH 02/15] Fixed problem with call to boost::mpl::for_each --- src/ifcgeom/mapping/IfcCurveSegment.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index 61fc913c15..9d13729d8d 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -121,7 +121,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { static int NUM_SEGMENTS = 64; curve_segment_evaluator cse(length_unit_, inst->ParentCurve(), inst->SegmentStart(), inst->SegmentLength()); - boost::mpl::for_each>(cse); + boost::mpl::for_each>(std::ref(cse)); std::vector polygon; From 3e272efce419bfa9a437f02283cd3bd2a031c99a Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Fri, 22 Sep 2023 13:49:18 -0700 Subject: [PATCH 03/15] Fixes cast to IfcCurveSegment for access to ParentCurve --- src/ifcgeom/mapping/IfcCompositeCurve.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ifcgeom/mapping/IfcCompositeCurve.cpp b/src/ifcgeom/mapping/IfcCompositeCurve.cpp index 7ed02dfabb..136953322a 100644 --- a/src/ifcgeom/mapping/IfcCompositeCurve.cpp +++ b/src/ifcgeom/mapping/IfcCompositeCurve.cpp @@ -66,7 +66,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) { } #ifdef SCHEMA_HAS_IfcCurveSegment else if (segment->as()) { - auto crv = map(segment->as()->ParentCurve()); + auto crv = map(segment->as()->ParentCurve()); for (auto& s : taxonomy::cast(crv)->children) { loop->children.push_back(s); } From 6ae0cce186bb10ee298a2febfdb186a77b246144 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Fri, 22 Sep 2023 13:52:48 -0700 Subject: [PATCH 04/15] Implements mapping of IfcCurveSegment to taxonomy::geometry --- src/ifcgeom/mapping/IfcCompositeCurve.cpp | 2 +- src/ifcgeom/mapping/IfcCurveSegment.cpp | 209 ++++++++++++++++++++-- 2 files changed, 198 insertions(+), 13 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCompositeCurve.cpp b/src/ifcgeom/mapping/IfcCompositeCurve.cpp index 136953322a..1124cd083e 100644 --- a/src/ifcgeom/mapping/IfcCompositeCurve.cpp +++ b/src/ifcgeom/mapping/IfcCompositeCurve.cpp @@ -66,7 +66,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) { } #ifdef SCHEMA_HAS_IfcCurveSegment else if (segment->as()) { - auto crv = map(segment->as()->ParentCurve()); + auto crv = map(segment->as()); for (auto& s : taxonomy::cast(crv)->children) { loop->children.push_back(s); } diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index 9d13729d8d..1de6bfaf5f 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -33,6 +33,8 @@ typedef boost::mpl::vector< #ifdef SCHEMA_HAS_IfcClothoid , IfcSchema::IfcClothoid #endif + , IfcSchema::IfcPolyline + , IfcSchema::IfcCircle > curve_seg_types; class curve_segment_evaluator { @@ -65,32 +67,189 @@ public: // Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes void operator()(IfcSchema::IfcClothoid* c) { // @todo verify - auto L = start_ + length_; + auto sign = [](double v)->int{return v < 0 ? -1 : (0 < v ? 1 : 0); }; + auto sign_s = sign(start_); + auto sign_l = sign(length_); + double L = 0; + if (sign_s == 0) L = fabs(length_); + else if (sign_s == sign_l) L = fabs(start_ + length_); + else L = fabs(start_); + auto A = c->ClothoidConstant(); auto R = A * A / L; - auto RL = R * L; + auto RL = (A < 0 ? -1.0 : 1.0) * R * L; - eval_ = [RL](double u) { + auto position = c->Position(); + auto placement = position->as(); + auto ref_direction = placement->RefDirection(); + double theta = 0.0; // angle the circle's placement X-axis makes with respect to global X axis + if (ref_direction) + { + auto dr = ref_direction->DirectionRatios(); + auto dx = dr[0]; + auto dy = dr[1]; + theta = atan2(dy, dx); + } + + auto C = placement->Location(); + if (!C->as()) + { + throw std::runtime_error("Only IfcCartesianPoint is supported for center of IfcCircle"); + // @todo add support for other IfcPoint subtypes + } + auto Cx = C->as()->Coordinates()[0]; + auto Cy = C->as()->Coordinates()[1]; + + eval_ = [RL,Cx,Cy,theta](double u) { + // coordinate along clothoid is local coordinates auto xterm_1 = u; auto xterm_2 = std::pow(u, 5) / (40 * std::pow(RL, 2)); auto xterm_3 = std::pow(u, 9) / (3456 * std::pow(RL, 4)); auto xterm_4 = std::pow(u, 13) / (599040 * std::pow(RL, 6)); - auto x = xterm_1 - xterm_2 + xterm_3 - xterm_4; + auto xl = xterm_1 - xterm_2 + xterm_3 - xterm_4; auto yterm_1 = std::pow(u, 3) / (6 * RL); auto yterm_2 = std::pow(u, 7) / (336 * std::pow(RL, 3)); auto yterm_3 = std::pow(u, 11) / (42240 * std::pow(RL, 5)); auto yterm_4 = std::pow(u, 15) / (9676800 * std::pow(RL, 7)); - auto y = yterm_1 - yterm_2 + yterm_3 - yterm_4; + auto yl = yterm_1 - yterm_2 + yterm_3 - yterm_4; - return Eigen::Vector3d(x, y, 0.); + // transform point into clothoid's coodinate system + auto x = xl * cos(theta) - yl * sin(theta) + Cx; + auto y = xl * sin(theta) + yl * cos(theta) + Cy; + return Eigen::Vector3d(x, y, 0.0); }; } #endif - // Another IfcCurve subtype - void operator()(IfcSchema::IfcLine*) { - throw std::runtime_error("not implemented"); + void operator()(IfcSchema::IfcCircle* c) + { + auto R = c->Radius(); + + auto position = c->Position(); + auto placement = position->as(); + auto ref_direction = placement->RefDirection(); + double theta = 0.0; // angle the circle's placement X-axis makes with respect to global X axis + if (ref_direction) + { + auto dr = ref_direction->DirectionRatios(); + auto dx = dr[0]; + auto dy = dr[1]; + theta = atan2(dy, dx); + } + + // center of circle location + auto C = placement->Location(); + if (!C->as()) + { + throw std::runtime_error("Only IfcCartesianPoint is supported for center of IfcCircle"); + // @todo add support for other IfcPoint subtypes + } + auto Cx = C->as()->Coordinates()[0]; + auto Cy = C->as()->Coordinates()[1]; + + eval_ = [R, Cx, Cy, theta](double u) + { + auto angle = u / R; // angle subtended by arc length u + + // compute point on circle centered at (0,0) with x-axis horizontal and y-axis vertical + auto xl = R * cos(angle); + auto yl = R * sin(angle); + + // transform point into circle's coodinate system + auto x = xl * cos(theta) - yl * sin(theta) + Cx; + auto y = xl * sin(theta) + yl * cos(theta) + Cy; + return Eigen::Vector3d(x, y, 0.0); + }; + } + + void operator()(IfcSchema::IfcPolyline* pl) + { + struct Range + { + double u_start; + double u_end; + std::function compare; + bool operator<(const Range& r) const { return u_start < r.u_start; } + }; + using Function = std::function(double u)>; + std::map fns; + + auto p = pl->Points(); + if (p->size() < 2) + { + throw std::runtime_error("invalid polyline - must have at least 2 points"); // this should never happen, but just in case it does + } + + 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 iter = p->begin(); + auto end = p->end(); + auto last = std::prev(end); + auto p1 = *(iter++); + auto u = 0.0; + for (; iter != end; iter++) + { + auto p2 = *iter; + + auto p1x = p1->Coordinates()[0]; + auto p1y = p1->Coordinates()[1]; + + auto p2x = p2->Coordinates()[0]; + auto p2y = p2->Coordinates()[1]; + + auto dx = p2x - p1x; + auto dy = p2y - p1y; + auto l = sqrt(dx * dx + dy * dy); + if (l == 0.0) + { + // @todo use closeness tolerance instead of absolute 0.0 + throw std::runtime_error("invalid polyline - points must not be coincident"); + } + + dx /= l; + dy /= l; + + auto fn = [p1x, p1y, dx, dy](double u) { return std::make_pair(p1x + u * dx, p1y + u * dy); }; + + fns.insert(std::make_pair(Range{ u, u + l,iter == last ? end_compare : std_compare }, fn)); + + p1 = p2; + u = u + l; + } + + eval_ = [fns](double u) { + auto iter = std::find_if(fns.cbegin(), fns.cend(), [=](const auto& fn) + { + auto [u_start, u_end, compare] = fn.first; + return compare(u_start, u, u_end); + }); + + if (iter == fns.end()) throw std::runtime_error("invalid distance from start"); // this should never happen, but just in case it does + + auto [u_start, u_end, compare] = iter->first; + auto [x,y] = (iter->second)(u - u_start); // (u - u_start) is distance from start of this segment of the polyline + return Eigen::Vector3d(x, y, 0); + }; + } + + void operator()(IfcSchema::IfcLine* l) { + auto s = l->Pnt(); + auto c = s->Coordinates(); + auto v = l->Dir(); + auto dr = v->Orientation()->DirectionRatios(); + auto m = v->Magnitude(); + auto px = c[0]; + auto py = c[1]; + auto dx = dr[0] / m; + auto dy = dr[1] / m; + + eval_ = [px, py, dx, dy](double u) { + auto x = px + u * dx; + auto y = py + u * dy; + return Eigen::Vector3d(x, y, 0); + }; } // Take the boost::type value from mpl::for_each and test it against our curve instance @@ -118,6 +277,7 @@ public: taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { // @todo fixed number of segments or fixed interval? // @todo placement + // @todo figure out what to do with the zero length segments at the end of compound curves static int NUM_SEGMENTS = 64; curve_segment_evaluator cse(length_unit_, inst->ParentCurve(), inst->SegmentStart(), inst->SegmentLength()); @@ -125,9 +285,34 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { std::vector polygon; - for (int i = 0; i <= NUM_SEGMENTS; ++i) { - auto p = cse(cse.length() * i / NUM_SEGMENTS); - polygon.push_back(taxonomy::make(p(0), p(1), p(2))); + auto placement = inst->Placement(); + auto location = placement->Location(); + auto Cx = location->as()->Coordinates()[0]; + auto Cy = location->as()->Coordinates()[1]; + auto ref_dir = placement->as()->RefDirection(); + auto dx = ref_dir->DirectionRatios()[0]; + auto dy = ref_dir->DirectionRatios()[1]; + auto angle = atan2(dy, dx); + + auto cos_angle = cos(angle); + auto sin_angle = sin(angle); + + auto length = cse.length(); + if (0.001 < fabs(length)) + { + for (int i = 0; i <= NUM_SEGMENTS; ++i) { + auto u = length * i / NUM_SEGMENTS; + + auto p = cse(u); + auto xl = p(0); + auto yl = p(1); + auto z = p(2); + + auto x = xl * cos_angle - yl * sin_angle + Cx; + auto y = xl * sin_angle + yl * cos_angle + Cy; + + polygon.push_back(taxonomy::make(x,y,z)); + } } return polygon_from_points(polygon); From 4f2b083ffc2fd6fe26b5448cf54e84a45c9c1054 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Thu, 28 Sep 2023 11:51:26 -0700 Subject: [PATCH 05/15] Revises implementation of operator(IfcSchema::IfcClothoid*) to use trapezoid rule numeric integration. This is a more general approach than the Taylor Series approximation and is applicable to all the other spiral types. I stubbed out a functor for IfcSecondOrderPolynomialSpiral, but haven't tested it. The purpose is to show the concept. --- src/ifcgeom/mapping/IfcCurveSegment.cpp | 191 ++++++++++++++++++++---- 1 file changed, 161 insertions(+), 30 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index 1de6bfaf5f..cd92d229bb 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -28,10 +28,38 @@ using namespace ifcopenshell::geometry; #include #include +// @todo use std::numbers::pi when upgrading to C++ 20 +#define PI 3.1415926535897932384626433832795 + + +namespace +{ + // trapezoid rule integration + // @todo is there a well established math library we can use instead of + // creating our own integrator? + double integrate(double a, double b, unsigned n, std::function fn) + { + double area = 0; + double h = (b - a) / (n + 1); + for (auto i = 0; i < n; i++) + { + auto x1 = a + h * i; + auto x2 = a + h * (i + 1); + auto f1 = fn(x1); + auto f2 = fn(x2); + area += h * (f1 + f2) / 2.0; + } + return area; + } +} + typedef boost::mpl::vector< IfcSchema::IfcLine #ifdef SCHEMA_HAS_IfcClothoid , IfcSchema::IfcClothoid +#endif +#if defined SCHEMA_HAS_IfcSecondOrderPolynomialSpiral + , IfcSchema::IfcSecondOrderPolynomialSpiral #endif , IfcSchema::IfcPolyline , IfcSchema::IfcCircle @@ -63,24 +91,80 @@ public: length_ = *le->as() * length_unit; } -#ifdef SCHEMA_HAS_IfcClothoid - // Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes - void operator()(IfcSchema::IfcClothoid* c) { - // @todo verify - auto sign = [](double v)->int{return v < 0 ? -1 : (0 < v ? 1 : 0); }; - auto sign_s = sign(start_); - auto sign_l = sign(length_); +// Clothoid using Taylor Series approximation +//#ifdef SCHEMA_HAS_IfcClothoid +// // Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes +// void operator()(IfcSchema::IfcClothoid* c) { +// // @todo verify +// auto sign = [](double v)->int{return v < 0 ? -1 : (0 < v ? 1 : 0); }; +// auto sign_s = sign(start_); +// auto sign_l = sign(length_); +// double L = 0; +// if (sign_s == 0) L = fabs(length_); +// else if (sign_s == sign_l) L = fabs(start_ + length_); +// else L = fabs(start_); +// +// auto A = c->ClothoidConstant(); +// auto R = A * A / L; +// auto RL = (A < 0 ? -1.0 : 1.0) * R * L; +// +// auto position = c->Position(); +// auto placement = position->as(); +// auto ref_direction = placement->RefDirection(); +// double theta = 0.0; // angle the circle's placement X-axis makes with respect to global X axis +// if (ref_direction) +// { +// auto dr = ref_direction->DirectionRatios(); +// auto dx = dr[0]; +// auto dy = dr[1]; +// theta = atan2(dy, dx); +// } +// +// auto C = placement->Location(); +// if (!C->as()) +// { +// throw std::runtime_error("Only IfcCartesianPoint is supported for center of IfcCircle"); +// // @todo add support for other IfcPoint subtypes +// } +// auto Cx = C->as()->Coordinates()[0]; +// auto Cy = C->as()->Coordinates()[1]; +// +// eval_ = [RL,Cx,Cy,theta](double u) { +// // coordinate along clothoid is local coordinates +// auto xterm_1 = u; +// auto xterm_2 = std::pow(u, 5) / (40 * std::pow(RL, 2)); +// auto xterm_3 = std::pow(u, 9) / (3456 * std::pow(RL, 4)); +// auto xterm_4 = std::pow(u, 13) / (599040 * std::pow(RL, 6)); +// auto xl = xterm_1 - xterm_2 + xterm_3 - xterm_4; +// +// auto yterm_1 = std::pow(u, 3) / (6 * RL); +// auto yterm_2 = std::pow(u, 7) / (336 * std::pow(RL, 3)); +// auto yterm_3 = std::pow(u, 11) / (42240 * std::pow(RL, 5)); +// auto yterm_4 = std::pow(u, 15) / (9676800 * std::pow(RL, 7)); +// auto yl = yterm_1 - yterm_2 + yterm_3 - yterm_4; +// +// // transform point into clothoid's coodinate system +// auto x = xl * cos(theta) - yl * sin(theta) + Cx; +// auto y = xl * sin(theta) + yl * cos(theta) + Cy; +// return Eigen::Vector3d(x, y, 0.0); +// }; +// } +//#endif + + void set_spiral_functor(IfcSchema::IfcSpiral* 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 binary_sign = [](double v)->int {return v < 0 ? -1 : (0 < v ? 1 : 0); }; // returns -1, 0, or 1 + auto sign_s = binary_sign(start_); + auto sign_l = binary_sign(length_); double L = 0; - if (sign_s == 0) L = fabs(length_); - else if (sign_s == sign_l) L = fabs(start_ + length_); - else L = fabs(start_); + 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 A = c->ClothoidConstant(); - auto R = A * A / L; - auto RL = (A < 0 ? -1.0 : 1.0) * R * L; - - auto position = c->Position(); - auto placement = position->as(); + auto position = s->Position(); + auto placement = position->as(); // @todo Update, this could be IfcAxis2Placement2D or IfcAxisPlacement3D + if (!placement) { throw std::runtime_error("Only IfcAxis2Placement2D is supported right now"); } auto ref_direction = placement->RefDirection(); double theta = 0.0; // angle the circle's placement X-axis makes with respect to global X axis if (ref_direction) @@ -94,31 +178,78 @@ public: auto C = placement->Location(); if (!C->as()) { - throw std::runtime_error("Only IfcCartesianPoint is supported for center of IfcCircle"); + throw std::runtime_error("Only IfcCartesianPoint is supported right now"); // @todo add support for other IfcPoint subtypes } auto Cx = C->as()->Coordinates()[0]; auto Cy = C->as()->Coordinates()[1]; - eval_ = [RL,Cx,Cy,theta](double u) { - // coordinate along clothoid is local coordinates - auto xterm_1 = u; - auto xterm_2 = std::pow(u, 5) / (40 * std::pow(RL, 2)); - auto xterm_3 = std::pow(u, 9) / (3456 * std::pow(RL, 4)); - auto xterm_4 = std::pow(u, 13) / (599040 * std::pow(RL, 6)); - auto xl = xterm_1 - xterm_2 + xterm_3 - xterm_4; + eval_ = [L, Cx, Cy, theta, signX, fnX, signY, fnY](double u) { + // integration limits, integrate from a to b + // from 8.9.3.19.1, integration limits are 0.0 to u where u is a normalized parameter + auto a = 0.0; + auto b = fabs(u / L); - auto yterm_1 = std::pow(u, 3) / (6 * RL); - auto yterm_2 = std::pow(u, 7) / (336 * std::pow(RL, 3)); - auto yterm_3 = std::pow(u, 11) / (42240 * std::pow(RL, 5)); - auto yterm_4 = std::pow(u, 15) / (9676800 * std::pow(RL, 7)); - auto yl = yterm_1 - yterm_2 + yterm_3 - yterm_4; + auto n = 10; // use 10 steps in the numeric integration + + auto xl = signX(u)*integrate(a, b, n, fnX); + auto yl = signY(u)*integrate(a, b, n, fnY); // transform point into clothoid's coodinate system auto x = xl * cos(theta) - yl * sin(theta) + Cx; auto y = xl * sin(theta) + yl * cos(theta) + Cy; return Eigen::Vector3d(x, y, 0.0); - }; + }; + } + +// Clothoid using numerical integration +#ifdef SCHEMA_HAS_IfcClothoid +// Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes + void operator()(IfcSchema::IfcClothoid* c) { + + auto A = c->ClothoidConstant(); + + // the integration is for the +X, +Y quadrant - need to adjust the signs of the resulting X and Y values + // so that the results are in the correct quadrant. + // A > 0 and u > 0 -> +X, +Y + // A < 0 and u > 0 -> +X, -Y + // A > 0 and u < 0 -> -X, -Y + // A < 0 and u < 0 -> -X, +Y + // X depends only on u, Y depends on u and A. + auto sign = [](double v)->int {return v < 0 ? -1 : 1; }; // returns -1 or 1 + auto sign_x = [sign](double t) {return sign(t); }; + auto sign_y = [sign, A](double t) {return sign(t) == sign(A) ? 1.0 : -1.0; }; + auto fn_x = [A](double t)->double {return A * sqrt(PI) * cos(PI * A * t * t / (2 * fabs(A))); }; + auto fn_y = [A](double t)->double {return A * sqrt(PI) * sin(PI * A * t * t / (2 * fabs(A))); }; + + set_spiral_functor(c->as(), sign_x, fn_x, sign_y, fn_y); + } +#endif + +#ifdef SCHEMA_HAS_IfcSecondOrderPolynomialSpiral + void operator()(IfcSchema::IfcSecondOrderPolynomialSpiral* s) + { + // @todo verify - this is an example implementation of a different kind of spiral - lots of clean up needed + auto A0 = s->ConstantTerm(); + auto A1 = s->LinearTerm(); + auto A2 = s->QuadraticTerm(); + + auto theta = [A0, A1, A2](double t) + { + auto a0 = A0.has_value() ? t / A0.value() : 0.0; + auto a1 = A1.has_value() ? A1.value() * std::pow(t, 2) / (2 * fabs(std::pow(A1.value(), 3))) : 0.0; + auto a2 = std::pow(t, 3) / (3 * std::pow(A2, 3)); + return a0 + a1 + a2; + }; + + auto sign = [](double v)->int {return v < 0 ? -1 : 1; }; // returns -1 or 1 + auto sign_x = [sign](double t) {return sign(t); }; + auto sign_y = [sign](double t) {return sign(t); }; // @todo fix - not sure about sign_y yet, need to find some plots of this spiral + + auto fn_x = [theta](double t)->double {return cos(theta(t)); }; + auto fn_y = [theta](double t)->double {return sin(theta(t)); }; + + set_spiral_functor(s->as(), sign_x, fn_x, sign_y, fn_y); } #endif From 1564b53ef5aee95b12c1c45c31deb0a6ea16d6f2 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Thu, 28 Sep 2023 14:28:17 -0700 Subject: [PATCH 06/15] Minor revision to numerical integration code --- src/ifcgeom/mapping/IfcCurveSegment.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index cd92d229bb..727032895d 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -40,11 +40,11 @@ namespace double integrate(double a, double b, unsigned n, std::function fn) { double area = 0; - double h = (b - a) / (n + 1); - for (auto i = 0; i < n; i++) + double h = (b - a) / n; + for (auto i = 1; i <= n; i++) { - auto x1 = a + h * i; - auto x2 = a + h * (i + 1); + auto x1 = a + h * (i - 1); + auto x2 = a + h * i; auto f1 = fn(x1); auto f2 = fn(x2); area += h * (f1 + f2) / 2.0; From d4f543d2ad0bc5402d1e1cc62b94be53fdee1708 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Thu, 28 Sep 2023 14:58:21 -0700 Subject: [PATCH 07/15] Implemented placement handling with transformation matrix. --- src/ifcgeom/mapping/IfcCurveSegment.cpp | 24 +++++------------------- 1 file changed, 5 insertions(+), 19 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index 727032895d..7d81a29eed 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -416,17 +416,9 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { std::vector polygon; - auto placement = inst->Placement(); - auto location = placement->Location(); - auto Cx = location->as()->Coordinates()[0]; - auto Cy = location->as()->Coordinates()[1]; - auto ref_dir = placement->as()->RefDirection(); - auto dx = ref_dir->DirectionRatios()[0]; - auto dy = ref_dir->DirectionRatios()[1]; - auto angle = atan2(dy, dx); - - auto cos_angle = cos(angle); - auto sin_angle = sin(angle); + // @todo - for some reason this isn't working, the matrix gets all messed up + //const auto& transformation_matrix = taxonomy::cast(map(inst->Placement()))->ccomponents(); + auto transformation_matrix = taxonomy::cast(map(inst->Placement()))->ccomponents(); auto length = cse.length(); if (0.001 < fabs(length)) @@ -435,14 +427,8 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { auto u = length * i / NUM_SEGMENTS; auto p = cse(u); - auto xl = p(0); - auto yl = p(1); - auto z = p(2); - - auto x = xl * cos_angle - yl * sin_angle + Cx; - auto y = xl * sin_angle + yl * cos_angle + Cy; - - polygon.push_back(taxonomy::make(x,y,z)); + auto result = transformation_matrix * Eigen::Vector4d(p(0),p(1),p(2), 1.); + polygon.push_back(taxonomy::make(result(0),result(1),result(2))); } } From 47c319aa2d954c25774d61b0924bcaf87faa3a5c Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Fri, 29 Sep 2023 08:02:33 -0700 Subject: [PATCH 08/15] Adds a continuity check for IfcCurveSegment --- src/ifcgeom/mapping/IfcCompositeCurve.cpp | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCompositeCurve.cpp b/src/ifcgeom/mapping/IfcCompositeCurve.cpp index 1124cd083e..9f5221fe9f 100644 --- a/src/ifcgeom/mapping/IfcCompositeCurve.cpp +++ b/src/ifcgeom/mapping/IfcCompositeCurve.cpp @@ -67,9 +67,22 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) { #ifdef SCHEMA_HAS_IfcCurveSegment else if (segment->as()) { auto crv = map(segment->as()); - for (auto& s : taxonomy::cast(crv)->children) { - loop->children.push_back(s); + auto crv_as_loop = taxonomy::cast(crv); + + // The end of the previous segment must be at the same location as the start of this segment + // @todo - need to apply some tolerancing + if (!loop->children.empty() and !crv_as_loop->children.empty() + and + boost::get(loop->children.back()->end)->components() != boost::get(crv_as_loop->children.front()->start)->components()) + { + std::ostringstream os; + auto& prev = boost::get(loop->children.back()->end)->components(); + auto& next = boost::get(crv_as_loop->children.front()->start)->components(); + os << "Common points are not continuous: (" << prev.x() << ", " << prev.y() << ", " << prev.z() << ")" << " " << "(" << next.x() << ", " << next.y() << ", " << next.z() << ")" << std::endl; + Logger::Notice(os.str()); } + + loop->children.insert(loop->children.end(), crv_as_loop->children.begin(), crv_as_loop->children.end()); } #endif } From d9b534b66d3f1c5686cd31f22750ac76b79fc2b7 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Fri, 29 Sep 2023 08:23:20 -0700 Subject: [PATCH 09/15] Implemented position handling with transformation matrix --- src/ifcgeom/mapping/IfcCurveSegment.cpp | 85 +++++++------------------ 1 file changed, 24 insertions(+), 61 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index 7d81a29eed..88d02f527b 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -53,6 +53,7 @@ namespace } } +// types of entities that can be IfcCurveSegment.ParentCurve typedef boost::mpl::vector< IfcSchema::IfcLine #ifdef SCHEMA_HAS_IfcClothoid @@ -67,6 +68,7 @@ typedef boost::mpl::vector< class curve_segment_evaluator { private: + mapping* mapping_; double length_unit_; double start_; double length_; @@ -76,8 +78,9 @@ private: public: // First constructor, takes parameters from IfcCurveSegment - curve_segment_evaluator(double length_unit, IfcSchema::IfcCurve* curve, IfcSchema::IfcCurveMeasureSelect* st, IfcSchema::IfcCurveMeasureSelect* le) - : length_unit_(length_unit) + curve_segment_evaluator(mapping* mapping,double length_unit, IfcSchema::IfcCurve* curve, IfcSchema::IfcCurveMeasureSelect* st, IfcSchema::IfcCurveMeasureSelect* le) + : mapping_(mapping) + , length_unit_(length_unit) , curve_(curve) { // @todo in IFC4X3_ADD2 this needs to be length measure @@ -151,7 +154,7 @@ public: // } //#endif - void set_spiral_functor(IfcSchema::IfcSpiral* s, std::function signX,std::function fnX, std::function signY, std::function fnY) + void set_spiral_functor(mapping* mapping,IfcSchema::IfcSpiral* 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 binary_sign = [](double v)->int {return v < 0 ? -1 : (0 < v ? 1 : 0); }; // returns -1, 0, or 1 @@ -162,29 +165,10 @@ public: 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 position = s->Position(); - auto placement = position->as(); // @todo Update, this could be IfcAxis2Placement2D or IfcAxisPlacement3D - if (!placement) { throw std::runtime_error("Only IfcAxis2Placement2D is supported right now"); } - auto ref_direction = placement->RefDirection(); - double theta = 0.0; // angle the circle's placement X-axis makes with respect to global X axis - if (ref_direction) - { - auto dr = ref_direction->DirectionRatios(); - auto dx = dr[0]; - auto dy = dr[1]; - theta = atan2(dy, dx); - } + //const auto& transformation_matrix = taxonomy::cast(mapping->map(s->Position()))->ccomponents(); + auto transformation_matrix = taxonomy::cast(mapping->map(s->Position()))->ccomponents(); - auto C = placement->Location(); - if (!C->as()) - { - throw std::runtime_error("Only IfcCartesianPoint is supported right now"); - // @todo add support for other IfcPoint subtypes - } - auto Cx = C->as()->Coordinates()[0]; - auto Cy = C->as()->Coordinates()[1]; - - eval_ = [L, Cx, Cy, theta, signX, fnX, signY, fnY](double u) { + eval_ = [L, transformation_matrix, signX, fnX, signY, fnY](double u) { // integration limits, integrate from a to b // from 8.9.3.19.1, integration limits are 0.0 to u where u is a normalized parameter auto a = 0.0; @@ -192,13 +176,12 @@ public: auto n = 10; // use 10 steps in the numeric integration - auto xl = signX(u)*integrate(a, b, n, fnX); - auto yl = signY(u)*integrate(a, b, n, fnY); + auto x = signX(u)*integrate(a, b, n, fnX); + auto y = signY(u)*integrate(a, b, n, fnY); // transform point into clothoid's coodinate system - auto x = xl * cos(theta) - yl * sin(theta) + Cx; - auto y = xl * sin(theta) + yl * cos(theta) + Cy; - return Eigen::Vector3d(x, y, 0.0); + auto result = transformation_matrix * Eigen::Vector4d(x, y, 0.0, 1.0); + return Eigen::Vector3d(result(0),result(1),result(2)); }; } @@ -222,7 +205,7 @@ public: auto fn_x = [A](double t)->double {return A * sqrt(PI) * cos(PI * A * t * t / (2 * fabs(A))); }; auto fn_y = [A](double t)->double {return A * sqrt(PI) * sin(PI * A * t * t / (2 * fabs(A))); }; - set_spiral_functor(c->as(), sign_x, fn_x, sign_y, fn_y); + set_spiral_functor(mapping_,c->as(), sign_x, fn_x, sign_y, fn_y); } #endif @@ -249,48 +232,27 @@ public: auto fn_x = [theta](double t)->double {return cos(theta(t)); }; auto fn_y = [theta](double t)->double {return sin(theta(t)); }; - set_spiral_functor(s->as(), sign_x, fn_x, sign_y, fn_y); + set_spiral_functor(mapping_, s->as(), sign_x, fn_x, sign_y, fn_y); } #endif void operator()(IfcSchema::IfcCircle* c) { auto R = c->Radius(); + //const auto& transformation_matrix = taxonomy::cast(mapping_->map(c->Position()))->ccomponents(); + auto transformation_matrix = taxonomy::cast(mapping_->map(c->Position()))->ccomponents(); - auto position = c->Position(); - auto placement = position->as(); - auto ref_direction = placement->RefDirection(); - double theta = 0.0; // angle the circle's placement X-axis makes with respect to global X axis - if (ref_direction) - { - auto dr = ref_direction->DirectionRatios(); - auto dx = dr[0]; - auto dy = dr[1]; - theta = atan2(dy, dx); - } - - // center of circle location - auto C = placement->Location(); - if (!C->as()) - { - throw std::runtime_error("Only IfcCartesianPoint is supported for center of IfcCircle"); - // @todo add support for other IfcPoint subtypes - } - auto Cx = C->as()->Coordinates()[0]; - auto Cy = C->as()->Coordinates()[1]; - - eval_ = [R, Cx, Cy, theta](double u) + eval_ = [R, transformation_matrix](double u) { auto angle = u / R; // angle subtended by arc length u // compute point on circle centered at (0,0) with x-axis horizontal and y-axis vertical - auto xl = R * cos(angle); - auto yl = R * sin(angle); + auto x = R * cos(angle); + auto y = R * sin(angle); // transform point into circle's coodinate system - auto x = xl * cos(theta) - yl * sin(theta) + Cx; - auto y = xl * sin(theta) + yl * cos(theta) + Cy; - return Eigen::Vector3d(x, y, 0.0); + auto result = transformation_matrix * Eigen::Vector4d(x, y, 0.0, 1.0); + return Eigen::Vector3d(result(0), result(1), result(2)); }; } @@ -411,7 +373,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { // @todo figure out what to do with the zero length segments at the end of compound curves static int NUM_SEGMENTS = 64; - curve_segment_evaluator cse(length_unit_, inst->ParentCurve(), inst->SegmentStart(), inst->SegmentLength()); + curve_segment_evaluator cse(this, length_unit_, inst->ParentCurve(), inst->SegmentStart(), inst->SegmentLength()); boost::mpl::for_each>(std::ref(cse)); std::vector polygon; @@ -420,6 +382,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { //const auto& transformation_matrix = taxonomy::cast(map(inst->Placement()))->ccomponents(); auto transformation_matrix = taxonomy::cast(map(inst->Placement()))->ccomponents(); + // @todo - is there a better way to deal with tolerance and "nearly zero" values? auto length = cse.length(); if (0.001 < fabs(length)) { From f9d259a3d85177649992e0248fd2daa183a50093 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Fri, 29 Sep 2023 11:48:26 -0700 Subject: [PATCH 10/15] Removes redundant function evaluation in numeric integration --- src/ifcgeom/mapping/IfcCurveSegment.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index 88d02f527b..55d5edda80 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -41,13 +41,16 @@ namespace { double area = 0; double h = (b - a) / n; + auto x1 = a; + auto f1 = fn(x1); for (auto i = 1; i <= n; i++) { - auto x1 = a + h * (i - 1); auto x2 = a + h * i; - auto f1 = fn(x1); auto f2 = fn(x2); area += h * (f1 + f2) / 2.0; + + x1 = x2; + f1 = f2; } return area; } From 00d8ffeba09ede67168c3be2f4ad4d53fe30cc2c Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Fri, 29 Sep 2023 11:49:43 -0700 Subject: [PATCH 11/15] Reverts to Taylor Series approximation of Clothoid curve. If you compare the results generated by the numeric integration and the approximation, they are very different. The approximation compares will with other software so there is likely a problem with the equation of the clothoid --- src/ifcgeom/mapping/IfcCurveSegment.cpp | 160 ++++++++++++------------ 1 file changed, 80 insertions(+), 80 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index 55d5edda80..397ac5d2a3 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -98,64 +98,64 @@ public: } // Clothoid using Taylor Series approximation -//#ifdef SCHEMA_HAS_IfcClothoid -// // Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes -// void operator()(IfcSchema::IfcClothoid* c) { -// // @todo verify -// auto sign = [](double v)->int{return v < 0 ? -1 : (0 < v ? 1 : 0); }; -// auto sign_s = sign(start_); -// auto sign_l = sign(length_); -// double L = 0; -// if (sign_s == 0) L = fabs(length_); -// else if (sign_s == sign_l) L = fabs(start_ + length_); -// else L = fabs(start_); -// -// auto A = c->ClothoidConstant(); -// auto R = A * A / L; -// auto RL = (A < 0 ? -1.0 : 1.0) * R * L; -// -// auto position = c->Position(); -// auto placement = position->as(); -// auto ref_direction = placement->RefDirection(); -// double theta = 0.0; // angle the circle's placement X-axis makes with respect to global X axis -// if (ref_direction) -// { -// auto dr = ref_direction->DirectionRatios(); -// auto dx = dr[0]; -// auto dy = dr[1]; -// theta = atan2(dy, dx); -// } -// -// auto C = placement->Location(); -// if (!C->as()) -// { -// throw std::runtime_error("Only IfcCartesianPoint is supported for center of IfcCircle"); -// // @todo add support for other IfcPoint subtypes -// } -// auto Cx = C->as()->Coordinates()[0]; -// auto Cy = C->as()->Coordinates()[1]; -// -// eval_ = [RL,Cx,Cy,theta](double u) { -// // coordinate along clothoid is local coordinates -// auto xterm_1 = u; -// auto xterm_2 = std::pow(u, 5) / (40 * std::pow(RL, 2)); -// auto xterm_3 = std::pow(u, 9) / (3456 * std::pow(RL, 4)); -// auto xterm_4 = std::pow(u, 13) / (599040 * std::pow(RL, 6)); -// auto xl = xterm_1 - xterm_2 + xterm_3 - xterm_4; -// -// auto yterm_1 = std::pow(u, 3) / (6 * RL); -// auto yterm_2 = std::pow(u, 7) / (336 * std::pow(RL, 3)); -// auto yterm_3 = std::pow(u, 11) / (42240 * std::pow(RL, 5)); -// auto yterm_4 = std::pow(u, 15) / (9676800 * std::pow(RL, 7)); -// auto yl = yterm_1 - yterm_2 + yterm_3 - yterm_4; -// -// // transform point into clothoid's coodinate system -// auto x = xl * cos(theta) - yl * sin(theta) + Cx; -// auto y = xl * sin(theta) + yl * cos(theta) + Cy; -// return Eigen::Vector3d(x, y, 0.0); -// }; -// } -//#endif +#ifdef SCHEMA_HAS_IfcClothoid + // Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes + void operator()(IfcSchema::IfcClothoid* c) { + // @todo verify + auto sign = [](double v)->int{return v < 0 ? -1 : (0 < v ? 1 : 0); }; + auto sign_s = sign(start_); + auto sign_l = sign(length_); + double L = 0; + if (sign_s == 0) L = fabs(length_); + else if (sign_s == sign_l) L = fabs(start_ + length_); + else L = fabs(start_); + + auto A = c->ClothoidConstant(); + auto R = A * A / L; + auto RL = (A < 0 ? -1.0 : 1.0) * R * L; + + auto position = c->Position(); + auto placement = position->as(); + auto ref_direction = placement->RefDirection(); + double theta = 0.0; // angle the circle's placement X-axis makes with respect to global X axis + if (ref_direction) + { + auto dr = ref_direction->DirectionRatios(); + auto dx = dr[0]; + auto dy = dr[1]; + theta = atan2(dy, dx); + } + + auto C = placement->Location(); + if (!C->as()) + { + throw std::runtime_error("Only IfcCartesianPoint is supported for center of IfcCircle"); + // @todo add support for other IfcPoint subtypes + } + auto Cx = C->as()->Coordinates()[0]; + auto Cy = C->as()->Coordinates()[1]; + + eval_ = [RL,Cx,Cy,theta](double u) { + // coordinate along clothoid is local coordinates + auto xterm_1 = u; + auto xterm_2 = std::pow(u, 5) / (40 * std::pow(RL, 2)); + auto xterm_3 = std::pow(u, 9) / (3456 * std::pow(RL, 4)); + auto xterm_4 = std::pow(u, 13) / (599040 * std::pow(RL, 6)); + auto xl = xterm_1 - xterm_2 + xterm_3 - xterm_4; + + auto yterm_1 = std::pow(u, 3) / (6 * RL); + auto yterm_2 = std::pow(u, 7) / (336 * std::pow(RL, 3)); + auto yterm_3 = std::pow(u, 11) / (42240 * std::pow(RL, 5)); + auto yterm_4 = std::pow(u, 15) / (9676800 * std::pow(RL, 7)); + auto yl = yterm_1 - yterm_2 + yterm_3 - yterm_4; + + // transform point into clothoid's coodinate system + auto x = xl * cos(theta) - yl * sin(theta) + Cx; + auto y = xl * sin(theta) + yl * cos(theta) + Cy; + return Eigen::Vector3d(x, y, 0.0); + }; + } +#endif void set_spiral_functor(mapping* mapping,IfcSchema::IfcSpiral* s, std::function signX,std::function fnX, std::function signY, std::function fnY) { @@ -189,28 +189,28 @@ public: } // Clothoid using numerical integration -#ifdef SCHEMA_HAS_IfcClothoid -// Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes - void operator()(IfcSchema::IfcClothoid* c) { - - auto A = c->ClothoidConstant(); - - // the integration is for the +X, +Y quadrant - need to adjust the signs of the resulting X and Y values - // so that the results are in the correct quadrant. - // A > 0 and u > 0 -> +X, +Y - // A < 0 and u > 0 -> +X, -Y - // A > 0 and u < 0 -> -X, -Y - // A < 0 and u < 0 -> -X, +Y - // X depends only on u, Y depends on u and A. - auto sign = [](double v)->int {return v < 0 ? -1 : 1; }; // returns -1 or 1 - auto sign_x = [sign](double t) {return sign(t); }; - auto sign_y = [sign, A](double t) {return sign(t) == sign(A) ? 1.0 : -1.0; }; - auto fn_x = [A](double t)->double {return A * sqrt(PI) * cos(PI * A * t * t / (2 * fabs(A))); }; - auto fn_y = [A](double t)->double {return A * sqrt(PI) * sin(PI * A * t * t / (2 * fabs(A))); }; - - set_spiral_functor(mapping_,c->as(), sign_x, fn_x, sign_y, fn_y); - } -#endif +//#ifdef SCHEMA_HAS_IfcClothoid +//// Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes +// void operator()(IfcSchema::IfcClothoid* c) { +// +// auto A = c->ClothoidConstant(); +// +// // the integration is for the +X, +Y quadrant - need to adjust the signs of the resulting X and Y values +// // so that the results are in the correct quadrant. +// // A > 0 and u > 0 -> +X, +Y +// // A < 0 and u > 0 -> +X, -Y +// // A > 0 and u < 0 -> -X, -Y +// // A < 0 and u < 0 -> -X, +Y +// // X depends only on u, Y depends on u and A. +// auto sign = [](double v)->int {return v < 0 ? -1 : 1; }; // returns -1 or 1 +// auto sign_x = [sign](double t) {return sign(t); }; +// auto sign_y = [sign, A](double t) {return sign(t) == sign(A) ? 1.0 : -1.0; }; +// auto fn_x = [A](double t)->double {return A * sqrt(PI) * cos(PI * A * t * t / (2 * fabs(A))); }; +// auto fn_y = [A](double t)->double {return A * sqrt(PI) * sin(PI * A * t * t / (2 * fabs(A))); }; +// +// set_spiral_functor(mapping_,c->as(), sign_x, fn_x, sign_y, fn_y); +// } +//#endif #ifdef SCHEMA_HAS_IfcSecondOrderPolynomialSpiral void operator()(IfcSchema::IfcSecondOrderPolynomialSpiral* s) From c3aeae86616e7889fe006b79acb78eb3b3f000b2 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Fri, 29 Sep 2023 13:27:23 -0700 Subject: [PATCH 12/15] Implements recommendation to use map_impl(IfcPolyline*) --- src/ifcgeom/mapping/IfcCurveSegment.cpp | 27 ++++++++++--------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index 397ac5d2a3..ddc7f969b0 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -261,6 +261,8 @@ public: void operator()(IfcSchema::IfcPolyline* pl) { + auto points = taxonomy::cast(mapping_->map_impl(pl)); + struct Range { double u_start; @@ -271,29 +273,23 @@ public: using Function = std::function(double u)>; std::map fns; - auto p = pl->Points(); - if (p->size() < 2) - { - throw std::runtime_error("invalid polyline - must have at least 2 points"); // this should never happen, but just in case it does - } - 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 iter = p->begin(); - auto end = p->end(); + auto iter = points->children.begin(); + auto end = points->children.end(); auto last = std::prev(end); - auto p1 = *(iter++); auto u = 0.0; for (; iter != end; iter++) { - auto p2 = *iter; + auto edge(*iter); + auto& start_point = boost::get(edge->start); + auto p1x = start_point->components_->x(); + auto p1y = start_point->components_->y(); - auto p1x = p1->Coordinates()[0]; - auto p1y = p1->Coordinates()[1]; - - auto p2x = p2->Coordinates()[0]; - auto p2y = p2->Coordinates()[1]; + auto& end_point = boost::get(edge->end); + auto p2x = end_point->components_->x(); + auto p2y = end_point->components_->y(); auto dx = p2x - p1x; auto dy = p2y - p1y; @@ -311,7 +307,6 @@ public: fns.insert(std::make_pair(Range{ u, u + l,iter == last ? end_compare : std_compare }, fn)); - p1 = p2; u = u + l; } From 7b0e2c52ac0f4cee1260ae8252637a3f324dadfc Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Mon, 2 Oct 2023 16:30:39 -0700 Subject: [PATCH 13/15] Manually merged in aothms revisions for gradient curves and implicit_item taxonomy --- src/ifcgeom/AbstractKernel.h | 1 + src/ifcgeom/mapping/IfcCompositeCurve.cpp | 52 +- src/ifcgeom/mapping/IfcCurveSegment.cpp | 288 +-- src/ifcgeom/mapping/IfcGradientCurve.cpp | 76 + src/ifcgeom/mapping/bind_convert_impl.i | 4 +- src/ifcgeom/mapping/mapping.cpp | 22 +- src/ifcgeom/mapping/mapping.h | 2 + src/ifcgeom/mapping/mapping.i | 3 + src/ifcgeom/taxonomy.cpp | 120 +- src/ifcgeom/taxonomy.h | 2333 +++++++++++---------- 10 files changed, 1565 insertions(+), 1336 deletions(-) create mode 100644 src/ifcgeom/mapping/IfcGradientCurve.cpp diff --git a/src/ifcgeom/AbstractKernel.h b/src/ifcgeom/AbstractKernel.h index 34eb5bde8a..a2c99ef441 100644 --- a/src/ifcgeom/AbstractKernel.h +++ b/src/ifcgeom/AbstractKernel.h @@ -54,6 +54,7 @@ namespace ifcopenshell { namespace geometry { namespace kernels { virtual bool convert_impl(const taxonomy::surface_curve_sweep::ptr, IfcGeom::ConversionResults&) { throw std::runtime_error("Not implemented"); } virtual bool convert_impl(const taxonomy::loft::ptr, IfcGeom::ConversionResults&) { throw std::runtime_error("Not implemented"); } virtual bool convert_impl(const taxonomy::collection::ptr, IfcGeom::ConversionResults&); + virtual bool convert_impl(const taxonomy::piecewise_function::ptr item, IfcGeom::ConversionResults& cs) { return convert(item->evaluate(), cs); } /* virtual void set_offset(const std::array &p_offset); diff --git a/src/ifcgeom/mapping/IfcCompositeCurve.cpp b/src/ifcgeom/mapping/IfcCompositeCurve.cpp index 9f5221fe9f..008f6375ef 100644 --- a/src/ifcgeom/mapping/IfcCompositeCurve.cpp +++ b/src/ifcgeom/mapping/IfcCompositeCurve.cpp @@ -23,6 +23,7 @@ using namespace ifcopenshell::geometry; taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) { auto loop = taxonomy::make(); + auto pwf = taxonomy::make(); #ifdef SCHEMA_HAS_IfcSegment // 4x3 @@ -30,7 +31,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) { #else IfcSchema::IfcCompositeCurveSegment::list::ptr segments = inst->Segments(); #endif - + for (auto& segment : *segments) { if (segment->as() && segment->as()->ParentCurve()->as()) { Logger::Notice("Infinite IfcLine used as ParentCurve of segment, treating as a segment", segment); @@ -47,14 +48,16 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) { e->orientation_2.reset(segment->as()->SameSense()); loop->children.push_back(e); - } else if (segment->as()) { + } + else if (segment->as()) { auto crv = map(segment->as()->ParentCurve()); if (crv) { if (crv->kind() == taxonomy::EDGE) { auto ecrv = taxonomy::cast(crv); ecrv->orientation_2.reset(segment->as()->SameSense()); loop->children.push_back(ecrv); - } else if (crv->kind() == taxonomy::LOOP) { + } + else if (crv->kind() == taxonomy::LOOP) { if (!segment->as()->SameSense()) { crv->reverse(); } @@ -66,31 +69,30 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* inst) { } #ifdef SCHEMA_HAS_IfcCurveSegment else if (segment->as()) { + // @todo check that we don't get a mixture of implicit and explicit definitions auto crv = map(segment->as()); - auto crv_as_loop = taxonomy::cast(crv); - - // The end of the previous segment must be at the same location as the start of this segment - // @todo - need to apply some tolerancing - if (!loop->children.empty() and !crv_as_loop->children.empty() - and - boost::get(loop->children.back()->end)->components() != boost::get(crv_as_loop->children.front()->start)->components()) - { - std::ostringstream os; - auto& prev = boost::get(loop->children.back()->end)->components(); - auto& next = boost::get(crv_as_loop->children.front()->start)->components(); - os << "Common points are not continuous: (" << prev.x() << ", " << prev.y() << ", " << prev.z() << ")" << " " << "(" << next.x() << ", " << next.y() << ", " << next.z() << ")" << std::endl; - Logger::Notice(os.str()); + if (crv->kind() == taxonomy::LOOP) { + for (auto& s : taxonomy::cast(crv)->children) { + loop->children.push_back(s); + } + } + else if (crv->kind() == taxonomy::PIECEWISE_FUNCTION) { + auto seg = taxonomy::cast(crv); + pwf->spans.insert(pwf->spans.end(), seg->spans.begin(), seg->spans.end()); } - - loop->children.insert(loop->children.end(), crv_as_loop->children.begin(), crv_as_loop->children.end()); } #endif } - aggregate_of_instance::ptr profile = inst->data().getInverse(&IfcSchema::IfcProfileDef::Class(), -1); - const bool force_close = profile && profile->size() > 0; - loop->closed = force_close; - return loop; + if (pwf->spans.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; + return loop; + } + else { + return pwf; + } } /* @@ -110,7 +112,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* l, TopoDS_Wi TopTools_ListOfShape converted_segments; - + for (auto it = segments->begin(); it != segments->end(); ++it) { if (!(*it)->declaration().is(IfcSchema::IfcCompositeCurveSegment::Class())) { @@ -121,7 +123,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* l, TopoDS_Wi IfcSchema::IfcCurve* curve = ((IfcSchema::IfcCompositeCurveSegment*)(*it))->ParentCurve(); // The type of ParentCurve is IfcCurve, but the documentation says: - // ParentCurve: The *bounded curve* which defines the geometry of the segment. + // ParentCurve: The *bounded curve* which defines the geometry of the segment. // At least let's exclude IfcLine as an infinite linear segment // definitely does not make any sense. TopoDS_Wire segment; @@ -178,4 +180,4 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCompositeCurve* l, TopoDS_Wi return true; } -*/ +*/ \ No newline at end of file diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index ddc7f969b0..a0c78d5a78 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -27,36 +27,11 @@ using namespace ifcopenshell::geometry; #include #include +#include // @todo use std::numbers::pi when upgrading to C++ 20 -#define PI 3.1415926535897932384626433832795 +static const double PI = boost::math::constants::pi(); - -namespace -{ - // trapezoid rule integration - // @todo is there a well established math library we can use instead of - // creating our own integrator? - double integrate(double a, double b, unsigned n, std::function fn) - { - double area = 0; - double h = (b - a) / n; - auto x1 = a; - auto f1 = fn(x1); - for (auto i = 1; i <= n; i++) - { - auto x2 = a + h * i; - auto f2 = fn(x2); - area += h * (f1 + f2) / 2.0; - - x1 = x2; - f1 = f2; - } - return area; - } -} - -// types of entities that can be IfcCurveSegment.ParentCurve typedef boost::mpl::vector< IfcSchema::IfcLine #ifdef SCHEMA_HAS_IfcClothoid @@ -69,40 +44,83 @@ typedef boost::mpl::vector< , IfcSchema::IfcCircle > curve_seg_types; +enum segment_type_t { + ST_HORIZONTAL, ST_VERTICAL, ST_CANT +}; + class curve_segment_evaluator { private: mapping* mapping_; double length_unit_; double start_; double length_; + segment_type_t segment_type_; IfcSchema::IfcCurve* curve_; - std::optional> eval_; + std::optional> eval_; public: // First constructor, takes parameters from IfcCurveSegment - curve_segment_evaluator(mapping* mapping,double length_unit, IfcSchema::IfcCurve* curve, IfcSchema::IfcCurveMeasureSelect* st, IfcSchema::IfcCurveMeasureSelect* le) + 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 - + 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; } -// Clothoid using Taylor Series approximation + void set_spiral_functor(mapping* mapping_,IfcSchema::IfcSpiral* 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 binary_sign = [](double v)->int {return v < 0 ? -1 : (0 < v ? 1 : 0); }; // returns -1, 0, or 1 + 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 + + //const auto& transformation_matrix = taxonomy::cast(mapping_->map(s->Position()))->ccomponents(); + auto transformation_matrix = taxonomy::cast(mapping_->map(s->Position()))->ccomponents(); + + eval_ = [L, transformation_matrix, signX, fnX, signY, fnY](double u) { + using boost::math::quadrature::trapezoidal; + + // integration limits, integrate from a to b + // from 8.9.3.19.1, integration limits are 0.0 to u where u is a normalized parameter + auto a = 0.0; + auto b = fabs(u / L); + + // @todo where to plug this in? + // auto n = 10; // use 10 steps in the numeric integration + + auto x = signX(u) * trapezoidal(fnX, a, b); + auto y = signY(u) * trapezoidal(fnY, a, b); + + // transform point into spiral's coodinate system + auto result = transformation_matrix * Eigen::Vector4d(x, y, 0.0, 1.0); + Eigen::VectorXd vec(4); + vec << result(0), result(1), 0.0, 1.0; + return vec; + }; + } + + + // Clothoid using Taylor Series approximation #ifdef SCHEMA_HAS_IfcClothoid // Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes void operator()(IfcSchema::IfcClothoid* c) { // @todo verify - auto sign = [](double v)->int{return v < 0 ? -1 : (0 < v ? 1 : 0); }; + auto sign = [](double v)->int {return v < 0 ? -1 : (0 < v ? 1 : 0); }; auto sign_s = sign(start_); auto sign_l = sign(length_); double L = 0; @@ -114,81 +132,33 @@ public: auto R = A * A / L; auto RL = (A < 0 ? -1.0 : 1.0) * R * L; - auto position = c->Position(); - auto placement = position->as(); - auto ref_direction = placement->RefDirection(); - double theta = 0.0; // angle the circle's placement X-axis makes with respect to global X axis - if (ref_direction) - { - auto dr = ref_direction->DirectionRatios(); - auto dx = dr[0]; - auto dy = dr[1]; - theta = atan2(dy, dx); - } + //const auto& transformation_matrix = taxonomy::cast(mapping_->map(c->Position()))->ccomponents(); + auto transformation_matrix = taxonomy::cast(mapping_->map(c->Position()))->ccomponents(); - auto C = placement->Location(); - if (!C->as()) - { - throw std::runtime_error("Only IfcCartesianPoint is supported for center of IfcCircle"); - // @todo add support for other IfcPoint subtypes - } - auto Cx = C->as()->Coordinates()[0]; - auto Cy = C->as()->Coordinates()[1]; - - eval_ = [RL,Cx,Cy,theta](double u) { + eval_ = [RL, transformation_matrix](double u) { // coordinate along clothoid is local coordinates auto xterm_1 = u; auto xterm_2 = std::pow(u, 5) / (40 * std::pow(RL, 2)); auto xterm_3 = std::pow(u, 9) / (3456 * std::pow(RL, 4)); auto xterm_4 = std::pow(u, 13) / (599040 * std::pow(RL, 6)); - auto xl = xterm_1 - xterm_2 + xterm_3 - xterm_4; + auto x = xterm_1 - xterm_2 + xterm_3 - xterm_4; auto yterm_1 = std::pow(u, 3) / (6 * RL); auto yterm_2 = std::pow(u, 7) / (336 * std::pow(RL, 3)); auto yterm_3 = std::pow(u, 11) / (42240 * std::pow(RL, 5)); auto yterm_4 = std::pow(u, 15) / (9676800 * std::pow(RL, 7)); - auto yl = yterm_1 - yterm_2 + yterm_3 - yterm_4; - - // transform point into clothoid's coodinate system - auto x = xl * cos(theta) - yl * sin(theta) + Cx; - auto y = xl * sin(theta) + yl * cos(theta) + Cy; - return Eigen::Vector3d(x, y, 0.0); - }; - } -#endif - - void set_spiral_functor(mapping* mapping,IfcSchema::IfcSpiral* 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 binary_sign = [](double v)->int {return v < 0 ? -1 : (0 < v ? 1 : 0); }; // returns -1, 0, or 1 - 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 - - //const auto& transformation_matrix = taxonomy::cast(mapping->map(s->Position()))->ccomponents(); - auto transformation_matrix = taxonomy::cast(mapping->map(s->Position()))->ccomponents(); - - eval_ = [L, transformation_matrix, signX, fnX, signY, fnY](double u) { - // integration limits, integrate from a to b - // from 8.9.3.19.1, integration limits are 0.0 to u where u is a normalized parameter - auto a = 0.0; - auto b = fabs(u / L); - - auto n = 10; // use 10 steps in the numeric integration - - auto x = signX(u)*integrate(a, b, n, fnX); - auto y = signY(u)*integrate(a, b, n, fnY); + auto y = yterm_1 - yterm_2 + yterm_3 - yterm_4; // transform point into clothoid's coodinate system auto result = transformation_matrix * Eigen::Vector4d(x, y, 0.0, 1.0); - return Eigen::Vector3d(result(0),result(1),result(2)); + Eigen::VectorXd vec(4); + vec << result(0), result(1), 0.0, 1.0; + return vec; }; } +#endif -// Clothoid using numerical integration + // Clothoid using numerical integration //#ifdef SCHEMA_HAS_IfcClothoid //// Then initialize Function(double) -> Vector3, by means of IfcCurve subtypes // void operator()(IfcSchema::IfcClothoid* c) { @@ -235,7 +205,7 @@ public: auto fn_x = [theta](double t)->double {return cos(theta(t)); }; auto fn_y = [theta](double t)->double {return sin(theta(t)); }; - set_spiral_functor(mapping_, s->as(), sign_x, fn_x, sign_y, fn_y); + set_spiral_functor(mapping_,s->as(), sign_x, fn_x, sign_y, fn_y); } #endif @@ -255,14 +225,14 @@ public: // transform point into circle's coodinate system auto result = transformation_matrix * Eigen::Vector4d(x, y, 0.0, 1.0); - return Eigen::Vector3d(result(0), result(1), result(2)); + Eigen::VectorXd vec(4); + vec << result(0), result(1), 0.0, 1.0; + return vec; }; } void operator()(IfcSchema::IfcPolyline* pl) { - auto points = taxonomy::cast(mapping_->map_impl(pl)); - struct Range { double u_start; @@ -273,23 +243,29 @@ public: using Function = std::function(double u)>; std::map fns; - 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 p = pl->Points(); + if (p->size() < 2) + { + throw std::runtime_error("invalid polyline - must have at least 2 points"); // this should never happen, but just in case it does + } - auto iter = points->children.begin(); - auto end = points->children.end(); + 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 iter = p->begin(); + auto end = p->end(); auto last = std::prev(end); + auto p1 = *(iter++); auto u = 0.0; for (; iter != end; iter++) { - auto edge(*iter); - auto& start_point = boost::get(edge->start); - auto p1x = start_point->components_->x(); - auto p1y = start_point->components_->y(); + auto p2 = *iter; - auto& end_point = boost::get(edge->end); - auto p2x = end_point->components_->x(); - auto p2y = end_point->components_->y(); + auto p1x = p1->Coordinates()[0]; + auto p1y = p1->Coordinates()[1]; + + auto p2x = p2->Coordinates()[0]; + auto p2y = p2->Coordinates()[1]; auto dx = p2x - p1x; auto dy = p2y - p1y; @@ -307,6 +283,7 @@ public: fns.insert(std::make_pair(Range{ u, u + l,iter == last ? end_compare : std_compare }, fn)); + p1 = p2; u = u + l; } @@ -316,12 +293,14 @@ public: auto [u_start, u_end, compare] = fn.first; return compare(u_start, u, u_end); }); - + if (iter == fns.end()) throw std::runtime_error("invalid distance from start"); // this should never happen, but just in case it does - + auto [u_start, u_end, compare] = iter->first; - auto [x,y] = (iter->second)(u - u_start); // (u - u_start) is distance from start of this segment of the polyline - return Eigen::Vector3d(x, y, 0); + auto [x, y] = (iter->second)(u - u_start); // (u - u_start) is distance from start of this segment of the polyline + Eigen::VectorXd vec(4); + vec << x, y, 0.0, 1.0; + return vec; }; } @@ -336,11 +315,27 @@ public: auto dx = dr[0] / m; auto dy = dr[1] / m; - eval_ = [px, py, dx, dy](double u) { - auto x = px + u * dx; - auto y = py + u * dy; - return Eigen::Vector3d(x, y, 0); - }; + if (segment_type_ == ST_HORIZONTAL) { + + eval_ = [px, py, dx, dy](double u) { + auto x = px + u * dx; + auto y = py + u * dy; + Eigen::VectorXd vec(4); + vec << x, y, 0.0, 1.0; + return vec; + }; + + } + else if (segment_type_ == ST_VERTICAL) { + + eval_ = [py, dy](double u) { + auto z = py + u * dy; + Eigen::VectorXd vec(4); + vec << 0.0, 0.0, z, 1.0; + return vec; + }; + + } } // Take the boost::type value from mpl::for_each and test it against our curve instance @@ -352,10 +347,11 @@ public: } // Then, with function populated based on IfcCurve subtype, we can evaluate to points - Eigen::Vector3d operator()(double u) { + Eigen::VectorXd operator()(double u) { if (eval_) { return (*eval_)((u + start_) * length_unit_); - } else { + } + else { throw std::runtime_error(curve_->declaration().name() + " not implemented"); } } @@ -363,30 +359,73 @@ public: double length() const { return length_; } + + const std::optional>& evaluation_function() const { + return eval_; + } }; taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { - // @todo fixed number of segments or fixed interval? - // @todo placement // @todo figure out what to do with the zero length segments at the end of compound curves - static int NUM_SEGMENTS = 64; - curve_segment_evaluator cse(this, length_unit_, inst->ParentCurve(), inst->SegmentStart(), inst->SegmentLength()); + bool is_horizontal = false; + bool is_vertical = false; + bool is_cant = false; + + { + aggregate_of_instance::ptr segment_owners = inst->data().getInverse(&IfcSchema::IfcCompositeCurve::Class(), 0); + if (segment_owners) { + for (auto& cc : *segment_owners) { + if (cc->as()) { + is_cant = true; + } + else if (cc->as()) { + is_vertical = true; + } + else { + is_horizontal = true; + } + } + } + } + + if ((is_horizontal + is_vertical + is_cant) != 1) { + // We have to choose the correct functor based on usage. We can't + // support multiple, because we don't know the caller at this point. + return nullptr; + } + + 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()); boost::mpl::for_each>(std::ref(cse)); - std::vector polygon; + auto fn = *cse.evaluation_function(); + auto length = fabs(cse.length()); // @todo - for some reason this isn't working, the matrix gets all messed up //const auto& transformation_matrix = taxonomy::cast(map(inst->Placement()))->ccomponents(); auto transformation_matrix = taxonomy::cast(map(inst->Placement()))->ccomponents(); - // @todo - is there a better way to deal with tolerance and "nearly zero" values? + auto fn_transformed = [fn, transformation_matrix](double u) { + return transformation_matrix * fn(u); + }; + + // @todo it might be suboptimal that we no longer have the spans now + auto pwf = taxonomy::make(); + pwf->spans.push_back({ length, fn_transformed }); + return pwf; + + /* + static int NUM_SEGMENTS = 64; + std::vector polygon; + auto length = cse.length(); if (0.001 < fabs(length)) { for (int i = 0; i <= NUM_SEGMENTS; ++i) { auto u = length * i / NUM_SEGMENTS; - + auto p = cse(u); auto result = transformation_matrix * Eigen::Vector4d(p(0),p(1),p(2), 1.); polygon.push_back(taxonomy::make(result(0),result(1),result(2))); @@ -394,6 +433,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { } return polygon_from_points(polygon); + */ } #endif \ No newline at end of file diff --git a/src/ifcgeom/mapping/IfcGradientCurve.cpp b/src/ifcgeom/mapping/IfcGradientCurve.cpp new file mode 100644 index 0000000000..9901c829bb --- /dev/null +++ b/src/ifcgeom/mapping/IfcGradientCurve.cpp @@ -0,0 +1,76 @@ +/******************************************************************************** + * * + * This file is part of IfcOpenShell. * + * * + * IfcOpenShell is free software: you can redistribute it and/or modify * + * it under the terms of the Lesser GNU General Public License as published by * + * the Free Software Foundation, either version 3.0 of the License, or * + * (at your option) any later version. * + * * + * IfcOpenShell is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * Lesser GNU General Public License for more details. * + * * + * You should have received a copy of the Lesser GNU General Public License * + * along with this program. If not, see . * + * * + ********************************************************************************/ + +#include "mapping.h" +#define mapping POSTFIX_SCHEMA(mapping) +using namespace ifcopenshell::geometry; + +#ifdef SCHEMA_HAS_IfcGradientCurve + +taxonomy::ptr mapping::map_impl(const IfcSchema::IfcGradientCurve* inst) { + auto horizontal = taxonomy::cast(map(inst->BaseCurve())); + auto vertical = taxonomy::make(); + + auto segments = inst->Segments(); + + for (auto& segment : *segments) { + if (segment->as()) { + // @todo check that we don't get a mixture of implicit and explicit definitions + auto crv = map(segment->as()); + if (crv->kind() == taxonomy::PIECEWISE_FUNCTION) { + auto seg = taxonomy::cast(crv); + vertical->spans.insert(vertical->spans.end(), seg->spans.begin(), seg->spans.end()); + } else { + Logger::Error("Unsupported"); + return nullptr; + } + } else { + Logger::Error("Unsupported"); + return nullptr; + } + } + + // @todo does this really make sense? + auto composition = [horizontal, vertical](double u) { + auto xy = horizontal->evaluate(u); + auto z = vertical->evaluate(u); + Eigen::VectorXd vec(3); + vec << xy(0), xy(1), z(0); + return vec; + }; + + // @todo where do we get the startdistalong from @civilx64's code? + std::array both = { horizontal , vertical }; + double min_length = std::numeric_limits::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(); + pwf->spans.push_back({ min_length, composition }); + return pwf; +} + +#endif \ No newline at end of file diff --git a/src/ifcgeom/mapping/bind_convert_impl.i b/src/ifcgeom/mapping/bind_convert_impl.i index e9c6f3f6e2..229b7487ed 100644 --- a/src/ifcgeom/mapping/bind_convert_impl.i +++ b/src/ifcgeom/mapping/bind_convert_impl.i @@ -5,7 +5,7 @@ #define BIND(T) \ if (inst->as()) { \ try { \ - taxonomy::ptr item = map_impl(inst->as()); \ + item = map_impl(inst->as()); \ if (item != nullptr) { \ item->instance = inst; \ try { \ @@ -24,11 +24,9 @@ } else {\ Logger::Message(Logger::LOG_ERROR,"Failed to convert:", inst);\ } \ - return item; \ } catch (const std::exception& e) { \ Logger::Message(Logger::LOG_ERROR, std::string(e.what()) + "\nFailed to convert:", inst); \ } \ - return nullptr; \ } #include "mapping.i" diff --git a/src/ifcgeom/mapping/mapping.cpp b/src/ifcgeom/mapping/mapping.cpp index 7a98125191..2cd9b05730 100644 --- a/src/ifcgeom/mapping/mapping.cpp +++ b/src/ifcgeom/mapping/mapping.cpp @@ -490,10 +490,26 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcStyledItem* inst) { taxonomy::ptr mapping::map(const IfcBaseInterface* inst) { - // std::wcout << inst->data().toString().c_str() << std::endl; + auto iden = inst->as()->identity(); + auto it = cache_.find(iden); + if (it != cache_.end()) { + return it->second; + } + taxonomy::ptr item = nullptr; + + // @todo we should check whether there is a notice performance impact on the large sequence + // of if-statements and whether a switch on e.g inst->declaration()->index_in_schema() + // isn't more efficient (which would disable inheritance though). + #include "bind_convert_impl.i" - Logger::Message(Logger::LOG_ERROR, "No operation defined for:", inst); - return nullptr; + + if (item) { + cache_.insert({ iden, item }); + } + else { + Logger::Message(Logger::LOG_ERROR, "No operation defined for:", inst); + } + return item; } namespace { diff --git a/src/ifcgeom/mapping/mapping.h b/src/ifcgeom/mapping/mapping.h index 248fc4aca2..fa825788fc 100644 --- a/src/ifcgeom/mapping/mapping.h +++ b/src/ifcgeom/mapping/mapping.h @@ -22,6 +22,8 @@ namespace geometry { double length_unit_, angle_unit_; std::string length_unit_name_; + std::map cache_; + const IfcParse::declaration* placement_rel_to_type_; const IfcUtil::IfcBaseEntity* placement_rel_to_instance_; diff --git a/src/ifcgeom/mapping/mapping.i b/src/ifcgeom/mapping/mapping.i index 435366a0c0..57de49a16d 100644 --- a/src/ifcgeom/mapping/mapping.i +++ b/src/ifcgeom/mapping/mapping.i @@ -113,6 +113,9 @@ BIND(IfcEdge); BIND(IfcEdgeLoop); BIND(IfcPolyline); BIND(IfcPolyLoop); +#ifdef SCHEMA_HAS_IfcGradientCurve +BIND(IfcGradientCurve); +#endif BIND(IfcCompositeCurve); BIND(IfcTrimmedCurve); BIND(IfcArbitraryOpenProfileDef); diff --git a/src/ifcgeom/taxonomy.cpp b/src/ifcgeom/taxonomy.cpp index 495209136c..a040ab749f 100644 --- a/src/ifcgeom/taxonomy.cpp +++ b/src/ifcgeom/taxonomy.cpp @@ -1,4 +1,5 @@ #include "taxonomy.h" +#include "profile_helper.h" using namespace ifcopenshell::geometry::taxonomy; @@ -23,9 +24,11 @@ namespace { bool compare(const eigen_base& t, const eigen_base& u) { if (t.components_ == nullptr && u.components_ == nullptr) { return false; - } else if (t.components_ == nullptr && u.components_ != nullptr) { + } + else if (t.components_ == nullptr && u.components_ != nullptr) { return true; - } else if (t.components_ != nullptr && u.components_ == nullptr) { + } + else if (t.components_ != nullptr && u.components_ == nullptr) { return false; } @@ -86,11 +89,14 @@ namespace { int less_to_order_optional(const boost::optional& a, const boost::optional& b) { if (a && b) { return less_to_order(*a, *b); - } else if (!a && !b) { + } + else if (!a && !b) { return 0; - } else if (a) { + } + else if (a) { return 1; - } else { + } + else { return -1; } } @@ -100,7 +106,8 @@ namespace { if (a.which() == 0) { a_lt_b = compare(*boost::get(a), *boost::get(b)); b_lt_a = compare(*boost::get(b), *boost::get(a)); - } else { + } + else { a_lt_b = std::less()(boost::get(a), boost::get(b)); b_lt_a = std::less()(boost::get(b), boost::get(a)); } @@ -144,7 +151,11 @@ namespace { bool compare(const surface_curve_sweep&, const surface_curve_sweep&) { throw std::runtime_error("not implemented"); } - + + bool compare(const piecewise_function&, const piecewise_function&) { + throw std::runtime_error("not implemented"); + } + bool compare(const style& a, const style& b) { const int order[5] = { less_to_order(a.name, b.name), @@ -166,7 +177,8 @@ namespace { auto A = static_cast*>(a); auto B = static_cast*>(b); return compare(*A, *B); - } else { + } + else { return dispatch_comparison::dispatch(a, b); } } @@ -223,23 +235,28 @@ namespace { if (!a_has_basis) { // Finally, equality return false; - } else { + } + else { return less(a.basis, b.basis); } - } else { + } + else { return a_has_basis < b_has_basis; } - } else { + } + else { return end_state == -1; } - } else { + } + else { return start_state == -1; } - } else { + } + else { return std::tie(a.orientation, a_which_start, a_which_end) < std::tie(b.orientation, b_which_start, b_which_end); @@ -262,7 +279,8 @@ namespace { } // Vectors equal, compare matrix (in case of mapped items). return compare(*a.matrix, *b.matrix); - } else { + } + else { return a.children.size() < b.children.size(); } } @@ -314,10 +332,10 @@ ifcopenshell::geometry::taxonomy::solid::ptr ifcopenshell::geometry::create_box( shell->children.push_back(face); std::array points{ - taxonomy::make(x+0, y+0, z+ 0), - taxonomy::make(x+0, y+dy, z+ 0), - taxonomy::make(x+0, y+dy, z+dz), - taxonomy::make(x+0, y+0, z+dz) + taxonomy::make(x + 0, y + 0, z + 0), + taxonomy::make(x + 0, y + dy, z + 0), + taxonomy::make(x + 0, y + dy, z + dz), + taxonomy::make(x + 0, y + 0, z + dz) }; loop->children.push_back(make(points[0], points[1])); @@ -335,10 +353,10 @@ ifcopenshell::geometry::taxonomy::solid::ptr ifcopenshell::geometry::create_box( shell->children.push_back(face); std::array points{ - taxonomy::make(x+dx, y+0, z+ 0), - taxonomy::make(x+dx, y+0, z+dz), - taxonomy::make(x+dx, y+dy, z+dz), - taxonomy::make(x+dx, y+dy, z+ 0) + taxonomy::make(x + dx, y + 0, z + 0), + taxonomy::make(x + dx, y + 0, z + dz), + taxonomy::make(x + dx, y + dy, z + dz), + taxonomy::make(x + dx, y + dy, z + 0) }; loop->children.push_back(make(points[0], points[1])); @@ -356,10 +374,10 @@ ifcopenshell::geometry::taxonomy::solid::ptr ifcopenshell::geometry::create_box( shell->children.push_back(face); std::array points{ - taxonomy::make(x+0, y+0, z+ 0), - taxonomy::make(x+0, y+0, z+dz), - taxonomy::make(x+dx, y+0, z+dz), - taxonomy::make(x+dx, y+0, z+ 0) + taxonomy::make(x + 0, y + 0, z + 0), + taxonomy::make(x + 0, y + 0, z + dz), + taxonomy::make(x + dx, y + 0, z + dz), + taxonomy::make(x + dx, y + 0, z + 0) }; loop->children.push_back(make(points[0], points[1])); @@ -377,10 +395,10 @@ ifcopenshell::geometry::taxonomy::solid::ptr ifcopenshell::geometry::create_box( shell->children.push_back(face); std::array points{ - taxonomy::make(x+ 0, y+dy, z+ 0), - taxonomy::make(x+dx, y+dy, z+ 0), - taxonomy::make(x+dx, y+dy, z+dz), - taxonomy::make(x+ 0, y+dy, z+dz) + taxonomy::make(x + 0, y + dy, z + 0), + taxonomy::make(x + dx, y + dy, z + 0), + taxonomy::make(x + dx, y + dy, z + dz), + taxonomy::make(x + 0, y + dy, z + dz) }; loop->children.push_back(make(points[0], points[1])); @@ -398,10 +416,10 @@ ifcopenshell::geometry::taxonomy::solid::ptr ifcopenshell::geometry::create_box( shell->children.push_back(face); std::array points{ - taxonomy::make(x+ 0, y+ 0, z+0), - taxonomy::make(x+dx, y+ 0, z+0), - taxonomy::make(x+dx, y+dy, z+0), - taxonomy::make(x+ 0, y+dy, z+0) + taxonomy::make(x + 0, y + 0, z + 0), + taxonomy::make(x + dx, y + 0, z + 0), + taxonomy::make(x + dx, y + dy, z + 0), + taxonomy::make(x + 0, y + dy, z + 0) }; loop->children.push_back(make(points[0], points[1])); @@ -419,10 +437,10 @@ ifcopenshell::geometry::taxonomy::solid::ptr ifcopenshell::geometry::create_box( shell->children.push_back(face); std::array points{ - taxonomy::make(x+ 0, y+ 0, z+dz), - taxonomy::make(x+ 0, y+dy, z+dz), - taxonomy::make(x+dx, y+dy, z+dz), - taxonomy::make(x+dx, y+ 0, z+dz) + taxonomy::make(x + 0, y + 0, z + dz), + taxonomy::make(x + 0, y + dy, z + dz), + taxonomy::make(x + dx, y + dy, z + dz), + taxonomy::make(x + dx, y + 0, z + dz) }; loop->children.push_back(make(points[0], points[1])); @@ -434,11 +452,31 @@ ifcopenshell::geometry::taxonomy::solid::ptr ifcopenshell::geometry::create_box( return solid; } +ifcopenshell::geometry::taxonomy::item::ptr ifcopenshell::geometry::taxonomy::piecewise_function::evaluate() const { + // @todo configure resolution + //double length = std::accumulate(spans.begin(), spans.end(), 0.0); // don't know why this doesn't compile + double length = 0.0; + for (auto& s : spans) + length += s.first; + + static const double resolution = 0.5; + std::vector polygon; + + int num_steps = std::ceil(length / resolution); + for (int i = 0; i < num_steps; ++i) { + auto u = resolution * i; + auto p = evaluate(u); + polygon.push_back(taxonomy::make(p(0), p(1), p(2))); + } + + return polygon_from_points(polygon); +} + ifcopenshell::geometry::taxonomy::collection::ptr ifcopenshell::geometry::flatten(taxonomy::collection::ptr deep) { auto flat = make(); ifcopenshell::geometry::visit(deep, [&flat](taxonomy::ptr i) { flat->children.push_back(taxonomy::cast(clone(i))); - }); + }); return flat; } @@ -446,10 +484,10 @@ const std::string& ifcopenshell::geometry::taxonomy::kind_to_string(kinds k) { using namespace std::string_literals; static std::string values[] = { - "matrix4"s, "point3"s, "direction3"s, "line"s, "circle"s, "ellipse"s, "bspline_curve"s, "offset_curve"s, "plane"s, "cylinder"s, "bspline_surface"s, "edge"s, "loop"s, "face"s, "shell"s, "solid"s, "loft"s, "extrusion"s, "revolve"s, "surface_curve_sweep"s, "node"s, "collection"s, "boolean_result"s + "matrix4"s, "point3"s, "direction3"s, "line"s, "circle"s, "ellipse"s, "bspline_curve"s, "offset_curve"s, "plane"s, "cylinder"s, "bspline_surface"s, "edge"s, "loop"s, "face"s, "shell"s, "solid"s, "loft"s, "extrusion"s, "revolve"s, "surface_curve_sweep"s, "node"s, "collection"s, "boolean_result"s, "piecewise_function"s, "colour"s, "style"s, }; return values[k]; } -std::atomic_uint32_t item::counter_(0); +std::atomic_uint32_t item::counter_(0); \ No newline at end of file diff --git a/src/ifcgeom/taxonomy.h b/src/ifcgeom/taxonomy.h index 92dfcf9b75..63cd39665f 100644 --- a/src/ifcgeom/taxonomy.h +++ b/src/ifcgeom/taxonomy.h @@ -27,41 +27,41 @@ namespace ifcopenshell { -namespace geometry { + namespace geometry { -namespace taxonomy { + namespace taxonomy { #ifdef TAXONOMY_USE_SHARED_PTR - template - T clone(T& t) { - return t; - } - template - std::shared_ptr cast(const std::shared_ptr& u); - template - std::shared_ptr dcast(const std::shared_ptr& u); + template + T clone(T& t) { + return t; + } + template + std::shared_ptr cast(const std::shared_ptr& u); + template + std::shared_ptr dcast(const std::shared_ptr& u); #endif #ifdef TAXONOMY_USE_UNIQUE_PTR - // untested currently - template - T clone(T& t) { - return t->clone_(); - } - template - T* cast(const std::unique_ptr& u); - template - T* dcast(const std::unique_ptr& u); + // untested currently + template + T clone(T& t) { + return t->clone_(); + } + template + T* cast(const std::unique_ptr& u); + template + T* dcast(const std::unique_ptr& u); #endif #ifdef TAXONOMY_USE_NAKED_PTR - // untested currently - template - T clone(T& t) { - return t->clone_(); - } - template - T* cast(const U*& u); - template - T* dcast(const U*& u); + // untested currently + template + T clone(T& t) { + return t->clone_(); + } + template + T* cast(const U*& u); + template + T* dcast(const U*& u); #endif #ifdef TAXONOMY_USE_SHARED_PTR @@ -80,1164 +80,1217 @@ typedef item* ptr; \ typedef item const* ptr; #endif -class topology_error : public std::runtime_error { -public: - topology_error() : std::runtime_error("Generic topology error") {} - topology_error(const char* const s) : std::runtime_error(s) {} -}; + class topology_error : public std::runtime_error { + public: + topology_error() : std::runtime_error("Generic topology error") {} + topology_error(const char* const s) : std::runtime_error(s) {} + }; -enum kinds { MATRIX4, POINT3, DIRECTION3, LINE, CIRCLE, ELLIPSE, BSPLINE_CURVE, OFFSET_CURVE, PLANE, CYLINDER, BSPLINE_SURFACE, EDGE, LOOP, FACE, SHELL, SOLID, LOFT, EXTRUSION, REVOLVE, SURFACE_CURVE_SWEEP, NODE, COLLECTION, BOOLEAN_RESULT, COLOUR, STYLE }; + enum kinds { MATRIX4, POINT3, DIRECTION3, LINE, CIRCLE, ELLIPSE, BSPLINE_CURVE, OFFSET_CURVE, PLANE, CYLINDER, BSPLINE_SURFACE, EDGE, LOOP, FACE, SHELL, SOLID, LOFT, EXTRUSION, REVOLVE, SURFACE_CURVE_SWEEP, NODE, COLLECTION, BOOLEAN_RESULT, PIECEWISE_FUNCTION, COLOUR, STYLE }; -const std::string& kind_to_string(kinds k); + const std::string& kind_to_string(kinds k); -struct item { -private: - uint32_t identity_; - static std::atomic_uint32_t counter_; - mutable size_t computed_hash_; -public: - DECLARE_PTR(item) + struct item { + private: + uint32_t identity_; + static std::atomic_uint32_t counter_; + mutable size_t computed_hash_; + public: + DECLARE_PTR(item) - const IfcUtil::IfcBaseInterface* instance; + const IfcUtil::IfcBaseInterface* instance; - boost::optional orientation; + boost::optional orientation; - virtual item* clone_() const = 0; - virtual kinds kind() const = 0; - virtual void print(std::ostream&, int indent=0) const = 0; - virtual void reverse() { throw taxonomy::topology_error(); } - virtual size_t calc_hash() const = 0; - virtual size_t hash() const { - if (computed_hash_) { - return computed_hash_; - } - computed_hash_ = calc_hash(); - if (computed_hash_ == 0) { - computed_hash_++; - } - return computed_hash_; - } + virtual item* clone_() const = 0; + virtual kinds kind() const = 0; + virtual void print(std::ostream&, int indent = 0) const = 0; + virtual void reverse() { throw taxonomy::topology_error(); } + virtual size_t calc_hash() const = 0; + virtual size_t hash() const { + if (computed_hash_) { + return computed_hash_; + } + computed_hash_ = calc_hash(); + if (computed_hash_ == 0) { + computed_hash_++; + } + return computed_hash_; + } - item(const IfcUtil::IfcBaseInterface* instance = nullptr) : identity_(counter_++), instance(instance), computed_hash_(0) {} + item(const IfcUtil::IfcBaseInterface* instance = nullptr) : identity_(counter_++), instance(instance), computed_hash_(0) {} - virtual ~item() {} + virtual ~item() {} - uint32_t identity() const { return identity_; } -}; + uint32_t identity() const { return identity_; } + }; + struct implicit_item : public item { + DECLARE_PTR(implicit_item) + + virtual item::ptr evaluate() const = 0; + }; + + struct piecewise_function : public implicit_item { + DECLARE_PTR(piecewise_function) + + std::vector>> spans; + + void print(std::ostream& o, int indent = 0) const { + o << "piecewise_function" << std::endl; + } + + virtual piecewise_function* clone_() const { return new piecewise_function(*this); } + virtual kinds kind() const { return PIECEWISE_FUNCTION; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(PIECEWISE_FUNCTION), 0); + return boost::hash{}(v); + } + + virtual item::ptr evaluate() const; + + Eigen::VectorXd evaluate(double u) const { + // @todo optimize, assume monotonic evaluation and store last evaluated segment? + for (auto& [length, fn] : spans) { + if (u < length) { + return fn(u); + } + u -= length; + } + } + }; #ifdef TAXONOMY_USE_SHARED_PTR -typedef std::shared_ptr ptr; -typedef std::shared_ptr const_ptr; -template -std::shared_ptr make(Args&&... args) { - return std::make_shared(std::forward(args)...); -} + typedef std::shared_ptr ptr; + typedef std::shared_ptr const_ptr; + template + std::shared_ptr make(Args&&... args) { + return std::make_shared(std::forward(args)...); + } #endif #ifdef TAXONOMY_USE_UNIQUE_PTR -typedef std::uniqe_ptr ptr; -typedef std::uniqe_ptr ptr; -template -std::uniqe_ptr make(Args&&... args) { - return new T(std::forward(args)...)); -} + typedef std::uniqe_ptr ptr; + typedef std::uniqe_ptr ptr; + template + std::uniqe_ptr make(Args&&... args) { + return new T(std::forward(args)...)); + } #endif #ifdef TAXONOMY_USE_NAKED_PTR -typedef item* ptr; -typedef item const* ptr; -template -T* make(Args&&... args) { - return new T(std::forward(args)...)); -} + typedef item* ptr; + typedef item const* ptr; + template + T* make(Args&&... args) { + return new T(std::forward(args)...)); + } #endif -bool less(item::const_ptr, item::const_ptr); + bool less(item::const_ptr, item::const_ptr); -struct less_functor { - bool operator()(item::const_ptr a, item::const_ptr b) const { - return less(a, b); - } -}; + struct less_functor { + bool operator()(item::const_ptr a, item::const_ptr b) const { + return less(a, b); + } + }; -namespace { + namespace { - template - const T& eigen_defaults(); + template + const T& eigen_defaults(); - template <> - const Eigen::Vector3d& eigen_defaults() { - static Eigen::Vector3d identity = Eigen::Vector3d::Zero(); - return identity; - } + template <> + const Eigen::Vector3d& eigen_defaults() { + static Eigen::Vector3d identity = Eigen::Vector3d::Zero(); + return identity; + } - template <> - const Eigen::Matrix4d& eigen_defaults() { - static Eigen::Matrix4d identity = Eigen::Matrix4d::Identity(); - return identity; - } + template <> + const Eigen::Matrix4d& eigen_defaults() { + static Eigen::Matrix4d identity = Eigen::Matrix4d::Identity(); + return identity; + } -} - -template -struct eigen_base { - T* components_; - - eigen_base() { - components_ = nullptr; - } - - eigen_base(const eigen_base& other) { - this->components_ = other.components_ ? new T(*other.components_) : nullptr; - } - - eigen_base(const T& other) { - this->components_ = new T(other); - } - - eigen_base& operator=(const eigen_base& other) { - if (this != &other) { - this->components_ = other.components_ ? new T(*other.components_) : nullptr; - } - return *this; - } - - void print_impl(std::ostream& o, const std::string& class_name, int indent = 0) const { - o << std::string(indent, ' ') << class_name; - if (this->components_) { - int n = T::RowsAtCompileTime * T::ColsAtCompileTime; - for (size_t i = 0; i < n; ++i) { - o << " " << (*components_)(i); } - } - o << std::endl; - } - - virtual ~eigen_base() { - delete this->components_; - } - - const T& ccomponents() const { - if (this->components_) { - return *this->components_; - } else { - return eigen_defaults(); - } - } - - T& components() { - if (!this->components_) { - this->components_ = new T(eigen_defaults()); - } - return *this->components_; - } - - explicit operator bool() const { - return components_; - } - - uint32_t hash_components() const { - size_t h = std::hash{}(T::RowsAtCompileTime); - boost::hash_combine(h, std::hash{}(T::ColsAtCompileTime)); - if (components_) { - for (size_t i = 0; i < components_->size(); ++i) { - auto elem = *(components_->data() + i); - boost::hash_combine(h, std::hash()(elem)); - } - } - return h; - } -}; - -struct matrix4 : public item, public eigen_base { -private: - void init(const Eigen::Vector3d& o, const Eigen::Vector3d& z, const Eigen::Vector3d& x) { - auto X = x.normalized(); - auto Y = z.cross(x).normalized(); - auto Z = z.normalized(); - components_ = new Eigen::Matrix4d; - (*components_) << - X(0), Y(0), Z(0), o(0), - X(1), Y(1), Z(1), o(1), - X(2), Y(2), Z(2), o(2), - 0, 0, 0, 1.; - if (is_identity()) { - // @todo detect this earlier to save us the heapalloc. - delete components_; - components_ = nullptr; - tag = IDENTITY; - } - } -public: - DECLARE_PTR(matrix4) - - enum tag_t { - IDENTITY, AFFINE_WO_SCALE, AFFINE_W_UNIFORM_SCALE, AFFINE_W_NONUNIFORM_SCALE, OTHER - }; - tag_t tag; - - matrix4() : eigen_base(), tag(IDENTITY) {} - matrix4(const Eigen::Matrix4d& c) : eigen_base(c), tag(OTHER) {} - matrix4(const Eigen::Vector3d& o, const Eigen::Vector3d& z, const Eigen::Vector3d& x) : tag(AFFINE_WO_SCALE) { - init(o, z, x); - } - matrix4(const Eigen::Vector3d& o, const Eigen::Vector3d& z) : tag(AFFINE_WO_SCALE) { - auto x = Eigen::Vector3d(1, 0, 0); - auto y = z.cross(x); - if (y.squaredNorm() < 1.e-7) { - x = Eigen::Vector3d(0, 0, 1); - } - init(o, z, x); - } - - bool is_identity() const { - return !components_ || components_->isIdentity(); - } - - void print(std::ostream& o, int indent = 0) const { - print_impl(o, "matrix4", indent); - } - - virtual matrix4* clone_() const { return new matrix4(*this); } - virtual kinds kind() const { return MATRIX4; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(MATRIX4), hash_components()); - return boost::hash{}(v); - } - - Eigen::Vector3d translation_part() const { return ccomponents().col(3).head<3>(); } -}; - -struct colour : public item, public eigen_base { - DECLARE_PTR(colour) - - void print(std::ostream& o, int indent = 0) const { - print_impl(o, "colour", indent); - } - - virtual colour* clone_() const { return new colour(*this); } - virtual kinds kind() const { return COLOUR; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(COLOUR), hash_components()); - return boost::hash{}(v); - } - - colour() : eigen_base() {} - colour(double r, double g, double b) { components() << r, g, b; } - - const double& r() const { return ccomponents()[0]; } - const double& g() const { return ccomponents()[1]; } - const double& b() const { return ccomponents()[2]; } -}; - -struct style : public item { - DECLARE_PTR(style) - - std::string name; - colour diffuse; - colour specular; - double specularity, transparency; - - void print(std::ostream& o, int indent = 0) const { - o << std::string(indent, ' ') << "style" << std::endl; - o << std::string(indent, ' ') << " " << "name" << (name) << std::endl; - if (diffuse.components_) { - o << std::string(indent, ' ') << " " << "diffuse" << (name) << std::endl; - diffuse.print(o, indent + 5 + 7); - } - if (specular.components_) { - o << std::string(indent, ' ') << " " << "specular" << (name) << std::endl; - specular.print(o, indent + 5 + 8); - } - // @todo - } - - virtual style* clone_() const { return new style(*this); } - virtual kinds kind() const { return STYLE; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(STYLE), name, diffuse.hash(), specular.hash(), specularity, transparency); - return boost::hash{}(v); - } - - // @todo equality implementation based on values? - bool operator==(const style& other) const { return instance == other.instance; } - - style() : specularity(std::numeric_limits::quiet_NaN()), transparency(std::numeric_limits::quiet_NaN()) {} - style(const std::string& name) : name(name), specularity(std::numeric_limits::quiet_NaN()), transparency(std::numeric_limits::quiet_NaN()) {} - - bool has_specularity() const { - return !std::isnan(specularity); - } - - bool has_transparency() const { - return !std::isnan(transparency); - } -}; - -struct geom_item : public item { - DECLARE_PTR(geom_item) - - style::ptr surface_style; - matrix4::ptr matrix; - - geom_item(const IfcUtil::IfcBaseClass* instance = nullptr) : item(instance), surface_style(nullptr) {} - geom_item(const IfcUtil::IfcBaseClass* instance, matrix4::ptr m) : item(instance), surface_style(nullptr), matrix(m) {} - geom_item(matrix4::ptr m) : surface_style(nullptr), matrix(m) {} -}; - -// @todo make 4d for easier multiplication -template -struct cartesian_base : public item, public eigen_base { - cartesian_base() : eigen_base() {} - cartesian_base(double x, double y, double z = 0.) : eigen_base(Eigen::Vector3d(x, y, z)) {} -}; - -struct point3 : public cartesian_base<3> { - DECLARE_PTR(point3) - - virtual point3* clone_() const { return new point3(*this); } - virtual kinds kind() const { return POINT3; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(POINT3), hash_components()); - return boost::hash{}(v); - } - - void print(std::ostream& o, int indent = 0) const { - print_impl(o, "point3", indent); - } - - point3() : cartesian_base() {} - point3(double x, double y, double z = 0.) : cartesian_base(x, y, z) {} -}; - -struct direction3 : public cartesian_base<3> { - DECLARE_PTR(direction3) - - virtual direction3* clone_() const { return new direction3(*this); } - virtual kinds kind() const { return DIRECTION3; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(DIRECTION3), hash_components()); - return boost::hash{}(v); - } - - void print(std::ostream& o, int indent = 0) const { - print_impl(o, "direction3", indent); - } - - direction3() : cartesian_base() {} - direction3(double x, double y, double z = 0.) : cartesian_base(x, y, z) {} -}; - -struct curve : public geom_item { - void print_impl(std::ostream& o, const std::string& classname, int indent = 0) const { - o << std::string(indent, ' ') << classname << std::endl; - this->matrix->print(o, indent + 4); - } -}; - -struct line : public curve { - DECLARE_PTR(line) - - virtual line* clone_() const { return new line(*this); } - virtual kinds kind() const { return LINE; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(LINE), matrix->hash_components()); - return boost::hash{}(v); - } - - void print(std::ostream& o, int indent = 0) const { - print_impl(o, "line", indent); - } -}; - -struct circle : public curve { - DECLARE_PTR(circle) - - double radius; - - virtual circle* clone_() const { return new circle(*this); } - virtual kinds kind() const { return CIRCLE; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(CIRCLE), matrix->hash_components(), radius); - return boost::hash{}(v); - } - - void print(std::ostream& o, int indent = 0) const { - print_impl(o, "circle", indent); - } - - static circle::ptr from_3_points(const Eigen::Vector3d& p1, const Eigen::Vector3d& p2, const Eigen::Vector3d& p3) { - Eigen::Vector3d t = p2 - p1; - Eigen::Vector3d u = p3 - p1; - Eigen::Vector3d v = p3 - p2; - - auto norm = t.cross(u); - auto mag = norm.dot(norm); - - auto iwsl2 = 1. / (2. * mag); - auto tt = t.dot(t); - auto uu = u.dot(u); - - auto orig = p1 + (u * tt * u.dot(v) - t * uu * t.dot(v)) * iwsl2; - - if (!orig.array().isNaN().any()) { - auto radius = std::sqrt(tt * uu * v.dot(v) * iwsl2 * 0.5f); - auto ax = norm / std::sqrt(mag); - - - auto c = make(); - c->radius = radius; - c->matrix = taxonomy::make(orig, ax); - return c; - } - - return nullptr; - } -}; - -struct ellipse : public circle { - DECLARE_PTR(ellipse) - - double radius2; - - virtual ellipse* clone_() const { return new ellipse(*this); } - virtual kinds kind() const { return ELLIPSE; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(ELLIPSE), matrix->hash_components(), radius, radius2); - return boost::hash{}(v); - } - - void print(std::ostream& o, int indent = 0) const { - print_impl(o, "ellipse", indent); - } -}; - -struct bspline_curve : public curve { - DECLARE_PTR(bspline_curve) - - virtual bspline_curve* clone_() const { return new bspline_curve(*this); } - virtual kinds kind() const { return BSPLINE_CURVE; } - - virtual size_t calc_hash() const { - size_t h = std::hash{}(BSPLINE_CURVE); - for (auto& x : control_points) { - boost::hash_combine(h, x->hash()); - } - for (auto& x : multiplicities) { - boost::hash_combine(h, std::hash{}(x)); - } - for (auto& x : knots) { - boost::hash_combine(h, std::hash{}(x)); - } - if (weights) { - for (auto& x : *weights) { - boost::hash_combine(h, std::hash{}(x)); - } - } - boost::hash_combine(h, std::hash{}(degree)); - return h; - } - - std::vector control_points; - std::vector multiplicities; - std::vector knots; - boost::optional> weights; - int degree; - - void print(std::ostream& o, int indent = 0) const { - o << std::string(indent, ' ') << "bspline curve" << std::endl; - } -}; - -struct offset_curve : public curve { - DECLARE_PTR(offset_curve) - - direction3::ptr reference; - double offset; - item::ptr basis; - - virtual offset_curve* clone_() const { return new offset_curve(*this); } - virtual kinds kind() const { return OFFSET_CURVE; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(OFFSET_CURVE), reference->hash(), offset, basis ? basis->hash() : size_t(0)); - return boost::hash{}(v); - } - - void print(std::ostream& o, int indent = 0) const { - o << std::string(indent, ' ') << "offset_curve" << std::endl; - } -}; - -struct trimmed_curve : public item { - DECLARE_PTR(trimmed_curve) - - // @todo The copy constructor of point3 within the variant fails on the avx instruction - // on the default gcc in Ubuntu 18.04 and a recent AMD Ryzen. Probably due to allignment. - boost::variant start, end; - - // @todo somehow account for the fact that curve in IFC can be trimmed curve, polyline and composite curve as well. - item::ptr basis; - - // @todo does this make sense? this is to accomodate for the fact that orientation is defined on both TrimmedCurve as well CompCurveSegment - boost::optional orientation_2; - - trimmed_curve() : basis(nullptr), orientation_2(true) {} - trimmed_curve(const point3::ptr& a, const point3::ptr& b) : start(a), end(b), basis(nullptr) {} - - virtual void reverse() { - // std::swap(start, end); - orientation = !orientation; - } - - void print(std::ostream& o, int indent = 0) const { - o << std::string(indent, ' ') << "trimmed_curve" << std::endl; - if (basis) { - basis->print(o, indent + 4); - } - - const boost::variant * const start_end[2] = { &start, &end }; - for (int i = 0; i < 2; ++i) { - o << std::string(indent + 4, ' ') << (i == 0 ? "start" : "end") << std::endl; - if (start_end[i]->which() == 0) { - boost::get(*start_end[i])->print(o, indent + 4); - } else if (start_end[i]->which() == 1) { - o << std::string(indent + 4, ' ') << "parameter " << boost::get(*start_end[i]) << std::endl; - } - } - - if (this->instance) { - o << std::string(indent, ' ') << this->instance->data().toString() << std::endl; - } - } -}; - -struct edge : public trimmed_curve { - DECLARE_PTR(edge) - - edge() : trimmed_curve() {} - edge(const point3::ptr& a, const point3::ptr& b) : trimmed_curve(a, b) {} - - // @todo how to express similarity between trimmed_curve and edge? - virtual edge* clone_() const { return new edge(*this); } - virtual kinds kind() const { return EDGE; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(EDGE), start, end, basis ? basis->hash() : size_t(0), orientation_2 ? *orientation_2 ? 2 : 1 : 0); - return boost::hash{}(v); - } -}; - -template -struct collection_base : public geom_item { - std::vector children; - - collection_base() {} - collection_base(const collection_base& other) { - std::transform(other.children.begin(), other.children.end(), std::back_inserter(children), [](typename T::ptr p) { return clone(p); }); - } - - /* - template - std::vector children_as() const { - std::vector ts; - ts.reserve(children.size()); - std::for_each(children.begin(), children.end(), [&ts](ptr i){ - auto v = dcast(i); - if (v) { - ts.push_back(v); - } - }); - return ts; - } - */ - - virtual void reverse() { - // @todo this needs to create copies of the children in case of shared_ptr - std::reverse(children.begin(), children.end()); - for (auto& child : children) { - child->reverse(); - } - } - - void print(std::ostream& o, int indent = 0) const { - o << std::string(indent, ' ') << kind_to_string(kind()) << std::endl; - if (!matrix->is_identity()) { - matrix->print(o, indent + 4); - } - for (auto& c : children) { - c->print(o, indent + 4); - } - } - - virtual ~collection_base() { + + template + struct eigen_base { + T* components_; + + eigen_base() { + components_ = nullptr; + } + + eigen_base(const eigen_base& other) { + this->components_ = other.components_ ? new T(*other.components_) : nullptr; + } + + eigen_base(const T& other) { + this->components_ = new T(other); + } + + eigen_base& operator=(const eigen_base& other) { + if (this != &other) { + this->components_ = other.components_ ? new T(*other.components_) : nullptr; + } + return *this; + } + + void print_impl(std::ostream& o, const std::string& class_name, int indent = 0) const { + o << std::string(indent, ' ') << class_name; + if (this->components_) { + int n = T::RowsAtCompileTime * T::ColsAtCompileTime; + for (size_t i = 0; i < n; ++i) { + o << " " << (*components_)(i); + } + } + o << std::endl; + } + + virtual ~eigen_base() { + delete this->components_; + } + + const T& ccomponents() const { + if (this->components_) { + return *this->components_; + } + else { + return eigen_defaults(); + } + } + + T& components() { + if (!this->components_) { + this->components_ = new T(eigen_defaults()); + } + return *this->components_; + } + + explicit operator bool() const { + return components_; + } + + uint32_t hash_components() const { + size_t h = std::hash{}(T::RowsAtCompileTime); + boost::hash_combine(h, std::hash{}(T::ColsAtCompileTime)); + if (components_) { + for (size_t i = 0; i < components_->size(); ++i) { + auto elem = *(components_->data() + i); + boost::hash_combine(h, std::hash()(elem)); + } + } + return h; + } + }; + + struct matrix4 : public item, public eigen_base { + private: + void init(const Eigen::Vector3d& o, const Eigen::Vector3d& z, const Eigen::Vector3d& x) { + auto X = x.normalized(); + auto Y = z.cross(x).normalized(); + auto Z = z.normalized(); + components_ = new Eigen::Matrix4d; + (*components_) << + X(0), Y(0), Z(0), o(0), + X(1), Y(1), Z(1), o(1), + X(2), Y(2), Z(2), o(2), + 0, 0, 0, 1.; + if (is_identity()) { + // @todo detect this earlier to save us the heapalloc. + delete components_; + components_ = nullptr; + tag = IDENTITY; + } + } + public: + DECLARE_PTR(matrix4) + + enum tag_t { + IDENTITY, AFFINE_WO_SCALE, AFFINE_W_UNIFORM_SCALE, AFFINE_W_NONUNIFORM_SCALE, OTHER + }; + tag_t tag; + + matrix4() : eigen_base(), tag(IDENTITY) {} + matrix4(const Eigen::Matrix4d& c) : eigen_base(c), tag(OTHER) {} + matrix4(const Eigen::Vector3d& o, const Eigen::Vector3d& z, const Eigen::Vector3d& x) : tag(AFFINE_WO_SCALE) { + init(o, z, x); + } + matrix4(const Eigen::Vector3d& o, const Eigen::Vector3d& z) : tag(AFFINE_WO_SCALE) { + auto x = Eigen::Vector3d(1, 0, 0); + auto y = z.cross(x); + if (y.squaredNorm() < 1.e-7) { + x = Eigen::Vector3d(0, 0, 1); + } + init(o, z, x); + } + + bool is_identity() const { + return !components_ || components_->isIdentity(); + } + + void print(std::ostream& o, int indent = 0) const { + print_impl(o, "matrix4", indent); + } + + virtual matrix4* clone_() const { return new matrix4(*this); } + virtual kinds kind() const { return MATRIX4; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(MATRIX4), hash_components()); + return boost::hash{}(v); + } + + Eigen::Vector3d translation_part() const { return ccomponents().col(3).head<3>(); } + }; + + struct colour : public item, public eigen_base { + DECLARE_PTR(colour) + + void print(std::ostream& o, int indent = 0) const { + print_impl(o, "colour", indent); + } + + virtual colour* clone_() const { return new colour(*this); } + virtual kinds kind() const { return COLOUR; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(COLOUR), hash_components()); + return boost::hash{}(v); + } + + colour() : eigen_base() {} + colour(double r, double g, double b) { components() << r, g, b; } + + const double& r() const { return ccomponents()[0]; } + const double& g() const { return ccomponents()[1]; } + const double& b() const { return ccomponents()[2]; } + }; + + struct style : public item { + DECLARE_PTR(style) + + std::string name; + colour diffuse; + colour specular; + double specularity, transparency; + + void print(std::ostream& o, int indent = 0) const { + o << std::string(indent, ' ') << "style" << std::endl; + o << std::string(indent, ' ') << " " << "name" << (name) << std::endl; + if (diffuse.components_) { + o << std::string(indent, ' ') << " " << "diffuse" << (name) << std::endl; + diffuse.print(o, indent + 5 + 7); + } + if (specular.components_) { + o << std::string(indent, ' ') << " " << "specular" << (name) << std::endl; + specular.print(o, indent + 5 + 8); + } + // @todo + } + + virtual style* clone_() const { return new style(*this); } + virtual kinds kind() const { return STYLE; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(STYLE), name, diffuse.hash(), specular.hash(), specularity, transparency); + return boost::hash{}(v); + } + + // @todo equality implementation based on values? + bool operator==(const style& other) const { return instance == other.instance; } + + style() : specularity(std::numeric_limits::quiet_NaN()), transparency(std::numeric_limits::quiet_NaN()) {} + style(const std::string& name) : name(name), specularity(std::numeric_limits::quiet_NaN()), transparency(std::numeric_limits::quiet_NaN()) {} + + bool has_specularity() const { + return !std::isnan(specularity); + } + + bool has_transparency() const { + return !std::isnan(transparency); + } + }; + + struct geom_item : public item { + DECLARE_PTR(geom_item) + + style::ptr surface_style; + matrix4::ptr matrix; + + geom_item(const IfcUtil::IfcBaseClass* instance = nullptr) : item(instance), surface_style(nullptr) {} + geom_item(const IfcUtil::IfcBaseClass* instance, matrix4::ptr m) : item(instance), surface_style(nullptr), matrix(m) {} + geom_item(matrix4::ptr m) : surface_style(nullptr), matrix(m) {} + }; + + // @todo make 4d for easier multiplication + template + struct cartesian_base : public item, public eigen_base { + cartesian_base() : eigen_base() {} + cartesian_base(double x, double y, double z = 0.) : eigen_base(Eigen::Vector3d(x, y, z)) {} + }; + + struct point3 : public cartesian_base<3> { + DECLARE_PTR(point3) + + virtual point3* clone_() const { return new point3(*this); } + virtual kinds kind() const { return POINT3; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(POINT3), hash_components()); + return boost::hash{}(v); + } + + void print(std::ostream& o, int indent = 0) const { + print_impl(o, "point3", indent); + } + + point3() : cartesian_base() {} + point3(double x, double y, double z = 0.) : cartesian_base(x, y, z) {} + }; + + struct direction3 : public cartesian_base<3> { + DECLARE_PTR(direction3) + + virtual direction3* clone_() const { return new direction3(*this); } + virtual kinds kind() const { return DIRECTION3; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(DIRECTION3), hash_components()); + return boost::hash{}(v); + } + + void print(std::ostream& o, int indent = 0) const { + print_impl(o, "direction3", indent); + } + + direction3() : cartesian_base() {} + direction3(double x, double y, double z = 0.) : cartesian_base(x, y, z) {} + }; + + struct curve : public geom_item { + void print_impl(std::ostream& o, const std::string& classname, int indent = 0) const { + o << std::string(indent, ' ') << classname << std::endl; + this->matrix->print(o, indent + 4); + } + }; + + struct line : public curve { + DECLARE_PTR(line) + + virtual line* clone_() const { return new line(*this); } + virtual kinds kind() const { return LINE; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(LINE), matrix->hash_components()); + return boost::hash{}(v); + } + + void print(std::ostream& o, int indent = 0) const { + print_impl(o, "line", indent); + } + }; + + struct circle : public curve { + DECLARE_PTR(circle) + + double radius; + + virtual circle* clone_() const { return new circle(*this); } + virtual kinds kind() const { return CIRCLE; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(CIRCLE), matrix->hash_components(), radius); + return boost::hash{}(v); + } + + void print(std::ostream& o, int indent = 0) const { + print_impl(o, "circle", indent); + } + + static circle::ptr from_3_points(const Eigen::Vector3d& p1, const Eigen::Vector3d& p2, const Eigen::Vector3d& p3) { + Eigen::Vector3d t = p2 - p1; + Eigen::Vector3d u = p3 - p1; + Eigen::Vector3d v = p3 - p2; + + auto norm = t.cross(u); + auto mag = norm.dot(norm); + + auto iwsl2 = 1. / (2. * mag); + auto tt = t.dot(t); + auto uu = u.dot(u); + + auto orig = p1 + (u * tt * u.dot(v) - t * uu * t.dot(v)) * iwsl2; + + if (!orig.array().isNaN().any()) { + auto radius = std::sqrt(tt * uu * v.dot(v) * iwsl2 * 0.5f); + auto ax = norm / std::sqrt(mag); + + + auto c = make(); + c->radius = radius; + c->matrix = taxonomy::make(orig, ax); + return c; + } + + return nullptr; + } + }; + + struct ellipse : public circle { + DECLARE_PTR(ellipse) + + double radius2; + + virtual ellipse* clone_() const { return new ellipse(*this); } + virtual kinds kind() const { return ELLIPSE; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(ELLIPSE), matrix->hash_components(), radius, radius2); + return boost::hash{}(v); + } + + void print(std::ostream& o, int indent = 0) const { + print_impl(o, "ellipse", indent); + } + }; + + struct bspline_curve : public curve { + DECLARE_PTR(bspline_curve) + + virtual bspline_curve* clone_() const { return new bspline_curve(*this); } + virtual kinds kind() const { return BSPLINE_CURVE; } + + virtual size_t calc_hash() const { + size_t h = std::hash{}(BSPLINE_CURVE); + for (auto& x : control_points) { + boost::hash_combine(h, x->hash()); + } + for (auto& x : multiplicities) { + boost::hash_combine(h, std::hash{}(x)); + } + for (auto& x : knots) { + boost::hash_combine(h, std::hash{}(x)); + } + if (weights) { + for (auto& x : *weights) { + boost::hash_combine(h, std::hash{}(x)); + } + } + boost::hash_combine(h, std::hash{}(degree)); + return h; + } + + std::vector control_points; + std::vector multiplicities; + std::vector knots; + boost::optional> weights; + int degree; + + void print(std::ostream& o, int indent = 0) const { + o << std::string(indent, ' ') << "bspline curve" << std::endl; + } + }; + + struct offset_curve : public curve { + DECLARE_PTR(offset_curve) + + direction3::ptr reference; + double offset; + item::ptr basis; + + virtual offset_curve* clone_() const { return new offset_curve(*this); } + virtual kinds kind() const { return OFFSET_CURVE; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(OFFSET_CURVE), reference->hash(), offset, basis ? basis->hash() : size_t(0)); + return boost::hash{}(v); + } + + void print(std::ostream& o, int indent = 0) const { + o << std::string(indent, ' ') << "offset_curve" << std::endl; + } + }; + + struct trimmed_curve : public item { + DECLARE_PTR(trimmed_curve) + + // @todo The copy constructor of point3 within the variant fails on the avx instruction + // on the default gcc in Ubuntu 18.04 and a recent AMD Ryzen. Probably due to allignment. + boost::variant start, end; + + // @todo somehow account for the fact that curve in IFC can be trimmed curve, polyline and composite curve as well. + item::ptr basis; + + // @todo does this make sense? this is to accomodate for the fact that orientation is defined on both TrimmedCurve as well CompCurveSegment + boost::optional orientation_2; + + trimmed_curve() : basis(nullptr), orientation_2(true) {} + trimmed_curve(const point3::ptr& a, const point3::ptr& b) : start(a), end(b), basis(nullptr) {} + + virtual void reverse() { + // std::swap(start, end); + orientation = !orientation; + } + + void print(std::ostream& o, int indent = 0) const { + o << std::string(indent, ' ') << "trimmed_curve" << std::endl; + if (basis) { + basis->print(o, indent + 4); + } + + const boost::variant* const start_end[2] = { &start, &end }; + for (int i = 0; i < 2; ++i) { + o << std::string(indent + 4, ' ') << (i == 0 ? "start" : "end") << std::endl; + if (start_end[i]->which() == 0) { + boost::get(*start_end[i])->print(o, indent + 4); + } + else if (start_end[i]->which() == 1) { + o << std::string(indent + 4, ' ') << "parameter " << boost::get(*start_end[i]) << std::endl; + } + } + + if (this->instance) { + o << std::string(indent, ' ') << this->instance->data().toString() << std::endl; + } + } + }; + + struct edge : public trimmed_curve { + DECLARE_PTR(edge) + + edge() : trimmed_curve() {} + edge(const point3::ptr& a, const point3::ptr& b) : trimmed_curve(a, b) {} + + // @todo how to express similarity between trimmed_curve and edge? + virtual edge* clone_() const { return new edge(*this); } + virtual kinds kind() const { return EDGE; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(EDGE), start, end, basis ? basis->hash() : size_t(0), orientation_2 ? *orientation_2 ? 2 : 1 : 0); + return boost::hash{}(v); + } + }; + + template + struct collection_base : public geom_item { + std::vector children; + + collection_base() {} + collection_base(const collection_base& other) { + std::transform(other.children.begin(), other.children.end(), std::back_inserter(children), [](typename T::ptr p) { return clone(p); }); + } + + /* + template + std::vector children_as() const { + std::vector ts; + ts.reserve(children.size()); + std::for_each(children.begin(), children.end(), [&ts](ptr i){ + auto v = dcast(i); + if (v) { + ts.push_back(v); + } + }); + return ts; + } + */ + + virtual void reverse() { + // @todo this needs to create copies of the children in case of shared_ptr + std::reverse(children.begin(), children.end()); + for (auto& child : children) { + child->reverse(); + } + } + + void print(std::ostream& o, int indent = 0) const { + o << std::string(indent, ' ') << kind_to_string(kind()) << std::endl; + if (!matrix->is_identity()) { + matrix->print(o, indent + 4); + } + for (auto& c : children) { + c->print(o, indent + 4); + } + } + + virtual ~collection_base() { #ifdef TAXONOMY_USE_NAKED_PTR - for (auto& c : children) { - delete c; - } + for (auto& c : children) { + delete c; + } #endif - } - - uint32_t hash_elements() const { - size_t h = 0; - for (auto& c : children) { - boost::hash_combine(h, c->hash()); - } - return h; - } -}; - -struct collection : public collection_base { - DECLARE_PTR(collection) - - virtual collection* clone_() const { return new collection(*this); } - virtual kinds kind() const { return COLLECTION; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(COLLECTION), hash_elements()); - return boost::hash{}(v); - } -}; - - -struct loop : public collection_base { - DECLARE_PTR(loop) - - boost::optional external, closed; - - bool is_polyhedron() const { - for (auto& e : children) { - if (e->basis != nullptr) { - if (e->basis->kind() != LINE) { - return false; } - } - } - return true; - } - virtual loop* clone_() const { return new loop(*this); } - virtual kinds kind() const { return LOOP; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(LOOP), hash_elements(), external ? *external ? 2 : 1 : 0, closed ? *closed ? 2 : 1 : 0); - return boost::hash{}(v); - } -}; - -struct face : public collection_base { - DECLARE_PTR(face) - - item::ptr basis; - - virtual face* clone_() const { return new face(*this); } - virtual kinds kind() const { return FACE; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(FACE), hash_elements(), basis ? basis->hash() : size_t(0)); - return boost::hash{}(v); - } -}; - -struct shell : public collection_base { - DECLARE_PTR(shell) - - boost::optional closed; - - virtual shell* clone_() const { return new shell(*this); } - virtual kinds kind() const { return SHELL; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(SHELL), hash_elements(), closed ? *closed ? 2 : 1 : 0); - return boost::hash{}(v); - } -}; - -struct solid : public collection_base { - DECLARE_PTR(solid) - - virtual solid* clone_() const { return new solid(*this); } - virtual kinds kind() const { return SOLID; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(SOLID), hash_elements()); - return boost::hash{}(v); - } -}; - -struct loft : public collection_base { - DECLARE_PTR(loft) - - item::ptr axis; - - virtual loft* clone_() const { return new loft(*this); } - virtual kinds kind() const { return LOFT; } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(LOFT), hash_elements(), axis ? axis->hash() : size_t(0)); - return boost::hash{}(v); - } -}; - -struct surface : public geom_item {}; - -struct plane : public surface { - DECLARE_PTR(plane) - - virtual plane* clone_() const { return new plane(*this); } - virtual kinds kind() const { return PLANE; } - - void print(std::ostream& o, int) const { - o << "not implemented"; - } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(PLANE), matrix->hash_components()); - return boost::hash{}(v); - } -}; - -struct cylinder : public surface { - DECLARE_PTR(cylinder) - - double radius; - - virtual cylinder* clone_() const { return new cylinder(*this); } - virtual kinds kind() const { return CYLINDER; } - - void print(std::ostream& o, int) const { - o << "not implemented"; - } - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(CYLINDER), matrix->hash_components()); - return boost::hash{}(v); - } -}; - -struct bspline_surface : public surface { - DECLARE_PTR(bspline_surface) - - virtual bspline_surface* clone_() const { return new bspline_surface(*this); } - virtual kinds kind() const { return BSPLINE_SURFACE; } - - virtual size_t calc_hash() const { - size_t h = std::hash{}(BSPLINE_SURFACE); - boost::hash_combine(h, std::hash{}(control_points.size())); - for (auto& xs : control_points) { - for (auto& x : xs) { - boost::hash_combine(h, x->hash()); - } - } - for (auto& xs : multiplicities) { - for (auto& x : xs) { - boost::hash_combine(h, std::hash{}(x)); - } - } - for (auto& xs : knots) { - for (auto& x : xs) { - boost::hash_combine(h, std::hash{}(x)); - } - } - if (weights) { - for (auto& xs : *weights) { - for (auto& x : xs) { - boost::hash_combine(h, std::hash{}(x)); + uint32_t hash_elements() const { + size_t h = 0; + for (auto& c : children) { + boost::hash_combine(h, c->hash()); + } + return h; } - } - } - boost::hash_combine(h, std::hash{}(degree[0])); - boost::hash_combine(h, std::hash{}(degree[1])); - return h; - } + }; - std::vector> control_points; - std::array, 2> multiplicities; - std::array, 2> knots; - boost::optional>> weights; - std::array degree; + struct collection : public collection_base { + DECLARE_PTR(collection) - void print(std::ostream& o, int) const { - o << "not implemented"; - } -}; + virtual collection* clone_() const { return new collection(*this); } + virtual kinds kind() const { return COLLECTION; } -struct sweep : public geom_item { - DECLARE_PTR(sweep) - - face::ptr basis; - - sweep(face::ptr b) : basis(b) {} - sweep(matrix4::ptr m, face::ptr b) : geom_item(m), basis(b) {} -}; - -struct extrusion : public sweep { - DECLARE_PTR(extrusion) - - direction3::ptr direction; - double depth; - - virtual extrusion* clone_() const { return new extrusion(*this); } - virtual kinds kind() const { return EXTRUSION; } - - extrusion(matrix4::ptr m, face::ptr basis, direction3::ptr dir, double d) : sweep(m, basis), direction(dir), depth(d) {} - - void print(std::ostream& o, int indent = 0) const { - o << std::string(indent, ' ') << "extrusion " << depth << std::endl; - direction->print(o, indent + 4); - basis->print(o, indent + 4); - } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(EXTRUSION), matrix->hash_components(), basis->calc_hash(), direction->hash_components(), depth); - return boost::hash{}(v); - } -}; - -struct revolve : public sweep { - DECLARE_PTR(revolve) - - point3::ptr axis_origin; - direction3::ptr direction; - boost::optional angle; - - virtual revolve* clone_() const { return new revolve(*this); } - virtual kinds kind() const { return REVOLVE; } - - revolve(matrix4::ptr m, face::ptr basis, point3::ptr pnt, direction3::ptr dir, const boost::optional& a) : sweep(m, basis), axis_origin(pnt), direction(dir), angle(a) {} - - void print(std::ostream& o, int indent = 0) const { - o << std::string(indent, ' ') << "revolve" << std::endl; - } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(REVOLVE), matrix->hash_components(), basis->calc_hash(), axis_origin->hash_components(), direction->hash_components(), angle ? *angle : 1000.); - return boost::hash{}(v); - } -}; - -struct surface_curve_sweep : public sweep { - DECLARE_PTR(surface_curve_sweep) - - item::ptr surface; - item::ptr curve; - - virtual surface_curve_sweep* clone_() const { return new surface_curve_sweep(*this); } - virtual kinds kind() const { return SURFACE_CURVE_SWEEP; } - - surface_curve_sweep(matrix4::ptr m, face::ptr basis, item::ptr surf, item::ptr crv) : sweep(m, basis), surface(surf), curve(crv) {} - - void print(std::ostream& o, int indent = 0) const { - o << std::string(indent, ' ') << "surface_curve_sweep" << std::endl; - } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(SURFACE_CURVE_SWEEP), matrix->hash_components(), basis->calc_hash(), surface->calc_hash(), curve->calc_hash()); - return boost::hash{}(v); - } -}; - -struct node : public item { - DECLARE_PTR(node) - - // std::map representations; - - virtual node* clone_() const { return new node(*this); } - virtual kinds kind() const { return NODE; } - - void print(std::ostream&, int = 0) const {} - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(NODE)); - return boost::hash{}(v); - } -}; - -struct boolean_result : public collection_base { - DECLARE_PTR(boolean_result) - - enum operation_t { - UNION, SUBTRACTION, INTERSECTION - }; - - virtual boolean_result* clone_() const { return new boolean_result(*this); } - virtual kinds kind() const { return BOOLEAN_RESULT; } - operation_t operation; - - static const std::string& operation_str(operation_t op) { - using namespace std::string_literals; - static std::string s[] = { "union"s, "subtraction"s, "intersection"s }; - return s[(size_t)op]; - } - - virtual size_t calc_hash() const { - auto v = std::make_tuple(static_cast(BOOLEAN_RESULT), hash_elements(), static_cast(operation)); - return boost::hash{}(v); - } -}; - -namespace impl { - typedef std::tuple KindsTuple; - typedef std::tuple CurvesTuple; - typedef std::tuple SurfacesTuple; -} - -struct type_by_kind { - template - using type = typename std::tuple_element::type; - - static const size_t max = std::tuple_size::value; -}; - -struct curves { - template - using type = typename std::tuple_element::type; - - static const size_t max = std::tuple_size::value; -}; - -struct surfaces { - template - using type = typename std::tuple_element::type; - - static const size_t max = std::tuple_size::value; -}; - -// Hacks around not wanting to use if constexpr -template -class loop_to_face_upgrade { -public: - loop_to_face_upgrade(taxonomy::ptr) {} - - operator bool() const { - return false; - } - - operator taxonomy::face::ptr() const { - throw taxonomy::topology_error(); - } - - operator typename T::ptr() const { - throw taxonomy::topology_error(); - } -}; - -template <> -class loop_to_face_upgrade { -private: - boost::optional face_; -public: - loop_to_face_upgrade(taxonomy::ptr item) { - auto loop = taxonomy::dcast(item); - if (loop) { - loop->external = true; - - face_ = taxonomy::make(); - (*face_)->instance = loop->instance; - (*face_)->matrix = loop->matrix; - (*face_)->children = { taxonomy::clone(loop) }; - } - } - - operator bool() const { - return face_.is_initialized(); - } - - operator taxonomy::face::ptr() const { - return *face_; - } -}; - - -#ifdef TAXONOMY_USE_SHARED_PTR -template -std::shared_ptr cast(const std::shared_ptr& u) { - loop_to_face_upgrade upg(u); - if (upg) { - return upg; - } - return std::static_pointer_cast(u); -} -template -std::shared_ptr dcast(const std::shared_ptr& u) { - loop_to_face_upgrade upg(u); - if (upg) { - return upg; - } - return std::dynamic_pointer_cast(u); -} -#endif -#ifdef TAXONOMY_USE_UNIQUE_PTR -template -T* cast(const std::unique_ptr& u) { - loop_to_face_upgrade upg(u); - if (upg) { - return upg; - } - return static_cast(&*u); -} -template -T* dcast(const std::unique_ptr& u) { - loop_to_face_upgrade upg(u); - if (upg) { - return upg; - } - return dynamic_cast(&*u); -} -#endif -#ifdef TAXONOMY_USE_NAKED_PTR -template -T* cast(const U*& u) { - loop_to_face_upgrade upg(u); - if (upg) { - return upg; - } - return std::static_cast(u); -} -template -T* dcast(const U*& u) { - loop_to_face_upgrade upg(u); - if (upg) { - return upg; - } - return std::dynamic_cast(u); -} -#endif - -} - - template - void visit(typename U::ptr deep, Fn fn) { - for (auto& i : deep->children) { - // @todo Sad... now that we have templated collection members, - // we can't generally use collection_base anymore as a cast target. - if (auto s = taxonomy::dcast(i)) { - visit(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit(s, fn); - } else { - fn(i); - } - } - } - - template - void visit_2(typename U::ptr c, const Fn& fn) { - static_assert(std::is_same::value, "@todo Only implemented for point3"); - for (auto& i : c->children) { - // @todo Sad... now that we have templated collection members, - // we can't generally use collection_base anymore as a cast target. - if (auto s = taxonomy::dcast(i)) { - visit_2(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit_2(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit_2(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit_2(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit_2(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit_2(s, fn); - } else if (auto s = taxonomy::dcast(i)) { - visit_2(s, fn); - } else if (auto pt = taxonomy::dcast(i)) { - fn(pt); - } else if (auto l = taxonomy::dcast(i)) { - // @todo maybe make edge a collection then as well? - if (l->start.which() == 0) { - fn(boost::get(l->start)); + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(COLLECTION), hash_elements()); + return boost::hash{}(v); } - if (l->end.which() == 0) { - fn(boost::get(l->end)); - } - } - } - } + }; - taxonomy::collection::ptr flatten(taxonomy::collection::ptr deep); - template - bool apply_predicate_to_collection(taxonomy::ptr i, Fn fn) { - if (i->kind() == taxonomy::COLLECTION) { - auto c = taxonomy::cast(i); - for (auto& child : c->children) { - if (apply_predicate_to_collection(child, fn)) { + struct loop : public collection_base { + DECLARE_PTR(loop) + + boost::optional external, closed; + + bool is_polyhedron() const { + for (auto& e : children) { + if (e->basis != nullptr) { + if (e->basis->kind() != LINE) { + return false; + } + } + } return true; } - } - } else { - return fn(i); - } - } - // @nb traverses nested collections - template - taxonomy::collection::ptr filter(taxonomy::collection::ptr collection, Fn fn) { - auto filtered = taxonomy::make(); - for (auto& child : collection->children) { - if (apply_predicate_to_collection(child, fn)) { - filtered->children.push_back(clone(child)); + virtual loop* clone_() const { return new loop(*this); } + virtual kinds kind() const { return LOOP; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(LOOP), hash_elements(), external ? *external ? 2 : 1 : 0, closed ? *closed ? 2 : 1 : 0); + return boost::hash{}(v); + } + }; + + struct face : public collection_base { + DECLARE_PTR(face) + + item::ptr basis; + + virtual face* clone_() const { return new face(*this); } + virtual kinds kind() const { return FACE; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(FACE), hash_elements(), basis ? basis->hash() : size_t(0)); + return boost::hash{}(v); + } + }; + + struct shell : public collection_base { + DECLARE_PTR(shell) + + boost::optional closed; + + virtual shell* clone_() const { return new shell(*this); } + virtual kinds kind() const { return SHELL; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(SHELL), hash_elements(), closed ? *closed ? 2 : 1 : 0); + return boost::hash{}(v); + } + }; + + struct solid : public collection_base { + DECLARE_PTR(solid) + + virtual solid* clone_() const { return new solid(*this); } + virtual kinds kind() const { return SOLID; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(SOLID), hash_elements()); + return boost::hash{}(v); + } + }; + + struct loft : public collection_base { + DECLARE_PTR(loft) + + item::ptr axis; + + virtual loft* clone_() const { return new loft(*this); } + virtual kinds kind() const { return LOFT; } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(LOFT), hash_elements(), axis ? axis->hash() : size_t(0)); + return boost::hash{}(v); + } + }; + + struct surface : public geom_item {}; + + struct plane : public surface { + DECLARE_PTR(plane) + + virtual plane* clone_() const { return new plane(*this); } + virtual kinds kind() const { return PLANE; } + + void print(std::ostream& o, int) const { + o << "not implemented"; + } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(PLANE), matrix->hash_components()); + return boost::hash{}(v); + } + }; + + struct cylinder : public surface { + DECLARE_PTR(cylinder) + + double radius; + + virtual cylinder* clone_() const { return new cylinder(*this); } + virtual kinds kind() const { return CYLINDER; } + + void print(std::ostream& o, int) const { + o << "not implemented"; + } + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(CYLINDER), matrix->hash_components()); + return boost::hash{}(v); + } + }; + + struct bspline_surface : public surface { + DECLARE_PTR(bspline_surface) + + virtual bspline_surface* clone_() const { return new bspline_surface(*this); } + virtual kinds kind() const { return BSPLINE_SURFACE; } + + virtual size_t calc_hash() const { + size_t h = std::hash{}(BSPLINE_SURFACE); + boost::hash_combine(h, std::hash{}(control_points.size())); + for (auto& xs : control_points) { + for (auto& x : xs) { + boost::hash_combine(h, x->hash()); + } + } + for (auto& xs : multiplicities) { + for (auto& x : xs) { + boost::hash_combine(h, std::hash{}(x)); + } + } + for (auto& xs : knots) { + for (auto& x : xs) { + boost::hash_combine(h, std::hash{}(x)); + } + } + if (weights) { + for (auto& xs : *weights) { + for (auto& x : xs) { + boost::hash_combine(h, std::hash{}(x)); + } + } + } + boost::hash_combine(h, std::hash{}(degree[0])); + boost::hash_combine(h, std::hash{}(degree[1])); + return h; + } + + std::vector> control_points; + std::array, 2> multiplicities; + std::array, 2> knots; + boost::optional>> weights; + std::array degree; + + void print(std::ostream& o, int) const { + o << "not implemented"; + } + }; + + struct sweep : public geom_item { + DECLARE_PTR(sweep) + + face::ptr basis; + + sweep(face::ptr b) : basis(b) {} + sweep(matrix4::ptr m, face::ptr b) : geom_item(m), basis(b) {} + }; + + struct extrusion : public sweep { + DECLARE_PTR(extrusion) + + direction3::ptr direction; + double depth; + + virtual extrusion* clone_() const { return new extrusion(*this); } + virtual kinds kind() const { return EXTRUSION; } + + extrusion(matrix4::ptr m, face::ptr basis, direction3::ptr dir, double d) : sweep(m, basis), direction(dir), depth(d) {} + + void print(std::ostream& o, int indent = 0) const { + o << std::string(indent, ' ') << "extrusion " << depth << std::endl; + direction->print(o, indent + 4); + basis->print(o, indent + 4); + } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(EXTRUSION), matrix->hash_components(), basis->calc_hash(), direction->hash_components(), depth); + return boost::hash{}(v); + } + }; + + struct revolve : public sweep { + DECLARE_PTR(revolve) + + point3::ptr axis_origin; + direction3::ptr direction; + boost::optional angle; + + virtual revolve* clone_() const { return new revolve(*this); } + virtual kinds kind() const { return REVOLVE; } + + revolve(matrix4::ptr m, face::ptr basis, point3::ptr pnt, direction3::ptr dir, const boost::optional& a) : sweep(m, basis), axis_origin(pnt), direction(dir), angle(a) {} + + void print(std::ostream& o, int indent = 0) const { + o << std::string(indent, ' ') << "revolve" << std::endl; + } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(REVOLVE), matrix->hash_components(), basis->calc_hash(), axis_origin->hash_components(), direction->hash_components(), angle ? *angle : 1000.); + return boost::hash{}(v); + } + }; + + struct surface_curve_sweep : public sweep { + DECLARE_PTR(surface_curve_sweep) + + item::ptr surface; + item::ptr curve; + + virtual surface_curve_sweep* clone_() const { return new surface_curve_sweep(*this); } + virtual kinds kind() const { return SURFACE_CURVE_SWEEP; } + + surface_curve_sweep(matrix4::ptr m, face::ptr basis, item::ptr surf, item::ptr crv) : sweep(m, basis), surface(surf), curve(crv) {} + + void print(std::ostream& o, int indent = 0) const { + o << std::string(indent, ' ') << "surface_curve_sweep" << std::endl; + } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(SURFACE_CURVE_SWEEP), matrix->hash_components(), basis->calc_hash(), surface->calc_hash(), curve->calc_hash()); + return boost::hash{}(v); + } + }; + + struct node : public item { + DECLARE_PTR(node) + + // std::map representations; + + virtual node* clone_() const { return new node(*this); } + virtual kinds kind() const { return NODE; } + + void print(std::ostream&, int = 0) const {} + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(NODE)); + return boost::hash{}(v); + } + }; + + struct boolean_result : public collection_base { + DECLARE_PTR(boolean_result) + + enum operation_t { + UNION, SUBTRACTION, INTERSECTION + }; + + virtual boolean_result* clone_() const { return new boolean_result(*this); } + virtual kinds kind() const { return BOOLEAN_RESULT; } + operation_t operation; + + static const std::string& operation_str(operation_t op) { + using namespace std::string_literals; + static std::string s[] = { "union"s, "subtraction"s, "intersection"s }; + return s[(size_t)op]; + } + + virtual size_t calc_hash() const { + auto v = std::make_tuple(static_cast(BOOLEAN_RESULT), hash_elements(), static_cast(operation)); + return boost::hash{}(v); + } + }; + + namespace impl { + typedef std::tuple KindsTuple; + typedef std::tuple CurvesTuple; + typedef std::tuple SurfacesTuple; + } + + struct type_by_kind { + template + using type = typename std::tuple_element::type; + + static const size_t max = std::tuple_size::value; + }; + + struct curves { + template + using type = typename std::tuple_element::type; + + static const size_t max = std::tuple_size::value; + }; + + struct surfaces { + template + using type = typename std::tuple_element::type; + + static const size_t max = std::tuple_size::value; + }; + + // Hacks around not wanting to use if constexpr + template + class loop_to_face_upgrade { + public: + loop_to_face_upgrade(taxonomy::ptr) {} + + operator bool() const { + return false; + } + + operator taxonomy::face::ptr() const { + throw taxonomy::topology_error(); + } + + operator typename T::ptr() const { + throw taxonomy::topology_error(); + } + }; + + template <> + class loop_to_face_upgrade { + private: + boost::optional face_; + public: + loop_to_face_upgrade(taxonomy::ptr item) { + auto loop = taxonomy::dcast(item); + if (loop) { + loop->external = true; + + face_ = taxonomy::make(); + (*face_)->instance = loop->instance; + (*face_)->matrix = loop->matrix; + (*face_)->children = { taxonomy::clone(loop) }; + } + } + + operator bool() const { + return face_.is_initialized(); + } + + operator taxonomy::face::ptr() const { + return *face_; + } + }; + + +#ifdef TAXONOMY_USE_SHARED_PTR + template + std::shared_ptr cast(const std::shared_ptr& u) { + loop_to_face_upgrade upg(u); + if (upg) { + return upg; + } + return std::static_pointer_cast(u); + } + template + std::shared_ptr dcast(const std::shared_ptr& u) { + loop_to_face_upgrade upg(u); + if (upg) { + return upg; + } + return std::dynamic_pointer_cast(u); } - } - if (filtered->children.empty()) { -#ifdef TAXONOMY_USE_NAKED_PTR - delete filtered; #endif - return nullptr; - } - return filtered; - } - - // @nb traverses nested collections - template - taxonomy::collection::ptr filter_in_place(taxonomy::collection::ptr collection, Fn fn) { - for (auto it = --collection->children.end(); it >= collection->children.begin(); --it) { - if (!apply_predicate_to_collection(*it, fn)) { -#ifdef TAXONOMY_USE_NAKED_PTR - delete *it; +#ifdef TAXONOMY_USE_UNIQUE_PTR + template + T* cast(const std::unique_ptr& u) { + loop_to_face_upgrade upg(u); + if (upg) { + return upg; + } + return static_cast(&*u); + } + template + T* dcast(const std::unique_ptr& u) { + loop_to_face_upgrade upg(u); + if (upg) { + return upg; + } + return dynamic_cast(&*u); + } #endif - collection->children.erase(it); +#ifdef TAXONOMY_USE_NAKED_PTR + template + T* cast(const U*& u) { + loop_to_face_upgrade upg(u); + if (upg) { + return upg; + } + return std::static_cast(u); + } + template + T* dcast(const U*& u) { + loop_to_face_upgrade upg(u); + if (upg) { + return upg; + } + return std::dynamic_cast(u); + } +#endif + + } + + template + void visit(typename U::ptr deep, Fn fn) { + for (auto& i : deep->children) { + // @todo Sad... now that we have templated collection members, + // we can't generally use collection_base anymore as a cast target. + if (auto s = taxonomy::dcast(i)) { + visit(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit(s, fn); + } + else { + fn(i); + } } } - return collection; + + template + void visit_2(typename U::ptr c, const Fn& fn) { + static_assert(std::is_same::value, "@todo Only implemented for point3"); + for (auto& i : c->children) { + // @todo Sad... now that we have templated collection members, + // we can't generally use collection_base anymore as a cast target. + if (auto s = taxonomy::dcast(i)) { + visit_2(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit_2(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit_2(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit_2(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit_2(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit_2(s, fn); + } + else if (auto s = taxonomy::dcast(i)) { + visit_2(s, fn); + } + else if (auto pt = taxonomy::dcast(i)) { + fn(pt); + } + else if (auto l = taxonomy::dcast(i)) { + // @todo maybe make edge a collection then as well? + if (l->start.which() == 0) { + fn(boost::get(l->start)); + } + if (l->end.which() == 0) { + fn(boost::get(l->end)); + } + } + } + } + + taxonomy::collection::ptr flatten(taxonomy::collection::ptr deep); + + template + bool apply_predicate_to_collection(taxonomy::ptr i, Fn fn) { + if (i->kind() == taxonomy::COLLECTION) { + auto c = taxonomy::cast(i); + for (auto& child : c->children) { + if (apply_predicate_to_collection(child, fn)) { + return true; + } + } + } + else { + return fn(i); + } + } + + // @nb traverses nested collections + template + taxonomy::collection::ptr filter(taxonomy::collection::ptr collection, Fn fn) { + auto filtered = taxonomy::make(); + for (auto& child : collection->children) { + if (apply_predicate_to_collection(child, fn)) { + filtered->children.push_back(clone(child)); + } + } + if (filtered->children.empty()) { +#ifdef TAXONOMY_USE_NAKED_PTR + delete filtered; +#endif + return nullptr; + } + return filtered; + } + + // @nb traverses nested collections + template + taxonomy::collection::ptr filter_in_place(taxonomy::collection::ptr collection, Fn fn) { + for (auto it = --collection->children.end(); it >= collection->children.begin(); --it) { + if (!apply_predicate_to_collection(*it, fn)) { +#ifdef TAXONOMY_USE_NAKED_PTR + delete* it; +#endif + collection->children.erase(it); + } + } + return collection; + } + + taxonomy::solid::ptr create_box(double dx, double dy, double dz); + taxonomy::solid::ptr create_box(double x, double y, double z, double dx, double dy, double dz); + + struct layerset_information { + std::vector thicknesses; + std::vector layers; + std::vector styles; + }; + + enum connection_type { + ATPATH, + ATSTART, + ATEND, + NOTDEFINED + }; + + typedef std::tuple endpoint_connection; } - taxonomy::solid::ptr create_box(double dx, double dy, double dz); - taxonomy::solid::ptr create_box(double x, double y, double z, double dx, double dy, double dz); - - struct layerset_information { - std::vector thicknesses; - std::vector layers; - std::vector styles; - }; - - enum connection_type { - ATPATH, - ATSTART, - ATEND, - NOTDEFINED - }; - - typedef std::tuple endpoint_connection; -} - } #endif \ No newline at end of file From a5c07d5e246fc4db887108e5a585eda60b4984e6 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Wed, 4 Oct 2023 13:25:31 -0700 Subject: [PATCH 14/15] Adds IfcPolynoialCurve for vertical curves --- src/ifcgeom/mapping/IfcCurveSegment.cpp | 49 +++++++++++++++++++++++-- 1 file changed, 45 insertions(+), 4 deletions(-) diff --git a/src/ifcgeom/mapping/IfcCurveSegment.cpp b/src/ifcgeom/mapping/IfcCurveSegment.cpp index a0c78d5a78..b2198fdb55 100644 --- a/src/ifcgeom/mapping/IfcCurveSegment.cpp +++ b/src/ifcgeom/mapping/IfcCurveSegment.cpp @@ -42,6 +42,7 @@ typedef boost::mpl::vector< #endif , IfcSchema::IfcPolyline , IfcSchema::IfcCircle + , IfcSchema::IfcPolynomialCurve > curve_seg_types; enum segment_type_t { @@ -338,6 +339,38 @@ public: } } + void operator()(IfcSchema::IfcPolynomialCurve* p) { + + if (segment_type_ == ST_HORIZONTAL) { + auto coeffX = p->CoefficientsX(); + auto coeffY = p->CoefficientsY(); + eval_ = [coeffX,coeffY](double u) { + + Eigen::VectorXd vec(4); + vec << 0.0, 0.0, 0.0, 1.0; + return vec; + }; + + } + else if (segment_type_ == ST_VERTICAL) { + auto coeffY = p->CoefficientsY(); + + eval_ = [coeffY](double u) { + const auto& coeffs = coeffY.get(); + auto exp = coeffs.size() - 1; + auto z = 0.0; + for (auto c : coeffs) + { + z += c * pow(u, exp--); + } + Eigen::VectorXd vec(4); + vec << 0.0, 0.0, z, 1.0; + return vec; + }; + + } + } + // Take the boost::type value from mpl::for_each and test it against our curve instance template void operator()(boost::type) { @@ -399,16 +432,24 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcCurveSegment* inst) { curve_segment_evaluator cse(this,length_unit_, segment_type, inst->ParentCurve(), inst->SegmentStart(), inst->SegmentLength()); boost::mpl::for_each>(std::ref(cse)); - - auto fn = *cse.evaluation_function(); + + 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()); // @todo - for some reason this isn't working, the matrix gets all messed up //const auto& transformation_matrix = taxonomy::cast(map(inst->Placement()))->ccomponents(); auto transformation_matrix = taxonomy::cast(map(inst->Placement()))->ccomponents(); - auto fn_transformed = [fn, transformation_matrix](double u) { - return transformation_matrix * fn(u); + auto fn_transformed = [fn, transformation_matrix](double u)->Eigen::VectorXd { + auto result = fn(u); + Eigen::Vector4d v(result.x(), result.y(), result.z(), 1.0); + // return transformation_matrix * fn(u); + auto r = transformation_matrix * v; + Eigen::VectorXd d(4); + d << r(0), r(1), r(2), r(3); + return d; }; // @todo it might be suboptimal that we no longer have the spans now From b9c234457d51dc5ad981cdd7b10738acbe29cc13 Mon Sep 17 00:00:00 2001 From: Richard Brice <37087370+RickBrice@users.noreply.github.com> Date: Wed, 4 Oct 2023 13:25:58 -0700 Subject: [PATCH 15/15] Fixes minor issues with mapping code --- src/ifcgeom/mapping/IfcGradientCurve.cpp | 6 +++--- src/ifcgeom/mapping/mapping.i | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ifcgeom/mapping/IfcGradientCurve.cpp b/src/ifcgeom/mapping/IfcGradientCurve.cpp index 9901c829bb..6794c1731d 100644 --- a/src/ifcgeom/mapping/IfcGradientCurve.cpp +++ b/src/ifcgeom/mapping/IfcGradientCurve.cpp @@ -47,11 +47,11 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcGradientCurve* inst) { } // @todo does this really make sense? - auto composition = [horizontal, vertical](double u) { + auto composition = [horizontal, vertical](double u)->Eigen::VectorXd { auto xy = horizontal->evaluate(u); auto z = vertical->evaluate(u); Eigen::VectorXd vec(3); - vec << xy(0), xy(1), z(0); + vec << xy(0), xy(1), z(2); return vec; }; @@ -69,7 +69,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcGradientCurve* inst) { } auto pwf = taxonomy::make(); - pwf->spans.push_back({ min_length, composition }); + pwf->spans.emplace_back( min_length, composition ); return pwf; } diff --git a/src/ifcgeom/mapping/mapping.i b/src/ifcgeom/mapping/mapping.i index 57de49a16d..5a1dc62572 100644 --- a/src/ifcgeom/mapping/mapping.i +++ b/src/ifcgeom/mapping/mapping.i @@ -113,10 +113,10 @@ BIND(IfcEdge); BIND(IfcEdgeLoop); BIND(IfcPolyline); BIND(IfcPolyLoop); +BIND(IfcCompositeCurve); #ifdef SCHEMA_HAS_IfcGradientCurve BIND(IfcGradientCurve); #endif -BIND(IfcCompositeCurve); BIND(IfcTrimmedCurve); BIND(IfcArbitraryOpenProfileDef); #ifdef SCHEMA_HAS_IfcIndexedPolyCurve