If two lines are defined in homogeneous coordinates, we can find the intersection point as cross product of two lines.
Check
Let's assume we have $l_1= [a,b,c]^T$ and $l_2 = [d,e,f]^T$, then traditionally to find the point that belong to the line, we need to solve the system of linear equations. $$ \begin{cases} ax + by + c = 0 \\ dx + ey + f = 0 \\ \end{cases} $$ taken $x$ from first equation results in $x = \frac{-c-by}{a}$ then substituting x into second equation, we will have $$ \begin{array} {lcl} \frac{d}{a}(-c-by)+ey+f = 0\\ -\frac{d}{a} - \frac{db}{a}y + ey +f = 0 \\ (e-\frac{db}{a}y = \frac{dc}{a}- f) \\ y = \frac{dc-fa}{ea-db} \end{array} $$
Substituting $y$ back to find $x$, we need to perform a sequence of simplifications: $$ \begin{array} {lcl} x &=& -\frac{c}{a}- \frac{b}{a}(\frac{dc-fa}{ea-db}) \\ &=& \frac{1}{a}\left(\frac{-c(ae-bd)-b(cd-af)}{ae-bd}\right) \\ &=& \frac{1}{a(ae-bd)}\left(-cae+cbd-bcd+baf\right) \\ &=& \frac{-cae+baf)}{a(ae-bd)} = \frac{a(bf-ec)}{a(ae-bd)} = \frac{bf-ec}{ae-bd} \end{array} $$
So solving the system of linear equations, we have: $$ \begin{cases} x= \frac{bf-ec}{ae-bd} \\ y = \frac{dc-fa}{ea-db} \end{cases} $$
Showing that by decomposing the cross product we will get the same results $$ \begin{array} {lcl} l_1 \times l_2 &=& \begin{bmatrix} i & j & k \\ a & b & c \\ d & e & f \end{bmatrix} = i \begin{bmatrix}b & c \\ e &f \end{bmatrix} - j \begin{bmatrix}a & c \\ d & f \end{bmatrix} + k \begin{bmatrix}a & b \\ d & e \end{bmatrix} \\ &=& (bf-ec)i - j(af-cd) + k(ae-bd) \\ &=& (bf-ec)i + (cd-af)j + k(ae-bd) \end{array} $$ Rewriting it a little bit leads to $$ \begin{bmatrix}i\\ j\\ k \end{bmatrix} = \begin{bmatrix} bf-ec \\ cd-af \\ ae-bd \end{bmatrix} = \begin{bmatrix} \frac{bf-ec}{ae-bd} \\ \frac{cd-af}{ae-bd} \\ 1 \end{bmatrix} = \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} $$
Checked
## Checking this with l1:= y=3 and l2 := y=x
import numpy as np
## data
A = np.matrix([[0,1], [-1,1]])
b = np.array([3,0])
l1 = np.array([0,1,-3])
l2 = np.array([-1,1,0])
## solving
sol_sys = np.linalg.solve(A,b)
print("Solution sys: ", sol_sys)
sol_cross = np.cross(l1,l2)
sol_cross = sol_cross /sol_cross[2]
print("Solution cross: ", sol_cross)
import matplotlib.pyplot as plt
fig=plt.figure(1)
x = np.arange(-2, 6, 0.1)
n = x.shape[0]
y1 = [3]*n ## generates vector of n elements all has value 3
y2 = x
plt.plot(x, y1, 'r')
plt.plot(x, y2, 'b');
plt.grid(linestyle='--', linewidth=1);
Solution sys: [3. 3.] Solution cross: [3. 3. 1.]
## Checking with parallel lines y = 3 and y = -1
## data
A = np.matrix([[-1,1], [-1,1]])
b = np.array([3,-1])
l1 = np.array([0,1,-3])
l2 = np.array([0,1,1])
## solving
# sol_sys = np.linalg.solve(A,b) ---> crashes as expected
# print("Solution sys: ", sol_sys)
sol_cross = np.cross(l1,l2) ## --> solution exists, but point to the infinity
# sol_cross = sol_cross /sol_cross[2] --> impossible, division by zero
print("Solution cross: ", sol_cross)
Solution cross: [4 0 0]
p1 = np.array([3, 1, 1])
p2 = np.array([-4, 5, 1])
l = np.cross(p1,p2)
print(l)
[-4 -7 19]