slopeless.py

# Some slopeless formulations for working with line segments
# Author: @avonmoll

# There are two functions, each with two versions, one which computes the slope
# and one which circumvents direct calculation of the slope.

import numpy as np

def closestPointOnLineToPoint(start, end, point):
    """Find the point on the line that is closest to the given 
    point. The line is given as the x and y coordinates of 
    the endpoints (i.e. the line is actually a line segment).
    """
    # unpack points into their x- and y-components
    x1, y1 = start
    x2, y2 = end
    m = (y2 - y1) / (x2 - x1)  # slope, potentially numerically unstable
    
    # unpack point of interest
    xp, yp = point
    
    # closest point on the line to this point must form a line that is perpendicular
    # to the original line segment (start, end). Mathematically, this means that the
    # dot product is zero. Dot product of p = (p_x, p_y) and q = (q_x, q_y) is given 
    # as p_x * q_x + p_y * q_y.
    
    # Now let us paramaterize the coordinates of a point, P, on the line as:
    # x = x1 + t
    # y = y1 + m * t
    # Thus this t is like a change in x. The use of slope here ensures the point lies on
    # on the line. Now we want dot product of the line from this point to our point of 
    # interest and the line from start to end to be zero:
    # (P - point) dot (end - start)
    # (x1 + t - x1) * (x2 - x1) + (y1 + m * t - y1) * (y2 - y1) = 0
    # Solve for t:
    t = -(x1 + y1 * m - xp - yp * m) / (1 + m**2)
    
    # Now compute x and y coordinate of closest point using our parameterization
    xt = x1 + t
    yt = y1 + m * t
    pt = np.array([xt, yt])
    
    # Some fancy footwork is needed to move the point back onto the line segment if it
    # lies before or beyond the start or end, respectively...

    return pt 


def closestPointOnLineToPoint_slopeless(start, end, point):
    """Slopeless formulation of the above
    """
    
    # unpack each 2D point into x- and y- components
    x1, y1 = start
    x2, y2 = end
    x3, y3 = point
    
    # Now let us paramaterize the line segment as start + t * (end - start) where t is
    # any scalar value. Note this exression is a vector expression. We want the line
    # segment from this point to our point of interest to be perpendicular to the segment
    # (end - start):
    # start + t *(end - start) - point _|_ (end - start)
    # which implies dot product is zero:
    # (x1 + t * (x2 - x1) - xp) * (x2 - x1) + (y1 + t * (y2 - y1) - yp) * (y2 - y1) = 0
    # We have one equation, and one unknown: t. Solving for t gives:
    t = (-x2*x3 + x2*x1 + x1*x3 - x1**2 - y2*y3 + y2*y1 + y1*y3 - y1**2) / \
        (-(x2 - x1)**2 - (y2 - y1)**2)
    # The above is "numerically stable" in the sense that the denominaor is only 
    # zero iff start and end are the same point.
        
    # Now, if t > 1 then the closest point lies past the end of the line segment so we
    # should push it back to 1. Similarly, if t < 0 then closest point lies before the 
    # start of the line segment so we should push it up to 0 (to be actually ON the 
    # line segment).
    t = min(max(0, t), 1)
    
    return start + (end - start) * t


def intersectionOfTwoLines(line1, line2):
    """Find intersection of two lines which are each given as the start and end point
    of a segment of the line.
    """

    # unpack lines into coordinates of each point
    [x1, y1], [x2, y2] = line1
    [x3, y3], [x4, y4] = line2

    # First compute slope
    m12 = (y2 - y1) / (x2 - x1)  # slope of line1
    m34 = (y4 - y3) / (x4 - x3)  # slope of line2, either could be infinite

    # Intersection satisfies both of these lines' point-slope formulae:
    # (y - y1) = m12 * (x - x1)
    # (y - y3) = m34 * (x - x3)
    # Two equations, two unknowns, solve for x and y:
    x = (y3 - y1 + m12 * x1 - m23 * x3) / (m12 - m23)
    y = y1 + (x - x1) * m12

    return np.array([x, y])

def intersectionOfTwoLines_slopeless(line1, line2):
    """Slopeless formulation of the above
    """

    # unpack lines into coordinates of each point
    [x1, y1], [x2, y2] = line1
    [x3, y3], [x4, y4] = line2

    # Now we employ a parameterization of a point on the first line as p + t * (q - p)
    # (as before) where p and q are the endpoints ofthe segment and t is a scalar number.
    # This point must also lie on the second line, which we parameterize as:
    # r + u * (s - r) where r and s are the endpoints and u is a scalar number which, 
    # like t, is saying "how far along this line segment is the point?". So now we
    # have two equations for the intercept point. The coordinates of the intercept point
    # are the two unknowns. We can solve algebraically, but we will end up with 
    # the same issue as before: numerical instabality when y1 = y2 or some such
    # condition. However, we can solve as a system of equations using linear 
    # algebra. To do so, we need to invert a 2x2 matrix. Sparing the details here...

    idet = 1/((x2 - x1) * (y3 - y4) - (y2 - y1) * (x3 - x4))
    iA = np.array([[y3 - y4, x4 - x3],
                   [y1 - y2, x2 - x1]])
    b = np.array([[x3 - x1],
                  [y3 - y1]])
    t, u = idet * iA.dot(b)

    # Now we know t AND u, but we only really need either one to compute the intercept
    # point:
    x_intercept = x1 + t * (x2 - x1)  # or x3 + u * (x4 - x3)
    y_intercept = y1 + t * (y2 - y1)  # or y3 + u * (y4 - y3)

    return np.array([x_intercept, y_intercept])