Add IObject::distance method

Computes the distance to the first intersection point
on a object regardless if it is entry/exit/glancing
track.
Refs #26690

Framework/Geometry/inc/MantidGeometry/Instrument/Container.h

@@ -57,6 +57,9 @@ public:
   int interceptSurface(Geometry::Track &t) const override {
     return m_shape->interceptSurface(t);
   }
+  double distance(const Geometry::Track &t) const override {
+    return m_shape->distance(t);
+  }
   double solidAngle(const Kernel::V3D &observer) const override {
     return m_shape->solidAngle(observer);
   }

Framework/Geometry/inc/MantidGeometry/Objects/CSGObject.h

@@ -129,7 +129,8 @@ public:
   void write(std::ostream &) const; ///< MCNPX output
 
   // INTERSECTION
-  int interceptSurface(Geometry::Track &) const override;
+  int interceptSurface(Geometry::Track &track) const override;
+  double distance(const Track &track) const override;
 
   // Solid angle - uses triangleSolidAngle unless many (>30000) triangles
   double solidAngle(const Kernel::V3D &observer) const override;

Framework/Geometry/inc/MantidGeometry/Objects/IObject.h

@@ -53,6 +53,7 @@ public:
   virtual int getName() const = 0;
 
   virtual int interceptSurface(Geometry::Track &) const = 0;
+  virtual double distance(const Geometry::Track &) const = 0;
   // Solid angle
   virtual double solidAngle(const Kernel::V3D &observer) const = 0;
   // Solid angle with a scaling of the object

Framework/Geometry/inc/MantidGeometry/Objects/MeshObject.h

@@ -87,6 +87,7 @@ public:
 
   // INTERSECTION
   int interceptSurface(Geometry::Track &) const override;
+  double distance(const Track &track) const override;
 
   // Solid angle - uses triangleSolidAngle unless many (>30000) triangles
   double solidAngle(const Kernel::V3D &observer) const override;

Framework/Geometry/inc/MantidGeometry/Objects/MeshObject2D.h

@@ -49,6 +49,7 @@ public:
       const Kernel::V3D &point) const override; ///< Check if a point is inside
   bool isOnSide(const Kernel::V3D &) const override;
   int interceptSurface(Geometry::Track &ut) const override;
+  double distance(const Geometry::Track &ut) const override;
   MeshObject2D *clone() const override;
   MeshObject2D *
   cloneWithMaterial(const Kernel::Material &material) const override;

Framework/Geometry/src/Objects/CSGObject.cpp

