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.