1
$\begingroup$

I have to figure out how to get the path of points that make up a geodesic on a triangular mesh:

  1. 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?

  2. How do we find the geodesic along the surface of the mesh, from one known point to another known point?

For example, Triangular mesh shows one such path in red.

We know which triangles are touching each other.


I need to figure out the following:

  1. The point at which the hypothetical red line intersects the edge of the triangle.

  2. Assuming it intersects the triangle, what is the new direction based on this path?

$\endgroup$
4
  • 1
    $\begingroup$ Your model is usually called a triangle mesh I think. Your mesh has strange properties, but it may be worth mentioning that name so that people knowledgeable about the topic can get into your problem faster. For your question 1, I guess it depends how the direction is given. For question 2, I'm afraid you'll have to deal with the levi civita connection... I'm a newbie on the topic, so I'm not too sure. $\endgroup$
    – N.Bach
    May 23, 2017 at 16:42
  • $\begingroup$ Actually, forget about levi civita, you're only moving a point around so I guess you don't need to worry about this... $\endgroup$
    – N.Bach
    May 23, 2017 at 16:43
  • $\begingroup$ See stackoverflow.com/questions/25814549/… $\endgroup$
    – lhf
    Jan 28, 2020 at 7:18
  • $\begingroup$ You can probably implement a pseudogedesics algorithm, e.g. the straightest geodesics by K. Polthier would come to mind dl.acm.org/doi/abs/10.1145/1185657.1185664 $\endgroup$
    – teodron
    Jan 20, 2021 at 8:22

2 Answers 2

5
$\begingroup$

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.)

Triangle coordinate system

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 )$$

$\endgroup$
9
  • $\begingroup$ @TheGreatDuck: Good question. The "triangle coordinates" are plane coordinates, with origin at $\vec{v}_i$, and $u$ axis towards $\vec{v}_j$; i.e. $\vec{v}_i$ is at $(0, 0)$, and $\vec{v}_j$ at $(\lVert\vec{v}_j - \vec{v}_i\rVert, 0)$. The $v$ axis is defined as perpendicular to $u$ axis, and is on the half-space towards the third point. Is this sufficient to clarify it, or should I expand this part in my answer? $\endgroup$ May 24, 2017 at 18:01
  • $\begingroup$ @TheGreatDuck: I edited the answer to explain the "triangle coordinates" and their calculation in detail; let me know if there is something odd or unclear (me often fail English). There is probably an easier way to transform triangle coordinates between order permutations directly (without computing the full 3D coordinates first), but I haven't bothered to check. $\endgroup$ May 24, 2017 at 19:18
  • 1
    $\begingroup$ @TheGreatDuck: Do not give up; the hardest part is getting your brain to grok the ideas here. The math will follow. For myself, sketching and playing with tangible models definitely helps. I created a printable tetrahedron model, where each vertex is numbered (0 to 3), and each face has the $u$,$v$ axes marked (with $u$ axis starting at the smallest numbered vertex, pointing to the next vertex in counterclockwise order). Print it, cut the shape, fold it, apply glue to the three tabs, and you get a real 3D object. $\endgroup$ Jun 5, 2017 at 12:13
  • $\begingroup$ @TheGreatDuck: No worries! I find basic vector algebra extremely useful. I'm definitely not a mathematician, and strictly use math as a tool; and vector algebra in 2 and 3 dimensions has been extremely useful tool for me across a big swath of problems I've solved. Because of that, I do very warmly recommend you pick it up; it really is useful in solving a number of problems. Also, I personally don't like juggling too many things at once, so if that is one too many thing to attack right now, just make a mental note to come back to it later on: it is good to choose where to concentrate on. $\endgroup$ Jun 5, 2017 at 16:25
  • 1
    $\begingroup$ @TheGreatDuck: By "d is in the plane of the triangle", I mean that the dot product between vector d and the surface normal of the triangle must be 0. $\endgroup$ Jun 5, 2017 at 16:49
3
$\begingroup$

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: walk over a tetrahedron

$\endgroup$
5
  • $\begingroup$ @TheGreatDuck: I think I might be able to edit your question in a way that would make it easier for others to understand. (I think I understand you because I wondered about basically the same question years ago -- I was thinking about procedurally generated planetoids, starting with an icosahedron, with each face subdivided into four triangles, recursively, with a radius perturbation function applied at each step. The problem was, how would a vehicle travel on the surface.) Do you mind? You can also email me directly; my address is at my homepage, linked to from my profile here. $\endgroup$ Jun 8, 2017 at 14:34
  • $\begingroup$ @TheGreatDuck: Also, I only right now noticed that your triangles are right triangles. (I do not believe they are precisely that, as typical meshes tend to have triangles that differ from right triangles very slightly, unless the mesh is exactly flat.) My solutions work for general triangular meshes. However, if you do have only right triangles, I'm sure that could be leveraged to make the solution even simpler -- for one, the shorted edges make natural planar axes. $\endgroup$ Jun 8, 2017 at 15:10
  • $\begingroup$ @TheGreatDuck: Okay; I just wanted to know that you are not opposed to editing your question in principle, that you'd be willing to consider it on its merits. Just don't make a hasty decision; I shall suggest a complete rewrite, so that the mathematicians over here see the problem the way (I think you and) I do. Also, I think an illustration is necessary, so the edit will require some effort on my part too. $\endgroup$ Jun 8, 2017 at 18:14
  • $\begingroup$ @TheGreatDuck: I tried to distill the question into two sub-questions I can parse from it. Feel free to edit it in any way you want, but do try to keep to the core of the question. You could even consider splitting the second part into a separate question. (For the first part I have shown an example program, above. For the second part, I think the path can be approximated by any 3D plane intersecting the object, leading to a 2D walker to perceive walking in a straight line. This does not take into account any local features in the mesh object, of course.) $\endgroup$ Jun 8, 2017 at 19:46
  • $\begingroup$ @TheGreatDuck: Excellent! Do feel free to edit it to reflect your question. (I am not invested in the question, or my wording, at all; if you can rewrite it to express your question better, absolutely you should!) What I would like to see, is the alternate solution you mentioned. Note that "triangle mesh" is the correct underlying term, because walking from one triangle to another is exactly that: walking on a triangle mesh. $\endgroup$ Jun 8, 2017 at 20:11

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .