I have to figure out how to get the path of points that make up a geodesic on a triangular mesh:
If we know the position and initial direction of a 2D walker restricted to the surface of the mesh, how do we find the path it will take?
How do we find the geodesic along the surface of the mesh, from one known point to another known point?
For example,
shows one such path in red.
We know which triangles are touching each other.
I need to figure out the following:
The point at which the hypothetical red line intersects the edge of the triangle.
Assuming it intersects the triangle, what is the new direction based on this path?
Let's say the vertices of the current triangle are $\vec{v}_1 = ( x_1 , y_1 , z_1 )$, $\vec{v}_2 = ( x_2 , y_2 , z_2 )$, and $\vec{v}_3 = ( x_3 , y_3 , z_3 )$, the point in the triangle is $\vec{p} = ( x_p , y_p , z_p )$, and the direction to move from point $\vec{p}$ is $\vec{d} = ( x_d , y_d , z_d )$.
Point $\vec{p}$ must be within the triangle plane, $$\left ( \vec{p} - \vec{v}_i \right ) \cdot \vec{n} = 0$$ where $i$ is $1$, $2$, or $3$ -- if it is true for one, it is true for all three --, and $\vec{n}$ is the normal of the triangle plane, $$\vec{n} = \left ( \vec{v}_2 - \vec{v}_1 \right ) \times \left ( \vec{v}_3 - \vec{v}_1 \right )$$ Similarly, $\vec{d}$ must be within the triangle plane too, $$\vec{d} \cdot \vec{n} = 0$$
The line starting at $\vec{p}$ in direction $\vec{d}$ will intersect one or more of the edges, $$\vec{p} + r_{ij} \vec{d} = \vec{v}_i + t_{ij} \left ( \vec{v}_j - \vec{v}_i \right )$$ where $i,j$ is $1,2$ or $2,3$ or $3,1$. For each $i,j$ pair we have three equations, $$\begin{cases} x_p + r_{ij} x_d = x_i + t_{ij} ( x_j - x_i ) \\ y_p + r_{ij} y_d = y_i + t_{ij} ( y_j - y_i ) \\ z_p + r_{ij} z_d = z_i + t_{ij} ( z_j - z_i ) \end{cases}$$ You can use any pair above to solve $r_{ij}$ and $t_{ij}$. (I'd use the pair with the largest $\lvert c_d \rvert \lvert c_j - c_i \rvert$.) The edge that is intersected first, is the edge between vertices $i$ and $j$ where $r_{ij}$ reaches its smallest positive value.
Let's assume that the edge is shared with another triangle. If we rotate the other triangle around the shared edge, so that it becomes planar with the current triangle, we can continue the line from the intersection point in direction $\vec{d}$. Rotating the triangle and the direction back to the original orientation, the direction becomes $\vec{q}$.
If the surface normal for the current triangle is $\vec{n}_1$ and the new triangle $\vec{n}_2$, then the rotation axis unit vector $\hat{a}$ and angle $\varphi$ fulfill $$\begin{array}{c} \hat{a} = \frac{ \vec{n}_1 \times \vec{n}_2 }{ \lVert \vec{n}_1 \times \vec{n}_2 \rVert} \\ \sin\varphi = \lVert \left ( \frac{\vec{n}_1}{\lVert\vec{n}_1\rVert} \right ) \times \left ( \frac{\vec{n}_2}{\lVert\vec{n}_2\rVert} \right ) \lVert \\ \cos\varphi = \left ( \frac{\vec{n}_1}{\lVert\vec{n}_1\rVert} \right ) \cdot \left ( \frac{\vec{n}_2}{\lVert\vec{n}_2\rVert} \right ) \end{array}$$ You can use Rodrigues' rotation formula to rotate $\vec{d}$: $$\vec{q} = \vec{d} \cos\varphi + \left ( \hat{a} \times \vec{d} \right ) \sin\varphi + \hat{a} \left ( \hat{a} \cdot \vec{d} \right ) ( 1 - \cos\varphi )$$
Another option is to use planar coordinates, say $(u, v)$, within each triangle. (Because these coordinates are specific to each triangle, and even specific to how the vertices are labeled, I call these triangle coordinates.)
Origin is at triangle vertex $\vec{v}_j$, and the unit vector $\hat{e}_{ij}$ is $$\hat{e}_{ij} = \frac{ \vec{v}_j - \vec{v}_i }{ \lVert \vec{v}_j - \vec{v}_i \rVert }$$ If we consider a pair of triangles sharing the edge $\vec{v}_i - \vec{v}_j$, only the $v$ axis ($\hat{e}_{k}$) differs for the two triangles. In the first triangle, it is $$\hat{e}_k = \frac{ \vec{v}_k - \vec{v}_i - \hat{e}_{ij} \left ( \hat{e}_{ij} \cdot ( \vec{v}_k - \vec{v}_i ) \right ) }{\lVert \vec{v}_k - \vec{v}_i - \hat{e}_{ij} \left ( \hat{e}_{ij} \cdot ( \vec{v}_k - \vec{v}_i ) \right ) \rVert }$$ It is computed the exact same way for the second triangle, too, except that the $\vec{v}_k$ is the third vertex for the second triangle.
Since $\hat{e}_{ij}$ and $\hat{e}_k$ are unit vectors, $\lVert\hat{e}_{ij}\rVert = 1$ (and $\hat{e}_{ij} \cdot \hat{e}_{ij} = 1$), and $\lVert\hat{e}_{k}\rVert = 1$ (and $\hat{e}_k \cdot \hat{e}_k = 1$).
Both $\hat{e}_k$'s are perpendicular to $\hat{e}_{ij}$. If the two triangles are coplanar, and we use $\hat{e}_{k , 1}$ for $\hat{e}_k$ in the first triangle, and $\hat{e}_{k , 2}$ for $\hat{e}_k$ in the second triangle, then that means that $$\begin{cases} \hat{e}_{k , 1} \cdot \hat{e}_{ij} = 0 \\ \hat{e}_{k , 2} \cdot \hat{e}_{ij} = 0 \end{cases}$$
The key observation is this:
If the two triangles were coplanar, $\hat{e}_{k,2} = -\hat{e}_{k,1}$
If a 3D direction vector in the plane of the two triangles corresponds to $( u , v )$ in the first triangle, then $(u , -v )$ in the second triangle corresponds to the exact same 3D direction
If the two triangles are not coplanar, then direction $(u, v)$ in the first triangle corresponds to direction $(u, -v)$ in the second triangle in the "geodesic sense" (that is, if the two triangles were coplanar, the directions would be the same).
Any 3D point $\vec{p}$ on the triangle plane can be described using triangle coordinates $(u, v)$: $$\begin{cases} u = \left ( \vec{p} - \vec{v}_i \right ) \cdot \hat{e}_{ij} \\ v = \left ( \vec{p} - \vec{v}_i \right ) \cdot \hat{e}_{k} \end{cases} \iff \vec{p} = \vec{v}_i + u \, \hat{e}_{ij} + v \, \hat{e}_{k}$$ For the direction vector $\vec{d}$, we use $$\begin{cases} u = \vec{d} \cdot \hat{e}_{ij} \\ v = \vec{d} \cdot \hat{e}_{k} \end{cases} \iff \vec{d} = u \, \hat{e}_{ij} + v \, \hat{e}_{k}$$
Note that if $\vec{p}$ is on the plane, then $$\vec{p} \cdot \left ( \hat{e}_{ij} \times \hat{e}_k \right ) = 0$$ because the triangle normal $\vec{n}$ is parallel to $\hat{e}_ij \times \hat{e}_k$.
If you have solved $t_{ij}$ for the shared edge intersection point using the previous method, the intersection point is at $$\begin{cases} u = \frac{t_{ij}}{\lVert \vec{v}_j - \vec{v}_i \rVert} \\ v = 0 \end{cases}$$ by definition: $t_{ij} = 0$ if it is at $\vec{v}_i$, $1$ if at $\vec{v}_j$, with $t_{ij}$ linear with respect to location.
It is possible to solve $t_{ij}$ in the triangle coordinates directly. Essentially, you calculate the point and the direction in three orientations ($(i,j,k)$ is $(1,2,3)$, $(2,3,1)$, or $(3,1,2)$), and pick the one that yields the smallest $r \ge 0$. In each orientation, you calculate the $(u_0 , v_0)$ corresponding to the starting point, and $(u_\Delta , v_\Delta)$ corresponding to the direction, and if $v_\Delta \lt 0$, $$r = -\frac{v_0}{v_\Delta}$$
Note that if $v_\Delta \ge 0$ for some orientation $i, j, k$, or if $$u_0 + r u_\Delta \lt 0$$ or if $$u_0 + r u_\Delta \gt \lVert \vec{v}_j - \vec{v}_i \rVert$$ then this orientation is not valid. At least one orientation will be valid for a non-degenerate triangle (triangle with area greater than zero).
Choosing the $i, j, k$ that yields the smallest positive valid $r$ basically chooses the orientation where the chosen direction will intersect edge $i,j$ first.
The intersection in the chosen orientation $i, j, k$ will occur at coordinates $$\left ( u_0 + r u_\Delta , 0 \right )$$ and the "geodesically same direction" in the new triangle will be $$\left ( u_\Delta , -v_\Delta \right )$$