|
Dithering with just one neighbor doesn't sound very helpful - unless you're clever about how you visit the neighbors.
December 01, 1998
URL:http://www.drdobbs.com/a-balanced-dithering-technique/184403590
/* Riemersma dither
*
* This program reads in an uncompressed
* gray-scale image with one byte per
* pixel and a size of 256*256 pixels (no
* image header). It dithers the image and
* writes an output image in the same
* format.
*
* This program was tested with Borland C++
* 3.1 16-bit (DOS), compiled in large
* memory model. For other compilers, you
* may have to replace the calls to
* _fmalloc() and _ffree() with straight
* malloc() and free() calls.
*/
#include <malloc.h> /* for _fmalloc() and
_ffree() */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define WIDTH 256
#define HEIGHT 256
#define BLOCKSIZE 16384 // read in chunks
// of 16 kBytes
enum {
NONE,
UP,
LEFT,
DOWN,
RIGHT,
};
/* variables needed for the Riemersma
* dither algorithm */
static int cur_x=0, cur_y=0;
static int img_width=0, img_height=0;
static unsigned char *img_ptr;
#define SIZE 16 /* queue size: number of
* pixels remembered */
#define MAX 16 /* relative weight of
* youngest pixel in the
* queue, versus the oldest
* pixel */
static int weights[SIZE]; /* weights for
* the errors
* of recent
* pixels */
static void
init_weights(int a[],int size,int max)
{
double m = exp(log(max)/(size-1));
double v;
int i;
for (i=0, v=1.0; i<size; i++) {
a[i]=(int)(v+0.5); /* store rounded
* value */
v*=m; /* next value */
} /*for */
}
static void
dither_pixel(unsigned char *pixel)
{
static int error[SIZE]; /* queue with error
* values of recent
* pixels */
int i,pvalue,err;
for (i=0,err=0L; i<SIZE; i++)
err+=error[i]*weights[i];
pvalue=*pixel + err/MAX;
pvalue = (pvalue>=128) ? 255 : 0;
/* shift queue */
memmove(error, error+1,
(SIZE-1)*sizeof error[0]);
error[SIZE-1] = *pixel - pvalue;
*pixel=(unsigned char)pvalue;
}
static void move(int direction)
{
/* dither the current pixel */
if (cur_x >= 0 && cur_x < img_width &&
cur_y >= 0 && cur_y < img_height)
dither_pixel(img_ptr);
/* move to the next pixel */
switch (direction) {
case LEFT:
cur_x--;
img_ptr--;
break;
case RIGHT:
cur_x++;
img_ptr++;
break;
case UP:
cur_y--;
img_ptr-=img_width;
break;
case DOWN:
cur_y++;
img_ptr+=img_width;
break;
} /* switch */
}
void hilbert_level(int level,int direction)
{
if (level==1) {
switch (direction) {
case LEFT:
move(RIGHT);
move(DOWN);
move(LEFT);
break;
case RIGHT:
move(LEFT);
move(UP);
move(RIGHT);
break;
case UP:
move(DOWN);
move(RIGHT);
move(UP);
break;
case DOWN:
move(UP);
move(LEFT);
move(DOWN);
break;
} /* switch */
} else {
switch (direction) {
case LEFT:
hilbert_level(level-1,UP);
move(RIGHT);
hilbert_level(level-1,LEFT);
move(DOWN);
hilbert_level(level-1,LEFT);
move(LEFT);
hilbert_level(level-1,DOWN);
break;
case RIGHT:
hilbert_level(level-1,DOWN);
move(LEFT);
hilbert_level(level-1,RIGHT);
move(UP);
hilbert_level(level-1,RIGHT);
move(RIGHT);
hilbert_level(level-1,UP);
break;
case UP:
hilbert_level(level-1,LEFT);
move(DOWN);
hilbert_level(level-1,UP);
move(RIGHT);
hilbert_level(level-1,UP);
move(UP);
hilbert_level(level-1,RIGHT);
break;
case DOWN:
hilbert_level(level-1,RIGHT);
move(UP);
hilbert_level(level-1,DOWN);
move(LEFT);
hilbert_level(level-1,DOWN);
move(DOWN);
hilbert_level(level-1,LEFT);
break;
} /* switch */
} /* if */
}
int log2(int value)
{
int result=0;
while (value>1) {
value >>= 1;
result++;
} /*while */
return result;
}
void
Riemersma(void *image, int width,
int height)
{
int level,size;
/* determine the required order of the
* Hilbert curve */
size=max(width,height);
level=log2(size);
if ((1L << level) < size)
level++;
init_weights(weights,SIZE,MAX);
img_ptr=image;
img_width=width;
img_height=height;
cur_x=0;
cur_y=0;
if (level>0)
hilbert_level(level,UP);
move(NONE);
}
void main(int argc,char *argv[])
{
unsigned char *ptr;
unsigned char *image;
FILE *fp;
long filesize;
char filename[128];
/* check arguments */
if (argc<2 || argc>3) {
printf(
"Usage: riemer <input> [output]\n\n"
"Input and output files must be raw "
"gray-scale files with a size of "
"256*256\npixels; one byte per
"pixel.\n");
return;
} /* if */
if ((fp=fopen(argv[1],"rb"))==NULL) {
printf("File not found (%s)\n",
argv[1]);
return;
} /* if */
if ((image=_fmalloc((long)WIDTH*HEIGHT))
== NULL) {
printf("Insufficient memory\n");
return;
} /* if */
/* read in the file */
filesize=(long)WIDTH*HEIGHT;
ptr=image;
while (filesize>BLOCKSIZE) {
fread(ptr,1,BLOCKSIZE,fp);
filesize-=BLOCKSIZE;
ptr+=BLOCKSIZE;
} /* while */
if (filesize>0)
fread(ptr,1,(int)filesize,fp);
fclose(fp);
/* dither (replacing the original */
Riemersma(image,WIDTH,HEIGHT);
/* save the dithered file */
if (argc==3)
strcpy(filename,argv[2]);
else
strcpy(filename,"riemer.raw");
fp=fopen(filename,"wb");
if (fp==NULL) {
printf("Cannot create file \"%s\"\n",
filename);
return;
} /* if */
filesize=(long)WIDTH*HEIGHT;
ptr=image;
while (filesize>BLOCKSIZE) {
fwrite(ptr,1,BLOCKSIZE,fp);
filesize-=BLOCKSIZE;
ptr+=BLOCKSIZE;
} /* while */
if (filesize>0)
fwrite(ptr,1,(int)filesize,fp);
fclose(fp);
_ffree(image);
}
/* End of File */
enum {
UP,
LEFT,
DOWN,
RIGHT,
};
void hilbert(int level,int direction=UP)
{
if (level==1) {
switch (direction) {
/* move() could draw a line in... */
/* ...the indicated direction */
case LEFT:
move(RIGHT);
move(DOWN);
move(LEFT);
break;
case RIGHT:
move(LEFT);
move(UP);
move(RIGHT);
break;
case UP:
move(DOWN);
move(RIGHT);
move(UP);
break;
case DOWN:
move(UP);
move(LEFT);
move(DOWN);
break;
} /* switch */
} else {
switch (direction) {
case LEFT:
hilbert_level(level-1,UP);
move(RIGHT);
hilbert_level(level-1,LEFT);
move(DOWN);
hilbert_level(level-1,LEFT);
move(LEFT);
hilbert_level(level-1,DOWN);
break;
case RIGHT:
hilbert_level(level-1,DOWN);
move(LEFT);
hilbert_level(level-1,RIGHT);
move(UP);
hilbert_level(level-1,RIGHT);
move(RIGHT);
hilbert_level(level-1,UP);
break;
case UP:
hilbert_level(level-1,LEFT);
move(DOWN);
hilbert_level(level-1,UP);
move(RIGHT);
hilbert_level(level-1,UP);
move(UP);
hilbert_level(level-1,RIGHT);
break;
case DOWN:
hilbert_level(level-1,RIGHT);
move(UP);
hilbert_level(level-1,DOWN);
move(LEFT);
hilbert_level(level-1,DOWN);
move(DOWN);
hilbert_level(level-1,LEFT);
break;
} /* switch */
} /* if */
}
/* End of File */
Dithering with just one neighbor doesn't sound very helpful — unless you're clever about how you visit the neighbors.
Dithering is a technique for displaying a color or gray scale that is not available in the display system's palette. A dithering algorithm selects colors from the palette and assigns them to the pixels in a given region such that their average value is close to the desired color or gray scale. Human observers who are far enough from the image will not see the individual pixels, but will perceive their average color. Thus, through dithering, an image can display realistic colors, even with a restricted palette. The biggest costs of dithering are loss of image resolution and extra computation required to do the dithering. In this article, I will present a new algorithm that provides a good balance between the benefits and costs of dithering.
The two most common types of dithering algorithms are ordered dither and error diffusion dither. These are described in more detail in the sidebar, "Ordered Dither and Error Diffusion Dither." You may also wish to consult reference [1]. These algorithms both have their advantages and disadvantages. The advantage of error diffusion dither over ordered dither is that it can choose colors from any color palette. When you reduce a full color image to palette color (e.g., with 256 colors), you will want to create an optimal palette for the image for best quality. The error diffusion dither algorithm can use such a palette, but the ordered dither algorithm can dither only to colors from a special "uniform" palette (see "Ordered Dither and Error Diffusion Dither" sidebar). Ordered dither has the advantage that the dithering of one pixel does not influence the dithering of surrounding pixels; that is, the dither is "localized." This is especially useful when the dithered image is used in an animation sequence.
I will present a novel dithering algorithm, which I call Riemersma dither, that can reduce a gray-scale or color image to any color palette, and that restricts the influence of a dithered pixel to a small surrounding area. Riemersma dither thereby combines the advantages of ordered dither and error diffusion dither. Applications for the new algorithm include image-processing systems using quadtrees and animation file formats that support palette files (especially if those animation file formats use some kind of delta-frame encoding).
The first step towards a dither procedure that combines localized dither and optimized color mapping is to reduce the error diffusion algorithm to its bare essentials. Error diffusion works by compensating for the difference between a source pixel's color (the desired color) and its palette-mapped color as it appears in an output pixel. This compensation is accomplished by distributing the difference among a set of neighboring pixels.
I have reduced this set of neighboring pixels to a single pixel. Normally, the error diffusion algorithm spreads the error (the difference between the source pixel and the mapped pixel) both horizontally and vertically, but if you only have one pixel in which to carry the error, the algorithm implicitly becomes one-dimensional. This one-dimensional nature is relevant because the next step is to align the algorithm to a Hilbert curve. A Hilbert curve is a "space-filling curve." It is a kind of fractal that runs through a square grid and visits every grid point exactly once without intersecting itself. A Hilbert curve that fills a 4x4 grid is shown below.
By applying this unidimensional algorithm to a Hilbert curve, I have obtained a very simple, but effective way to dither a two-dimensional picture. See the sidebar, "The Hilbert Curve," for more information on Hilbert curves. The forwarded error mentioned previously holds the accumulated total of all quantization errors. This accumulated error is added to the value of a pixel prior to quantizing it. For example, when dithering a gray-level image to black and white, every pixel is mapped to either black (level 0) or white (level 255), and the error introduced by this mapping is in the range -127 to +128 for every pixel (this is the quantization error). The single error accumulator is the sum of the quantization errors of all pixels that were quantized so far.
The image in Figure 1 is the well-known "Lena" picture dithered with the Hilbert curve and the error propagation as described above. Although I have simplified the error diffusion algorithm, a change in a pixel near the start of the image may still affect, through the error accumulator, the quantized values of a large number of subsequent pixels.
Conventional error diffusion dithering algorithms accumulate the quantization error continuously throughout the entire image. While continuous accumulation simplifies the algorithms, it provides no clear benefit in terms of image quality. It would make more sense to restrict the accumulation of error to a region surrounding the pixel in question. This localization of accumulated error is simple in concept to accomplish. Instead of adding the accumulated error of all previous pixels to the current pixel before quantizing it, you just add the accumulated error of the last, say, 16 pixels. In other words, you keep a list of the quantization errors of a number of recent pixels. At every new pixel to quantize, you add to the pixel value all the quantization errors in the list. Then you quantize the pixel. This creates a new quantization error (for the pixel just quantized). This error is added to the list, while removing the oldest error entry from the list.
There is a pitfall in this scheme. When removing an entry from the quantization error list, you are effectively adding a bias of the negative value of the removed entry to the list. An example will make this clearer. Assume that you keep a list of the four most recent quantization errors. When dithering a gray surface (gray value 128) to black and white, the progression would be as shown in Table 1.
The following will help in interpreting the table:
On the sixth line, the leftmost value was dropped from the list. The accumulated error of the list is suddenly 127 higher than what it should have been. Dropping the oldest entry from the root creates a severe error in the dithered result.
In the quantization error list of the above example, all entries have the same weight. The solution to the problem of the "bias" introduced by dropping a value from the list is to give the old entries in the quantization error list a (much) smaller weight than the new entries. As an entry moves down from the top of the list towards the bottom, its weight gradually lowers. In my experience, an exponentially decaying curve works well.
The weight of every list item is:
where i is the position of the item in the list (starting at 0 and ending at list length -1), r is the ratio of the highest weight to the lowest weight, and b is the base value for the exponential curve. b is computed as:
where q represents the length of the list.
The results in Table 2 apply to the same example, but now the four values have a weight attached to them. The oldest entry has the weight of 1/4th of the newest entry.
The resulting dither pattern in Table 2 is what we would expect for the dither pattern for mid-level gray. Note that the quantization error list is much too short to achieve a decent dither quality. I regard an error list of 16 as a minimum. In my own implementation, I have also set the weight ratio of the oldest entry versus the newest entry to 1:16 (instead of to 1:4 as in Table 2).
A simple example program in C (Listing 1), which created the dithered "Lena" picture in Figure 1, implements the concepts of this article. It was tested with Borland C++ for DOS in large memory model, but I have tried to avoid non-portable constructs. The bulk of the code deals with the Hilbert curve generation. The implementation of the error propagation is straightforward.
Although this article focuses on dithering to black and white, the same principles apply to the case where intermediate gray levels are available. In a similar vein, the algorithm can be adapted to color images (and restricted color palettes) by keeping three quantization error lists for each of the red, green, and blue channels.
Dithering along a space-filling curve is a relatively new subject of study. I expect that the results of this simple dither can be adjusted and improved, for example:
Related to dithering are the topics of computing the best color map (creating an optimal palette) and of performing quantization as accurately as possible.
[1] Robert Ulichney. Digital Halftoning (MIT Press, 1987). This book presents an overview and discussion of various dithering and halftoning algorithms. It is becoming somewhat outdated: new error diffusion matrices like Burkes and Sierra are absent from the book, and the research on dithering using space-filling curves is also missing.
[2] Pieter Gosselink. "Gosselink Dither." Dr. Dobb's Journal, December 1994. Gosselink dither combines ordered dither and error diffusion by creating a large set of square tiles, one for each possible color in an RGB cube. Each tile is then dithered using a simplified error diffusion algorithm. The dithered tiles are used as ordered dither matrices, where the RGB indices of each pixel in the original picture select the tile to use. The memory requirements of Gosselink dither are steep.
[3] Douglas Voorhies. "Space-Filling Curves and a Measure Of Coherence." Graphic Gems II, edited by James Arvo (Academic Press, 1991), pp 26-30. This article discusses the coherence level of the Peano and the Hilbert curves. The coherence level is a measure of the compactness of the area occupied by pixels at sequential positions on a space-filling curve. The smaller and more circular shaped that this area is, the higher the coherence between the pixels on the curve.
Thiadmer Riemersma writes software toolkits for computer graphics and animation for his company, CompuPhase, in the Netherlands. Taking an engineering approach to software development, he tries to balance performance, visual quality, and feasibility of the implementation. You may contact Thiadmer via the company homepage at www.compuphase.com.
Dithering is a subject that could cover an entire book, and indeed such a book has been written. Digital Halftoning by Robert Ulichney [1] discusses and analyzes several well-known dithering techniques. Of all dithering algorithms, the two most common are ordered dither and error diffusion dither. Both of these algorithms run through the source image scan line by scan line; they do not employ a Hilbert curve or other space-filling curves.
For the following discussion I will assume that we are given a gray-scale image with 256 shades of gray, which are assigned intensity values from 0 to 255 (for black to white respectively), and that we are dithering the image to a bi-level, black-and-white image.
Ordered dither algorithms use a fixed matrix; two common approaches are "clustered dot ordered dither" and "dispersed dot ordered dither." The halftoning procedure used in newspapers is an example of clustered dot dither.
The ordered dither algorithm compares the gray level of a pixel to a level in a dither matrix. Depending on this comparison, the pixel is set to either black or white. The matrix contains threshold values between 0 and 255. The size of the matrix and the arrangement of the values have an important effect on the dither process. When a picture is larger than the dither matrix (and it usually is), the matrix is tiled over the image. Every pixel in the source image is compared to only one level in the dither matrix. Neighboring pixels have no effect on the dither result.
Error diffusion dither algorithms spread the error (caused by quantizing a pixel) over neighboring pixels; the best known example of error diffusion dither is the "Floyd Steinberg" dither. In error diffusion, the value of a pixel is compared to a single threshold of 50 percent gray (that is, when black is 0 and white is 255, each pixel is compared to 128). The pixel is set to either black or white, based on that comparison. After outputting the dithered pixel, the algorithm calculates the difference between the source pixel and the output pixel (a simple subtraction), and then it spreads this difference (the "error") over neighboring pixels. In other words, the algorithm changes the source image as it goes through it. The error is partitioned between pixels to the right of the current pixel and pixels below the current pixel. The "diffusion" part of the error diffusion algorithm is thus two-dimensional.
The ordered dither matrices have the advantage that the dither operation is a point process. The value to which a pixel is quantized depends only on the pixel color and its position. With error diffusion dithers, the dither result of any pixel also depends on the dither results of its neighboring pixels. This means that modifying one pixel in the source image can provoke changes in a large area of the dithered result.
The error diffusion dither algorithms have the advantage that they create images containing less low-frequency content, which makes the dithered pictures more detailed and more attractive (see Digital Halftoning, page 236). More important still is that error diffusion dither does not require a uniform color palette when dithering to a gray scale or to a reduced color palette. On computer displays, dithering is often used in combination with color quantization, which seeks to produce an optimized palette. The error diffusion algorithms can exploit an optimized (and probably non-uniform) palette. The ordered dither must use either a uniform palette or a subset of the optimized palette for dithering.
Especially in the case of (frame) animation, the "point process" attribute of ordered dither is important, because it can dramatically reduce the size of the delta between one frame and the next. On the other hand, for best image quality, the value of an optimal palette (and a dither that can exploit that full palette) should not be underestimated. Hence, the focus on new dither algorithms that combine these features.
The Riemersma dither algorithm presented here is a compromise: it is an area process instead of a point process, and the produced dither is coarser than that of the error diffusion algorithm. The number of neighboring pixels that are affected by the dither value of the current pixel can be selected; I regard 16 neighboring pixels as a minimum to get a respectable dither. In the main article, I said that the accumulated error list used in Riemersma dither should be at least 16 entries long. In the case of Riemersma dither (but not other dithering algorithms) these two statements really amount to the same thing: a respectable dither requires that a change in one pixel adds to the accumulated error that affects the dither value of at least 16 subsequent pixels along the Hilbert curve.
The Hilbert curve is a space-filling curve that visits every point in a square grid with a size of 2x2, 4x4, 8x8, 16x16, or any other power of two.
It was described by David Hilbert in 1891. I have heard that Niklaus Wirth published a recursive algorithm to produce Hilbert curves in Algorithms + Data Structures = Programs, Second Edition (Prentice Hall, 1976). I am unable to confirm this because the book is now out of print. I developed the algorithm presented here independently; however, if Wirth's book does indeed contain an implementation of the algorithm, I suspect that my approach is quite similar to his.
The Hilbert curve has applications in image processing, especially image compression and dithering. It has advantages in those operations where the coherence between neighboring pixels is important [3]. The Hilbert curve is also a special version of a quadtree; any image-processing function that benefits from the use of quadtrees may also use a Hilbert curve.
The basic elements of a Hilbert curve are what I call "cups" (a square with one open side) and "joins" (a vector that joins two cups). The "open" side of a cup can be top, bottom, left, or right. In addition, every cup has two endpoints, and each of these can be the "entry" point or the "exit" point. So, there are eight possible varieties of cups. In practice, a Hilbert curve uses only four types of cups. In a similar vein, a join has a direction: up, down, left, or right.
A first-order Hilbert curve is just a single cup:
It fills a 2x2 space. The second-order Hilbert curve replaces that cup by four (smaller) cups, which are linked together by three joins
The link between a cup and a join has been marked with a fat dot in the figure.
Figure 2 shows how to produce a second-order curve from any of the four possible cups. Read this figure from left to right. For example, to produce a second-order Hilbert curve from an upward pointing ("up") cup (top row of figure), first produce a left cup; attach a vertical join to its lower endpoint; attach the bottom end of the vertical join to the left endpoint of an up cup; attach a horizontal join to the right endpoint of the up cup; attach the right end of the horizontal join to the left endpoint of another up cup; attach a vertical join to the right endpoint of this up cup; and attach the top of this vertical join to the bottom endpoint of a right cup. Using this method, you can produce a Hilbert curve of the next-highest order from any existing Hilbert curve. Simply replace every cup in the first curve with four cups as outline by Figure 2.
Listing 2 presents a function that computes the Hilbert curve. Note that the curve is symmetrical around the vertical axis. If you were drawing Hilbert curves, it would therefore be sufficient to draw half of the Hilbert curve and flip it on its axis of symmetry.
The Hilbert curve is easy to generate. When applied over a digitized photograph or a ray-traced image, it makes better use of the coherence of neighboring pixels than the traditional scan-line based approach.
For those interested in a three-dimensional variant of the Hilbert curve (a curve that fills a cube), see http://math.uwaterloo.ca/~wgilbert/Research/HilbertCurve/HilbertCurve.html.
| Quantization error list | source pixel | adjusted pixel | quantized pixel | |||
| 0 | 0 | 0 | 0 | 128 | 128 | 255 |
| 0 | 0 | 0 | -127 | 128 | 1 | 0 |
| 0 | 0 | -127 | 128 | 128 | 176 | 255 |
| 0 | -127 | 128 | -127 | 128 | 30 | 0 |
| -127 | 128 | -127 | 128 | 128 | 195 | 255 | 128 | -127 | 128 | -127 | 128 | 63 | 0 |
| weights | ||||||
| 0.25 | 0.40 | 0.63 | 1.0 | |||
Terms of Service | Privacy Statement | Copyright © 2024 UBM Tech, All rights reserved.