A PIXEL ORDERING ALGORITHM
A shortcut for interactive development
This article contains the following executables> ALLEN.EXE
Norton T. Allen
Norton is a systems analyst/programmer with the Harvard University Atmospheric Research Project and can be reached at 202 Ridge St., Winchester, MA 01890.
Personal computers can generate and display some very interesting graphics. Unfortunately, some of the most interesting programs require hours to produce a single graphics screen. Mandelbrot sets and ray-tracing programs do extensive floating-point calculations for each pixel of the final display. The iterations can be painfully slow if the resulting image needs touching up. Often, a general idea of how the final image will look is enough to act on, but the standard linear approach to drawing the screen can mean a long wait just to find out what's going on in the lower right-hand corner.
For programs that produce graphics one pixel at a time, a better method of pixel ordering provides a shortcut during interactive development. By selecting pixels that are always evenly distributed on the screen, the general character of the final image can be seen early on. Subsequent calculations continuously improve the focus until full resolution is achieved. This makes it possible to jump in early and adjust colors, move objects, or zoom in on a particular area of interest. I have implemented this algorithm for viewing the Mandelbrot set, but it also applies as well to ray-tracing or any other pixel-by-pixel task.
Reversing Bit Order
Part of my job is writing software to receive data transmitted from high-altitude research balloons. Information from many different sensors is collected on the balloon, sorted into a standard data structure, and then broadcast as a serial bit stream. On the ground, the process is reversed. Bits are collected into bytes, and the bytes become members of structures. To help sort things out, we use a frame counter in the data stream. This is a 1- or 2-byte integer that increments every time it is transmitted. Once we successfully locked onto a data stream, but the frame counter readout showed no apparent pattern. By changing the display to binary, the problem became obvious -- the normal bit order (MSB - >LSB) was reversed (LSB - >MSB in our case).
In that case, bit-reversal was the root of the problem. When I began thinking about generating pixel-oriented graphics, however, bit-reversal was the key to a solution. When counting the binary integers from 0002 to 1112, the most significant bit on the left is not set until the halfway point at 1002. If the bits are reversed, the least significant bit is not set until halfway through, and all bit-reversed numbers up to that point are even. As Figure 1 shows, these bit-reversed numbers can be used to address a row of pixels. Only even pixels are addressed until Z=1002. At each step, the distribution of filled and unfilled pixels is as uniform as possible.
This is fine for addressing a single row of pixels, but both rows and columns need to be addressed. Two independent counters in a nested loop could be used, but this would not give the desired result. The inside loop would completely fill one row before the next evenly-spaced row was selected.
The following algorithm combines both indices/coordinates (X, Y) in a single counter, Z, by interleaving the bits of their binary representations in reverse order. The variable map (shown in Figure 2) keeps track of which bits correspond to which coordinate (X or Y). Bits in Z, which correspond to zero bits in the map are part of the X coordinate. Bits corresponding to one are part of the Y coordinate. The results are demonstrated in a 4 x 4 region in Figure 2.
This is close to the desired result, but the distribution at the halfway point (shown in the second pattern in Figure 2) could be more uniform. Half the rows are completely filled in, while the other half haven't been filled in at all. This problem can be seen in a 2 x 2 example. The algorithm described produces X, Y pairs in the order (0,0), (1,0), (0, 1), (1, 1), which can be represented:
--------
| 0 1 |
| 2 3 |
--------
A more desirable X, Y pair ordering is (0, 0), (1, 1), (0, 1), (1, 0), or:
--------
| 0 3 |
| 2 1 |
--------
The X coordinates for both orders are the same, but the Y coordinates are different. The following inset shows a truth table comparing the original order (Xz, Yz), and the new order (Xz, Y'z).
| X<SUB>z</SUB>
Y<SUB>'z</SUB>| 0 1
------------
Y<SUB>z</SUB> 0 | 0 1
1 | 1 0
The truth table shows that Y'z = Xz 
Yz. designates the exclusive-OR operation. Once Xz and Yz are calculated, it is easy to generate Y'z, and report the new point (Xz, Y'z). Figure 3 shows how this works in a 4 x 4 case.
There are a few rough edges dealing with arbitrary rectangular regions. First, the ranges may not be even powers of two. This is handled by rounding up to the next even power of two and then discarding any points outside of the rectangular region. Second, even after rounding, the regions may not be square. One coordinate may have more bits than the other. This is handled by storing the extra bits on the LSB end of map. If X has 2 bits and Y has 4 bits, map = 1010112. The exclusive-OR operation imposes an additional complication. The Y range must be at least as large as the X range. This is handled by reversing the roles of the two coordinates, if the X range is larger.
Compilation
I have implemented a Mandelbrot set example using grafix.lib--Kent Porter's graphics library (see DDJ, February 1989 through August 1989). This library is compatible with Microsoft C 5.1 and Turbo C 2.0. I compiled this project using Lattice C, Version 3.41, and some minor modifications were necessary. Lattice did not support mixed mode programming before Version 6.0. Grafix.lib contains explicit "far" designations. These designations are redundant when programming in the large model, and problematic when programming in the small model. Because my program is small, I removed the explicit declarations. I also ran into syntax difficulties using Kent's #define byte unsigned char, which vanished when I used typedef unsigned char byte.
On the assembly side, I modified the initializations to use the Lattice assembly macros to provide model-sensitive procedure definitions. I also followed Lattice's recommendation and saved ES when that register was used. Finally, Lattice function and data names are not prefaced with an underscore in assembly modules.
The Code
The file points.c (Listing One, page 116) contains two routines: init_points( ) and next_point( ). The X and Y limits are passed to init_points( ), which initializes the map variable. next_point( ) defines the coordinates of each successive pixel through pointer arguments, returning a non-zero value when the area has been completely filled.
Mandel.c (Listing Two, page 116) provides functions particular to the Mandelbrot example. The function set_range( ) creates a virtual coordinate system. This function is similar to grafix' setcoords( ), with two major differences: The function set_range( ) guarantees that the new coordinates preserve the aspect ratio of the screen. One inch on the screen should represent the same number of units in either the X or the Y direction. If one range is proportionately smaller than the other, the coordinates are extended to place that range in the center of the screen. The function setcoords( ) only provides functions to convert from virtual coordinates to device coordinates, but in this program we need to do the reverse. mandel_color( ) contains the Mandelbrot calculations to determine the color of a particular pixel.
Generate.c (Listing Three, page 116) provides the operational framework for the demonstration program and point.h (Listing Four, page 118) is the include file. generate( ) calls next_point( ) to step through the pixels, and calls mandel_color( ) to determine their color. menu( ) is called any time a key is pressed. menu( ) lets the user move and scale a rectangular box on the screen to identify a region for closer examination. The area in the rectangular box can be focused or zoomed on. Focusing confines subsequent calculations to the bounds of the box. This gives a sharper look at a small area, in the context of the larger screen. Zooming in rescales the virtual coordinates of the screen to the coordinates of the box.
This example demonstrates the advantages of this algorithm. A small region can be zoomed in without waiting for the completion of intermediate views. By displaying the full scale image at low resolution, the zooming can be accurately selected. Focusing can be used to choose from different regions. The beauty of this algorithm is that the same amount of detail is provided in the same amount of time, regardless of whether the image is zoomed or focused. If you zoom, a longer wait will produce even more detail.
XOR
To overlay a movable box on the graphics display, a significant addition was made to the library. I added the ability to draw pixels on the screen by calculating the exclusive-OR (XOR) of the current foreground color and the color already on the screen. XOR is very useful for three reasons. If you XOR a bit pattern with a pattern of ones you get the ones complement of the original bit pattern. Whether writing to a black or a white background, a change is visible. XORing a bit pattern with zeros has no effect. If a bit pattern is XORed with itself, the result is zero. XOR is associative. If bit pattern X is XORed with non-zero mask M, a result will be seen. If the result is XORed with M again, the result is (X M) M. This is equal to X (M M), which is equal to X 0, which is equal to X.
XOR allows an image to be put on and removed from the screen quickly without destroying the image being written over. This is very useful for interactive applications. Not only does it allow the cursor to move around a graphics screen, but it also enables experimentation with the placement of objects before committing them to an image.
This is so useful that it is implemented in the hardware of the EGA adapter. The write function register can be set to replace the stored value, or to perform OR, AND, or XOR functions between the stored value and the value being written. In this case, I want to move a rectangular box. The function to draw the box is draw_rect( ). This function calls draw_line( ), which calls draw_point( ). Adding a XOR mode to draw_point( ) allows the box to be moved. Listing Five, page 118, contains definitions to be included in grafix.h, and Listing Six page 118, contains wmode.c, which defines the function set_write_mode( ). This module should be compiled and included in grafix.lib by using the LIB grafix +wmode; command.
Two changes need to be made in drawpt.asm. Below line 14, where the vuport pointer is declared, add the line: EXTRN _write_mode: BYTE. Also change line 78 from: xor ah, ah to: mov ah, _write_mode. Assemble drawpt, then replace it in the library with the command: LIB grafix -+drawpt. A similar change can be made to hline.asm by adding the declaration after line 14, and changing line 106.
0ther Applications
To use this algorithm in a ray-tracing study, replace the functions in mandel.c with the ray-tracing code and change the calls in generate( ) accordingly. Another application of this algorithm is for demo programs, where one graphics image fades into another. Running the fade one pixel at a time is much too slow. Divide the image into small rectangles, use this algorithm to index the rectangles, and use a bitblt routine to copy the rectangles.
The pixel ordering presented here is closely related to dither matrices used for generating gray shading on monochrome monitors. For 16 shades of gray, the screen is tiled with imaginary 4 x 4 pixel blocks. Each pixel is numbered from 0 to 15, according to its position in the order. When the user wishes to color an area with a certain gray level, only those pixels with numbers less than that level are turned on.
Tradeoffs
This algorithm does much bit twiddling to calculate each set of coordinates. This adds CPU processing time. On the 4.77-MHz XT, this algorithm took 15 minutes longer to fill the screen than a simple nested loop. If the exact calculations for drawing an image are already known, it is preferable not to use this algorithm. When a lot of CPU time is being used, however, it's nice to know if the right picture is being generated. Seeing only the left quarter of the screen won't help, but one quarter resolution across the entire display is enough to see considerable detail.
_A PIXEL ORDERING ALGORITHM_ by Norton T. Allen
[LISTING ONE]<a name="0141_000b">
/* points.c Copyright (c) 1989 by Norton T. Allen */
#include "points.h"
int nbits(int i) {
int nb;
for (nb = 0; i != 0; nb++) i >>= 1;
return(nb);
}
static int focus_x, focus_y, focus_width, focus_height;
/* zeros in map correspond to x bits, 1's to y bits */
static long int z, map, zlim;
static int tbits, xor_y;
void init_points(int left, int top, int width, int height) {
int nbits_w, nbits_h;
long int mask;
focus_x = left;
focus_y = top;
focus_width = width;
focus_height = height;
nbits_w = nbits(focus_width);
nbits_h = nbits(focus_height);
tbits = nbits_w + nbits_h;
xor_y = nbits_w < nbits_h;
z = 0L;
zlim = 1L << tbits;
map = 0L;
for (mask = 1L; nbits_w || nbits_h; ) {
if (nbits_w > nbits_h || (nbits_w == nbits_h && xor_y)) {
nbits_w--;
mask <<= 1;
}
if (nbits_h > 0 &&
(nbits_h > nbits_w || (nbits_w == nbits_h && !xor_y))) {
nbits_h--;
map |= mask;
mask <<= 1;
}
}
}
int next_point(int *dx, int *dy) {
int x, y, n;
long int m, rz;
for (;;) {
if (z == zlim) return(1);
m = map; rz = z; n = tbits;
x = y = 0;
while (n-- > 0) {
if (m & 1) { /* this is a ybit */
y <<= 1;
if (rz & 1) y++;
} else {
x <<= 1;
if (rz & 1) x++;
}
m >>= 1; rz >>= 1;
}
z++;
if (xor_y) y = x ^ y;
else x = x ^ y;
if (x >= focus_width) continue;
if (y >= focus_height) continue;
break;
}
*dx = x + focus_x;
*dy = y + focus_y;
return(0);
}
<a name="0141_000c"><a name="0141_000c">
<a name="0141_000d">[LISTING TWO]<a name="0141_000d">
/* mandel.c will handle the actual Mandelbrot set calculations.
Copyright (c) 1989 by Norton T. Allen */
#include "grafix.h"
/* This is the aspect ratio calculated on my screen: */
#define ASPECT 0.739
static double xscale, yscale, xo, yo;
/* set_range() is similar to setcoords() except that it guarantees
equal x and y scales, taking the aspect ratio into account.
The axis with the smaller scale is adjust so the specified
range will be centered on the screen. Another reason for
not using setcoords() is that I need the inverse functions
mapping device coordinates to virtual coordinates.
*/
void set_range(double xmin, double ymin, double xmax, double ymax) {
double delta;
yscale = (ymax-ymin)/vp_height();
xscale = (xmax-xmin)/vp_width();
if (yscale*ASPECT > xscale) {
xscale = yscale*ASPECT;
delta = (xmin + xscale * vp_width() - xmax)/2.;
xmin -= delta;
xmax += delta;
} else {
yscale = xscale/ASPECT;
delta = (ymin + yscale * vp_height() - ymax)/2.;
ymin -= delta;
ymax += delta;
}
xo = xmin;
yo = ymin;
} /* ---------------------------------------------------------------- */
double vx (int dx) /* convert device x to virtual x */
{
return ((double) dx * xscale + xo);
} /* ---------------------------------------------------------------- */
double vy (int dy) /* convert device y to virtual y */
{
return ((double) dy * yscale + yo);
} /* ---------------------------------------------------------------- */
/* Membership in the Mandelbrot set is determined by an
iterative function. Given the complex starting
coordinate C the function is:
Z(0) = 0.
Z(n+1) = Z(n)^2 + C
A point is deemed to be in the set if NLOOP iterations fails
to produce a point with absolute value greater than LIMIT.
Points within the set are colored black. Points outside the
set are colored based on how many iterations passed before
exceeding LIMIT. Miriad other schemes are possible. Other
values for NLOOP and LIMIT are probably desirable, depending
on how deep you go.
*/
#define NLOOP 100
#define LIMIT 10000.
#define NCOLORS 15
int mandel_color(double cx, double cy) {
int k;
double zx, zy, zx2, zy2;
zx = zy = 0.;
for (k = 0; k < NLOOP; k++) {
zx2 = zx*zx;
zy2 = zy*zy;
if (zx2+zy2 > LIMIT) break;
zy = 2*zx*zy + cy;
zx = zx2 - zy2 + cx;
}
if (k < NLOOP) return((k % NCOLORS)+1);
return(0);
}
<a name="0141_000e"><a name="0141_000e">
<a name="0141_000f">[LISTING THREE]<a name="0141_000f">
/* generate.c Copyright (c) 1989 by Norton T. Allen */
#include <stdio.h>
#include <stdlib.h>
#include <dos.h>
#include "points.h"
#include "grafix.h"
/* This structure defines a rectangular box which we can move around the
screen using the cursor keys. fbox.on is TRUE if the box is
currently displayed. fbox.mode takes the values:
0 Not active
1 Cursor keys move the whole box
2 Cursor keys change box's size
*/
struct bx {
int x, y, dx, dy, on, mode;
} fbox = {300, 160, 50, 40, 0, 0};
void flip_box(void) {
set_write_mode(WM_XOR);
set_color1(15);
draw_rect(fbox.x, fbox.y, fbox.dx, fbox.dy);
set_write_mode(WM_REPLACE);
fbox.on = !fbox.on;
}
void help(void) {
printf("'Esc' to exit\n");
printf("? for this message\n");
printf("m to move the focus box.(use cursor keys)\n");
printf("s to size the focus box.(use cursor keys)\n");
printf("f to focus on the focus area\n");
printf("z to zoom in on the focus area\n");
getch();
}
static int dm = 1; /* How many pixels to move with each cursor step */
/* check_box makes sure the new box meets the following requirements:
1. It isn't off the screen.
2. It's not smaller then 2 pixels in either dimension
Check box turns off the old box, but leaves the new box off also;
menu() will turn it on when there's no more keyboard input.
*/
void check_box(struct bx *nb) {
if (nb->x < 0 || nb->y < 0 ||
nb->x + nb->dx > vp_width() ||
nb->y + nb->dy > vp_height() ||
nb->dx < 2 || nb->dy < 2 ||
(nb->x == fbox.x && nb->y == fbox.y &&
nb->dx == fbox.dx && nb->dy == fbox.dy))
return;
if (fbox.on) flip_box();
fbox = *nb;
fbox.on = 0;
}
/* These are scan codes for cursor keys: */
#define EX_UP 72
#define EX_DOWN 80
#define EX_RIGHT 77
#define EX_LEFT 75
/* Menu supports the following keys:
ESC Exit
M suspend calculations and put box on the screen in
mode 1 for 'moving'.
S As with M, but mode 2 for 'sizing'.
C Remove box and continue calculations as before
F Focus on boxed region.
Z Zoom in on boxed region.
Cursor keys Move or size box
0-9 Change step size for cursor keys
*/
void menu(void) {
int c;
struct bx new_box;
double xmin, ymin, xmax, ymax;
for (;;) {
while (kbhit()) {
c = getch();
switch (c) {
case '\033':
exit(0);
case 'c':
case 'C':
fbox.mode = 0;
break;
case 'm':
case 'M':
fbox.mode = 1;
break;
case 's':
case 'S':
fbox.mode = 2;
break;
case 'f': /* Focus on boxed region */
case 'F':
if (fbox.mode == 0) break;
init_points(fbox.x, fbox.y, fbox.dx, fbox.dy);
fbox.mode = 0;
break;
case 'z': /* Zoom in on boxed region */
case 'Z':
if (fbox.mode == 0) break;
xmin = vx(fbox.x);
ymin = vy(fbox.y+fbox.dy);
ymax = vy(fbox.y);
xmax = vx(fbox.x+fbox.dx);
set_range(xmin, ymin, xmax, ymax);
init_points(0, 0, vp_width(), vp_height());
pc_textmode(); /* Kluge to clear screen */
init_video(EGA);
fbox.mode = fbox.on = 0;
break;
case '1':
case '2':
case '3':
case '4':
case '5':
case '6':
case '7':
case '8':
case '9':
dm = c-'0';
break;
case 0:
c = getch();
new_box = fbox;
if (fbox.mode == 1) { /* moving the box */
switch (c) {
case EX_UP:
new_box.y -= dm;
break;
case EX_DOWN:
new_box.y += dm;
break;
case EX_RIGHT:
new_box.x += dm;
break;
case EX_LEFT:
new_box.x -= dm;
break;
default: break;
}
} else if (fbox.mode == 2) {
switch (c) {
case EX_UP:
new_box.dy -= dm;
break;
case EX_DOWN:
new_box.dy += dm;
break;
case EX_RIGHT:
new_box.dx += dm;
break;
case EX_LEFT:
new_box.dx -= dm;
break;
default: break;
}
}
check_box(&new_box);
break;
default:
break;
}
}
if (fbox.mode == 0) break;
else if (!fbox.on) flip_box();
}
if (fbox.on) flip_box();
}
/* Generate cycles through all the pixels, possibly starting over as
dictated by menu(). Generate never returns.
*/
void generate(void) {
int c, x, y;
double fx, fy;
for (;;) {
while (next_point(&x, &y) == 0) {
fx = vx(x);
fy = vy(y);
c = mandel_color(fx, fy);
set_color1(c);
draw_point(x, y);
if (kbhit()) menu();
}
menu();
}
}
void main(int argc, char **argv) {
if (init_video (EGA)) {
init_points(0, 0, vp_width(), vp_height());
set_range(-2., -.95, .75, .95);
generate();
} else printf("Cannot select graphics mode");
}
<a name="0141_0010"><a name="0141_0010">
<a name="0141_0011">[LISTING FOUR]<a name="0141_0011">
/* points.h include file for pixel ordering program. */
void init_points(int left, int top, int width, int height);
int next_point(int *dx, int *dy);
int mandel_color(double cx, double cy);
void set_range(double xmin, double ymin, double xmax, double ymax);
double vx (int dx); /* convert device x to virtual x */
double vy (int dy); /* convert device y to virtual y */
<a name="0141_0012"><a name="0141_0012">
<a name="0141_0013">
[LISTING FIVE]<a name="0141_0013">
To be included in grafix.h
/* Added by Norton Allen */
/* --------------------- */
void set_write_mode(int mode);
#define WM_REPLACE 0
#define WM_XOR 0x18
<a name="0141_0014"><a name="0141_0014">
<a name="0141_0015">
[LISTING SIX]<a name="0141_0015">
/* wmode.c for DDJ grafix library Copyright (c) 1989 by Norton T. Allen */
#include <dos.h>
#include "grafix.h"
int write_mode = WM_REPLACE;
void set_write_mode(int mode) {
write_mode = mode;
}