@@ -43,6 +43,10 @@ using namespace Mantid::Geometry;
 using namespace Mantid::Kernel;
 
 namespace {
+
+/// A shift to add/subtract to a point to test if it is an entry/exit point
+constexpr double VALID_INTERCEPT_POINT_SHIFT{2.5e-05};
+
 /**
  * Find the solid angle of a triangle defined by vectors a,b,c from point
  *"observer"
@@ -1094,16 +1098,15 @@ int CSGObject::procString(const std::string &Line) {
 
 /**
  * Given a track, fill the track with valid section
- * @param UT :: Initial track
+ * @param track :: Initial track
  * @return Number of segments added
  */
-int CSGObject::interceptSurface(Geometry::Track &UT) const {
-  int originalCount = UT.count(); // Number of intersections original track
+int CSGObject::interceptSurface(Geometry::Track &track) const {
+  int originalCount = track.count(); // Number of intersections original track
   // Loop over all the surfaces.
-  LineIntersectVisit LI(UT.startPoint(), UT.direction());
-  std::vector<const Surface *>::const_iterator vc;
-  for (vc = m_SurList.begin(); vc != m_SurList.end(); ++vc) {
-    (*vc)->acceptVisitor(LI);
+  LineIntersectVisit LI(track.startPoint(), track.direction());
+  for (auto &surface : m_SurList) {
+    surface->acceptVisitor(LI);
   }
   const auto &IPoints(LI.getPoints());
   const auto &dPoints(LI.getDistance());
@@ -1114,14 +1117,37 @@ int CSGObject::interceptSurface(Geometry::Track &UT) const {
     if (*ditr > 0.0) // only interested in forward going points
     {
       // Is the point and enterance/exit Point
-      const TrackDirection flag = calcValidType(*iitr, UT.direction());
+      const TrackDirection flag = calcValidType(*iitr, track.direction());
       if (flag != TrackDirection::INVALID)
-        UT.addPoint(flag, *iitr, *this);
+        track.addPoint(flag, *iitr, *this);
     }
   }
-  UT.buildLink();
+  track.buildLink();
   // Return number of track segments added
-  return (UT.count() - originalCount);
+  return (track.count() - originalCount);
+}
+
+/**
+ * Compute the distance to the first point of intersection with the surface
+ * @param track Track defining start/direction
+ * @return The distance to the object
+ * @throws std::runtime_error if no intersection was found
+ */
+double CSGObject::distance(const Geometry::Track &track) const {
+  LineIntersectVisit LI(track.startPoint(), track.direction());
+  for (auto &surface : m_SurList) {
+    surface->acceptVisitor(LI);
+  }
+  const auto &distances(LI.getDistance());
+  if (!distances.empty()) {
+    return std::abs(
+        *std::min_element(std::begin(distances), std::end(distances)));
+  } else {
+    std::ostringstream os;
+    os << "Unable to find intersection with object with track starting at "
+       << track.startPoint() << " in direction " << track.direction() << "\n";
+    throw std::runtime_error(os.str());
+  }
 }
 
 /**
@@ -1134,11 +1160,9 @@ int CSGObject::interceptSurface(Geometry::Track &UT) const {
  */
 TrackDirection CSGObject::calcValidType(const Kernel::V3D &point,
                                         const Kernel::V3D &uVec) const {
-  const Kernel::V3D shift(uVec * Kernel::Tolerance * 25.0);
-  const Kernel::V3D testA(point - shift);
-  const Kernel::V3D testB(point + shift);
-  const int flagA = isValid(testA);
-  const int flagB = isValid(testB);
+  const Kernel::V3D shift(uVec * VALID_INTERCEPT_POINT_SHIFT);
+  const int flagA = isValid(point - shift);
+  const int flagB = isValid(point + shift);
   if (!(flagA ^ flagB))
     return Geometry::TrackDirection::INVALID;
   return (flagA) ? Geometry::TrackDirection::LEAVING

Framework/Geometry/src/Objects/MeshObject.cpp

@@ -143,7 +143,6 @@ bool MeshObject::isOnSide(const Kernel::V3D &point) const {
  * @return Number of segments added
  */
 int MeshObject::interceptSurface(Geometry::Track &UT) const {
-
   int originalCount = UT.count(); // Number of intersections original track
   BoundingBox bb = getBoundingBox();
   if (!bb.doesLineIntersect(UT)) {
@@ -167,6 +166,28 @@ int MeshObject::interceptSurface(Geometry::Track &UT) const {
   return UT.count() - originalCount;
 }
 
+/**
+ * Compute the distance to the first point of intersection with the surface
+ * @param track Track defining start/direction
+ * @return The distance to the object
+ * @throws std::runtime_error if no intersection was found
+ */
+double MeshObject::distance(const Track &track) const {
+  Kernel::V3D vertex1, vertex2, vertex3, intersection;
+  TrackDirection unused;
+  for (size_t i = 0; getTriangle(i, vertex1, vertex2, vertex3); ++i) {
+    if (MeshObjectCommon::rayIntersectsTriangle(
+            track.startPoint(), track.direction(), vertex1, vertex2, vertex3,
+            intersection, unused)) {
+      return track.startPoint().distance(intersection);
+    }
+  }
+  std::ostringstream os;
+  os << "Unable to find intersection with object with track starting at "
+     << track.startPoint() << " in direction " << track.direction() << "\n";
+  throw std::runtime_error(os.str());
+}
+
 /**
  * Get intersection points and their in out directions on the given ray
  * @param start :: Start point of ray

Framework/Geometry/src/Objects/MeshObject2D.cpp

@@ -232,13 +232,13 @@ bool MeshObject2D::isOnSide(const Kernel::V3D &point) const {
  * @return Number of segments added
  */
 int MeshObject2D::interceptSurface(Geometry::Track &ut) const {
-
-  int originalCount = ut.count(); // Number of intersections original track
+  const int originalCount =
+      ut.count(); // Number of intersections original track
 
   const auto &norm = m_planeParameters.normal;
-  auto t = -(ut.startPoint().scalar_prod(norm) +
-             m_planeParameters.p0.scalar_prod(norm)) /
-           ut.direction().scalar_prod(norm);
+  const auto t = -(ut.startPoint().scalar_prod(norm) +
+                   m_planeParameters.p0.scalar_prod(norm)) /
+                 ut.direction().scalar_prod(norm);
 
   // Intersects infinite plane. No point evaluating individual segements if this
   // fails
@@ -260,6 +260,38 @@ int MeshObject2D::interceptSurface(Geometry::Track &ut) const {
   return ut.count() - originalCount;
 }
 
+/**
+ * Compute the distance to the first point of intersection with the mesh
+ * @param track Track defining start/direction
+ * @return The distance to the object
+ * @throws std::runtime_error if no intersection was found
+ */
+double MeshObject2D::distance(const Track &ut) const {
+  const auto &norm = m_planeParameters.normal;
+  const auto t = -(ut.startPoint().scalar_prod(norm) +
+                   m_planeParameters.p0.scalar_prod(norm)) /
+                 ut.direction().scalar_prod(norm);
+
+  // Intersects infinite plane. No point evaluating individual segements if this
+  // fails
+  if (t >= 0) {
+    Kernel::V3D intersection;
+    TrackDirection entryExit;
+    for (size_t i = 0; i < m_vertices.size(); i += 3) {
+      if (MeshObjectCommon::rayIntersectsTriangle(
+              ut.startPoint(), ut.direction(), m_vertices[i], m_vertices[i + 1],
+              m_vertices[i + 2], intersection, entryExit)) {
+        // All vertices on plane. So only one triangle intersection possible
+        return intersection.distance(ut.startPoint());
+      }
+    }
+  }
+  std::ostringstream os;
+  os << "Unable to find intersection with object with track starting at "
+     << ut.startPoint() << " in direction " << ut.direction() << "\n";
+  throw std::runtime_error(os.str());
+}
+
 MeshObject2D *MeshObject2D::clone() const {
   return new MeshObject2D(this->m_triangles, this->m_vertices,
                           this->m_material);

Framework/Geometry/test/CSGObjectTest.h

@@ -423,6 +423,24 @@ public:
     checkTrackIntercept(track, expectedResults);
   }
 
+  void testDistanceWithIntersectionReturnsResult() {
+    auto geom_obj = createCappedCylinder();
+    V3D dir(0., 1., 0.);
+    dir.normalize();
+    Track track(V3D(0, 0, 0), dir);
+
+    TS_ASSERT_DELTA(3.0, geom_obj->distance(track), 1e-08)
+  }
+
+  void testDistanceWithoutIntersectionThrows() {
+    auto geom_obj = createCappedCylinder();
+    V3D dir(-1., 0., 0.);
+    dir.normalize();
+    Track track(V3D(-10, 0, 0), dir);
+
+    TS_ASSERT_THROWS(geom_obj->distance(track), std::runtime_error)
+  }
+
   void testTrackTwoIsolatedCubes()
   /**
   Test a track going through an object

Framework/Geometry/test/MeshObjectTest.h

@@ -384,6 +384,24 @@ public:
     checkTrackIntercept(std::move(geom_obj), track, expectedResults);
   }
 
+  void testDistanceWithIntersectionReturnsResult() {
+    auto geom_obj = createCube(3);
+    V3D dir(0., 1., 0.);
+    dir.normalize();
+    Track track(V3D(0, 0, 0), dir);
+
+    TS_ASSERT_DELTA(3.0, geom_obj->distance(track), 1e-08)
+  }
+
+  void testDistanceWithoutIntersectionThrows() {
+    auto geom_obj = createCube(3);
+    V3D dir(-1., 0., 0.);
+    dir.normalize();
+    Track track(V3D(-10, 0, 0), dir);
+
+    TS_ASSERT_THROWS(geom_obj->distance(track), std::runtime_error)
+  }
+
   void testTrackTwoIsolatedCubes()
   /**
   Test a track going through two objects