Reading farther on the Linked Wikipedia page, we have two cases:
Case 1: the triangle has an angle $\ge120^{\circ}$: the Fermat point is the obtuse-angled vertex.
Case 2: the triangle does not have an angle $\ge120^{\circ}$: the Fermat point is the first isogonic center. Following the link to this page, we find that point has trilinear coordinates $\csc(A+120^{\circ}):\csc(B+120^{\circ}):\csc(C+120^{\circ})$. If, like me, you are not familiar with trilinear coordinates, follow the link to its wikipedia page, where we find a method to convert trilinear coordinates to cartesian:
Taking vertex $C$ as the origin, find vectors $\vec{A}$ and $\vec{B}$ corresponding to the vertexes $A$ and $B$ respectively. Given side lengths $a, b, c$ corresponding in the usual way to the sides of the triangle respectively opposite vertices $A, B, C$, a trilinear coordinate $x:y:z$ is converted to a vector in cartesian coordinates by $$\vec{P}=\frac{ax}{ax+by+cz}\vec{A}+\frac{by}{ax+by+cz}\vec{B}.$$
In short, you have a piecewise function that checks for angle size and, if greater than 120 degrees, returns the $\vec{P}$, plus a correction $\vec{C}$ representing the location of $C$ if it is not convenient to have $C$ at the origin. You could embed the computation for $\vec{A}$ and $\vec{B}$ into the formula to have a single function result. That computation is necessary only if you don't already have coordinates for the vertices, and is elementary trig.
The equation of the plane is given by
$$\left|\begin{matrix}x-A_x & y -A_y & z-A_ z \\ B_x-A_x & B_y -A_y & B_z-A_ z\\ C_x-A_x & C_y -A_y & C_z-A_ z\end{matrix}\right|=0.$$
Thus:
$$z=A_z+\frac{(B_x-A_x)(C_z-A_z)-(C_x-A_x)(B_z-A_z)}{(B_x-A_x)(C_y-A_y)-(C_x-A_x)(B_y-A_y)}(y-A_y)\\-\frac{(B_y-A_y)(C_z-A_z)-(C_y-A_y)(B_z-A_z)}{(B_x-A_x)(C_y-A_y)-(C_x-A_x)(B_y-A_y)}(x-A_x)$$
Best Answer
Although others have already answered the stated question, I'd like to give some advice to those doing similar to what OP is doing.
Precalculate the unit normal vector to the triangle, and store it.
If the triangle vertices are $\vec{v}_0$, $\vec{v}_1$, and $\vec{v}_2$, then the triangle normal is usually defined as $$\vec{n} = \bigl(\vec{v}_1 - \vec{v}_0\bigr)\times\bigl(\vec{v}_2 - \vec{v}_0\bigr)$$ where $\times$ denotes vector cross product. This is scaled to unit length, $$\hat{n} = \frac{\vec{n}}{\left\lVert\vec{n}\right\rVert} = \frac{\vec{n}}{\sqrt{\vec{n} \cdot \vec{n}}}$$
This is best precalculated and stored, because the square root operation tends to be slow, and it is just one vector per triangle.
A point is "above" the triangle if it is in the same halfspace as the triange normal vector, with respect to the plane formed by the triangle.
(If you extend the triangle plane to infinity, it splits space into two halves. If you look at the triangle sideways, the point is "above" the triangle if it is in the same side as the normal is; and "below" if it is on the other side. If the point is on the triangle plane exactly, I recommend you use a clear and unambigous term for it; I use "exactly on the triangle" myself. I'm sure mathematicians have more precise terms, but for those of us who only use math as a tool for solving problems, making sure the terms are unambiguous suffices.)
This means that point $\vec{p}$ has signed distance $$d = \bigl(\vec{p} - \vec{v}_0\bigr) \cdot \hat{n}$$ from the plane of the triangle. (Because all three triangle vertices are on the plane, you can use $\vec{v}_1$ or $\vec{v}_2$ in above; you'll still get the exact same result, within rounding error.)
If $\hat{n}$ is not an unit vector, i.e. $\left\lVert\hat{n}\right\rVert\ne 1$, the distance is in units of the normal vector length.
When you know the signed distance $d$ of the point $\vec{p}$ from the plane, you can trivially project the point to the plane: $$\vec{p}_\text{plane} = \vec{p} - d \hat{n}$$ If $\hat{n}$ is not an unit vector, i.e. $\left\lVert\hat{n}\right\rVert\ne1$, the right side is actually $\vec{p} - \hat{n}(d / \lVert\hat{n}\rVert)$, and you need a "slow" square root operation here. Because you probably test many points per each triangle, it is well worth precalculating and storing the $\hat{n}$ beforehand.
To determine if $\vec{p}_\text{plane}$ is within the triangle, find out the barycentric coordinates $(u, v)$ corresponding to $\vec{p}_\text{plane}$ with respect to that triangle.
To do this, check which axes have the smallest values in magnitude in $\hat{n}$. Their order does not matter. (This essentially picks the axis along which to look at the triangle, with the triangle being as large as possible. We do this to minimize numerical errors.) There are three possibilities: $x$ and $y$; $x$ and $z$; and $y$ and $z$. The solution is obtained using the coordinates on those two axes.
The solution for $\vec{p}_\text{plane} = (x, y, z)$ with $\hat{n}$ having $x$ and $y$ components closest to zero, is $$\begin{cases} \displaystyle u = \frac{ (x - x_0)(y_2 - y_0) - (y - y_0)(x_2 - x_0)}{ (x - x_2)(y_1 - y_0) - (y - y_2)(x_1 - x_0) } \\ \displaystyle v = \frac{ (y - y_0)(x_1 - x_0) - (x - x_0)(y_1 - y_0)}{ x_0 (y_1 - y_2) + x_1 (y_2 - y_0) + x_2 (y_0 - y_1) } \\ \end{cases}$$ The solution for the other axis pairs you'll find just by permuting the coordinates ($y$, $y_0$, $y_1$, $y_2$ to $z$, $z_0$, $z_1$, $z_2$ for $x$ and $z$ axes; and also $x$, $x_0$, $x_1$, $x_2$ to $y$, $y_0$, $y_1$, $y_2$ for $y$ and $z$).
Also note that you'll want to calculate the divisor for $u$ first. It will be zero at $\vec{v}_2$, where $u = 0$, so when the divisor is very close to zero, just use $u = 0$.
In barycentric coordinates, $\vec{v}_0$ is at $(0, 0)$; $\vec{v}_1$ is at $(1,0)$, and $\vec{v}_2$ is at $(0,1)$. The entire triangle is $0 \le u, v, u + v \le 1$. The edges are at $u = 0$, $v = 0$, and $u + v = 1$.
In practice, you define the triangle inclusion test as $$\Bigl(u \ge -\epsilon\Bigr) \text{ and } \Bigl(v \ge -\epsilon\Bigr) \text{ and } \Bigl(u \le 1+\epsilon\Bigr) \text{ and } \Bigl(v \le 1+\epsilon\Bigr) \text{ and } \Bigl(u + v \le 1+\epsilon\Bigr)$$ where $\epsilon$ represents the precision of your numerical coordinates. It is the largest positive number you consider equal to zero. (Similarly, if $D$ is the divisor for $u$ earlier, if $-\epsilon \le D \le +\epsilon$, use $u = 0$.)
To linearly interpolate any value $c$ within the triangle, with $c = c_0$ at $\vec{v}_0$, $c_1$ at $\vec{v}_1$, and $c_2$ at $\vec{v}_2$, you can use $$c(u,v) = (1 - v)\bigl((1-u) c_0 + u c_1\bigr) + v c_2$$
All together, this means that you need only 15 (scalar) multiplications and less than thirty additions or subtractions, to find out if a point is within the triangle, or within a triangular prism extending from the triangle along its normal, including the barycentric coordinates needed for texture or color interpolation. (Color interpolation itself is just four multiplications per component, so typically 12 additional multiplications.) Nifty, eh?