What are Voxels?

A hobbyist level overview of voxels, voxel space and Phong shading

By Geoffrey Clements, Chelmsford, Massachusetts

About the author

Geoffrey Clements can be reached on internet at clementsg@gw1.hanscom.af.mil.

Introduction

The medical industry is studying ways of viewing medical imagery on the computer. Doctors are using computer graphics to plan operations before going into the operating room. In fact, entire operations are dry run on the computer, before the doctor lifts a knife.

To create the data for a medical image, a radiologist will take a series of computer aided tomography (CAT) scan images at different depths. The images are lined up in order of depth to form a three dimensional representation of the object being scanned. Each data point in this three dimensional cube is called a voxel. The cube of voxels are called a voxel space. Think of a voxel as a measure of the density of matter at that point. You can also think of a voxel as a “volume pixel”.

In this article we’ll draw a sphere that appears three dimensional using the techniques used to draw medical imagery. I chose a sphere to make the problem simpler to understand. But first we’ll check out how a medical image is generated.

A viewing window defines the orientation of the voxel space, and what part of the voxel space is in view. This viewing window corresponds to the computer screen. The image is created by following a ray from the viewing window through the voxels and summing the contributions of the voxels the ray passes through. The ray is perpendicular to the viewing window. Figure 1 shows two voxels spaces. The first has the viewing window on the front face. The second has the viewing window on the top left front corner. The arrow is the ray followed through the voxel space.

Figure 1 Voxel Space and Viewing Windows

The contribution of each voxel is calculated from a transparency and a color. The transparency, sometimes called opacity, is calculated from the intensities associated with each voxel. Usually the intensity is just scaled to yield an opacity. In our example, inside the sphere is opaque, and outside the sphere is transparent. Figure 2 shows the voxel space we will be using and the sphere opacity.

Figure 2 Voxel Space

To characterize the viewing environment we’ll define several vectors. First a normalized vector pointing toward the viewing window. Next we define normal vectors at each voxel. Finally we need a normal vector pointing toward each light source.

The color of a particular pixel is calculated using the intensities of the voxels and these vectors defining the viewing environment. We’ll be using the Phong shading model to calculate the color.

Phong Shading

The Phong shading model is used in order to give objects in an image their three dimensional quality. We’ll go through the Phong shading equation by describing each piece. I use the words “shade” and “color” interchangeably. They mean the same thing. We start with:

Shade is the color of the voxel we are calculating. Ip is the intensity of a light source. ka is the ambient light reflection coefficient. The ambient light reflection coefficient, as you might expect, just sets an overall light level for the image. A gray disk is drawn if we use this equation to calculate the color.

Next, we add some depth cueing, which will start to give the image a 3D look.

kd1 and kd2 are depth cueing coefficients. We set kd2 = 1 in the example program. k is the distance from the viewing plane along the ray. This adds a small 3D effect. A gray disk is drawn that gets slightly lighter as we move from the edge of the sphere into the center.

Now we add the effect of diffuse and specular reflected light. (This step is a big jump, but these effects are what makes Phong shading so good.)

Figure 3 shows the difference between diffuse and specularly reflected light. Diffuse reflection is light reflected in all directions. This is caused by the roughness of the reflecting surface. Specular reflection, on the other hand, is light that is reflected in only one direction. Specular reflection is caused by the smoothness of the reflecting surface.

Figure 3 Reflected light

Here is the Phong shading equation with the effect of diffuse and specularly reflected light added.

d is the diffuse reflection coefficient, and s is the specular reflection coefficient. N is a normal vector to the voxel whose shade we are calculating, V is a vector pointing toward the viewer, L is a vector pointing at the light source and H is the normalized sum of V and L. N•L is the dot product of N and L.

We’ve come to the point where we can set our viewing window and define the vector space. Figure 4 shows the vector space. There are two L vectors shown in Figure 4. We’ll use two light sources to make the image a little more interesting. To add another light source we only need to add the effect of diffuse and specular reflection to the shade calculated above.

