Hp 50g Math question: Intersection of 2 triangulars

[peacecalc]

Hello all,

it's my first post here:

1) It is a good idea to make a new kategory, labeled "mathematics"? In my opinion: yes.

2) My problem six points are given in 3D and three of this points define the trinangular ONE and the other three the triangular TWO. I want for my hp 50g an input of six vectors and as an output an empty list (no intersection set), a list with two vectors (if the two vectors have the same value, then the intersecting set is a point, if they are different then the set is a piece of a line), if the list contains more then two vectors the intersecting set is a aera.

That is the plan, my brutal force solution looks for me clumsy and fanciless. And it is a solution which is far from a "fast" solution.

Does anybody know more about that kind of problems?

sincerely
peacecalc

P.S.: Edit: corrected headline

[David Hayden]

I'm doing this off the top of my head, so it will probably have problems.

Let's call the triangles A and B. First, see if any of lines in A intersects any of the lines in B. If they do then the triangles intersect. Otherwise they might still intersect, but it would be the special case where one triangle in enclosed within the other.

If they do intersect, the intersection could have between 3 and 6 sides I believe.

[Marcel Samek]

(12-20-2013 04:37 PM)peacecalc Wrote:  Hello all,

it's my first post here:

1) It is a good idea to make a new kategory, labeled "mathematics"? In my opinion: yes.

2) My problem six points are given in 3D and three of this points define the trinangular ONE and the other three the triangular TWO. I want for my hp 50g an input of six vectors and as an output an empty list (no intersection set), a list with two vectors (if the two vectors have the same value, then the intersecting set is a point, if they are different then the set is a piece of a line), if the list contains more then two vectors the intersecting set is a aera.

That is the plan, my brutal force solution looks for me clumsy and fanciless. And it is a solution which is far from a "fast" solution.

Does anybody know more about that kind of problems?

sincerely
peacecalc

There are several cases:

1) They do not intersect
2) They are not coplanar and intersect on a line
3) They are coplanar and intersect to create a plane
4) They are coplanar and intersect on an edge to create a line segment

For each triangle, you can determine the normal by taking the cross product of the two vectors defined by the three points. If the normals are identical, that means that the planes are parallel. If they are parallel and coplanar, they do not intersect; if they are parallel and non-coplanar, they might intersect to create an area. If the normals are not identical, then they might intersect on a line.

For each of the possible intersection cases you will probably need to do a separate algorithm.

For the non-coplanar test, simply find the line that is defined by the intersection of the two planes defined by the triangles. Then test to see if that line intersects each triangle. The segment of the line that intersects both triangles defines the intersection of the two triangles.

The coplanar test is the most complex because the intersection area is polygon whose vertices are defined by either the intersection of the triangle edges, or by the vertices of one triangle that fall inside the other. You need to construct that polygon. Poke around the internet and you will find a variety of approaches.

All that said, it is a tricky problem because it is sensitive to numerical precision if the triangles are of odd shapes, disparate sizes, far apart, small angles, etc.

[Han]

The equation of a line is:
\[ a \cdot (x-x_0) + b\cdot (y-y_0) = 0 \]
and of a plane is:
\[ a \cdot (x-x_0) + b\cdot (y-y_0) + c\cdot (z-z_0) = 0 \]
and generalizes in the natural way to higher dimension. What is both interesting and useful is that it is basically the dot product of a normal (i.e. perpendicular) vector \( \langle a,b \rangle \) with the vector \( \langle x-x_0, y-y_0 \rangle \), and the dot product is 0 since the vectors are perpendicular. The vector \( \langle a,b \rangle \) can be easily obtained from a simple slope calculation through 2 points. In 3D, you would need to create 2 vectors from a set of 3 points, and take a cross product to come up with a vector normal to the plane containing the three points. Here's the 2D version:

Since we now have a dot-product interpretation for a linear equation, it's quite simple to determine if two triangles intersect, and even how they intersect. Let the six points be
\[ P_1:=(x_1, y_1), \quad P_2:=(x_2, y_2), \quad \dotsm,
\quad P_6:=(x_6, y_6) \]
and let \(L_{i,j}\) be the line through points \( (x_i,y_i)\) and \((x_j,y_j)\) with \(i<j\). Define an associated function:
\[ L_{i,j}(X,Y) := -(y_j-y_i) \cdot (X-x_i) + (x_j-x_i) \cdot (Y-y_i) \]
Without loss of generality, let the first triangle be the one created by the lines \( L_{1,2}, L_{2,3}, L_{1,3} \). Then for each line, we can define the "in" side of that line to be the sign of \(L_{i,j}\) evaluated at the remaining point. For example, \( L_{1,2}(P_3) \) will either be negative or positive, or zero. If zero, then \(P_3\) lies on \(L_{1,2} \) -- not possible if we have a valid triangle. If negative, then \(P_3\) is on the opposite side of the normal vector, and positive if on the same side. Apply this to all three lines and it's easy to determine the "in"-side of a that triangle (based purely on the signs of the valuations of each line). So, to continue, if \( L_{i,j}(P_4) \) has the same corresponding signs as for \(P_1\), \(P_2\), and \(P_3\), then \(P_4\) lies on the "in"-side of that triangle. If \( L_{i,j}(P_4) = 0 \) then \(P_4\) lies on the edge defined by \(L_{i,j}\)... and hopefully you get the idea.

