Why do we need to homogenise these vectors?

Question

I was looking for a solution for finding the intersection point of two lines. I am aware that this can be done by finding their vector product.

I stumbled upon this example here:

Numpy and line intersections

def get_intersect(a1, a2, b1, b2):
    s = np.vstack([a1,a2,b1,b2])        # s for stacked
    h = np.hstack((s, np.ones((4, 1)))) # h for homogeneous
    l1 = np.cross(h[0], h[1])           # get first line
    l2 = np.cross(h[2], h[3])           # get second line
    x, y, z = np.cross(l1, l2)          # point of intersection
    if z == 0:                          # lines are parallel
        return (float('inf'), float('inf'))
    return (x/z, y/z)

I have gone through the example and used it in a few scenarios and it seems to work pretty well. However, there are three things I don't quite get:

1. Why does the vector need to be homogeneous (the part where we fill a column with ones)?
2. How does the homogeneous solution differ compared to a non-homogeneous solution (if at all)?
3. How come we only check the result for parallelism along the Z-axis and not X and Y as well?

I feel like I am missing something very obvious but I can't wrap my head around on what it is.

Answer 1

As a slight correction to the answers above, using 3D cross products to compute a 2D line intersection is pretty much the least efficient possible. Cross products do not SIMD optimise at all (unless you go to the lengths of using an SOA data structure, but even than it's an extremely wasteful approach). The best approach is to use good old simultaneous equations.

Starting with our input points that define the line:

line1: [(x1, y1), (x2, y2)]
line2: [(x3, y3), (x4, y4)]

Compute some direction vectors:

// 1st direction
u1 = x2 - x1
v1 = y2 - y1
D1 = [u1, v1]

// 2nd direction
u2 = x4 - x3
v2 = y4 - y3
D2 = [u2, v2]

Now lets reformulate the line equations as rays, and come up with an equation for any point along those rays:

// coords of a point 'd1' distance along the 1st ray
Px = x1 + d1*u1
Py = y1 + d1*v1

// coords of a point 'd2' distance along the 2nd ray
Px = x3 + d2*u2
Py = y3 + d2*v2

Let's assume the lines intersect, which means both P's will be the same, allowing us to state:

x1 + d1*u1 = x3 + d2*u2
y1 + d1*v1 = y3 + d2*v2

I wont go through each step, but re-arrange both equations in terms of d1, and we get:

 d1 = x3 + d2*u2 - x1
      ---------------
            u1

 d1 = y3 + d2*v2 - y1
      ---------------
            v1

Now we have two equations for d1, so lets do another simultaneous equation to get a value for d2:

x3 + d2*u2 - x1     y3 + d2*v2 - y1
---------------  =  ---------------
      u1                  v1

Rearrange to isolate d2:

 d2 = u1(y3 - y1) - v1(x3 - x1)
      -------------------------
           v1*u2 - u1*v2

If (v1*u2 - u1*v2) happens to be zero, there is no solution for this equation (let's call that the determinant, because it is what it is!). If the determinant is non zero, then simply compute d2 using the above equation, and pump back into one of our earlier equations to find the point value:

Px = x3 + d2*u2
Py = y3 + d2*v2 

Some untested C++ code:

bool computeIntersectPoint(
  float x1, float y1,
  float x2, float y2,
  float x3, float y3,
  float x4, float y4,
  float outPoint[2])
{
  // compute direction vectors
  // in some cases, it can be worth 
  // storing the lines as rays as an 
  // optimisation. (avoids 4 subs)
  const float u1 = x2 - x1;
  const float v1 = y2 - x1;
  const float u2 = x4 - x3;
  const float v2 = y4 - x3;

  // check to see if we have a solution
  // 1 mul, 1 fmsub
  float determinant = v1*u2 - u1*v1; 
  if(determinant == 0)
    return false;

  // 2 sub, 1 mul, 1 fmsub
  float temp = u1 * (y3 - y1) - v1 * (x3 - x1);

  // 1 div
  float intersectDistance = temp / determinant;

  // 2 fma
  outPoint[0] = intersectDistance * u2 + x3;
  outPoint[1] = intersectDistance * v2 + y3;

  return true;
}

Quick desmos proof: https://www.desmos.com/calculator/gtlmnmzn6l

It's worth at this point comparing the instruction count between the two approaches here. A 3d cross product requires 3 mult, and 3 fmsub instructions (or 6 mul + 3 sub if FMA is not available). Since we have 3 of those, we are up to: 9 mul and 9 fmsub ops. Add in the 2 divisions, and we get:

9 mul
9 fmsub
2 div

The approach I've posted would require:

1 div
6 sub
4 fma
2 mul

Although you could save 4 of those subs if you can get away with storing the lines as rays.

Answer 2

For reference, the equation from Wikipedia:

Let a1 = (x1, y1), a2 = (x2, y2), b1 = (x3, y3), b2 = (x4, y4).

---

Computational interpretation

Observe the first two cross-products in the linked answer:

l1 = np.cross(h[0], h[1]) = (x1, y1, 1) ^ (x2, y2, 1)
   = (y1 - y2, x2 - x1, x1*y2 - x2*y1)

l2 = np.cross(h[2], h[3]) = (x3, y3, 1) ^ (x4, y4, 1)
   = (y3 - y4, x4 - x3, x3*y4 - x4*y3)

These two lines are all it takes to calculate the 6 different terms in the equation above. And the last cross-product:

x, y, z = np.cross(l1, l2)
      x = (x2 - x1) * (x3*y4 - x4*y3) - (x4 - x3) * (x1*y2 - x2*y1)
-->   y = (y3 - y4) * (x1*y2 - x2*y1) - (y1 - x2) * (x3*y4 - x4*y3)
      z = (y1 - y2) * (x4 - y3) - (y3 - y4) * (x2 - x1)

These numbers are exactly equal to the numerators and denominator in Wikipedia's equation.

A fairly complex expression like this would take dozens of FPU instructions to compute term-by-term.

Homogenizing the vectors allows for this cross-product method, which can be optimized to just a handful of SIMD instructions – much more efficient.

---

Geometrical interpretation

Suppose you treat the homogenized vectors as points in 3D space, and join each pair with the origin to make two triangles:

All 4 points lie on the plane Z = 1 (gray).

The line L (green) is the intersection between the planes of the two triangles (blue + red), and passes through the origin O and the desired point of intersection P (yellow).

The normal of a triangle is given by the cross-product of its side vectors. In this case the side vectors are simply given by the 4 points, as the other point is the origin. In the code the normals are given by l1 and l2.

One definition of a plane is that all lines which lie in it must be perpendicular to its normal. Since the line L lies in the planes of both triangles, it must be perpendicular to l1 and l2, i.e. its direction is given by np.cross(l1, l2).

Homogenization allows for a clever final step which uses similar triangles to compute P:

if z == 0:                          # lines are parallel
    return (float('inf'), float('inf'))
return (x/z, y/z)                   # Px = x / z, Py = y / z