One very simple approach is to use a linear model. Let $I(x,y)$ be the pixel intensity of the pixel at coordinates $(x,y)$. We can fit to a linear model:
$$I(x,y) = \alpha x + \beta y + \gamma,$$
where $\alpha,\beta,\gamma$ are parameters to be determined. If we are given the pixel intensity at each corner of the triangle, then we know the value of $I(x,y)$ for each corner point $(x,y)$. Since a triangle has 3 corners, this gives us 3 linear equations in 3 unknowns. We can solve this system of linear equations to find $\alpha,\beta,\gamma$, and then use this linear model to fill in the value of all other points.
In particular, if the $i$th corner is at $(x_i,y_i)$ and has pixel intensity (value) $v_i$, then we get a system of three equations:
I(x_1,y_1) &= v_1\\
I(x_2,y_2) &= v_2\\
I(x_3,y_3) &= v_3
Expanding the linear model, this becomes
\alpha x_1 + \beta y_1 + \gamma &= v_1\\
\alpha x_2 + \beta y_2 + \gamma &= v_2\\
\alpha x_3 + \beta y_3 + \gamma &= v_3
Here we have three unknowns ($\alpha,\beta,\gamma$) and the rest are known constants. So, we can solve this system of linear equations for $\alpha,\beta,\gamma$. This involves inverting a $3 \times 3$ matrix and then multiplying it by a $3$-vector, so it should be relatively efficient.
Once you have derived the parameters $\alpha,\beta,\gamma$, you can now evaluate the value $I(x,y)$ at any other point $(x,y)$ in the triangle (or on its interior). This lets you interpolate the value at each fill-in point.
The approach I described above works for black-and-white images. If you have a color image, you can interpolate independently for each of the 3 channels (R, G, and B; you might get better results if you work in a transformed color space, e.g., HSV, CIELAB, CIECAM02).
I recommend you take a look at barycentric coordinates, as they simplify the formulas for interpolation in a triangular mesh.
If you are familiar with interpolation for image processing, bilinear interpolation is a slight generalization of this approach, where you add a $\delta xy$ term to the linear model. However, bilinear interpolation isn't applicable here, because it requires 4 points: 3 points isn't enough.
It's possible to come up with more sophisticated interpolation models. For instance, you could fit a linear model to more than 3 points (e.g., taking the corner points of all adjacent triangles as well), and use ordinary least squares linear regression to find a best-fit linear model, weighting the corners of the adjacent triangles lower than the corners of the target triangle. You could also build generalizations of bicubic interpolation, etc.
More generally, I recommend that you do a bit of research on "tetrahedral interpolation", as I believe that is the catch-phrase for the kind of problem you are asking about (or the generalization of it to 3 dimensions).