Computational geometry, tetrahedron signed volume

Question

I'm not sure if this is the right place to ask, but here goes...

Short version: I'm trying to compute the orientation of a triangle on a plane, formed by the intersection of 3 edges, without explicitly computing the intersection points.

Long version: I need to triangulate a PSLG on a triangle in 3D. The vertices of the PSLG are defined by the intersections of line segments with the plane through the triangle, and are guaranteed to lie within the triangle. Assuming I had the intersection points, I could project to 2D and use a point-line-side (or triangle signed area) test to determine the orientation of a triangle between any 3 intersection points.

The problem is I can't explicitly compute the intersection points because of the floating-point error that accumulates when I find the line-plane intersection. To figure out if the line segments strike the triangle in the first place, I'm using some freely available robust geometric predicates, which give the sign of the volume of a tetrahedron, or equivalently which side of a plane a point lies on. I can determine if the line segment endpoints are on opposite sides of the plane through the triangle, then form tetrahedra between the line segment and each edge of the triangle to determine whether the intersection point lies within the triangle.

Since I can't explicitly compute the intersection points, I'm wondering if there is a way to express the same 2D orient calculation in 3D using only the original points. If there are 3 edges striking the triangle that gives me 9 points in total to play with. Assuming what I'm asking is even possible (using only the 3D orient tests), then I'm guessing that I'll need to form some subset of all the possible tetrahedra between those 9 points. I'm having difficultly even visualizing this, let alone distilling it into a formula or code. I can't even google this because I don't know what the industry standard terminology might be for this type of problem.

Any ideas how to proceed with this? Thanks. Perhaps I should ask MathOverflow as well...

EDIT: After reading some of the comments, one thing that occurs to me... Perhaps if I could fit non-overlapping tetrahedra between the 3 line segments, then the orientation of any one of those that crossed the plane would be the answer I'm looking for. Other than when the edges enclose a simple triangular prism, I'm not sure this sub-problem is solvable either.

EDIT: The requested image.

Answer 1

I would write symbolic vector equations, you know, with dot and cross products, to find the normal of the intersection triangle. Then, the sign of the dot product of this normal with the initial triangle one gives the orientation. So finally you can express this in a form sign(F(p1,...,p9)), where p1 to p9 are your points and F() is an ugly formula including dot and cross products of differences (pi-pj). Don't know if this can be done simpler, but this general approach does the job.

Answer 2

As I understand it, you have three lines intersecting the plane, and you want to calculate the orientation of the triangle formed by the intersection points, without calculating the intersection points themselves?

If so: you have a plane

N·(x - x0) = 0

and six points...

l1a, l1b, l2a, l2b, l3a, l3b

...forming three lines

l1 = l1a + t(l1b - l1a)
l2 = l2a + u(l2b - l2a)
l3 = l3a + v(l3b - l3a)

The intersection points of these lines to the plane occur at specific values of t, u, v, which I'll call t_i, u_i, v_i

N·(l1a + ti(l1b - l1a) - x0) = 0

      N·(x0 - l1a)
ti =  ----------------
      N·(l1b - l1a)
(similarly for ui, vi)

Then the specific points of intersection are

intersect1 = l1a + ti(l1b - l1a)
intersect2 = l2a + ui(l2b - l2a)
intersect3 = l3a + vi(l3b - l3a)

Finally, the orientation of your triangle is

orientation = direction of (intersect2 - intersect1)x(intersect3 - intersect1)

(x is cross-product) Work backwards plugging the values, and you'll have an equation for orientation based only on N, x₀, and your six points.

Answer 3

Let's call your triangle vertices T[0], T[1], T[2], and the first line segment's endpoints are L[0] and L[1], the second is L[2] and L[3], and the third is L[4] and L[5]. I imagine you want a function

int Orient(Pt3 T[3], Pt3 L[6]); // index L by L[2*i+j], i=0..2, j=0..1

which returns 1 if the intersections have the same orientation as the triangle, and -1 otherwise. The result should be symmetric under interchange of j values, antisymmetric under interchange of i values and T indices. As long as you can compute a quantity with these symmetries, that's all you need.

Let's try

Sign(Product( Orient3D(T[i],T[i+1],L[2*i+0],L[2*i+1]) * -Orient3D(T[i],T[i+1],L[2*i+1],L[2*i+0]) ), i=0..2))

where the product should be taken over cyclic permutations of the indices (modulo 3). I believe this has all the symmetry properties required. Orient3D is Shewchuk's 4-point plane orientation test, which I assume you're using.

Answer 4

I am answering this on both MO & SO, expanding the comments I made on MO.

My sense is that no computational trick with signed tetrahedra volumes will avoid the precision issues that are your main concern. This is because, if you have tightly twisted segments, the orientation of the triangle depends on the precise positioning of the cutting plane.
[image removed; see below]
In the above example, the upper plane crosses the segments in the order (a,b,c) [ccw from above]: (red,blue,green), while the lower plane crosses in the reverse order (c,b,a): (green,blue,red). The height of the cutting plane could be determined by your last bit of precision.

Consequently, I think it makes sense to just go ahead and compute the points of intersection in the cutting plane, using enough precision to make the computation exact. If your segment endpoints coordinates and plane coefficients have L bits of precision, then there is just a small constant-factor increase needed. Although I am not certain of precisely what that factor is, it is small--perhaps 4. You will not need e.g., L² bits, because the computation is solving linear equations. So there will not be an explosion in the precision required to compute this exactly.

Good luck!

(I was prevented from posting the clarifying image because I don't have the reputation. See the MO answer instead.)

Edit: Do see the MO answer, but here's the image: