diff --git a/GeoModelCore/GeoModelHelpers/GeoModelHelpers/TransformToStringConverter.h b/GeoModelCore/GeoModelHelpers/GeoModelHelpers/TransformToStringConverter.h
index 5e15a85a0a572d2913ce0f816baf5e5c8a5404b9..cb912668e5596ffb0b3db9e37c784805ff3bb428 100644
--- a/GeoModelCore/GeoModelHelpers/GeoModelHelpers/TransformToStringConverter.h
+++ b/GeoModelCore/GeoModelHelpers/GeoModelHelpers/TransformToStringConverter.h
@@ -28,6 +28,10 @@ namespace GeoTrf{
std::string toString(const Transform3D& trans, bool useCoordAngles = false, int precision = 4);
std::string toString(const EulerAngles& angles, int precision = 4);
std::string toString(const CoordEulerAngles& angles, int precision = 4);
+
+ std::string toString(const Line3D& line, int precision = 4);
+ std::string toString(const Plane3D& plane, int precision = 4);
+
}
#endif
\ No newline at end of file
diff --git a/GeoModelCore/GeoModelHelpers/src/TransformToStringConverter.cxx b/GeoModelCore/GeoModelHelpers/src/TransformToStringConverter.cxx
index 13f5c885b0f6c14603af2a36abc9f360d774ee41..cef0fb0d9bec1dbebdc8fcc3fe4e9962747055f1 100644
--- a/GeoModelCore/GeoModelHelpers/src/TransformToStringConverter.cxx
+++ b/GeoModelCore/GeoModelHelpers/src/TransformToStringConverter.cxx
@@ -56,4 +56,15 @@ namespace GeoTrf{
ostr <<"gamma: "<<angles.gamma * toDeg;
return ostr.str();
}
+ std::string toString(const Line3D& line, int precision) {
+ std::stringstream ostr{};
+ ostr<<"line --- "<<toString(line.position(), precision)<<" + "<<toString(line.direction(), precision);
+ return ostr.str();
+ }
+ std::string toString(const Plane3D& plane, int precision) {
+ std::stringstream ostr{};
+ ostr << std::setiosflags(std::ios::fixed) << std::setprecision(precision);
+ ostr<<"plane -- normal: "<<toString(plane.normal(), 4)<<" - C: "<<plane.offSet();
+ return ostr.str();
+ }
}
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDefinitions.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDefinitions.h
index efe8538f8215b0d0feabacff3ed5d855742cfe8e..d4a85fa6fe12ae022e1204e89d89b84a0f38f801 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDefinitions.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDefinitions.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
#ifndef GEOMODELKERNEL_GEODEFINITIONS_H
@@ -14,64 +14,115 @@
#else
#include <Eigen/Dense>
#endif
+#include <optional>
namespace GeoTrf {
using Rotation3D = Eigen::Quaternion<double>;
using Translation3D = Eigen::Translation<double, 3>;
using AngleAxis3D = Eigen::AngleAxisd;
using Transform3D = Eigen::Affine3d;
- template<size_t N> using VectorN = Eigen::Matrix<double, N, 1>;
+ template<int N> using VectorN = Eigen::Matrix<double, N, 1>;
using Vector3D = VectorN<3>;
using Vector2D = VectorN<2>;
using RotationMatrix3D = Eigen::Matrix<double, 3, 3>;
using Diagonal3DMatrix = Eigen::DiagonalMatrix <double, 3>;
- /// check if we can use Eigen Hyperplane for this
- struct Plane3D {
- protected:
- double a_{0.};
- double b_{0.};
- double c_{1.};
- double d_{0.};
-
- public:
- /**
- * Default constructor - creates plane z=0. */
- Plane3D() = default;
- /**
- * Constructor from four numbers - creates plane a*x+b*y+c*z+d=0. */
- Plane3D(double a1, double b1, double c1, double d1)
- : a_(a1), b_(b1), c_(c1), d_(d1)
- {}
+ /** @brief Returns a unit vector from an arbitrary vector */
+ template <int N>
+ VectorN<N> unit(const VectorN<N>& a);
+ /** @brief Helper struct describing a Line in 3D space */
+ class Plane3D;
+ class Line3D{
+ private:
+ Vector3D m_pos{Vector3D::Zero()};
+ Vector3D m_dir{Vector3D::Zero()};
+ public:
+ /** @brief Default empty constructor */
+ Line3D() = default;
+ /*** @brief Constructor taking a point on the line and the line direction
+ * @param pos: Point on the line
+ * @param dir: Direction vector (Needs to be normalized) */
+ Line3D(const Vector3D& pos, const Vector3D& dir):
+ m_pos{pos},
+ m_dir{dir} {}
+ /** @brief Returns the position of the line */
+ const Vector3D& position() const {
+ return m_pos;
+ }
+ /** @brief Returns the direction of the line */
+ const Vector3D& direction() const {
+ return m_dir;
+ }
+ /** @brief Travels a distance x along the line */
+ GeoTrf::Vector3D travel(const double x) const {
+ return position() + x * direction();
+ }
+ /** @brief Returns the distance of the point along the line */
+ double distance(const Vector3D& v) const {
+ return direction().dot(v-position());
+ }
+ /** @brief Returns the intersection point with another line
+ * @param other: Other line to intersect */
+ std::optional<Vector3D> intersect(const Line3D& other) const;
+ /** @brief the intersection point with a plane */
+ std::optional<Vector3D> intersect(const Plane3D& other) const;
+ };
+ /** @brief Helper struct describing a Plane in 3D dimensional space */
+ class Plane3D {
+ private:
+ Vector3D m_normal{Vector3D::Zero()};
+ double m_offSet{0.};
+ public:
+ /**
+ * Default constructor - creates plane z=0. */
+ Plane3D() = default;
+ /**
+ * Constructor from four numbers - creates plane a*x+b*y+c*z+d=0. */
+ Plane3D(double a1, double b1, double c1, double d1):
+ m_normal{a1,b1,c1},
+ m_offSet{d1} {}
- /**
- * Constructor from normal and point. */
- Plane3D(const Vector3D& n, const Vector3D& p)
- : a_(n.x()), b_(n.y()), c_(n.z()), d_(-n.dot(p))
- {}
+ /*** Constructor from normal and point. */
+ Plane3D(const Vector3D& n, const Vector3D& p):
+ m_normal{n},
+ m_offSet{-n.dot(p)} {}
+
+ /** @brief Constructor from two direction vectors and one point in the plane
+ * @param dir1: First direction along the plane
+ * @param dir2: Second direction vector along the plane
+ * @param point: Arbitrary point inside the plane */
+ Plane3D(const Vector3D& dir1,
+ const Vector3D& dir2,
+ const Vector3D& point):
+ Plane3D{unit(dir1.cross(dir2)), point} {}
/**
* Returns the a-coefficient in the plane equation: a*x+b*y+c*z+d=0. */
- double a() const { return a_; }
+ double a() const { return m_normal.x(); }
/**
* Returns the b-coefficient in the plane equation: a*x+b*y+c*z+d=0. */
- double b() const { return b_; }
+ double b() const { return m_normal.y(); }
/**
* Returns the c-coefficient in the plane equation: a*x+b*y+c*z+d=0. */
- double c() const { return c_; }
+ double c() const { return m_normal.z(); }
/**
* Returns the free member of the plane equation: a*x+b*y+c*z+d=0. */
- double d() const { return d_; }
+ double d() const { return m_offSet; }
- /**
- * Returns normal. */
- Vector3D normal() const { return Vector3D(a_,b_,c_); }
+ /** @brief Returns the plane normal vector */
+ const Vector3D& normal() const { return m_normal; }
+ /** @brief Returns the offset or the d of the plane equation*/
+ double offSet() const { return m_offSet; }
+ /** @brief Returns the intersection of the plane with a line */
+ std::optional<Vector3D> intersect(const Line3D& other) const;
+ /** @brief Returns whether a point is inside a plane */
+ bool inside(const Vector3D& point) const;
};
class Translate3D : public Transform3D {
public:
- Translate3D(const GeoTrf::Vector3D& v):
+ Translate3D(const Vector3D& v):
Translate3D{v.x(),v.y(),v.z()}{}
Translate3D(double x, double y, double z)
@@ -264,5 +315,5 @@ namespace GeoTrf {
virtual ~GeoTransformRT() = default;
};
}
-
+#include "GeoModelKernel/GeoDefinitions.icc"
#endif // GEOMODELKERNEL_GEODEFINITIONS_H
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDefinitions.icc b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDefinitions.icc
new file mode 100644
index 0000000000000000000000000000000000000000..90c2d6f9b68c15af57a654e99392159440ba51b3
--- /dev/null
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDefinitions.icc
@@ -0,0 +1,38 @@
+/*
+ Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
+*/
+#ifndef GEOMODELKERNEL_GEODEFINITIONS_ICC
+#define GEOMODELKERNEL_GEODEFINITIONS_ICC
+
+#include <climits>
+namespace GeoTrf{
+ template <int N>
+ VectorN<N> unit(const VectorN<N>& a) {
+ return a / std::sqrt(a.dot(a));
+ }
+ inline std::optional<Vector3D> Line3D::intersect(const Line3D& other) const {
+ const double dirDots = direction().dot(other.direction());
+ const double divisor = 1. - dirDots*dirDots;
+ if (std::abs(divisor) < std::numeric_limits<double>::epsilon()){
+ return std::nullopt;
+ }
+ const Vector3D AminusB = other.position() - position();
+ const double lambda = (AminusB.dot(direction()) - AminusB.dot(other.direction()) * dirDots) / divisor;
+ return std::make_optional<Vector3D>(position()+ lambda*direction());
+ }
+ inline std::optional<Vector3D> Line3D::intersect(const Plane3D& other) const {
+ const double normDot = other.normal().dot(direction());
+ if (std::abs(normDot) < std::numeric_limits<double>::epsilon()) {
+ return std::nullopt;
+ }
+ const double lambda = (other.offSet() - position().dot(other.normal())) / normDot;
+ return std::make_optional<Vector3D>(position() + lambda * direction());
+ }
+ inline std::optional<Vector3D> Plane3D::intersect(const Line3D& other) const {
+ return other.intersect(*this);
+ }
+ inline bool Plane3D::inside(const Vector3D& point) const {
+ return (normal().dot(point) - offSet()) < std::numeric_limits<double>::epsilon();
+ }
+}
+#endif
\ No newline at end of file
diff --git a/GeoModelCore/GeoModelKernel/tests/testPlaneAndLine.cxx b/GeoModelCore/GeoModelKernel/tests/testPlaneAndLine.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..46fd96360aaeec77557f58184cba6a9b1221c18d
--- /dev/null
+++ b/GeoModelCore/GeoModelKernel/tests/testPlaneAndLine.cxx
@@ -0,0 +1,82 @@
+/*
+ Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+ */
+#include "GeoModelKernel/GeoDefinitions.h"
+#include "GeoModelKernel/Units.h"
+
+#include <stdlib.h>
+#include <iostream>
+
+using namespace GeoTrf;
+
+constexpr double tolerance = 1.e-9;
+std::ostream& operator<<( std::ostream& ostr, const Vector3D& v) {
+ ostr<<"("<<v.x()<<", "<<v.y()<<", "<<v.z()<<")";
+ return ostr;
+}
+std::ostream & operator<<(std::ostream& ostr, const Line3D& l) {
+ ostr<<"line "<<l.position()<<" + "<<l.direction();
+ return ostr;
+}
+std::ostream& operator<<(std::ostream& ostr, const Plane3D& p) {
+ ostr<<"Plane with normal: "<<p.normal()<<" & offset: "<<p.offSet();
+ return ostr;
+}
+
+#define ERROR(msg) \
+ { \
+ std::cerr<<__FILE__<<":"<<__LINE__<<" - "<<msg<<std::endl; \
+ return EXIT_FAILURE; \
+ } \
+
+int main() {
+
+
+ const Line3D lineA{ Vector3D{10.,0.,0.}, Vector3D::UnitX()};
+ const Line3D lineB{ Vector3D{30.,0.,0.}, Vector3D::UnitY()};
+
+ const std::optional<Vector3D> isectAtoB = lineA.intersect(lineB);
+
+ if (!isectAtoB) ERROR("The two lines A "<<lineA<<" & B "<<lineB<<" are not intersecting");
+
+ const std::optional<Vector3D> isectBtoA = lineB.intersect(lineA);
+ if (!isectAtoB) ERROR("The two lines B "<<lineB<<" & A "<<lineA<<" are not intersecting");
+ /// Check that the two lines intersect in the same point
+ const GeoTrf::Vector3D iSectDiff = (*isectAtoB) - (*isectBtoA);
+ if (iSectDiff.dot(iSectDiff) > tolerance)
+ ERROR("Intersections ended up at two different points "<<(*isectAtoB)<<" vs. "<<(*isectBtoA));
+ const GeoTrf::Vector3D iSectToExpect = (*isectAtoB) - Vector3D{30.,0.,0.};
+ if (iSectToExpect.dot(iSectToExpect) > tolerance)
+ ERROR("Intersections ended up at two different points "<<(*isectAtoB)<<" vs. "<<lineB.position());
+
+ /// Now make it a bit more difficult
+ for (const double dist : {-10., 10., 45. ,100., 666., 1789.}) {
+ for (double angle =0; angle<=180.*GeoModelKernelUnits::deg; angle+=15.*GeoModelKernelUnits::deg){
+ const GeoTrf::Vector3D dir{0.,std::cos(angle), std::sin(angle)};
+ const Line3D test{lineA.travel(dist) - dist*dir, dir};
+
+ auto isectWithTest = lineA.intersect(test);
+ if (!isectWithTest) ERROR("Intersection of "<<lineA<<" with "<<test<<" failed.");
+ const Vector3D backTest = (*isectWithTest) - lineA.travel(lineA.distance(*isectWithTest));
+ if (backTest.dot(backTest) > tolerance)
+ ERROR("Intersection with "<<lineA<<" recuperated a different travel distance: "
+ <<lineA.distance(*isectWithTest)<<" vs. "<<dist);
+ if (std::abs(dist - test.distance(*isectWithTest))> tolerance)
+ ERROR("Intersection with "<<test<<" recuperated a different travel distance: "
+ <<test.distance(*isectWithTest)<<" vs. "<<dist);
+ }
+ }
+ Plane3D plane{lineA.direction(), lineB.direction(), lineA.position()};
+
+ if (!plane.inside(lineA.position()))
+ ERROR("The point "<<lineA.position()<<" defining the plane "<<plane<<" is not inside");
+
+ Line3D testLine{300.*plane.normal(), plane.normal()};
+ auto planeIsect = testLine.intersect(plane);
+ if(!planeIsect || !plane.inside(*planeIsect))
+ ERROR("The intersection point "<<planeIsect.value_or(Vector3D::Zero())<<" is outside plane"
+ <<plane);
+
+ return EXIT_SUCCESS;
+}
+