1
$\begingroup$

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.

$\endgroup$
2
  • 2
    $\begingroup$ 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. $\endgroup$
    – njuffa
    Dec 21, 2022 at 8:28
  • $\begingroup$ 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. $\endgroup$
    – njuffa
    Dec 21, 2022 at 10:42

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$ Dec 21, 2022 at 17:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.