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.
That makes
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]
which gives us sign changes for (100,0)->(0,0) and (100,100)->(0,100).
Hence the line cuts the bottom edge and the top edge respectively.
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)
and fill that in in the line equation:
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
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