If all three points \(P_4, P_5, P_6\) all do not lie on the "outside" then either their vertices are in the interior or possibly boundary of the first triangle. The area would be a simple determinant calculation. The only remaining cases are one or two points are on the "out"-side. Since inside and outside are all determined by the valuations of a line \(L_{i,j}\), it's pretty easy to find the intersection points as you clearly know which lines are involved (one is \(L_{i,j}\) and the other consists of one point inside the triangle and one outside the triangle).

[Han]

The 3D case is more involved, but not that much different in principle from the 2D case. If no one posts a solution for the 3D case by tomorrow afternoon (or one different from mine), then I will post mine to share.

[peacecalc]

Wow,

I thought this a minor problem, but your answers show that it worth to think a little bit about that topic.

My solution will work in some instances:

1) testing wether the two planes which are defined by the 2 triangulars are truly parallel and don't intersect each other. That can be easly done by the cross product (don't forget I use the hp 50g) and using the Hesse-form for planes:

\[ \vec{n}_1 = \vec{a} \times \vec{b} \] first triangular (trig) \[ \vec{a},\vec{b} \] are defining the two sides of trig 1

\[ \vec{n}_2 = \vec{c} \times \vec{d} \] second trig \[\vec{c},\vec{d}\] are defining the two sides of trig 2

\[ (\vec{x} - \vec{t}_i )\cdot \vec{n}_i = 0 \]

When the two planes are not intersecting the procedure ends, otherwise: (the planes could be identical
or intersect with a line)

2) Now I decide to calculate the the radii-vectors (of the circumcircles) of the two trigs. If they obey the inequality:

\[ \| \vec{r}_2 - \vec{r}_1 \| > \| \vec{r}_1 \| + \| \vec{r}_2 \| \]

then the trigs cannot intersect and the procedure ends (I do this testing because I have two hulls build
with a lot of triangles and a lot of them are not intersecting. Let's say hull 1 has 15 and hull 2 has 20 trigs,
I have to test 300 trigs, that is a challenge for my hp 50g), otherwise

3) If the plane intersect with a line, I consider the the vectors \[\vec{c},\vec{d}\] and \[\vec{d}-\vec{c}\]
as the direction-vectors of three line equalities in analytical geometry:

i. e.: \[ g_1 ~ : ~ \vec{x} = \vec{t} + \gamma\vec{c} \] and so on. Vector \[ \vec{t} \] is the
position-vector of trig 2. With the Hesse form I can find out where are the intersection points with plane 1.
Every point can be expressed as a vector beginning on the position-vector of trig 1. And now this vector can be
expressed as a linear-combination of the vectors \[ \vec{a},\vec{b} \] (which describe the trig 1) \[ \vec{s}=\alpha\cdot\vec{a}+\beta\cdot\vec{b}\] and with the two parameters \[ \alpha, \beta\] it can be decided wether that intersecting point is inside from trig 1 or not.


4) If the plane 1 is identical to plane 2 (We know this from 1) we have the similar work to do as in 3). But we
intersect then lines with lines and not lines with a plane, like in 3).

Statisticaly I expect more trigs, which are not intersecting, but I've the
impression part 3) or 4) is also a lot of work for my calculator. Maybe somebody else knows
a more sharp heuristic argument (as in 2) formulated) for accelerating the program.
My Program has to work a lot with the cross and dot product of the hp 50g.


@Han, a special thank for you, you spend a lot of time formulating your answer and solution. I'm curious about your 3-D solution.

Sincerely
peacecalc

[Han]

For the 3D case, you have a lot of the right calculations already. I will use \( \vec{x}_i \) to mean both a vector as well as a point.