Figure 4 Vector space

The axes are slightly unusual, but the x and y axes correspond to the Macintosh coordinate system. N, H, L1, L2, and V are drawn to show their relative directions. This space corresponds to a viewer looking into the front face, (through the viewing window). The two light sources are pointing down from the top-left front corner of the cube, and into the right side. These vectors must be set before we start.

We can pre-calculate most of the shade equation before going into the loop that calculates the pixel color. This speeds up drawing the image.

Here are the definitions for V, L1, L2, H, and N.

r is the radius of the sphere, i, j, k are a point on the surface of the sphere. In our example, the center of the sphere is at (cv, cv, cv).

We can generate equations for N•L and N•H.

At this point you may be asking yourself, “How do you sum the contributions from the voxels?” Here’s how the example program does it. The indexes i and j step through the pixels on the viewing window. The index k corresponds to the ray going back into the voxels. As k varies, we check to see if point (i, j, k) is inside or outside the sphere. If it is outside the sphere, return a color of zero. If k gets to the back of the voxel space without hitting the sphere, draw a black pixel at (i, j) on the screen. If we hit the sphere, calculate a color for the sphere surface and draw that color at (i, j) on the screen, then go on to the next pixel. There is no need to process any farther inside the sphere, because there is no way the light can get there. Normally the effects of all voxels the ray passes through are used to calculate a color, but I have chosen to make the problem (and program) simpler to understand.

The Example Program

The code was developed under Think C. The MacTraps library is the only library that needs to be included. Turn on the Native floating-point format and Generate 68881 instruction switches in the Compiler settings screen if you have a floating point coprocessor. If not, just turn on the Native floating-point format switch. The program takes about a minute to generate a 128x128 pixel image on a Mac IIci with 68881 instructions on, and three minutes with it off.

The example program uses a standard Macintosh event loop shell and main routine. All of the interesting processing is done by DoColor() called by the Init() routine. Init() starts by initializing the the Macintosh managers and setting up the menus. If Color Quickdraw is not available, the program quietly exits. If Color Quickdraw is available, a window is opened and sized to our voxel data. The size of the volume and sphere are set with the defines:

/* size of the voxel data */
#define volSize  64
/* half the size of voxel data */
#define halfVolSize32
/* the radius of the sphere */
#define sphere_r 30
/* sphere_r*sphere_r */
#define volumeMag900.0
/* sqrt(3.0)*sphere_r */
#define sqrt3r   51.9615
/* sphere_r*sqrt(6+2*sqrt(3)) */
#define rsqrt6   92.2913

Notice that some of the constants for the N•L and N•H equation are defined. This speeds up processing. Because this program is calculation intensive, start with a small volume and increase it later when you have the effect you like.

A grayscale palette is loaded and attached to the window using SetPalette(). The palette is a 128 shade grayscale 'pltt' resource created in ResEdit. An offscreen drawing port is used because we only want to calculate the pixel colors once. The offscreen port is set up using the GWorld calls defined in Volume VI of “Inside Macintosh”.

This brings us to the drawing section. The indexes i and j cycle through all the pixels of the viewing window. The k index is the ray moving back through the voxels. All of the work of deciding the color of a pixel is done in the DoColor() routine. Once a non-zero color is returned we move on to the next pixel in the viewing plane. Inside DoColor() the CalcVolumeData() routine calculates whether or not i, j, k is inside the sphere or not. If it is, we calculate a shade. If not, return a RGBColor of zero.

Once the drawing to the offscreen port is done, we set the current port to the onscreen window and exit Init(). At this point the screen is still blank. When WaitNextEvent() receives an update event we use CopyBits() to copy the offscreen bit map onto the screen.

Stuff To Try

Use a small volume to start off. Start with 64x64x64 voxel set. The following are a couple of other sets of defines to try to get various sized spheres. Remember the bigger the voxel space the more time you have to get coffee while the program runs.

