/ Functional Programming

Solving a Problem: Segment Intersection

This post was written with purpose to demonstrate how the programming is done, from problem statement to working code. Also to learn MathJax and refresh knowledge of F# and analytic geometry.

Problem: Create function that checks whether two line segments intersect. The segments are defined by their ends { (x11, y11), (x12, y12) } and { (x21, y21), (x22, y22) }.

Solution:

Let's find point where two segment intersect by first finding a point where lines that contain these segments intersect and then check whether this point belongs to both these segments' bounded rectangles.

The equation

$$A x + B y + C = 0 \tag{1}$$

where A, B, and C are constant and A and B are not both zero, defines a line in two-dimensional space. If C != 0, the equation can be rewritten as

$$a x + b y + 1 = 0 \tag{2}$$

If it is known that \((x_1, y_1)\) and \((x_2, y_2)\) are two points on the line (2), the a and b can be obtained by solving the following system of equations:

$$a x_1 + b y_1 + 1 = 0\\a x_2 + b y_2 + 1 = 0 \tag{3}$$

get a from the first and put it to the second, to get b:

$$b = \frac{x_2 - x_1} {x_1 y_2 - x_2 y_1} \tag{4}$$

and do the same for a:

$$a = - \frac{ y_2 - y_1 } { x_1 y_2 - x_2 y_1 } \tag{5}$$

From (2), (4), (5):

$$( y_2 - y_1 ) x - ( x_2 - x_1 ) y - ( x_1 y_2 - x_2 y_1 ) = 0 \tag{6}$$

So the A, B, and C of the line equation (1) for line going through two points (x1, y1) and (x2, y2) can be calculated as follows:

$$A = y_2 - y_1, B = - (x_2 - x_1), C = - (x_1 y_2 - x_2 y_1) \tag{7}$$

Let's now assume that there are two lines defined by their equations:

$$a_1 x + b_1 y + c_1 = 0\\
a_2 x + b_2 y + c_2 = 0$$

The same method can be used to find intersection point coordinates:

$$x = - \frac{b_2 c_1 - b_1 c_2}{a_1 b_2 - a_2 b_1}\\
y = \frac{a_2 c_1 - a_1 c_2}{a_1 b_2 - a_2 b_1} \tag{8}$$

If \(a_1 b_2 - a_2 b_1 = 0\), then there is either no solutions if lines are parallel or there are infinite number of solution because the lines are equivalent.

The segments intersect when the lines intersection point is within the boundaries of their bounded rectangles.

The algorithm for finding intersection point of two segments is:

  1. Get coefficients for linear equation for segment 1
  2. Get coefficients for linear equation for segment 2
  3. If they are not parallel, get the intersection point
  4. Check if the intersection point is within the boundaries of segments

There are several commands that can be derived from this steps:

  • lineEquation - returns coefficients a, b, c for two points (x1, y1) and (x2, y2)
  • intersectLines - returns intersection point for a1, b1, c1, a2, b2, c2 or no value if there is no such point
  • segmentBounds - returns x, y, width, height for x1, y1, x2, y2
  • contains - returns whether point x, y is within area specified by location and size

Although it is already possible to start implementing the functions, thinking ahead and define a few data types can make the result cleaner.

The algorithm above can be coded into function as follows (in F#):

type point = double * double
type area = (double * double) * (double * double)
type line = double * double * double
type segment = point * point

let intersectSegments segment1 segment2 =
    let line1 = lineEquation segment1
    let line2 = lineEquation segment2
    match intersectLines line1 line2
    | None -> None
    | Some point ->
        let ok1 = segmentBounds segment1
                  |> contains p
        let ok2 = segmentBounds segment2
                  |> contains p
        match ok1, ok2
        | true, true -> Some(point)
        | _ -> None

Then code the implementation of the functions lineEquation, intersectLines, segmentBounds, and contains.

let lineEquation segment =
    let (x1, y1), (x2, y2) = segment
    (y2 - y1, - (x2 - x1), - (x1 * y2 - x2 * y1))

let intersectLines line1 line2 =
    let (a1, b1, c1), (a2, b2, c2) = line1, line2
    let d = a1 * b2 - a2 * b1
    match Math.Abs(double(d)) < 1e-10 with
    | true -> None
    | false ->
        let x = - (b2 * c1 - b1 * c2) / d
        let y = (a2 * c1 - a1 * c2) / d
        Some (x, y)

let segmentBounds segment =
    let (x1, y1), (x2, y2) = segment
    ((x1, y1), (x2 - x1, y2 - y1))

let contains point area =
    let location, size = area
    let x, y = location
    let width, height = size
    let px, py = point
    px > x && px < (x + width) && py > y && py < (y + height)

The first iteration of the code is not inteded to be correct, it is only to build confidence that the data structures are right and the data flows are clear.

What tests the code is um... tests!

let check value message =
    if !value then
        printfn "%s: FAIL" message

let lineEquationTest1 () =
    let a, b, c = lineEquation ((0.0, 0.0), (1.0, 1.0))
    check (Math.Abs(a - 1.0) < 1e-10) "a equals 1"
    check (Math.Abs(b - (-1.0)) < 1e-10) "b equals 1"
    check (Math.Abs(c) < 1e-10) "c equals 0"

let intersectLinesTest1 () =
    let p = intersectLines (1.0, -1.0, 0.0) (1.0, 1.0, 0.0)
    match p with
    | None -> printfn "Intersection check: FAIL"
    | Some (x, y) -> check (Math.Abs(x) < 1e-10) "x equals 0"
                     check (Math.Abs(y) < 1e-10) "y equals 0"

let segmentBoundsTest1 () =
    let (x, y), (w, h) = segmentBounds ((1.0, 1.0), (0.0, 0.0))
    check (Math.Abs(x) < 1e-10) "x equals 0"
    check (Math.Abs(y) < 1e-10) "y equals 0"
    check (Math.Abs(w - 1.0) < 1e-10) "w equals 1"
    check (Math.Abs(h - 1.0) < 1e-10) "h equals 1"

Tests reveal that segmentBounds constructs the bounding rectangle incorrectly. The corrected implementation is:

let segmentBounds segment =
    let (x1, y1), (x2, y2) = segment
    let x = Math.Min(double(x1), double(x2))
    let y = Math.Min(double(y1), double(y2))
    let w = Math.Max(double(x1), double(x2)) - x
    let h = Math.Max(double(y1), double(y2)) - y
    ((x, y), (w, h))

Finally, the intersectSegments returns correct intersection point:

    let s1 = ((10.0, 0.0), (0.0, 10.0))
    let s2 = ((9.0, 0.0), (0.0, 11.0))
    let p = intersectSegments s1 s2
    printfn "%A" p

prints "Some (4.5, 5.5)".

Alex Netkachov

Alex Netkachov

Alex Netkachov is a Senior Software Developer, currently working in Central London on new generation of energy trading solutions for brokers, traders and exchanges.

Read More

Why not to stay updated if the subject is interesting? Join Telegram channel Alex@Net or follow alex_at_net on Twitter. Or just, use the comments form below.