# robustness of geometric predicates in Euclidean vs homogeneous coordinates

The signed volume of the triangle formed by the points $$p, q, r$$ in the plane is defined to be $$\text{volume}(p, q, r) \equiv \det\left[\begin{matrix}q_1 - p_1 & r_1 - p_1 \\ q_2 - p_2 & r_2 - p_2\end{matrix}\right].\tag{1} \label{eq:1}$$ Whether the signed volume is positive, negative, or zero tells you if the input points are oriented in a way that obeys the right-hand rule, the left-hand rule, or if they're coplanar. The signed volume is a key part of point-in-polygon testing and in convex hull algorithms. A mathematically equivalent definition of the signed volume predicate is $$\text{volume}(p, q, r) = \det\left[\begin{matrix} 1 & 1 & 1 \\ p_1 & q_1 & r_1 \\ p_2 & q_2 & r_2\end{matrix}\right].\tag{2} \label{eq:2}$$ The two definitions are equivalent by performing a sequence of column operations on the matrix. You can think of the last equation as using the homogeneous or projective coordinates of the input points.

Computing the signed volume naively in floating-point arithmetic is a bad idea because underflow can give you the wrong answer for nearly co-linear points. This paper by Lutz Kettner and others gives some examples; in particular, taking

$$p = \left(\begin{matrix}0.5 \\ 0.5\end{matrix}\right), \quad q = \left(\begin{matrix}12 \\ 12\end{matrix}\right), \quad r = \left(\begin{matrix} 24 \\ 24\end{matrix}\right) \tag{3} \label{eq:3}$$

and perturbing $$p$$ by increments of $$2^{-53}$$ in each coordinate can give positively-oriented points that are wrongly classified as co-linear or negatively-oriented, and vice versa. There are remedies for these robustness problems but that's not what I'm interested in.

When I compute the signed volume of the points in equation \eqref{eq:3} using the Cartesian form \eqref{eq:1} and perturb $$p$$ in a grid by increments of an ulp, I get the expected floating-point failures and lack of robustness shown in the Kettner et al. paper. But when I use the homogeneous form \eqref{eq:2}, I don't -- all the signs are correct! What is a triple $$p, q, r$$ of co-linear points that can make even the homogeneous form fail? A recommendation for a search strategy would be welcome if an explicit example isn't obvious.

• Floating-point arithmetic is not associative, i.e. the order of operations matters, and this can make computation more or less numerically robust. Are you assuming a particular way of computing the 3x3 determinant, e.g. as a dot-product from the determinants of 2x2 sub-matrices, using the matrix concatenation "shortcut", or LU decomposition? Likewise, are you allowing or disallowing the use of FMA (fused multiply-add)? Off the cuff, it seems to me that unless the determinant can be computed as a correctly rounded result (max err < 0.5 ulp) one can make any of the alternatives I mentioned fail. Dec 21, 2022 at 8:28
• When I use a simple direct method of computing the 3x3 determinant I can provoke failures (+ vs -) easily by fuzzing $p_x$, $p_y$ by $c$ ulp, $c \in [-128,127]$. My actual code (KISS64 is my PRNG): m[1][0]= uint64_as_double(double_as_uint64(m[0][2])+((KISS64%256)-128)); m[2][0]= uint64_as_double(double_as_uint64(m[0][1])+((KISS64%256)-128)); With a more robust implementation of the 3x3 determinant I don't see such failures. Dec 21, 2022 at 10:42

I suggest trying the same search procedure with collinear points that either:

a. don't have their two coordinates equal, or

b. all have their coordinates strictly smaller than 1 in modulus.

I suspect that what fixes the problem is the fact that pivoting in the LU factorization (which is used to compute the determinant) always puts the largest value first, and this somewhat makes things better. With (a) the largest value does not depend on perturbations, so the pivoting permutation is fixed, and with (b) the first step of Gaussian elimination recovers your formula (1), so the two formulas should have the same behavior.

I am assuming that your linear algebra library uses LU to compute a 3x3 determinant; if it uses something different such as a hardcoded formula then this won't apply. Knowing which algorithm is used is useful.

• I'm using Eigen and, from reading the source code, it looks like they use a hand-coded formula for determinants up to 4x4 and LU with partial pivoting for larger matrices. But you make a good point about LU with partial pivoting acting as a kind of stabilization, I hadn't thought of that. Dec 21, 2022 at 17:47

I wrote a fuzzer as per njuffa's suggestion in the comments. In each trial, I generated two points with integer coordinates between -64 and +64 and then made a third point along the line that connects them. Here's one triple that gives interesting results: $$x_0 = \left(\begin{matrix}34 \\ -18\end{matrix}\right), \quad x_1 = \left(\begin{matrix}-61 \\ -43\end{matrix}\right), \quad x_2 = \left(\begin{matrix}-156 \\ -68\end{matrix}\right)$$ The top row of the figure below shows the results of evaluating $$\text{signed area}(x_0 + \delta x, x_1, x_2)$$ in the first column, $$\text{signed area}(x_1, x_2, x_0 + \delta x)$$ in the second, and $$\text{signed area}(x_2, x_0 + \delta x, x_1)$$ in the third for small values of $$\delta x$$. Mathematically, these are all supposed to be equal, but in floating-point arithmetic they are not because of e.g. loss of associativity. The remaining rows of the plot are generated from points in this JSON file. There's some more code and complete commentary in this notebook.

I don't have a compelling answer for why points with identical coordinates should work ok by computing the signed area in homogeneous coordinates but poorly in Cartesian coordinates.

Shameless plug: the notebook that I used to generate these examples is a demonstration for the library predicates which I wrote for doing robust geometric predicates.