/* 1 */

#define volSize  128
#define halfVolSize64
#define sphere_r 60
#define volumeMag3600.0
#define sqrt3r   103.9230
#define rsqrt6   184.5827
and 

#define volSize  256
#define halfVolSize128
#define sphere_r 120
#define volumeMag14400.0
#define sqrt3r   207.8461
#define rsqrt6   369.1654
The defines:

/* these constants define the Phong shading */
/* ambient reflection coefficient */
#define ambientReflCoef 0.1
/* depth cueing coefficient */
#define depthCueCoef 1.0
/* diffuse reflection coefficient */
#define diffReflCoef 2.0
/* specular reflection coefficient */
#define specReflCoef 3.0
/* first light source intensity */
#define light  0.6
/* second light source intensity */
#define light2   1.2
/* coefficient to approx highlight */
#define highlightCoef11

set the constants for the Phong shading. You can play around with these to change the shading effects in the displayed image. But be careful. The value of shade should fall between 0.0 and 1.0. If shade is greater than one, the color will roll over from white to black, and the image will appear with black blotches in the middle of an area that should be white.

A major performance improvement can be made by replacing the floating-point math with suitable integer arithmetic. Some improvement could be made by calling SetCPixel() from a pointer rather than leaving it to the trap dispatcher. Or, the code for the functions could be inserted into Init() to eliminate the overhead of the function calls.

Drawing in 3D is not hard; it just takes some math know-how and a good computer.

References

Computer Graphics: Principles and Practice, 2nd ed., by J. D. Foley, A. Van Dam, S. K. Feiner, and J. F. Hughes (Addison-Wessley, 1990)

Marc Levoy, “Display of Surfaces from Volume Data,” IEEE Computer Graphics and Applications, May 1988 pp. 29-37

Code Listing
#include <Palettes.h>
#include <SANE.h>
#include <QDOffscreen.h>

/* size of the voxel data */
#define volSize  128
/* half the size of voxel data */
#define halfVolSize64
/* the radius of the sphere */
#define sphere_r 60
/* sphere_r*sphere_r */
#define volumeMag3600.0
/* sqrt(3.0)*sphere_r */
#define sqrt3r   103.9230
/* sphere_r*sqrt(6+2*sqrt(3)) */
#define rsqrt6   184.5827

/* resource numbers for the window, palette and menus */
#define windowRscr 128
#define paletteRscr 128

#define appleID 128
#define appleM 1
#define appleAbout 1

#define fileID 129
#define fileM 2
#define fileQuit 1

#define editID 130
#define editM 3
#define editUndo 1
#define editCut 3
#define editCopy 4
#define editPaste 5
#define editClear 6

#define sleepTicks 30

#define aboutDialog 128

/* these constants define the Phong shading */
/* ambient reflection coefficient */
#define ambientReflCoef 0.1
/* depth cueing coefficient */
#define depthCueCoef 1.0
/* diffuse reflection coefficient */
#define diffReflCoef 5.0
/* specular reflection coefficient */
#define specReflCoef 5.0
/* first light source intensity */
#define light  1.0
/* coefficient to approx highlight */
#define highlightCoef30

char aChar;
WindowPtr currentWindow;
MenuHandle myMenus[editM+1];
Rect dragRect, growRect;
long newSize;
Boolean doneFlag;
EventRecord event;
WindowPtr whichWindow;
RGBColor pixColor;
short i, j, k;
PaletteHandle palH;
DialogPtr dPtr;
short doneDlg;
OSErr err;
SysEnvRec envRec;

Rect copyRect;
GWorldPtr wallyWorld;
GDHandle savedDevice;
CGrafPtr savedPort;
 
double PowerOfN (double x, short r) {
 double ans;
 
 ans = 1.0;
 while (r- > 0) ans *= x;
 return ans;
}

double fx, fy, fz;

