Voxels
Volume Number: 9
Issue Number: 3
Column Tag: C Workshop
What are Voxels?
A hobbyist level overview of voxels, voxel space and Phong shading
By Geoffrey Clements, Chelmsford, Massachusetts
Note: Source code files accompanying article are located on MacTech CD-ROM orsource code disks.
About the author
Geoffrey Clements can be reached on internet at clementsg@gw1.hanscom.af.mil.
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 halfVolSize 32
/* the radius of the sphere */
#define sphere_r 30
/* sphere_r*sphere_r */
#define volumeMag 900.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 halfVolSize 64
#define sphere_r 60
#define volumeMag 3600.0
#define sqrt3r 103.9230
#define rsqrt6 184.5827
and
#define volSize 256
#define halfVolSize 128
#define sphere_r 120
#define volumeMag 14400.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 highlightCoef 11
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
#include
#include
/* size of the voxel data */
#define volSize 128
/* half the size of voxel data */
#define halfVolSize 64
/* the radius of the sphere */
#define sphere_r 60
/* sphere_r*sphere_r */
#define volumeMag 3600.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 highlightCoef 30
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 (©Rect, 0, 0, volSize-1,
volSize-1);
LocalToGlobal (©Rect.top);
LocalToGlobal (©Rect.bottom);
err = NewGWorld (&wallyWorld, 0, ©Rect,
NULL, NULL, 0);
GlobalToLocal (©Rect.top);
GlobalToLocal (©Rect.bottom);
if (err != noErr)
doneFlag = 1;
else {
SetGWorld (wallyWorld, NULL);
if (LockPixels (wallyWorld->portPixMap)) {
/* draw off screen here */
for(i=0;i
for (j=0;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, ©Rect,
& 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);
}