diff --git a/Common/DataModel/vtkBox.cxx b/Common/DataModel/vtkBox.cxx
index 083da4aeef75a3ab4cf476ae1707153fd228779a..cd8c9cbf92f3977be9f22f894fa83b4e7fe9955c 100644
--- a/Common/DataModel/vtkBox.cxx
+++ b/Common/DataModel/vtkBox.cxx
@@ -20,6 +20,7 @@
 
 #include <algorithm> // for sorting
 #include <cassert>
+#include <limits> // for IntersectWithInfiniteLine
 #include <vector> // for IntersectWithPlane
 
 vtkStandardNewMacro(vtkBox);
@@ -319,7 +320,7 @@ void vtkBox::EvaluateGradient(double x[3], double n[3])
 // coordinate along line. (Notes: the intersection ray dir[3] is NOT
 // normalized.  Valid intersections will only occur between 0<=t<=1.)
 char vtkBox::IntersectBox(
-  double bounds[6], const double origin[3], double dir[3], double coord[3], double& t)
+  const double bounds[6], const double origin[3], const double dir[3], double coord[3], double& t)
 {
   bool inside = true;
   char quadrant[3];
@@ -523,6 +524,68 @@ int vtkBox::IntersectWithLine(const double bounds[6], const double p1[3], const
   return 1;
 }
 
+//----------------------------------------------------------------------------
+bool vtkBox::IntersectWithInfiniteLine(const double bounds[6], const double p1[3],
+  const double p2[3], double& t1, double& t2, double x1[3], double x2[3], int& plane1, int& plane2)
+{
+  plane1 = -1;
+  plane2 = -1;
+  t1 = -std::numeric_limits<double>::infinity();
+  t2 = std::numeric_limits<double>::infinity();
+
+  for (int j = 0; j < 3; j++)
+  {
+    for (int k = 0; k < 2; k++)
+    {
+      // Compute distances of p1 and p2 from the plane along the plane normal
+      int i = 2 * j + k;
+      double t =
+        std::abs(bounds[i] - p1[j]) < VTK_DBL_MIN ? 0.0 : (bounds[i] - p1[j]) / (p2[j] - p1[j]);
+      double xface = p1[(j + 1) % 3] + t * (p2[(j + 1) % 3] - p1[(j + 1) % 3]),
+             yface = p1[(j + 2) % 3] + t * (p2[(j + 2) % 3] - p1[(j + 2) % 3]);
+      // if (xface, yface) is inside the current face
+      if (xface >= bounds[(2 * j + 2) % 6] && xface <= bounds[(2 * j + 3) % 6] &&
+        yface >= bounds[(2 * j + 4) % 6] && yface <= bounds[(2 * j + 5) % 6])
+      {
+        if (t1 == -std::numeric_limits<double>::infinity())
+        {
+          t1 = t;
+          plane1 = 2 * j + k;
+        }
+        else if (t >= t1)
+        {
+          t2 = t;
+          plane2 = 2 * j + k;
+          break;
+        }
+        else
+        {
+          t2 = t1;
+          t1 = t;
+          plane2 = plane1;
+          plane1 = 2 * j + k;
+          break;
+        }
+      }
+    }
+  }
+
+  if (x1)
+  {
+    x1[0] = p1[0] + t1 * (p2[0] - p1[0]);
+    x1[1] = p1[1] + t1 * (p2[1] - p1[1]);
+    x1[2] = p1[2] + t1 * (p2[2] - p1[2]);
+  }
+  if (x2)
+  {
+    x2[0] = p1[0] + t2 * (p2[0] - p1[0]);
+    x2[1] = p1[1] + t2 * (p2[1] - p1[1]);
+    x2[2] = p1[2] + t2 * (p2[2] - p1[2]);
+  }
+
+  return t1 != -std::numeric_limits<double>::infinity();
+}
+
 //----------------------------------------------------------------------------
 vtkTypeBool vtkBox::IntersectWithPlane(double bounds[6], double origin[3], double normal[3])
 {
diff --git a/Common/DataModel/vtkBox.h b/Common/DataModel/vtkBox.h
index 97ce8e65e596fa45e92a36a2fde48295d77c86b6..e5cdfe418b10d15bf3321c282738860555e66396 100644
--- a/Common/DataModel/vtkBox.h
+++ b/Common/DataModel/vtkBox.h
@@ -96,8 +96,8 @@ public:
    * dir[3] is NOT normalized.  Valid intersections will only occur between
    * 0<=t<=1.)
    */
-  static char IntersectBox(
-    double bounds[6], const double origin[3], double dir[3], double coord[3], double& t);
+  static char IntersectBox(const double bounds[6], const double origin[3], const double dir[3],
+    double coord[3], double& t);
 
   /**
    * Intersect a line with the box.  Give the endpoints of the line in
@@ -114,6 +114,19 @@ public:
   static int IntersectWithLine(const double bounds[6], const double p1[3], const double p2[3],
     double& t1, double& t2, double x1[3], double x2[3], int& plane1, int& plane2);
 
+  /**
+   * Same method as vtkBox::IntersectWithLine, except that t1 and t2 can be outside of [0,1].
+   * t1 is the distance of x1 to p1 in parametric coordinates, and t2 is the distance of x2 to p1
+   * in parametric coordinates as well.
+   * In vtkBox::IntersectWithInLine, it is assumed that [p1,p2] is a segment, here, it is
+   * assumed that it is a line with no ends.
+   * t1 <= t2, which means that x1 is always "before" x2 on the line parameterized by [p1,p2].
+   * x1 and x2 can be set to nullptr without crash.
+   */
+  static bool IntersectWithInfiniteLine(const double bounds[6], const double p1[3],
+    const double p2[3], double& t1, double& t2, double x1[3], double x2[3], int& plane1,
+    int& plane2);
+
   /**
    * Plane intersection with the box. The plane is infinite in extent and
    * defined by an origin and normal. The function indicates whether the