short CalcVolumeData (short i, short j, short k) {
 long x, y, z;
 
 fx = -(double)(i - halfVolSize);
 fy = -(double)(j - halfVolSize);
 fz = -(double)(k - halfVolSize);
 if ((fx * fx + fy * fy + fz * fz) <= volumeMag) 
 return 1;
 else return 0;
}

void DoColor (short i, short j, short k,
 RGBColor *RGBVal) {
 double n_dot_h, n_dot_l;
 double n_dot_h2, n_dot_l2, shade;
 unsigned short color;
 
 if (CalcVolumeData (i, j, k)) {
 n_dot_l = (fx + fy + fz)/sqrt3r;
 n_dot_h = (fx + fy + 2.7321*fz)/rsqrt6;
 shade = light*ambientReflCoef+
 (light/((double)(k)+depthCueCoef)
 *(diffReflCoef*n_dot_l+specReflCoef
 *PowerOfN (n_dot_h, highlightCoef)));

 /* second light source */
 n_dot_l2 = -fx/sphere_r;
 n_dot_h2 = (-fx + fz)/(1.4142*sphere_r);
 shade +=  light/((double)(k)+depthCueCoef)
 *(diffReflCoef*n_dot_l2+specReflCoef
 *PowerOfN (n_dot_h2, highlightCoef));

 color = (unsigned short)(shade * 65534.0);
 

RGBVal->red = color;
 RGBVal->green = color;
 RGBVal->blue = color;
 }
 else {
 RGBVal->red = 0;
 RGBVal->green = 0;
 RGBVal->blue = 0;
 }
}

void OpenWindow (void) {
 currentWindow =  (WindowPtr)GetNewCWindow(
 windowRscr, NULL, (Ptr)-1);
 SetPort(currentWindow);
 SizeWindow(currentWindow, volSize + 25,
 volSize + 25, 1);
 SetWTitle(currentWindow, &”\pVol3D”);
 ShowWindow(currentWindow);
}

void Init (void) {
 short i, j, k;

 InitGraf(&thePort);
 InitFonts ();
 FlushEvents (everyEvent, 0);
 InitWindows ();
 InitMenus ();
 TEInit ();
 InitDialogs (NULL);

 myMenus[appleM] = GetMenu(appleID);
 AddResMenu(myMenus[appleM], ‘DRVR’);

 myMenus[fileM] = GetMenu(fileID);
 myMenus[editM] = GetMenu(editID);

 for (i=appleM;i<=editM;i++)
 InsertMenu(myMenus[i], 0);

 DrawMenuBar ();

 SetRect(&dragRect, 30, 20,
 screenBits.bounds.right - 10,
 screenBits.bounds.bottom - 30);
 SetRect(&growRect, 50, 50,
 screenBits.bounds.right - 20,
 screenBits.bounds.bottom - 50);

 doneFlag = 0;
 err = SysEnvirons(1, &envRec);
 if (!envRec.hasColorQD) doneFlag = 1;
 else {
 OpenWindow ();
 palH = GetNewPalette (paletteRscr);
 if (palH == NULL) {
 doneFlag = 1;
 }
 else {
 SetPalette (currentWindow, palH, 1);
 }
 
 /* set up the offscreen drawing port */
 GetGWorld (&savedPort, &savedDevice);
 SetRect (&copyRect, 0, 0, volSize-1,
 volSize-1);
 LocalToGlobal (&copyRect.top);
 LocalToGlobal (&copyRect.bottom);
 err = NewGWorld (&wallyWorld, 0, &copyRect,
 NULL, NULL, 0);
 GlobalToLocal (&copyRect.top);
 GlobalToLocal (&copyRect.bottom);

 if (err != noErr)
 doneFlag = 1;
 else {
 SetGWorld (wallyWorld, NULL);
 if (LockPixels (wallyWorld->portPixMap)) {
 /* draw off screen here */
 for(i=0;i<volSize;i++) 
 for (j=0;j<volSize;j++) {
 k = 0;
 do {
 DoColor(i, j, k, &pixColor);
 k++;
 } while ((pixColor.red == 0)
 & (k < volSize));
 SetCPixel (i, j, &pixColor);
 }
 UnlockPixels (wallyWorld->portPixMap);
 }
 else doneFlag = 1;
 
/* the drawing is done set the current port back to the display window 
*/
 }
 SetGWorld (savedPort, savedDevice);
 }
}

