I realized that a solution using barycentric coordinates for walking the triangular mesh would be much simpler.
Let's assume we have vertices
$$\begin{array}{l}\vec{v}_1 = ( x_1 , y_1 , z_1 ) \\
\vec{v}_2 = ( x_2 , y_2 , z_2 ) \\
\vec{v}_3 = ( x_3 , y_3 , z_3 ) \end{array}$$
defining a triangle.
We can define any position within the triangle using barycentric coordinates $(u , v)$, with
$$\begin{cases}
0 \le u \le 1 \\
0 \le v \le 1 \\
0 \le u + v \le 1 \end{cases}$$
so that point $\vec{p}$ on the triangle is
$$\vec{p} = (x , y , z ) = \vec{v}_1 + u \left ( \vec{v}_2 - \vec{v}_1 \right ) + v \left ( \vec{v}_3 - \vec{v}_1 \right )$$
There are three, mathematically equivalent ways to compute $(u, v)$, because we have three equations (one for each dimension in $\vec{p}$), but only two unknowns. Numerically, the one with the largest divisor in magnitude should yield most stable result:
$$\begin{cases}
u = \frac{ x ( y_2 - y_0 ) + x_0 ( y - y_2 ) + x_2 ( y_0 - y )
}{ x_0 ( y_1 - y_2 ) + x_1 ( y_2 - y_0 ) + x_2 ( y_0 - y_1 ) } \\
v = \frac{ x ( y_0 - y_1 ) + x_0 ( y_1 - y ) + x_1 ( y - y_0 )
}{ x_0 ( y_1 - y_2 ) + x_1 ( y_2 - y_0 ) + x_2 ( y_0 - y_1 ) }
\end{cases}$$
$$\iff \begin{cases}
u = \frac{ x ( z_2 - z_0 ) + x_0 ( z - z_2 ) + x_2 ( z_0 - z )
}{ x_0 ( z_1 - z_2 ) + x_1 ( z_2 - z_0 ) + x_2 ( z_0 - z_1 ) } \\
v = \frac{ x ( z_0 - z_1 ) + x_0 ( z_1 - z ) + x_1 ( z - z_0 )
}{ x_0 ( z_1 - z_2 ) + x_1 ( z_2 - z_0 ) + x_2 ( z_0 - z_1 ) } \\
\end{cases}$$
$$\iff \begin{cases}
u = \frac{ y ( z_2 - z_0 ) + y_0 ( z - z_2 ) + y_2 ( z_0 - z )
}{ y0 ( z_1 - z_2 ) + y_1 ( z_2 - z_0 ) + y_2 ( z_0 - z_1 ) } \\
v = \frac{ y ( z_0 - z_1 ) + y_0 ( z_1 - z ) + y1 ( z - z_0 )
}{ y0 ( z_1 - z_2 ) + y_1 ( z_2 - z_0 ) + y_2 ( z_0 - z_1 ) }
\end{cases}$$
Similarly, we can describe any direction on the plane using $(du, dv)$:
$$\vec{d} = ( dx, dy, dz ) = du \left ( \vec{v}_2 - \vec{v}_1 \right )
+ dv \left ( \vec{v}_3 - \vec{v}_1 \right )$$
Again, there are three mathematically equivalent solutions for $(du, dv)$ (assuming $\vec{d}$ is parallel to the plane, i.e. $\vec{v}_1 + \vec{d}$ is on the plane of the triangle):
$$\begin{cases}
du = \frac{ dx ( y_2 - y_0 ) + dy ( x_0 - x_2 )
}{ x_0 ( y_1 - y_2 ) + x_1 ( y_2 - y_0 ) + x_2 ( y_0 - y_1 ) } \\
dv = \frac{ dx ( y_0 - y_1 ) + dy ( x_1 - x_0 )
}{ x_0 ( y_1 - y_2 ) + x_1 ( y_2 - y_0 ) + x_2 ( y_0 - y_1 ) } \\
\end{cases}$$
$$\iff \begin{cases}
du = \frac{ dx ( z_2 - z_0 ) + dz ( x_0 - x_2 )
}{ x_0 ( z_1 - z_2 ) + x_1 ( z_2 - z_0 ) + x_2 ( z_0 - z_1 ) } \\
dv = \frac{ dx ( z_0 - z_1 ) + dz ( x_1 - x_0 )
}{ x_0 ( z_1 - z_2 ) + x_1 ( z_2 - z_0 ) + x_2 ( z_0 - z_1 ) } \\
\end{cases}$$
$$\iff \begin{cases}
du = \frac{ dy ( z_2 - z_0 ) + dz ( y_0 - y_2 )
}{ y0 ( z_1 - z_2 ) + y_1 ( z_2 - z_0 ) + y_2 ( z_0 - z_1 ) } \\
dv = \frac{ dy ( z_0 - z_1 ) + dz ( y_1 - y_0 )
}{ y0 ( z_1 - z_2 ) + y_1 ( z_2 - z_0 ) + y_2 ( z_0 - z_1 ) }
\end{cases}$$
Using the barycentric coordinates $(u, v)$ and direction $(du, dv)$, we can trivially calculate the (signed) distance in that direction to each of the triangle edges:
Signed distance to edge between $\vec{v}_1$ and $\vec{v}_2$ ($v = 0$)is
$$d_{12} = - \frac{v}{dv}$$
signed distance to edge between $\vec{v}_2$ and $\vec{v}_3$ ($u + v = 1$) is
$$d_{23} = \frac{1 - u - v}{du + dv}$$
and signed distance to edge between $\vec{v}_1$ and $\vec{v}_3$ ($u = 0$) is
$$d_{31} = - \frac{u}{dv}$$
If the signed distance is zero, it means $(u, v)$ is on that specific edge. If the signed distance is negative, it means $(du, dv)$ is pointed away from that edge. Thus, the smallest positive signed distance determines which edge is intersected first.
When the walker arrives at an edge, we convert the $(du, dv)$ direction to a three-dimensional direction vector $\vec{d}$:
$$\vec{d} = du \left ( \vec{v}_2 - \vec{v}_1 \right )
+ dv \left ( \vec{v}_3 - \vec{v}_1 \right )$$
We can then rotate $\vec{d}$ in three dimensions, so that it is parallel to the new triangle face. To do this, we need the unit normal vectors for the old face, as well as the new face.
We can compute the normal from the vertex coordinates trivially, using vector cross product:
$$\hat{n} = \frac{ \left ( \vec{v}_2 - \vec{v}_1 \right ) \times \left ( \vec{v}_3 - \vec{v}_1 \right )}{\lVert \left ( \vec{v}_2 - \vec{v}_1 \right ) \times \left ( \vec{v}_3 - \vec{v}_1 \right ) \rVert}$$
If the unit normal vector in the old face is $\hat{n}_o$, and the unit normal vector in the new face is $\hat{n}_n$, then we need to calculate
$$\begin{array}{l}
\vec{a} = \hat{n}_o \times \hat{n}_n \\
\sin\varphi = \lVert \vec{a} \rVert \\
\cos\varphi = \hat{n}_o \cdot \hat{n}_n \\
\hat{a} = \frac{\vec{a}}{\sin\varphi} \end{array}$$
so that we can apply Rodrigues' rotation formula to rotate the walking direction:
$$\vec{d}' = \vec{d} \cos\varphi - \left ( \vec{d} \times \hat{a} \right ) \sin\varphi + \hat{a} \left ( \vec{d} \cdot \hat{a} \right ) ( 1 - \cos\varphi )$$
Numerically, we'll want to ensure $\vec{d}'$ is parallel to the new face, i.e.
$$\vec{d} = \vec{d}' - \hat{n}_n \left ( \frac{ \vec{d}' \cdot \hat{n}_n }{ \vec{d}' \cdot \vec{d}' } \right )$$
and that the position we are at, $\vec{p}$, is also on the plane of the face, to stop numerical errors from accumulating.
Here is an example C program, tetra.c, that walks along the faces of a tetrahedron:
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#define SQRT_HALF 0.7071067811865475244008443621048490392850
/*
* 3D vector operations
*/
typedef struct {
double x;
double y;
double z;
} vec3d;
static vec3d vec3d_def(const double x, const double y, const double z)
{
const vec3d result = { x, y, z };
return result;
}
static vec3d vec3d_sub(const vec3d v0, const vec3d v1)
{
const vec3d result = { v0.x - v1.x, v0.y - v1.y, v0.z - v1.z };
return result;
}
static double vec3d_dot(const vec3d v0, const vec3d v1)
{
return v0.x * v1.x + v0.y * v1.y + v0.z * v1.z;
}
static vec3d vec3d_cross(const vec3d v0, const vec3d v1)
{
const vec3d result = { v0.y * v1.z - v0.z * v1.y,
v0.z * v1.x - v0.x * v1.z,
v0.x * v1.y - v0.y * v1.x };
return result;
}
static double vec3d_len(const vec3d v)
{
return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
}
static vec3d vec3d_unit(const vec3d v)
{
const double vv = v.x*v.x + v.y*v.y + v.z*v.z;
if (vv > 0.0) {
const double d = sqrt(vv);
const vec3d result = { v.x / d, v.y / d, v.z / d };
return result;
} else {
const vec3d result = { 0.0, 0.0, 0.0 };
return result;
}
}
static vec3d vec3d_perpendicular(const vec3d v, const vec3d p)
{
const double pp = p.x*p.x + p.y*p.y + p.z*p.z;
if (pp > 0.0) {
const double c = (v.x*p.x + v.y*p.y + v.z*p.z) / pp;
const vec3d result = { v.x - c*p.x,
v.y - c*p.y,
v.z - c*p.z };
return result;
} else
return v;
}
/*
* 2D vectors
*/
typedef struct {
double x;
double y;
} vec2d;
/*
* Triangular face stuff
*/
typedef struct triface triface;
struct triface {
/* Three vertices, preferably in counter-clockwise order */
vec3d vertex[3];
/* Triangle plane equation is p.x*normal.x + p.y*normal.y + p.z*normal.z = distance */
vec3d normal;
double distance;
/* Neighbor faces sharing edges [0],[1]; [1],[2]; and [2],[0]. */
triface *neighbor[3];
};
/* Planar coordinates:
0 <= p.x <= 1
0 <= p.y <= 1
0 <= p.x + p.y <= 1
p.y == 0 on the edge between vertices 0 and 1
p.x == 0 on the edge between vertices 2 and 0
p.x+p.y == 1 on the edge between vertices 1 and 2.
*/
/* Return the number of the edge (0, 1, 2) we intersect at (p + (*ds)*d),
or -1 if none. */
static int planar_edge(const vec2d p, const vec2d d, double *const ds, const double eps)
{
double d0 = (d.y != 0.0) ? -p.y / d.y : -1.0;
double d1 = (d.x + d.y != 0.0) ? (1.0 - p.x - p.y) / (d.x + d.y) : -1.0;
double d2 = (d.x != 0.0) ? -p.x / d.x : -1.0;
if (d0 >= eps && (d0 <= d1 || d1 <= eps) && (d0 <= d2 || d2 <= eps)) {
if (ds)
*ds = d0;
return 0;
} else
if (d1 >= eps && (d1 <= d0 || d0 <= eps) && (d1 <= d2 || d2 <= eps)) {
if (ds)
*ds = d1;
return 1;
} else
if (d2 > eps && (d2 <= d0 || d0 <= eps) && (d2 <= d1 || d1 <= eps)) {
if (ds)
*ds = d2;
return 2;
} else
return -1;
}
/* Rodrigues' rotation: rotates v by theta around unit vector a,
with cosa = cos(theta) and sina = sin(theta). */
static vec3d vec3d_rotate(const vec3d v, const vec3d a,
const double cosa, const double sina)
{
const vec3d vxa = vec3d_cross(v, a);
const double ca = (1 - cosa) * vec3d_dot(v, a);
const vec3d result = { v.x * cosa + a.x * ca - vxa.x * sina,
v.y * cosa + a.y * ca - vxa.y * sina,
v.z * cosa + a.z * ca - vxa.z * sina };
return result;
}
static vec2d point_to_planar(const vec3d p, const vec3d v[3])
{
const double d1 = v[0].x*(v[1].y - v[2].y) + v[1].x*(v[2].y - v[0].y) + v[2].x*(v[0].y - v[1].y),
d2 = v[0].x*(v[1].z - v[2].z) + v[1].x*(v[2].z - v[0].z) + v[2].x*(v[0].z - v[1].z),
d3 = v[0].y*(v[1].z - v[2].z) + v[1].y*(v[2].z - v[0].z) + v[2].y*(v[0].z - v[1].z);
const double a1 = fabs(d1),
a2 = fabs(d2),
a3 = fabs(d3);
if (a1 >= a2 && a1 >= a3) {
const vec2d result = { ( p.x*(v[2].y-v[0].y) + v[0].x*(p.y-v[2].y) + v[2].x*(v[0].y-p.y) ) / d1,
( p.x*(v[0].y-v[1].y) + v[0].x*(v[1].y-p.y) + v[1].x*(p.y-v[0].y) ) / d1 };
return result;
} else
if (a2 >= a1 && a2 >= a3) {
const vec2d result = { ( p.x*(v[2].z-v[0].z) + v[0].x*(p.z-v[2].z) + v[2].x*(v[0].z-p.z) ) / d2,
( p.x*(v[0].z-v[1].z) + v[0].x*(v[1].z-p.z) + v[1].x*(p.z-v[0].z) ) / d2 };
return result;
} else {
const vec2d result = { ( p.y*(v[2].z-v[0].z) + v[0].y*(p.z-v[2].z) + v[2].y*(v[0].z-p.z) ) / d3,
( p.y*(v[0].z-v[1].z) + v[0].y*(v[1].z-p.z) + v[1].y*(p.z-v[0].z) ) / d3 };
return result;
}
}
static vec3d planar_to_point(const vec2d p, const vec3d v[3])
{
const vec3d result = { v[0].x + p.x*(v[1].x - v[0].x) + p.y*(v[2].x - v[0].x),
v[0].y + p.x*(v[1].y - v[0].y) + p.y*(v[2].y - v[0].y),
v[0].z + p.x*(v[1].z - v[0].z) + p.y*(v[2].z - v[0].z) };
return result;
}
static vec2d direction_to_planar(const vec3d d, const vec3d v[3])
{
const double d1 = v[0].x*(v[1].y-v[2].y) + v[1].x*(v[2].y-v[0].y) + v[2].x*(v[0].y-v[1].y),
d2 = v[0].x*(v[1].z-v[2].z) + v[1].x*(v[2].z-v[0].z) + v[2].x*(v[0].z-v[1].z),
d3 = v[0].y*(v[1].z-v[2].z) + v[1].y*(v[2].z-v[0].z) + v[2].y*(v[0].z-v[1].z);
const double a1 = fabs(d1),
a2 = fabs(d2),
a3 = fabs(d3);
if (a1 >= a2 && a1 >= a3) {
const vec2d result = { ( d.x*(v[2].y-v[0].y) + d.y*(v[0].x-v[2].x) ) / d1,
( d.x*(v[0].y-v[1].y) + d.y*(v[1].x-v[0].x) ) / d1 };
return result;
} else
if (a2 >= a1 && a2 >= a3) {
const vec2d result = { ( d.x*(v[2].z-v[0].z) + d.z*(v[0].x-v[2].x) ) / d2,
( d.x*(v[0].z-v[1].z) + d.z*(v[1].x-v[0].x) ) / d2 };
return result;
} else {
const vec2d result = { ( d.y*(v[2].z-v[0].z) + d.z*(v[0].y-v[2].y) ) / d3,
( d.y*(v[0].z-v[1].z) + d.z*(v[1].y-v[0].y) ) / d3 };
return result;
}
}
static vec3d planar_to_direction(const vec2d d, const vec3d v[3])
{
const vec3d result = { d.x*(v[1].x - v[0].x) + d.y*(v[2].x - v[0].x),
d.x*(v[1].y - v[0].y) + d.y*(v[2].y - v[0].y),
d.x*(v[1].z - v[0].z) + d.y*(v[2].z - v[0].z) };
return result;
}
int main(int argc, char *argv[])
{
vec3d v[4];
triface t[4], *face, *next;
vec3d p3, d3, a;
vec2d p, d;
double ds, sina, cosa, dn;
size_t i, j;
int k;
long n;
char dummy;
/* Tetrahedron. */
v[0] = vec3d_def( -1.0, 0.0, -SQRT_HALF );
v[1] = vec3d_def( +1.0, 0.0, -SQRT_HALF );
v[2] = vec3d_def( 0.0, -1.0, SQRT_HALF );
v[3] = vec3d_def( 0.0, +1.0, SQRT_HALF );
t[0].vertex[0] = v[0];
t[0].vertex[1] = v[1];
t[0].vertex[2] = v[2];
t[0].neighbor[0] = t + 1; /* Edge v[0],v[1] */
t[0].neighbor[1] = t + 3; /* Edge v[1],v[2] */
t[0].neighbor[2] = t + 2; /* Edge v[2],v[0] */
t[1].vertex[0] = v[0];
t[1].vertex[1] = v[3];
t[1].vertex[2] = v[1];
t[1].neighbor[0] = t + 2; /* Edge v[0],v[3] */
t[1].neighbor[1] = t + 3; /* Edge v[3],v[1] */
t[1].neighbor[2] = t + 0; /* Edge v[1],v[0] */
t[2].vertex[0] = v[0];
t[2].vertex[1] = v[2];
t[2].vertex[2] = v[3];
t[2].neighbor[0] = t + 0; /* Edge v[0],v[2] */
t[2].neighbor[1] = t + 3; /* Edge v[2],v[3] */
t[2].neighbor[2] = t + 1; /* Edge v[3],v[0] */
t[3].vertex[0] = v[1];
t[3].vertex[1] = v[3];
t[3].vertex[2] = v[2];
t[3].neighbor[0] = t + 1; /* Edge v[1],v[3] */
t[3].neighbor[1] = t + 2; /* Edge v[3],v[2] */
t[3].neighbor[2] = t + 0; /* Edge v[2],v[1] */
/* Calculate unit normals for each face,
pointing "outwards". */
for (i = 0; i < 4; i++) {
const vec3d n = vec3d_unit( vec3d_cross( vec3d_sub(t[i].vertex[1], t[i].vertex[0]),
vec3d_sub(t[i].vertex[2], t[i].vertex[0]) ) );
t[i].normal = n;
t[i].distance = vec3d_dot(n, t[i].vertex[0]);
}
if (argc != 6) {
fprintf(stderr, "\n");
fprintf(stderr, "Usage: %s\n", argv[0]);
fprintf(stderr, " %s u v du dv steps\n", argv[0]);
fprintf(stderr, "\n");
fprintf(stderr, "Tetrahedron:\n");
for (i = 0; i < 4; i++) {
fprintf(stderr, " Face %zu: %12.9f x %+12.9f y %+12.9f z = %.9f\n",
i + 1,
t[i].normal.x,
t[i].normal.y,
t[i].normal.z,
t[i].distance);
for (j = 0; j < 3; j++)
fprintf(stderr, " Vertex %zu: %12.9f %12.9f %12.9f\n",
j + 1,
t[i].vertex[j].x,
t[i].vertex[j].y,
t[i].vertex[j].z);
}
fprintf(stderr, "\n");
return EXIT_FAILURE;
}
if (sscanf(argv[1], " %lf %c", &p.x, &dummy) != 1 || p.x < 0.0 || p.x > 1.0) {
fprintf(stderr, "%s: Invalid initial u coordinate.\n", argv[1]);
return EXIT_FAILURE;
}
if (sscanf(argv[2], " %lf %c", &p.y, &dummy) != 1 || p.y < 0.0 || p.y > 1.0) {
fprintf(stderr, "%s: Invalid initial v coordinate.\n", argv[2]);
return EXIT_FAILURE;
}
if (p.x + p.y > 1.0) {
fprintf(stderr, "%s %s: Invalid initial u and v coordinates. Sum must not exceed 1.\n", argv[1], argv[2]);
return EXIT_FAILURE;
}
if (sscanf(argv[3], " %lf %c", &d.x, &dummy) != 1) {
fprintf(stderr, "%s: Invalid du (u component of initial direction).\n", argv[3]);
return EXIT_FAILURE;
}
if (sscanf(argv[4], " %lf %c", &d.y, &dummy) != 1) {
fprintf(stderr, "%s: Invalid dv (v component of initial direction).\n", argv[4]);
return EXIT_FAILURE;
}
if (d.x*d.x + d.y*d.y < 0.0001) {
fprintf(stderr, "%s %s: Direction vector (du and dv) is too short.\n", argv[3], argv[4]);
return EXIT_FAILURE;
}
if (sscanf(argv[5], " %ld %c", &n, &dummy) != 1 || n < 1L) {
fprintf(stderr, "%s: Invalid number of steps.\n", argv[5]);
return EXIT_FAILURE;
}
/* Initial face is face 0. Could make it too a parameter.. */
face = t+0;
p3 = planar_to_point(p, face->vertex);
d3 = vec3d_unit(planar_to_direction(d, face->vertex));
printf("%12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f\n",
p3.x, p3.y, p3.z, d3.x, d3.y, d3.z, p.x, p.y, d.x, d.y);
fflush(stdout);
while (n-->0L) {
/* Which edge will we intersect next? */
k = planar_edge(p, d, &ds, 0.0000000005);
if (k < 0)
break;
/* Advance to the edge. */
p.x += ds*d.x;
p.y += ds*d.y;
/* Location in 3D. */
p3 = planar_to_point(p, face->vertex);
/* Direction in 3D. */
d3 = vec3d_unit(planar_to_direction(d, face->vertex));
/* Output point and direction. */
printf("%12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f\n",
p3.x, p3.y, p3.z, d3.x, d3.y, d3.z, p.x, p.y, d.x, d.y);
fflush(stdout);
/* Do we have a next face? */
next = face->neighbor[k];
if (!next)
break;
/* Rotation axis for direction vector. */
a = vec3d_cross(face->normal, next->normal);
sina = vec3d_len(a);
a = vec3d_unit(a);
cosa = vec3d_dot(face->normal, next->normal);
/* Rotate direction around "bend". */
d3 = vec3d_rotate(d3, a, cosa, sina);
/* We switch to the new face. */
face = next;
/* Remove non-planar component from direction vector,
and renormalize to unit length for numerical stability. */
d3 = vec3d_unit(vec3d_perpendicular(d3, face->normal));
if (d3.x*d3.x + d3.y*d3.y + d3.z*d3.z <= 0.0)
break;
/* Ensure p3 is on the new face. */
dn = vec3d_dot(p3, face->normal);
if (dn != face->distance && dn != 0.0 && face->distance != 0.0) {
p3.x = p3.x * face->distance / dn;
p3.y = p3.y * face->distance / dn;
p3.z = p3.z * face->distance / dn;
}
/* Switch back to planar coordinates, but this time
using the current triangular face. */
p = point_to_planar(p3, face->vertex);
d = direction_to_planar(d3, face->vertex);
/* Output point and direction. */
printf("%12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f %12.9f\n",
p3.x, p3.y, p3.z, d3.x, d3.y, d3.z, p.x, p.y, d.x, d.y);
fflush(stdout);
}
return EXIT_SUCCESS;
}
In Linux and Macs, you can compile it using
gcc -Wall -O2 tetra.c -lm -o tetra
Run it without arguments to see the usage.
The edge detection is not terribly numerically stable. We could make it better by noting the edge normals are $(0, -1)$, $(1/\sqrt{2}, 1/\sqrt{2})$, and $(-1, 0)$, and only considering the edges with normals in the same half-plane as $(du, dv)$.
As an example, I ran
./tetra 0.1 0.3 31 71 100 > out
which I plot in Gnuplot using
splot "out" u 1:2:3 notitle w lines
which shows a nice, tetrahedral thread spool:
shows one such path in red.
