Return to Answer

The Abstract-Is-Better way

Below we talk about n-dimensional space. I'll use "n-dimensional hyper-xyz" to denote something that would be called xyz in 3 dimensions. For example, a 2-dimensional hyper-block is a retangle, etc).

Say you want to divide an n-dimensional hyper-block with a n-1 dimensional hyper-plane. The hyper-plane has a normal unit vector N and its points X satisfy:

N·X + d = 0

where N·X is a scalar: the dot product of the vectors N and X, and d is a scalar giving the distance from the origin to the hyper-plane in units of N: if d is negative then the hyper-plane just lays in the direction of -N from the origin.

The n-dimensional hyper-block can be thought of as being spanned by a set of n orthogonal "base vectors" B=(v0, v1, v2, ...) and has 2n hyper-planes as borders and 2^n corners. Lets denote these corners with an n-dimensional boolean vector (aka, C=(0,0,1,1,0,1,...)) where the i-th boolean gives us whether that corner is a part of the lower or higher hyper-plane that is perpendicular to the i-th base vector v_i.

Note the exact coordinates of corner C is now given by B·C, aka corner (0,0,0...) is given by O + O + O + ... = O (where O is the origin) and (1,1,1,...) is given by v0 + v1 + v2 + ....

We have to run over all corners and see on which side of the hyper-plane they are on; that is, if the distance from the corner to the hyper-plane requires adding or subtracting (a fraction of) N.

We can determine this distance k by solving the line equation for B·C + kN (aka, how many times do we have to add the normal N to the corner to reach the hyper-plane?):

N·(B·C + kN) + d = 0

and since N is a unit vector N·kN = k N·N = k * 1 = k, therefore

k = - N·(B·C) - d

If this k is negative then this Corner is on one side, if k is positive is it on the other side and if k=0 then the hyper-plane goes through that corner and it doesn't matter which side you pick, so lets say:

k' = 0 if k <= 0, and k' = 1 if k > 0.

and we have as many k' values as corners (one for each C): 2^n of them. You can draw a line between any two corner points that only differ in one bit of C: n edges from every corner, where then you count every edge twice; so there are 2^n * n / 2 = 2^(n-1) * n edges. However, you can also just run over all corners and then only toggle a single 1 (if there is one) to a 0 (for each 1 that that corner has); if the corresponding k' changes value then the corresponding edge is cut by your hyper-plane and you can calculate the intersection point.

Lets apply the above to an axis aligned rectangle (where one corner is the origin); in that case:

B = [(100,0), (0,100)]

Having two points on the line: P1 = (50,50) and P2 = (50,150), the line goes -say- from P1 to P2: P2-P1 = (p2_x - p1_x, p2_y - p1_y) and the normal to that line N = (-(p2_y - p1_y), p2_x - p1_x). And the distance to the origin is for example d = -N·P1. Working that out we get:

N = (-(150 - 50), 50 - 50) = (-100, 0) = (-1, 0) (normalized to length 1)
d = -(-1*50 + 0*50) = 50

Assuming B is column vector, lets make C a matrix that includes all corners at once:

C = [(0,0),(1,0),(0,1),(1,1)]

where the first element of each column vector tells us if b0 participates and the second element if b1 does.

B·C = [(0,0), (100,0), (0,100), (100,100)]

the four corners. And the corresponding k values:

K = -N·(B·C) - d = [-50,50,-50,50]

To calculate where, just parameterize the edge using the two corners that made the sign of k change - say C1 and C2 (e.g. (100,0) and (0,0)):

Edge: C1 + g(C2 - C1)

N·(C1 + g(C2 - C1)) + d = 0

to find g:

g = -(N·C1 + d) / (N·(C2 - C1))

The intersection point then is C1 + g(C2 - C1).

For example, with C1=(100,0) and C2=(0,0):

g = -((-1,0)·(100,0) + 50)/((-1,0)·((0,0)-(100,0))) =
  = -(-100 + 50)/((-1,0)·(-100,0)) =
  = 50/100 = 1/2

And thus the intersection is at (100,0) + 1/2(-100,0) = (50,0). While for C1=(100,100) and C2=(0,100) we get (50,100).

