Ray intersecting a quad mesh I am trying to solve the math behind rendering a quad-mesh surface. MatLab for instance can take a regularly spaced (x,y) grid with arbitrary third-dimension (z) values, treat each four neighbouring coordinates as 3d quadrilaterals, and render each of them as a facet of the surface. 

I'm trying to do this with ray tracing which would be easy if the quads were planar (all corners on the same plane), but since every coordinate can be arbitrary most quads will likely not be on the same plane, a so-called "skew-quadrilateral" (you can see how most quads in the above MatLab image "twist and turn"). Does anyone know the steps or pseudocode to intersect a ray with a skew quadrilateral, I've tried looking around but couldn't find it? 
I don't even know if such a thing could be possible since you could have three corners of a quadrilateral with z-coordinates of 0 and suddenly the last corner with a very high z-coordinate, implying some sort of bending or breaking of the surface. I mean what would the following grid even look like, perhaps four quadrilaterals that start on the same level and bend towards a high point in the middle?
(0,0,0) (0,0,0) (0,0,0)
(0,0,0) (0,0,5) (0,0,0)   ===>    ?
(0,0,0) (0,0,0) (0,0,0)

I suppose the main difference could be that MatLab does not use ray tracing to render each quadrilateral, but rather some type of 3d to 2d reprojection. If so, does anyone know the usual non-raytracing way of rendering a quad mesh surface?
On the other hand I have read somewhere that quadrilateral surface meshes are indeed used in 3d movie or game developing where ray tracing is certainly used, so it should be possible. 
Links
Here are some relevant quotes and links that describe how MatLab treats the surface:
http://se.mathworks.com/help/matlab/visualize/representing-a-matrix-as-a-surface.html
"The plot is formed by joining adjacent points with straight lines." 
"Mesh plots are wire-frame surfaces that color only the lines connecting the defining points. Surface plots display both the connecting lines and the faces of the surface in color." 
"...generates a colored, faceted view of the surface and displays it in a 3-D view. Ordinarily, the facets are quadrilaterals, each of which is a constant color..."
http://se.mathworks.com/help/matlab/ref/surf.html
"Each point in the rectangular grid can be thought of as connected to its four nearest neighbors."
"This defines a mesh of quadrilaterals or a quad-mesh."
 A: In practice, raytracers tend to decompose surfaces into triangular meshes first, because the math needed is so much simpler.
However, you asked how to calculate the intersection between a line and a quadrilateral in 3D, so here goes.
Let's define your ray using an unit vector $\hat{n}$ (unit referring to unit length, $\lvert\hat{n}\rvert=1$) that passes through some point $\vec{p}_0$,
$$\vec{p}_{RAY}(t) = \vec{p}_0 + t\,\hat{n}$$
The quadrilateral is a bit more complicated. Let's say the four corners are $\vec{p}_1$, $\vec{p}_2$, $\vec{p}_3$, and $\vec{p}_1$, with $\vec{p}_1$ and $\vec{p}_4$ diagonally across from each others. We can trace the surface using two variables, $0 \le u, v \le 1$, where $u=0,v=0$ refers to $\vec{p}_1$, $u=1,v=0$ to $\vec{p}_2$, $u=0,v=1$ to $\vec{p}_3$, and $u=1,v=1$ to $\vec{p}_4$, using bilinear interpolation of the coordinates:
$$\vec{p}_{QUAD}(u,v) = \left ( \vec{p}_1 (1-u) + u \vec{p}_2 \right ) (1-v) + v \left ( \vec{p}_3 (1-u) + u \vec{p}_4 \right )$$
which is, after rearranging the terms, the same as
$$\vec{p}_{QUAD}(u,v) = \vec{p}_1 + u \; v \; ( \vec{p}_4 - \vec{p}_3 - \vec{p}_2 + \vec{p}_1 ) + u \; ( \vec{p}_2 - \vec{p}_1 ) + v \; ( \vec{p}_3 - \vec{p}_1 )$$
The intersection is, of course,
$$\vec{p}_{RAY}(t) = \vec{p}_{QUAD}(u,v)$$
In three dimensions, that is actually three equations with three unknowns,
$$\begin{cases}
x_0 + t n_x = x_1 + u \; v \; ( x_4 - x_3 - x_2 + x_1 ) + u \; ( x_2 - x_1 ) + v \; ( x_3 - x_1 ) \\
y_0 + t n_y = y_1 + u \; v \; ( y_4 - y_3 - y_2 + y_1 ) + u \; ( y_2 - y_1 ) + v \; ( y_3 - y_1 ) \\
z_0 + t n_z = z_1 + u \; v \; ( z_4 - z_3 - z_2 + z_1 ) + u \; ( z_2 - z_1 ) + v \; ( z_3 - z_1 )
\end{cases}$$
This can be solved, but the solution contains dozens of terms, and is therefore terribly slow to compute. (I asked Maple for the exact solution. There are two (i.e., $(t_1,u_1,v_1)$ and $(t_2,u_2,v_2)$, but as they don't fit in one screenful, I decided they are way too long to reproduce here.)
Biquadratic and bicubic Bézier surfaces can be solved the exact same way, it's just that there are more terms (up to $u^3$ and $v^3$, and 9 (biquadratic) 16 (bicubic) coordinates per dimension), and thus the result is even more complicated and slow to compute.
A: There are many questions, here. I'll take stab at a few of them.
The easiest way to render a set of triangles is a z-buffer algorithm. You paint each triangle into a pixel array, and you keep track of z-depth at each pixel. Pixels that are in front over-write pixels that are behind. Most graphics packages (like OpenGL or DirectX) will do this for you; all you have to do is give them the list of triangles.
Another option is the so-called "painter's algorithm". You sort the triangles by z-depth, and then paint them in order. Again, the front triangles over-write the back ones.
Ray-tracing is typically the slowest algorithm, and most people don't use it unless they want some very sophisticated graphical effects, like inter-object reflection, or refraction through a translucent material. Ray tracing is not typically used in games. It's used in movie animation, where image quality is more important than speed of computation.
How to break a quadrilateral into a triangle: there are two choices. Choose the one that produces the least deviation between the triangles and the original quad.
Intersecting  a ray with a "skew" quadrilateral. The quadrilateral is actually a hyperbolic paraboloid, which is a quadric surface (of degree 2). So, in general, any ray will intersect the quad in two places. The points of intersection can be found by solving a quadratic equation where the variable is the parameter value on the ray. That's what's described in the other answer. I'm surprised that it's so complicated. I would expect that you can eliminate the surface parameter values, leaving just a quadratic in the ray parameter. But, I haven't tried to do this, so maybe it's harder than I think.
I think decomposing into triangles is the right approach. If you want the graph to look more smooth and less "tesselated", you can vary the shade across each triangle. Instead of the color being constant on the triangle, you can linearly vary the color (Gouraud shading) or linearly vary the surface normal (Phong shading). Or if you want a "quad" appearance as opposed to a "triangulated" appearance, just make sure you use the same color to paint the two triangles that came from one original quad.
