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 or source code disks.

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.

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 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

/* 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;
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;
unsigned short color;

if (CalcVolumeData (i, j, k)) {
n_dot_l = (fx + fy + fz)/sqrt3r;
n_dot_h = (fx + fy + 2.7321*fz)/rsqrt6;
(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);
*(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 ();
TEInit ();
InitDialogs (NULL);

for (i=appleM;i<=editM;i++)

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);
}
}

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

void CleanUp (void) {

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

Str255 daName;
short daErr;

case appleID:
else {
daErr = OpenDeskAcc(daName);
if (currentWindow)
SetPort (currentWindow);
}
break;
case fileID:
case fileQuit:
CleanUp ();
break;
}
break;
}
}

void DoEvent (void) {

switch (event.what) {
case mouseDown:
switch (FindWindow(event.where,
&whichWindow)) {
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,
if (BitAnd (event.modifiers, cmdKey))
break;
case activateEvt:
if (BitAnd(event.modifiers, activeFlag))
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);
}
```

Community Search:
MacTech Search:

Latest Forum Discussions

Fresh From the Land Down Under – The Tou...
After a two week hiatus, we are back with another episode of The TouchArcade Show. Eli is fresh off his trip to Australia, which according to him is very similar to America but more upside down. Also kangaroos all over. Other topics this week... | Read more »
TouchArcade Game of the Week: ‘Dungeon T...
I’m a little conflicted on this week’s pick. Pretty much everyone knows the legend of Dungeon Raid, the match-3 RPG hybrid that took the world by storm way back in 2011. Everyone at the time was obsessed with it, but for whatever reason the... | Read more »
Hello gentle readers, and welcome to the SwitchArcade Round-Up for July 19th, 2024. In today’s article, we finish up the week with the unusual appearance of a review. I’ve spent my time with Hot Lap Racing, and I’m ready to give my verdict. After... | Read more »
Draknek Interview: Alan Hazelden on Thin...
Ever since I played my first release from Draknek & Friends years ago, I knew I wanted to sit down with Alan Hazelden and chat about the team, puzzle games, and much more. | Read more »
The Latest ‘Marvel Snap’ OTA Update Buff...
I don’t know about all of you, my fellow Marvel Snap (Free) players, but these days when I see a balance update I find myself clenching my… teeth and bracing for the impact to my decks. They’ve been pretty spicy of late, after all. How will the... | Read more »
‘Honkai Star Rail’ Version 2.4 “Finest D...
HoYoverse just announced the Honkai Star Rail (Free) version 2.4 “Finest Duel Under the Pristine Blue" update alongside a surprising collaboration. Honkai Star Rail 2.4 follows the 2.3 “Farewell, Penacony" update. Read about that here. | Read more »
‘Vampire Survivors+’ on Apple Arcade Wil...
Earlier this month, Apple revealed that poncle’s excellent Vampire Survivors+ () would be heading to Apple Arcade as a new App Store Great. I reached out to poncle to check in on the DLC for Vampire Survivors+ because only the first two DLCs were... | Read more »
Homerun Clash 2: Legends Derby opens for...
Since launching in 2018, Homerun Clash has performed admirably for HAEGIN, racking up 12 million players all eager to prove they could be the next baseball champions. Well, the title will soon be up for grabs again, as Homerun Clash 2: Legends... | Read more »
‘Neverness to Everness’ Is a Free To Pla...
Perfect World Games and Hotta Studio (Tower of Fantasy) announced a new free to play open world RPG in the form of Neverness to Everness a few days ago (via Gematsu). Neverness to Everness has an urban setting, and the two reveal trailers for it... | Read more »
Meditative Puzzler ‘Ouros’ Coming to iOS...
Ouros is a mediative puzzle game from developer Michael Kamm that launched on PC just a couple of months back, and today it has been revealed that the title is now heading to iOS and Android devices next month. Which is good news I say because this... | Read more »

Price Scanner via MacPrices.net

Amazon is still selling 16-inch MacBook Pros...
Prime Day in July is over, but Amazon is still selling 16-inch Apple MacBook Pros for \$500-\$600 off MSRP. Shipping is free. These are the lowest prices available this weekend for new 16″ Apple... Read more
Walmart continues to sell clearance 13-inch M...
Walmart continues to offer clearance, but new, Apple 13″ M1 MacBook Airs (8GB RAM, 256GB SSD) online for \$699, \$300 off original MSRP, in Space Gray, Silver, and Gold colors. These are new MacBooks... Read more
Apple is offering steep discounts, up to \$600...
Apple has standard-configuration 16″ M3 Max MacBook Pros available, Certified Refurbished, starting at \$2969 and ranging up to \$600 off MSRP. Each model features a new outer case, shipping is free,... Read more
Save up to \$480 with these 14-inch M3 Pro/M3...
Apple has 14″ M3 Pro and M3 Max MacBook Pros in stock today and available, Certified Refurbished, starting at \$1699 and ranging up to \$480 off MSRP. Each model features a new outer case, shipping is... Read more
Amazon has clearance 9th-generation WiFi iPad...
Amazon has Apple’s 9th generation 10.2″ WiFi iPads on sale for \$80-\$100 off MSRP, starting only \$249. Their prices are the lowest available for new iPads anywhere: – 10″ 64GB WiFi iPad (Space Gray or... Read more
Apple is offering a \$50 discount on 2nd-gener...
Apple has Certified Refurbished White and Midnight HomePods available for \$249, Certified Refurbished. That’s \$50 off MSRP and the lowest price currently available for a full-size Apple HomePod today... Read more
The latest MacBook Pro sale at Amazon: 16-inc...
Amazon is offering instant discounts on 16″ M3 Pro and 16″ M3 Max MacBook Pros ranging up to \$400 off MSRP as part of their early July 4th sale. Shipping is free. These are the lowest prices... Read more
14-inch M3 Pro MacBook Pros with 36GB of RAM...
B&H Photo has 14″ M3 Pro MacBook Pros with 36GB of RAM and 512GB or 1TB SSDs in stock today and on sale for \$200 off Apple’s MSRP, each including free 1-2 day shipping: – 14″ M3 Pro MacBook Pro (... Read more
14-inch M3 MacBook Pros with 16GB of RAM on s...
B&H Photo has 14″ M3 MacBook Pros with 16GB of RAM and 512GB or 1TB SSDs in stock today and on sale for \$150-\$200 off Apple’s MSRP, each including free 1-2 day shipping: – 14″ M3 MacBook Pro (... Read more
Amazon is offering \$170-\$200 discounts on new...
Amazon is offering a \$170-\$200 discount on every configuration and color of Apple’s M3-powered 15″ MacBook Airs. Prices start at \$1129 for models with 8GB of RAM and 256GB of storage: – 15″ M3... Read more

Jobs Board

*Apple* Systems Engineer - Chenega Corporati...
…LLC,** a **Chenega Professional Services** ' company, is looking for a ** Apple Systems Engineer** to support the Information Technology Operations and Maintenance Read more
Solutions Engineer - *Apple* - SHI (United...
**Job Summary** An Apple Solution Engineer's primary role is tosupport SHI customers in their efforts to select, deploy, and manage Apple operating systems and Read more
*Apple* / Mac Administrator - JAMF Pro - Ame...
Amentum is seeking an ** Apple / Mac Administrator - JAMF Pro** to provide support with the Apple Ecosystem to include hardware and software to join our team and Read more
Operations Associate - *Apple* Blossom Mall...
Operations Associate - Apple Blossom Mall Location:Winchester, VA, United States (https://jobs.jcp.com/jobs/location/191170/winchester-va-united-states) - Apple Read more
Cashier - *Apple* Blossom Mall - JCPenney (...
Cashier - Apple Blossom Mall Location:Winchester, VA, United States (https://jobs.jcp.com/jobs/location/191170/winchester-va-united-states) - Apple Blossom Mall Read more