Barycentric coordinates seem like a reasonable way to go here. To review, the barycentric coordinates of a point $\mathbf p=(x,y)$ relative to a triangle with vertices $\mathbf v_0=(x_0,y_0)$, $\mathbf v_1=(x_1,y_1)$, $\mathbf v_2=(x_2,y_2)$ are a triple $\mathbf\lambda=[\lambda_0:\lambda_1:\lambda_2]$ such that $\mathbf p=\lambda_0\mathbf v_0+\lambda_1\mathbf v_1+\lambda_2\mathbf v_2$ and $\lambda_0+\lambda_1+\lambda_2=1$. One way to think of these coordinates is as the masses that must be be placed at the vertices of the triangle so that the center of mass is at $\mathbf p$. We can remove the normalization condition and work with unnormalized coordinates, in which case $\mathbf p={\lambda_0\mathbf v_0+\lambda_1\mathbf v_1+\lambda_2\mathbf v_2\over\lambda_0+\lambda_1+\lambda_2}$. These coordinates are homogeneous—the important part is the ratios among them, and two sets of coordinates that are non-zero scalar multiples of each other represent the same point.
The features of (normalized) barycentric coordinates that I’ll make use of are the following:
- Barycentric and Cartesian coordinates are related by a linear transformation.
- If a point is interior to or on the triangle, then all of its barycentric coordinates lie in the closed interval $[0,1]$.
- A zero coordinate indicates that the point lies on the side line opposite the corresponding vertex. A negative coordinate indicates that the point lies on the opposite side of that side line as the vertex. E.g., if the first coordinate is negative, then the point lies on the opposite side of the line through $v_1$ and $v_2$ as does $v_0$.
The Cartesian-to-barycentric mapping can be computed by expressing the definition of barycentric coordinates as the matrix equation $$\begin{bmatrix}x_0&x_1&x_2\\y_0&y_1&y_2\\1&1&1\end{bmatrix}\begin{bmatrix}\lambda_0\\\lambda_1\\\lambda_2\end{bmatrix}=\begin{bmatrix}x\\y\\1\end{bmatrix}.$$ The triangle is not degenerate, so that the matrix on the left is nonsingular and thus the barycentric coordinates are found by multiplying $(x,y,1)^T$ by the inverse of that matrix. This $3\times3$ inversion can be reduced to $2\times2$ by rewriting the coordinate equation as $$\begin{bmatrix}\mathbf v_0-\mathbf v_2 & \mathbf v_1-\mathbf v_2\end{bmatrix}\begin{bmatrix}\lambda_0\\\lambda_1\end{bmatrix}=\mathbf p-\mathbf v_2.$$ The third coordinate can then be computed as $\lambda_2=1-\lambda_0-\lambda_1$.
Another way to avoid a matrix inversion is to use the fact that the barycentric coordinate ratios are equal to the ratios of signed areas of the triangles formed by the point and pairs of vertices of the reference triangle. These areas can be computed as $$\frac12\begin{vmatrix}x_{i+1} & x_{i+2} & x \\ y_{i+1} & y_{i+2} & y \\ 1&1&1 \end{vmatrix}=\frac12(x_{i+1},y_{i+1},1)\times(x_{i+2},y_{i+2},1)\cdot(x,y,1)$$ with indices wrapping around mod 3. Recalling that the components of the product of a matrix and column vector are the dot products of the vector with the rows of the matrix, the area computations can be written as $$\frac12\begin{bmatrix}(\tilde{\mathbf v}_1\times\tilde{\mathbf v}_2)^T \\ (\tilde{\mathbf v}_2\times\tilde{\mathbf v}_0)^T \\ (\tilde{\mathbf v}_0\times\tilde{\mathbf v}_1)^T \end{bmatrix}\tilde{\mathbf p}.$$ (The tildes denote the homogeneous coordinates obtained by appending a $1$ to the Cartesian coordinate pair.) These coordinates can then be normalized by dividing by the area of the reference triangle, so the normalized barycentric coordinates of $\mathbf p=(x,y)$ are $$\mathbf\lambda=\frac1{\det\begin{bmatrix}\tilde{\mathbf v}_0&\tilde{\mathbf v}_1&\tilde{\mathbf v}_2\end{bmatrix}}\begin{bmatrix}(\tilde{\mathbf v}_1\times\tilde{\mathbf v}_2)^T \\ (\tilde{\mathbf v}_2\times\tilde{\mathbf v}_0)^T \\ (\tilde{\mathbf v}_0\times\tilde{\mathbf v}_1)^T \end{bmatrix}\tilde{\mathbf p}.$$ (You can also work with unnormalized coordinates, but you’ll have to use a different interval for determining when a point is in/on the triangle: the coordinates must be bounded by $0$ and $\lambda_1+\lambda_2+\lambda_3$, but the latter might be negative.)
Let the barycentric coordinates of the line segment’s end points be $\mathbf P=[\lambda_1:\lambda_2:\lambda_3]$ and $\mathbf Q=[\mu_1:\mu_2:\mu_3]$. If either of these points is in/on the triangle, i.e., $0\le\lambda_1,\lambda_2,\lambda_3\le1$ or $0\le\mu_1,\mu_2,\mu_3\le1$, then we’re done: the line segment obviously intersects the triangle.
Otherwise, we have two exterior points. Thanks to the linearity of Cartesian-to-barycentric mapping, the line through $\mathbf P$ and $\mathbf Q$ can be parameterized as $(1-t)\,\mathbf P+t\,\mathbf Q$ in both coordinate systems, with $t\in[0,1]$ giving the line segment. Values of $t$ for which a barycentric coordinate of the above linear combination is $0$ are the intersections with the corresponding side lines. These values of $t$ are easily found to be $t_i={\lambda_i\over\lambda_i-\mu_i}$, unless $\lambda_i=\mu_i$ in which case the line is parallel to that side of the triangle. The line segment intersects the triangle if any of these side line crossings occurs on the triangle. In addition, for a line segment with two exterior end points to intersect the triangle, there must be at least two side line crossings, regardless of where those crossing points are.
So, compare the corresponding barycentric coordinates of $\mathbf P$ and $\mathbf Q$. If fewer than two pairs have opposite signs, then the segment doesn’t intersect the triangle and we’re done. Otherwise, for each such coordinate pair $\lambda_i$ and $\mu_i$, compute the intersection point ${\mu_i\over\mu_i-\lambda_i}\mathbf P+{\lambda_i\over\lambda_i-\mu_i}\mathbf Q$. (It might be more convenient to use $1-{\lambda_i\over\lambda_i-\mu_i}$ instead for the coefficient of $\mathbf P$.) If any of these intersection points lies on the triangle, then the segment intersects the triangle.
The above computations were for points in $\mathbb R^2$, but you’re working in $\mathbb R^3$. Affine transformations preserve lines, their intersections, and relative segment lengths along a line, so the problem can be made two-dimensional via a parallel projection onto any convenient plane that’s not orthogonal to the plane of the triangle. Projecting onto one of the coordinate planes is simplest since that just involves dropping one of the Cartesian coordinates. If you’re using the $3\times3$ coordinate transformation matrix, projecting onto $z=1$ is convenient since that sets the third coordinate to $1$ as required to use the matrix.
It’s also possible to avoid an explicit projection if the triangle’s plane doesn’t pass through the origin. In that case, the areas of triangles on the plane are proportional to the volumes of tetrahedra defined by the triangle vertices and the origin, so we can use triple products of the points directly to compute barycentric coordinates. That is, we can compute $$\mathbf\lambda=\begin{bmatrix}(\mathbf v_1\times\mathbf v_2)^T\\(\mathbf v_2\times\mathbf v_0)^T\\(\mathbf v_0\times\mathbf v_1)^T\end{bmatrix}\mathbf p$$ using the raw coordinate triples of these points. If the triangle’s plane goes through the origin, shift all of the points by a small distance in a direction not parallel to the plane.
This solution is fundamentally similar to the one given by Nominal Animal. However, using barycentric coordinates allows all of the triangle’s sides to be treated uniformly. Also, if you think about it a bit, you’ll see that these computations are a series of “left-of” tests, unrolled.
The cross-product form of the Cartesian-to-barycentric mapping prompts me to toss this out as well, though I doubt that it’s as efficient as other solutions that have been presented. It takes advantage of the fact that in $\mathbb{RP}^2$, the line through two points and the intersection of two lines can both be computed by taking cross products of homogeneous coordinate vectors. Thus, the coordinate mapping matrix can be seen as having the three side lines of the triangle as its rows.
The barycentric coordinates of the three side line intersections can be computed “in bulk” as $B(\tilde{\mathbf p}\times\tilde{\mathbf q})_\times B^T$, where $B$ is the coordinate transformation matrix from above and the central matrix is the “cross-product matrix” defined as follows: for any vector $\mathbf w=(w_1,w_2,w_3)T$, $$\mathbf w_\times=\begin{bmatrix}0&-w_3&w_2\\w_3&0&-w_1\\-w_2&w_1&0\end{bmatrix}.$$ As before, you will need to check if any of these points are on the triangle and you will also have to check that viable candidates for the triangle intersections fall between $\mathbf p$ and $\mathbf q$. Moreover, these coordinates will be unnormalized, so you’ll have to take that into account when making these tests. You will also have to deal with points at infinity, which will be produced when the segment is parallel to one of the sides of the triangle. In barycentric form, the coordinates of points at infinity sum to zero.