How to compare points in homogeneous coordinates universally? Problem
Given a plane, I have lines in the canonical form (described by 3D vectors). Using a cross product, I can find an intersection of any two lines. I want to check, whether two pairs of lines intersect at the same point, taking noisiness of the data into account. This means that I should compare them with respect to a chosen threshold
$$
p_1 = \ell_1 \times \ell_2 \\
p_2 = \ell_3 \times \ell_4 \\
\left( l_1 \nparallel l_2 \right) \land \left( l_3 \nparallel l_4 \right) \implies
\left( p_1^z \neq 0 \right) \land \left( p_2^z \neq 0 \right) \implies
\left\lVert \frac{p_1}{p_1^z} - \frac{p_2}{p_2^z} \right\rVert < \varepsilon?
$$
If two given lines are not parallel, their intersection point has a non-zero third coordinate. So, I can divide by it to get ordinary 2D coordinates of the point and compare it with another one in the usual way.
Otherwise, if the lines are parallel, the third coordinate of their intersection point is close to zero (meaning it's the point at infinity), so I cannot divide by it and should compare the directions of two points at infinity
$$
p_1 = \ell_1 \times \ell_2 \\
p_2 = \ell_3 \times \ell_4 \\
\left( l_1 \parallel l_2 \right) \land \left( l_3 \parallel l_4 \right) \implies
\left( p_1^z \approx 0 \approx p_2^z \right) \implies
\left\lVert
\frac{p_1}{\left\lVert p_1 \right\rVert} - \frac{p_2}{\left\lVert p_2 \right\rVert}
\right\rVert < \varepsilon?
$$
Question
Is there a universal way to compare points in a homogeneous coordinate system to avoid the separate check of the third coordinate?
 A: You can "normalize" as you've done: replace $p_1$ by $q_1 = \frac{p_1}{\|p_1\|}$, and similarly for $p_2$ to get a pair of unit vectors in the 2-sphere in 3-space. There's a small problem here, which is that even if $p_1$ and $p_2$ represent the same point, $q_1$ and $q_2$ might point in opposite directions (i.e., $p_1 = -q_2$).
So the next step is to compute $d = q_1 \cdot q_2$, and if $d < 0$, you negate $q_2$. (I'm going to assume you've done this, and now call the replaced vector $q_2$.)
Now you've got two vectors in the same hemisphere of $S^2$, so to measure how far apart they are, you can compute the angle between them:
$$
\theta = \arccos(q_1 \cdot q_2).
$$
(Noting that $\cos$ is an even function, you'll realize that the whole "negate $q_2$ isn't actually necessary...)
This value $\theta$ gives you a "distance" between the computed intersections.
Is $\theta$ invariant with respect to the various transformations that matter to you? You might, for instance, say "If I changed coordinates by translating everything by $(a,b)$ for any part of numbers $(a,b)$, then the 'distance' should be unchanged". If you want that, then you have to check whether $\theta$ has that property. Since I don't know what properties you wish to have, I can't really go any farther.
