Intersection of a Line Segment and a Rectangle

One easy approach is: check intersection with each rectangle line.

Here follow my lua code for that (Line line intersection)

function getVec2(x_,y_)
    vec2={}
    vec2.x=x_
    vec2.y=y_

    function vec2:dist() --set relative position
        return math.sqrt(vec2.x*vec2.x + vec2.y*vec2.y)
    end

    function vec2:sqrdist() --set relative position
        return vec2.x*vec2.x + vec2.y*vec2.y
    end

    function vec2.mulf(v,f  ) 
        return getVec2(v.x*f , v.y*f)
    end

    function vec2.minus(v1 ,v2  ) 
        return getVec2(v1.x-v2.x , v1.y-v2.y)
    end

    --cross product v × w to be vx wy − vy wx
    function vec2.cross(v1,v2  ) 
        return v1.x*v2.y - v1.y*v2.x
    end
    function vec2.dot(v1,v2  ) 
        return v1.x*v2.x + v1.y*v2.y
    end
    function vec2.add(v1,v2  ) 
        return getVec2(v1.x+v2.x , v1.y+v2.y)
    end

    function vec2.normalize(v1)
        local d = v1:dist() 
        if d> 0 then
            return getVec2(v1.x/d,v1.y/d)
        else
            return getVec2(0,0)
        end
    end

    function vec2.normal(v1)
        return getVec2(-v1.y ,v1.x)
    end

    function vec2.left(v1)
        return getVec2(-v1.y ,v1.x)
    end

    function vec2.right(v1)
        return getVec2(v1.y ,-v1.x)
    end
    return vec2 
end
function findIntersect(l1p1x,l1p1y, l1p2x,l1p2y, l2p1x,l2p1y, l2p2x,l2p2y)
    local p = getVec2(l1p1x,l1p1y)
    local p1 = getVec2(l1p2x,l1p2y)
    local q = getVec2(l2p1x,l2p1y)
    local q1 = getVec2(l2p2x,l2p2y)

    return  vec2findIntersect(p,p1,q,q1)

end

function vec2findIntersect(p,p1,q,q1)
--https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect

    local r = vec2.minus(p1 ,p  ) 
    local s = vec2.minus(q1 ,q  ) 

    --t = (q − p) × s / (r × s)
    --u = (q − p) × r / (r × s)

    local q_p = vec2.minus(q,p)
    local q_pdots = vec2.cross(q_p,s)
    local q_pdotr = vec2.cross(q_p,r)
    local rdots = vec2.cross(r,s)

    --If r × s = 0
    if rdots == 0 then
        return false,false,nil
    end

    local t = q_pdots / rdots
    local u = q_pdotr / rdots

    --if 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1, the two line segments meet at the point p + t r = q + u s.
    if ((t>= 0 ) and (t<=1) and (u>= 0 ) and (u<= 1 )) then
        --p + t  *r 
        local ret = vec2.mulf(r,t) 
        ret = vec2.add(p,ret  ) 
        return ret.x,ret.y,ret
    end
    return false,false,nil
end 

Technically, you have a line expressed as a series of points. Then, you express the rectangle as bounded constraints. Then you would create sub-lines from the original lines by the bounded constraints. Depending on your own algorithm, this may have different effects, but this is the general idea.

You can use the Line-Line intersection formula, Line-line intersection.

If you consider a class Line class that is made up of two Vector2s you can calculate, using the above formula, the intersection between them.

class Line {
    public Vector2 a;
    public Vector2 b;

    public Line(Vector2 a, Vector2 b) {
        this.a = a;
        this.b = b;
    }

    Vector2 findIntersection(Line a, Line b)
        float x1 = a.a.x;
        float y1 = a.a.y;
        float x2 = a.b.x;
        float y2 = a.b.y;

        float x3 = b.a.x;
        float y3 = b.a.y;
        float x4 = b.b.x;
        float y4 = b.b.y;

        float denominator = (x1 - x2)*(y3 - y4) - (y1 - y2)*(x3 - x4);
        if (denominator == 0)
            return null;

        float xNominator = (x1*y2 - y1*x2)*(x3 - x4) - (x1 - x2)*(x3*y4 - y3*x4);
        float yNominator = (x1*y2 - y1*x2)*(y3 - y4) - (y1 - y2)*(x3*y4 - y3*x4);

        float px = xNominator / denominator;
        float py = yNominator / denominator;

        return new Vector2D(px, py);
    }
}

The above maps onto the function on wikipedia where you've got two points on a line (Given two points on each line).

Then it's a matter of turning the rectangle into Line which is of course done by taking the different corners as the points on the line and call the above method for each line and corner-to-corner line. Somewhat ugly that can be solved like this;

Vector2 findIntersection(Line a, Vector2 c0, Vector2 c1, Vector2 c2, Vector2 c3) {
    Line l0 = new Line(c0, c1);
    Line l1 = new Line(c1, c2);
    Line l2 = new Line(c2, c3);
    Line l3 = new Line(c3, c0);

    Vector2 intersection = null;

    intersection = findIntersection(a, l0);
    if (intersection != null) return intersection;

    intersection = findIntersection(a, l1);
    if (intersection != null) return intersection;

    intersection = findIntersection(a, l2);
    if (intersection != null) return intersection;

    intersection = findIntersection(a, l3);
    if (intersection != null) return intersection;

    return null;

} 

An important thing to point out here is that the formula will give you an intersection even if that intersection isn't between a and b, so you likely will want to constrain the results to results that fall within that range (one way to do that would be to do an bounds-inclusive contains on the source rectangle).

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.

That makes

$$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]$$

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 \$C_1\$ and \$C_2\$ (e.g. (100,0) and (0,0)):

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

and fill that in in the line equation:

$$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.

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
Latest version