Value of determinant from a predicate

[DanielVandH]

Is it possible to get both the sign of a predicate, and the value? e.g. for

$$ \text{orient}(u, v, w) = \det \begin{bmatrix} u_1-w_1 & u_2-w_2 \ v_1-w_1 & v_2-w_2 \end{bmatrix}, $$

is it possible to get both $\text{sgn}(\text{orient}(u, v, w))$ and a reliable estimate for the value?

--

In case this might be a bit misguided, here's the context of what I actually want. Let $O(u, v, w)$ denote the orient predicate. In implementing a constrained Delaunay triangulation, a vector needs to be built of the form $d_i = O(v_0, vi, v{m-1})$, $i=0,\ldots,m-1$, and certain vertices are selected based on whether $di < d{i-1}$ and $di < d{i+1}$. One way to do this rigorously would be to just define a predicate for the sign of $di - d{i-1}$ and $di - d{i+1}$, but this might be very expensive to do many times.

[DanielVandH]

I have one possible idea here https://github.com/DanielVandH/ExactPredicates.jl/blob/withval/src/Codegen.jl for returning the computed values (assuming such an approach doesn't already exist that I just haven't found).

Judging from the code, though, it looks like different filters aren't necessarily looking at comparable values? Is there a way to compare them? e.g. in my code,

julia> using ExactPredicates
[ Info: Precompiling ExactPredicates [429591f6-91af-11e9-00e2-59fbe8cec110]

julia> p, q, r = (0.3, 0.391), (0.1, 0.991), (-5.0, 1.0)
((0.3, 0.391), (0.1, 0.991), (-5.0, 1.0))

julia> [
       ExactPredicates.orient_naive_val_dbg(p, q, r)
       ExactPredicates.orient_reference_val_dbg(p, q, r)
       ExactPredicates.orient_slow_val_dbg(p, q, r)
       ExactPredicates.orient_val_dbg(p, q, r)
       ]
4-element Vector{Tuple{Int64, ExactPredicates.Codegen.FilterEnum, Real}}:
 (1, ExactPredicates.Codegen.naive_flt, 3.0582)
 (1, ExactPredicates.Codegen.exact_flt, 496221320399100289854781446127249//162259276829213363391578010288128)
 (1, ExactPredicates.Codegen.interval_flt, [3.05819, 3.05821])
 (1, ExactPredicates.Codegen.fastfp_flt, 0x0000000000000400)

How can I ensure that these values all refer to orient(p, q, r) (which in this example is about 3.0582)?

[lairez]

The problem with “reliable estimates” is that they are not exact, so there is a chance that you cannot determine the sign of a difference of two reliable estimates. So output the computed value, not just the sign, is only one part of a solution.

For your specific case, I have a much simpler solution. As you know, $d_i = O(a, v_i, b)$ is the exterior product $(a-b)\wedge (v_i - b)$. So, by linearity, we have

$$ d_{i+1} - di = (a-b)\wedge(v{i+1} - v{i}) = O(b-a, v{i+1}-v_i)$$

Of course, you cannot directly call orient(b-a, v[i+1]-v[i]) because there will be round-off errors in the subtractions, that would not be “robust”. Good news: this sign predicate is already implemented as parallelorder(b, a, v[i], v[i+1]).

Honestly, I would just use that. I am afraid that the efforts you will put in a more complex solution might not pay off.

[DanielVandH]

That all makes sense - indeed the subtraction is a big problem here.

Thanks for your solution! Never used parallelorder before, but that's a very nice solution. I will use that.

Do you think it might still be of value to have #10? Obviously I won't use it anymore for my approach, but maybe there's still some use - don't really know any application of it currently though. If not, we can just close this issue / #10.

[lairez]

Do you think it might still be of value to have https://github.com/lairez/ExactPredicates.jl/pull/10?

It is not clear, sorry! There are several issues: the value output is not very well specified mathematically. Secondly, the functions will be type-unstable, so I doubt we can do efficient code with them.

[DanielVandH]

Hi again @lairez , I just wanted to follow up on something similar to this. I think the answer again is it can't be done within the scope of this package, but I was hoping you might have some insight. The context is segment splitting in mesh refinement.

The problem is to compute a point on a line $pq$, say $r = p + t(q-p)$ for $t \in (0, 1)$, and it is important that this computation be such that orient(p, q, r) = 0. It's not so uncommon that this is not the case, e.g. orient(p, q, r) could be like 1e-16. Do you know if there might be something in this package that I should reuse to help improve this computation, or maybe have any suggestions for handling this? "Improve" meaning perhaps make it much more frequent that orient(p, q, r) = 0 is true, I'm not sure if it's possible all the time.

One idea would be to slightly nudge r in the direction of the perpendicular bisector if orient(p, q, r) != 0, e.g.

r = p + t(q - p)
o = orient(p, q, r)
o == 0 && return nothing # no filtering
A2 = squared_triangle_area(p, q, r) # == 2orient(p, q, r)
L2 = (px - qx)^2 + (py - qy)^2
A2 <= 0.0 && return nothing 
a = sqrt(4A2/L2)
a = o == -1 ? a : -a 
rx += a * (qy - py)
ry += a * (px - qx)

but unfortunately this doesn't work as the computation of the area (= twice the orient) is not robust.

Thanks again.

[lairez]

You have a big problem here: the point r that you are looking for does not always exist. For example, if you consider the line segment between p = (0,1) and q = (1, 1.0000000000000002), then there is no point r, with Float64 coordinates that lies on the line segment defined by p and q (expect p and q themselves of course). So don't waste time trying to find an algorithm for computing r, it is a problem without a solution.