To test whether or not the second triangle is on the same plane as the first triangle, we only only need to do a very straightforward set of computations. If \(H: \vec{n} \cdot (\vec{x}-\vec{x}_1) = 0 \) is the equation of the plane \(H\) containing triangle one (where \( \vec{x}_1 \) is a vertex of triangle one), then we simply need to check whether
\[ \vec{n} \cdot (\vec{x}_i -\vec{x}_1) =0 \]
for \(i = 4, 5, 6 \). Within this check we can also deduce whether the second triangle even intersects the first. If the dot product is 0 for \(i = 4, 5, 6 \) then the second triangle is coplanar with the first, so we need only check for the position of \(\vec{x}_i\) relative to triangle one for \(i=4,5,6\). More on this later. If the dot product is non-zero in all three cases, and all of the same sign, then we know the two triangles do not intersect (as triangle two lies completely inside one of the half-spaces separated by \(H\)). That only leaves the case that two vertices of triangle two are one side of \(H\) while the third vertex is on the opposite side of \(H\). Without loss of generality, suppose the points \( \vec{x}_4 \) and \( \vec{x}_5 \) are on opposite sides of \( H \), and \( \vec{x}_6 \) is on the same side as \( \vec{x}_5 \). We can find a point of intersection between \( H \) and the line \( L \) through \( \vec{x}_4 \) and \( \vec{x}_5 \). This line \( L \) has the form
\[ L : t \cdot (\vec{x}_5 - \vec{x}_4) + \vec{x}_4, \quad \text{ where } \quad t\in [0,1]\subseteq \mathbb{R} \]
Substitute this into the equation for \( H \) to obtain an equation of only one variable: \( t \). Solve for \(t \) and we have the desired intersection point. If necessary, find the other intersection point using \( \vec{x}_4 \) and \( \vec{x}_6 \).

At this stage, we have basically reduced all cases down to \( H \) containing triangle one, and points that also lie on \( H\) -- either intersection points or actual vertices of triangle two. The most difficult step is to determine the position of a point the the plane \( H \) relative to triangle one.

To determine their positions relative to triangle one, you can use a linear system. Since we already have the three vertices from triangle one being coplanar, any other point in \( H \) must be a linear combination of the vertices of triangle one. The solution can even be normalized. That is, if \( \mathbf{A} \) is the matrix whose entries are the coordinates of the vertices (in column form) of triangle one, and \( \vec{b} \) is a point on \( H \), then the solution \(\vec{x} = \vec{\lambda} \) to \( \mathbf{A}\cdot \vec{x} = \vec{b} \) can also be made to satisfy
\[ \sum_{j=1}^3 \lambda_j = 1 \]
where \( \lambda_j \) is the \(j\)-th component of \( \vec{\lambda} \). In a sense, \( \lambda_j \) is a weight on the \(j\)-th vertex of triangle one and is a measure of how close \( \vec{b} \) is to that vertex. Moreover, the sign of \( \lambda_j \) determines the relative position of \( \vec{b} \) to vertex \( j \). (Look up Radon partitions; there are also lots of theorems regarding \(d+2\) points in a \(d\)-dimensional space.) That is, if \(\lambda_j\) is positive, then \( \vec{b} \) and vertex \( j \) are on the same side of the line passing through the remaining two vertices of triangle one. Thus, if \(\lambda_j \) are all positive, then \( \vec{b} \) lies in the relative interior of triangle one. If \( \lambda_j = 0 \) for a specific value of \( j \) then \( \vec{b} \) actually lies one an edge of triangle one.

This should give you more than enough to figure out the rest. As a friendly reminder, area of a triangle in \( \mathbb{R}^3 \) is merely half the magnitude of the cross product of the two vectors representing two sides of that triangle. So once you know the position of \( \vec{x}_4\), \(\vec{x}_5\), and \( \vec{x}_6 \) relative to triangle one, just partition into triangles and compute the area accordingly.

Edit: The worst case is of course when all six points are coplanar. In this case, you will actually have to repeat the last bit of linear algebra and switch the roles of triangle one and triangle two to reduce the complexity. It very well could be that triangle one is inside triangle two! However, the steps above would detect that by simply swapping triangle one out for triangle two.

[peacecalc]

Hello Han,

thank you for lucid answer!

Quote: The most difficult step is to determine the position of a point the the plane H relative to triangle one.

To determine their positions relative to triangle one, you can use a linear system. Since we already have the three vertices from triangle one being coplanar, any other point in H must be a linear combination of the vertices of triangle one. The solution can even be normalized. That is, if A is the matrix whose entries are the coordinates of the vertices (in column form) of triangle one, and b⃗ is a point on H, then the solution x⃗ =λ⃗ to A⋅x⃗ =b⃗ can also be made to satisfy
∑j=13λj=1

