I'm trying to find the intersection points of two triangles using the Tomas Moeller algorithm (https://web.stanford.edu/class/cs277/resources/papers/Moller1997b.pdf)
the intersection of plane1 and plane2 is a line, L = O + tD, where D = N1 x N2 is the direction of the line and $O$ is some point on it.
https://i.gyazo.com/89e717283d0a8784f9bc2dc3cf1cfebe.png
The paper says the intersection line is equal to $L = O + tD$ but I am unsure how to obtain this point. I have the direction of the line which is the cyan line from one of the vertex positions in the image, but it is not in the correct position. I am assuming this is because I am missing $O$ and it does not say how to calculate it.
The paper does mention that $O$ is not needed towards the end, but I am unsure how to get the correct line position as at the moment, my line direction is correct but is a normalised direction.
How do I calculate $O$ given that all other data is known, i.e. angles, vertices, distances
Edit:
Vector3 D = Vector3.Cross(N2, N1);
Vector3 a = Vector3.Dot(N2, triangle_wv0) * triangle_wv1;
Vector3 b = Vector3.Dot(N1, triangle_wv1) * triangle_wv0;
float c = Vector3.Dot(N2, triangle_wv0);
float d = Vector3.Dot(N1, triangle_wv1);
Vector3 O = (a - b) / (c - d);
$intersection: (21.8, 18.2, 34.5)$
$N1 = (0.0, 0.0, -1.0)$
$N2 = (0.0, 1.0, 0.0)$
$D = (-1.0, 0.0, 0.0)$
$triangle_wv0 = (21.8, 13.2, 31.0)$
$triangle_wv1 = (21.8, 25.4, 39.7)$
$triangle_wv2 = (36.8, 13.2, 31.0)$
$cutter_wv0 = (8.0, 20.0, 5.0)$
$cutter_wv1 = (8.0, 20.0, 80.0)$
$cutter_wv2 = (83.0, 20.0, 5.0)$
Edit 2: I think I found the problem and I now have a closer point to the line, but it doesn't seem to work in all cases
At the moment I am only drawing the red line to point $O$ to see if the point is on the intersection line and that does not seem to be correct.
https://i.gyazo.com/7238288c18f38457d13702da6357c691.png - Good, point is on intersection line
https://i.gyazo.com/29b478a873f2af1497f7e513942f2d37.png - Bad, point is not on intersection line
The calculated intersection for the triangles defined above is not correct
$intersection: (21.8, 18.2, 34.5)$
$O$ can be any point on the line of intersection of the planes. There are various ways to find one, which you can find with an Internet search. Given what you have to work with, though, there’s another direct computation that you can perform that’s not among those in the usual sources.
A general approach to finding a point on the intersection line is to find a line in one of the planes that’s guaranteed to intersect it and then compute the intersection of this second line with the other plane. A reasonably efficient way to do this is to use the Plücker matrix of the line: working in homogeneous coordinates, the line through points $\mathbf p$ and $\mathbf q$ can be represented by the matrix $\mathcal L = \mathbf p\mathbf q^T-\mathbf q\mathbf p^T$. (This is the 3-D analogue of representing a line through two points in 2-D by their cross product.) The intersection of this line with a plane $\mathbf\pi$ is simply $$\mathcal L\mathbf\pi = (\mathbf p\mathbf q^T-\mathbf q\mathbf p^T)\mathbf\pi = (\mathbf q^T\mathbf\pi)\mathbf p-(\mathbf p^T\mathbf\pi)\mathbf q.\tag1$$
In the course of applying Moeller’s algorithm you will have found a pair of triangle vertices $V_i$ and $V_j$ that lie on opposite sides of the intersection line. These are the points $\mathbf p$ and $\mathbf q$ (with a $1$ appended to make them homogeneous, of course). For $\mathbf\pi$ you’ll have the vector $[N^T;d]^T$, i.e., the concatenation of the parameters of the other triangle’s plane. Plugging these values into formula (1) gives $$\begin{align} ([V_j;1]\cdot[N;d])[V_i;1] - ([V_i;1]\cdot[N;d])[V_j;1] &= (N\cdot V_j+d)[V_i;1] - (N\cdot V_i+d)[V_j;1] \\ &= [(N\cdot V_j)V_i-(N\cdot V_i)V_j+d(V_i-V_j); N\cdot V_j-N\cdot V_i], \end{align}$$ which in inhomogeneous Cartesian coordinates is $$O = {(N\cdot V_j)V_i-(N\cdot V_i)V_j+d(V_i-V_j) \over N\cdot V_j-N\cdot V_i}.\tag2$$ This method isn’t quite as efficient as a related method that uses a different line and Cramer’s rule, but this way doesn’t require first finding the component of $D$ with the greatest absolute value and then picking out the other coordinates of $N_1$ and $N_2$.
Example: Using the vertices of the two triangles in the question, we compute $$\mathbf\pi_{tri} = [N_{tri};d_{tri}] = [0. : 0.580607 : -0.814184 : 17.5757] \\ \mathbf\pi_{cut} = [N_{cut};d_{cut}] = [0. : 1. : 0. : -20.],$$ so $$D = N_{tri}\times N_{cut} = (0.814184, 0., 0.)$$ normalized to $(1,0,0)$. Since we’re using the first triangle to compute $O$, we use the values from $\mathbf\pi_{cut}$ in formula (2): $$N\cdot V_j = (0., 1., 0.)\cdot(21.8, 25.4, 39.7) = 25.4 \\ N\cdot V_i = (0., 1., 0.)\cdot(21.8, 13.2, 31.0) = 13.2$$ and finally $$O = {(25.4-20.0)\,(21.8, 13.2, 31.0)-(13.2-20.0)\,(21.8, 25.4, 39.7) \over 25.4-13.2} = (21.8,20.0,35.8492).$$ (I’ve rolled the $d(V_i-V_j)$ term back into the first two terms in this computation.)