PS If N·(C2 - C1) = 0 we can't devide through it; this means that the hyper-plane goes through both corners and intersects the edge everywhere. In that case one should have picked the other possibility for when k=0 for at least one of the corners.

24 hours later

Online test case: https://wandbox.org/permlink/5MTAowhWLxUkQVK9 Latest version: https://github.com/CarloWood/machine-learning/blob/master/src/intersection_points.h

The Abstract-Is-Better way

Below we talk about \$n\$-dimensional space. I'll use "\$n\$-dimensional hyper-xyz" to denote something that would be called xyz in 3 dimensions. For example, a 2-dimensional hyper-block is a rectangle, etc).

Say you want to divide an \$n\$-dimensional hyper-block with a \$n-1\$ dimensional hyper-plane. The hyper-plane has a normal unit vector \$N\$ and its points \$X\$ satisfy:

$$N\cdot X + d = 0$$

where \$N\cdot X\$ is a scalar: the dot product of the vectors \$N\$ and \$X\$, and \$d\$ is a scalar giving the distance from the origin to the hyper-plane in units of \$N\$: if \$d\$ is negative then the hyper-plane just lays in the direction of \$-N\$ from the origin.

The \$n\$-dimensional hyper-block can be thought of as being spanned by a set of \$n\$ orthogonal "base vectors" \$B=(v_0, v_1, v_2, ...)\$ and has \$2n\$ hyper-planes as borders and \$2^n\$ corners. Let's denote these corners with an \$n\$-dimensional boolean vector (aka, \$C=(0,0,1,1,0,1,...)\$) where the \$i\$th boolean gives us whether that corner is a part of the lower or higher hyper-plane that is perpendicular to the \$i\$th base vector \$v_i\$.

Note the exact coordinates of corner \$C\$ is now given by \$B\cdot C\$, aka corner \$(0,0,0...)\$ is given by \$O + O + O + ... = O\$ (where \$O\$ is the origin) and \$(1,1,1,...)\$ is given by \$v0 + v1 + v2 + ....\$

We have to run over all corners and see on which side of the hyper-plane they are on; that is, if the distance from the corner to the hyper-plane requires adding or subtracting (a fraction of) \$N\$.

We can determine this distance \$k\$ by solving the line equation for \$B\cdot C + kN\$ (aka, how many times do we have to add the normal \$N\$ to the corner to reach the hyper-plane?):

$$N\cdot(B\cdot C + kN) + d = 0$$

and since \$N\$ is a unit vector \$N\cdot kN = k N\cdot N = k * 1 = k\$, therefore:

$$k = - N\cdot (B\cdot C) - d$$

If this \$k\$ is negative then this Corner is on one side, if \$k\$ is positive is it on the other side and if \$k=0\$ then the hyper-plane goes through that corner and it doesn't matter which side you pick, so let's say:

$$ \begin{align} k' &= 0 \text{ if } k <= 0 \text{ and}\\ k' &= 1 \text{ if } k > 0 \end{align} $$

and we have as many \$k'\$ values as corners (one for each \$C\$): \$2^n\$ of them. You can draw a line between any two corner points that only differ in one bit of \$C\$: \$n\$ edges from every corner, where then you count every edge twice; so there are \$2^n \frac{n}{2} = 2^{n-1} n\$ edges. However, you can also just run over all corners and then only toggle a single 1 (if there is one) to a 0 (for each 1 that that corner has); if the corresponding \$k'\$ changes value then the corresponding edge is cut by your hyper-plane and you can calculate the intersection point.

Let's apply the above to an axis aligned rectangle (where one corner is the origin); in that case:

\$B = [(100,0), (0,100)]\$

Having two points on the line: \$P_1 = (x_1,y_1)=(50,50)\$ and \$P_2 = (x_2,y_2)=(50,150)\$, the line goes -say- from \$P_1\$ to \$P_2\$: \$P_2-P_1 = (x_2 - x_1, y_2 - y_1)\$ and the normal to that line \$N = (-(y_2 - y_1), x_2 - x_1)\$. And the distance to the origin is for example \$d = -N\cdot P_1\$. Working that out we get:

$$ \begin{align} N &= (-(150 - 50), 50 - 50) = (-100, 0) = (-1, 0) \textit{ (normalized to length 1)} \\ d &= -(-1*50 + 0*50) = 50 \end{align} $$

Assuming \$B\$ is column vector, let's make \$C\$ a matrix that includes all corners at once:

$$C = [(0,0),(1,0),(0,1),(1,1)]$$

where the first element of each column vector tells us if \$b_0\$ participates and the second element if \$b_1\$ does.

$$B\cdot C = [(0,0), (100,0), (0,100), (100,100)]$$

the four corners. And the corresponding \$k\$ values:

$$K = -N\cdot (B\cdot C) - d = [-50,50,-50,50]$$

To calculate where, just parameterize the edge using the two corners that made the sign of \$k\$ change - say \$C_1\$ and \$C_2\$ (e.g. (100,0) and (0,0)):

Edge: \$C_1 + g(C_2 - C_1)\$

$$N\cdot (C_1 + g(C_2 - C_1)) + d = 0$$

to find \$g\$:

$$g = -\frac{N\cdot C_1 + d}{N\cdot (C_2 - C_1)}$$

The intersection point then is \$C_1 + g(C_2 - C_1)\$.

For example, with \$C_1=(100,0)\$ and \$C_2=(0,0)\$: $$ \begin{align} g &= -\frac{(-1,0)\cdot (100,0) + 50}{(-1,0)\cdot ((0,0)-(100,0))} \\ &= -\frac{-100 + 50}{(-1,0)\cdot (-100,0)} \\ &= \frac{50}{100} \\ &= \frac{1}{2} \\ \end{align} $$ And thus the intersection is at (100,0) + 1/2(-100,0) = (50,0). While for \$C_1=(100,100)\$ and \$C_2=(0,100)\$ we get (50,100).

PS If \$N\cdot (C_2 - C_1) = 0\$ we can't devide through it; this means that the hyper-plane goes through both corners and intersects the edge everywhere. In that case one should have picked the other possibility for when \$k=0\$ for at least one of the corners.

Online test case 
Latest version

24 hours later

I wrote an implementation:

#include <array>
#include <vector>
#include <concepts>
#include <stdexcept>
#include <iostream>

namespace intersections {

// An n-dimensional vector (a column vector with n elements).
template<std::floating_point FloatType, int n>
struct Vector
{
  std::array<FloatType, n> v_;                          // The elements of the vector.

  // Construct a zeroed vector.
  Vector() : v_{} { }

  // Initialize this Vector from an initializer list.
  Vector(std::initializer_list<FloatType> v)
  {
    if (v.size() != n)
      throw std::invalid_argument("Initializer list must have exactly n elements.");
    std::copy(v.begin(), v.end(), v_.begin());
  }

  // Element access.
  FloatType& operator[](int i) { return v_[i]; }
  FloatType operator[](int i) const { return v_[i]; }

  // Add v2.
  Vector& operator+=(Vector const& v2)
  {
    for (int i = 0; i < n; ++i)
      v_[i] += v2[i];
    return *this;
  }

  // Add v1 and v2.
  friend Vector operator+(Vector const& v1, Vector const& v2)
  {
    Vector result(v1);
    result += v2;
    return result;
  }

  // Subtract v2.
  Vector& operator-=(Vector const& v2)
  {
    for (int i = 0; i < n; ++i)
      v_[i] -= v2[i];
    return *this;
  }

  // Subtract v2 from v1.
  friend Vector operator-(Vector const& v1, Vector const& v2)
  {
    Vector result(v1);
    result -= v2;
    return result;
  }

  // Return the dot product of v1 and v2.
  friend FloatType operator*(Vector const& v1, Vector const& v2)
  {
    FloatType result = 0;
    for (int i = 0; i < n; ++i)
      result += v1[i] * v2[i];
    return result;
  }

  // Elementwise multiply with scalar g.
  friend Vector operator*(FloatType g, Vector const& v)
  {
    Vector result(v);
    for (int i = 0; i < n; ++i)
      result[i] *= g;
    return result;
  }