where λj is the j-th component of λ⃗ . In a sense, λj is a weight on the j-th vertex of triangle one and is a measure of how close b⃗ is to that vertex. Moreover, the sign of λj determines the relative position of b⃗ to vertex j. (Look up Radon partitions; there are also lots of theorems regarding d+2 points in a d-dimensional space.) That is, if λj is positive, then b⃗ and vertex j are on the same side of the line passing through the remaining two vertices of triangle one. Thus, if λj are all positive, then b⃗ lies in the relative interior of triangle one. If λj=0 for a specific value of j then b⃗ actually lies one an edge of triangle one.

That is maybe the solution to avoid a lot of "IF THEN" structures, which often slow down the execution speed of programs.

[SlideRule]

peacecalc

An interesting read ...
Triangle Triangle Intersection
...attached.

BEST!

SlideRule

[Han]

(12-22-2013 05:21 PM)SlideRule Wrote:  peacecalc

An interesting read ...
Triangle Triangle Intersection
...attached.

BEST!

SlideRule

Very curious if there is any similarity between the paper and what I had proposed. Too bad they didn't use \( \LaTeX\ \). Their notation borders on atrocious ;_;

[peacecalc]

Hello Han

the notation from this article is really bad. In a first overview they represent the problem in an such ugly way, so that is no wonder that they need a lot of pages for solving this.

When my program is ready, I will publish it here in that forum. The mathematics the problem needs is secondary II level (in germany, maybe highschool level in USA) in schools.
I never thought to work on university level. Although your solution Han, included a generalization for more dimensions as three.

[peacecalc]

Hello all,

The great program for solving this problem isn't ready. But the frame programe is ready and some sub programs which I want to present.
Some remarks: I use the decimal comma and not a point and the command \->NUM can be deleted if it is shure that all values are reals (otherwise the if conditions will not work correct). The programs are all "hard coded", I don't use the solve abilities of the HP 50g in hope that my solutions are faster and I can format the output in a way I want. If your are interested in the mathematics behind these programs, let me know.

The first one: "NONUL":

Its function is to find a non-zero value in a 3-dim vector.

Input: stack 1: a vector
Output: stack 2: value
stack 1: index

This program is used in program "CUTLL" avoiding a possible division by zero..

Code:


%%HP: T(3)A(R)F(,);
\<< OBJ\-> DROP DUP
  IF 0, \=/
       THEN 3,
       ELSE UNROT DUP
              IF 0, \=/
                       THEN 2,
            ELSE UNROT DUP
                  IF 0, \=/
                       THEN 1,
                       ELSE 0,
                       END
            END
       END 4 ROLLD 4 ROLLD DROP2
\>>

If the vector is the null-vector index is zero (and the value, too).

The second is called: "CUTPL"

Its function is to find a interection between a line and a plane.

Input: stack 4: position vector of the line
stack 3: direction vector of the line
stack 2: position vector of the plane
stack 1: normal vector of the plane
Output: stack 1: { } if there is no intersection
stack 1: { [position vector] [direction vector] } if the
line lies in the plane
stack 1: single value (if that value is inserted in the line equation
you get the position vector of the intersecting point.)
Because this value is importent later on, this program stops
here with the calculation.

Code:


%%HP: T(3)A(R)F(,);
\<< \-> TG R TE N
  \<< TG TE - N DOT NEG R N DOT DUP \->NUM
    IF 0, ==
         THEN SWAP \->NUM
                 IF 0, ==
              THEN TG R 2, \->LIST
              ELSE { }  END 
           SWAP DROP
         ELSE /  END
  \>>
\>>

The third one is "CUTLL".

Its function is to find a interection between a line and a line.

Input: stack 4: position vector of the line 1
stack 3: direction vector of the line 1
stack 2: position vector of the line 2
stack 1: direction vector of the line 2
Output: stack 1: { } if there is no intersection
stack 1: { [position vector 1] [direction vector 1] } if the
lines are identical
stack 1: single value (if that value is inserted in the line equation
of line 1 you get the position vector of the intersecting point.)
Because this value is importent later on, this program stops
here with the calculation.

Code:


%%HP: T(3)A(R)F(,);
\<< \-> T1 R1 T2 R2
  \<< T2 T1 - R2 NEG CROSS 
         DUP R1 R2 NEG CROSS 
         DUPDUP ABS \->NUM
    IF 0, ==
    THEN ROT ABS \->NUM
             IF 0, ==
             THEN 3, DROPN T1 R1 2, \->LIST
             ELSE 3, DROPN { } END
    ELSE ROT CROSS ABS \->NUM
             IF 0, \=/
             THEN DROP2 { }
             ELSE NONUL ROT SWAP GET SWAP / END
    END
  \>>
\>>

Sincerely peacecalc