void DoAboutBox (void) {

 dPtr = GetNewDialog (aboutDialog, NULL,
 (Ptr)-1);
 do
 ModalDialog(NULL, &doneDlg);
 while (!doneDlg);
 DisposDialog(dPtr);
}

void CleanUp (void) {
 
 HideWindow (currentWindow);
 DisposeGWorld (wallyWorld);
 DisposePalette (palH);
 DisposeWindow (currentWindow);
 doneFlag = 1;
}

void DoCommand (long menuResult) {
 short menuID, menuItem;
 Str255 daName;
 short daErr;

 menuItem = LoWord (menuResult);
 menuID = HiWord (menuResult);

 switch (menuID) {
 case appleID: 
 if (menuItem == appleAbout) DoAboutBox ();
 else {
 GetItem(myMenus[appleM], menuItem, daName);
 daErr = OpenDeskAcc(daName);
 if (currentWindow)
 SetPort (currentWindow);
 }
 break;
 case fileID: 
 switch (menuItem) { 
 case fileQuit: 
 CleanUp ();
 break;
 }
 break;
 }
 HiliteMenu(0);
}

void DoEvent (void) {

 switch (event.what) {
 case mouseDown: 
 switch (FindWindow(event.where, 
 &whichWindow)) {
 case inMenuBar: 
 DoCommand(MenuSelect(event.where));
 break;
 case inSysWindow: 
 SystemClick(&event, whichWindow);
 break;
 case inDrag: 
 DragWindow(whichWindow, event.where,
 &dragRect);
 break;
 case inGrow: 
 newSize = GrowWindow(whichWindow,
 event.where, &growRect);
 SizeWindow(whichWindow, LoWord(newSize),
 HiWord(newSize), 1);
 InvalRect(&currentWindow->portRect);
 break;
 case inGoAway: 
 if (TrackGoAway(whichWindow,
 event.where)) CleanUp ();
 break;
 } /* case findwindow (...) */
 break;
 case keyDown:
 case autoKey: 
 aChar = (char)(BitAnd (event.message,
 charCodeMask));
 if (BitAnd (event.modifiers, cmdKey))
 DoCommand(MenuKey(aChar));
 break;
 case activateEvt: 
 if (BitAnd(event.modifiers, activeFlag))
 DisableItem(myMenus[editM], 0);
 else EnableItem(myMenus[editM], 0);
 break;

 case updateEvt: 
 BeginUpdate(currentWindow);
 EraseRect(&currentWindow->portRect);
 DrawGrowIcon(currentWindow);
 InsetRect (&currentWindow->portRect, 8, 8);
 OffsetRect (&currentWindow->portRect,
 -8, -8);

 if (LockPixels (wallyWorld->portPixMap)) {
 CopyBits(&wallyWorld->portPixMap,
 &currentWindow->portBits, &copyRect,
 &currentWindow->portRect, srcCopy, NULL);
 UnlockPixels (wallyWorld->portPixMap);
 }
 
 OffsetRect (&currentWindow->portRect, 8, 8);
 InsetRect (&currentWindow->portRect, -8, -8);
 EndUpdate(currentWindow);
 break;
 }
}

void main (void) {
 currentWindow = NULL;
 Init ();
 InitCursor ();

 do {
 if (WaitNextEvent (everyEvent, &event,
 sleepTicks, NULL)) DoEvent ();
 } while (!doneFlag);
}