  // For debugging purposes.
  void print_on(std::ostream& os) const
  {
    char const* sep = "";
    os << '(';
    for (int i = 0; i < n; ++i)
    {
      os << sep << v_[i];
      sep = ", ";
    }
    os << ')';
  }

  friend std::ostream& operator<<(std::ostream& os, Vector const& v)
  {
    v.print_on(os);
    return os;
  }
};

// An n-1 dimensional hyper-plane, orthogonal to a given normal unit vector N, with offset dN from the origin.
template<std::floating_point FloatType, int n>
struct HyperPlane
{
  using VectorType = Vector<FloatType, n>;

  VectorType N_;                                        // The unit normal of the hyper-plane.
  FloatType d_;                                         // The signed distance (in Normal vectors) from origin to hyper-plane.

  // Create a hyper-plane that satisfies N·X + d = 0.
  HyperPlane(VectorType const& N, FloatType d) : N_(N), d_(d) { }

  // Return the number of N_'s that have to be added to P to end up on this HyperPlane.
  // The sign tells you which "side" of the hyper-plane the point is on.
  FloatType distance(VectorType const& P) const
  {
    FloatType dist = -(N_ * P + d_);
    return dist;
  }

  // Return intersection of the line through C1 and C2 with this HyperPlane.
  VectorType intersection(VectorType const& C1, VectorType const& C2) const
  {
    // Let E be a line through C1 and C2: E: C1 + g(C2 - C1), where g parameterizes the points on E.
    // Fill that in in the line equation to find the intersection:
    // N·(C1 + g(C2 - C1)) + d = 0 --> N·C1 + d + g N·(C2 - C1) = 0 --> g = -(N·C1 + d) / N·(C2 - C1)
    VectorType diff = C2 - C1;
    FloatType g = -(N_ * C1 + d_) / (N_ * diff);
    return C1 + g * diff;
  }
};

template<std::floating_point FloatType, int n>
struct HyperBlock
{
  using VectorType = Vector<FloatType, n>;
  static constexpr int number_of_corners = 1 << n;

  std::array<VectorType, number_of_corners> C_;         // The 2^n corners of the hyper-block.

  // Construct an axis-aligned hyper-block from two adjecent corner vectors.
  HyperBlock(VectorType const& c1, VectorType const& c2)
  {
    VectorType base = c2 - c1;
    for (int ci = 0; ci < number_of_corners; ++ci)
    {
      VectorType c;
      for (int d = 0; d < n; ++d)
      {
        int bit = 1 << d;
        if ((ci & bit))
          c[d] += base[d];
      }
      C_[ci] = c1 + c;
    }
  }

  std::vector<VectorType> intersection_points(HyperPlane<FloatType, n> const& plane);
};

template<std::floating_point FloatType, int n>
std::vector<typename HyperBlock<FloatType, n>::VectorType> HyperBlock<FloatType, n>::intersection_points(HyperPlane<FloatType, n> const& plane)
{
  std::vector<VectorType> intersections;

  std::array<bool, number_of_corners> side;
  for (int ci = 0; ci < number_of_corners; ++ci)
    side[ci] = plane.distance(C_[ci]) <= 0;

  for (int ci = 0; ci < number_of_corners; ++ci)
  {
    for (int d = 0; d < n; ++d)
    {
      int bit = 1 << d;
      int ci2 = ci & ~bit;
      if (side[ci] != side[ci2])
      {
        // Found two corners on opposite sides of the hyper-plane.
        intersections.push_back(plane.intersection(C_[ci], C_[ci2]));
      }
    }
  }

  return intersections;
}

} // namespace intersections

Online test case: https://wandbox.org/permlink/5MTAowhWLxUkQVK9 Latest version: https://github.com/CarloWood/machine-learning/blob/master/src/intersection_points.h

PS If N·(C2 - C1) = 0 we can't devide through it; this means that the hyper-plane goes through both corners and intersects the edge everywhere. In that case one should have picked the other possibility for when k=0 for at least one of the corners.