For each pixel, we have to evaluate the discrete VRI, $$\Lout = \displaystyle \sum_{i=0}^n C_i\ \alpha_i\ \left(\ \prod_{j=0}^{i-1}\ (1-\alpha_j)\ \right)$$
Let $\ell(t)$ be the line from the eye through the pixel and let $\ell(t_i)$ be the 3D position on that ray in the middle of the $\ith$ interval, as shown below:
For $i = n$ we'll let $C_i\ \alpha_i = \Lin$ for the backlighting.
However, for $i \lt n$, we need to evaluate $C_i$ and $\alpha_i$.
Computing $\mathbf{\alpha_i}$
The step size, $\Delta s$, is a constant set by the user. A smaller $\Delta s$ will give more accurate renderings, but will take more time.
Since $\alpha_i = \tau_i\ \Delta s$, we need only to compute $\tau_i$.
The CT is provided as a volume consisting of CT numbers in Hounsfield units, and is usually converted to a 3D texture of attenuation coefficients, $\tau(x,y,z)$, at each of the discrete $(x,y,z)$ points in the CT volume.
An accurate conversion would require a knowledge of the parameters of the CT (such as voltage), but the usual process is to just linearly transform the CT numbers into the range [0,1].
To determine $\tau_i$, we need to look up $\tau( \ell(t_i) )$ in the 3D texture.
Since $\ell(t_i)$ is usually not at the integer coordinates of one of the texels, we ask the GPU to compute a tri-linear blend of the eight nearest texels. That's done in OpenGL with the usual
glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR ); glTexParameteri( GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR );
and the texture unit is set up with
glActiveTexture( GL_TEXTURE0 + VOLUME_TEXTURE_UNIT ); // VOLUME_TEXTURE_UNIT = ID of texture unit glBindTexture( GL_TEXTURE_3D, volumeTextureID ); // volumeTextureID = ID of texture
So computing $\alpha_i = \tau_i\ \Delta s$ is pretty straightforward.
Computing $\mathbf{C_i}$
We'll use the Phong illumination model to bounce light off of "surfaces" in the volume:
Then $$C_i = k_d\ N \cdot L + k_s\ (R \cdot V)^n$$
The direction, $L$, to the light and the direction, $V$, to the viewer are known. And $n$ is the shininess constant.
Ways to determine $k_d$ are discussed in the Transfer Functions lecture.
The surface normal, $N$, will be the defined using the gradient of the attenuation field, or $$\nabla \tau(x,y,z) = ( {\partial\: \tau \over \partial\: x}, {\partial\: \tau \over \partial\: y}, {\partial\: \tau \over \partial\: z} )$$
Since the gradient points in the direction of fastest increase of $\tau$, it will point in to a surface, since density increases when going from above a surface to below a surface.
Since the surface normal, $N$, should point out from a surface, we define $$N = - {\nabla \tau \over |\nabla \tau|}$$
Computation of the gradient is discussed in the Computing Gradients lecture.
For volume rendering, it's often sufficient to use only the diffuse term of the Phong model.