Why do we need to homogenise these vectors?

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.

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