Creates a tetrahedral tessellation from a colour palette via 3D Delaunay triangulation. Intended for use in arbitrary-palette colour image dithering.
Tetrapal is shorthand for tetrahedral palette. It is a utility that takes a set of points in colour space and forms a triangulated irregular network of tetrahedra. As a result, any colour can be represented as the weighted sum of a number of existing palette colours.
The main motivation behind this library is to enable an efficient implementation of colour image ordered dithering for irregular palettes. Typically, ordered dithering is only optimal when the palette contains colours that are equally distributed in colour space, i.e. that they form a regular grid. However, the algorithm can be modified to accommodate irregular or arbitrary palettes by representing input colours as the weighted sum of existing palette colours. Tetrapal provides a data structure that can determine the necessary weights much faster and with more precision than existing implementations.
| Typical Algorithm | Tetrapal Algorithm |
|---|---|
The idea to use a 3D triangulated irregular network as a means to dither colour images is not a new one. The earliest source I could find that describes such an implementation is the 1988 article "Using tetrahedrons for dithering color pictures" by Eduard Gröller and Werner Purgathofer1. However, references to this technique online are scarce. To my knowledge, this repository is the only public implementation of its kind.
Tetrapal* tetrapal_new(const float *points, const int size);Creates a new triangulation from an existing palette. The parameter *points should be a pointer to a buffer of 3D coordinates in the following format:
Where size is the number of points represented in the buffer. Internally, points are indexed according to the order they appear in the buffer, where the starting index is 0. If successful this function will return an opaque pointer to the Tetrapal data structure, otherwise it will return NULL.
Tetrapal expects coordinate values to be in the range 0.0 to 1.0. Values beyond this range will be clamped. Internally, coordinates are transformed and represented as integers with 16 bits of precision.
void tetrapal_free(Tetrapal* tetrapal);Triangulation data can safely be freed by passing the Tetrapal pointer to the above function.
int tetrapal_interpolate(const Tetrapal* tetrapal, const float point[3], int* indices, float* weights);Performs barycentric interpolation within a triangulation. This will return an int between 1 and 4 depending on the number of points contributing to the interpolant given by point[3]. The indices of the points will be written to the values at *indices, and their respective weights written to *weights. In the case where the number of points
The output arrays should be large enough to hold the maximum expected number of points, which is equivalent to tetrapal_element_size (please refer to the utility functions section).
This is the recommended interpolation function as it provides the best combination of stability, quality, and performance.
int tetrapal_natural_neighbour(const Tetrapal* tetrapal, const float point[3], int* indices, float* weights, const int size);In addition to standard barycentric interplation, Tetrapal is able to perform natural neighbour interpolation as well. This function returns an int specifiying the number of natural neighbours contributing to the interpolant given by point[3]. The size of the output arrays *indices and *weights should be given by the parameter size. It is possible for this function to fail, either due to a lack of available system memory or because the number of natural neighbours exceeds the size of the output arrays. In both cases, the function will return 0.
In theory the number of natural neighbours for a given input point can range from 1 to the total number of vertices in the triangulation, but in practice the maximum number is much lower for most point sets. However, to prevent failure it is advised that the output arrays are at least as large as the number of vertices in the triangulation.
It is recommended to use barycentric interpolation for most cases. Natural neighbour interpolation is much slower and the resulting dither quality may not be much better than barycentric interpolation, perceptually speaking.
int tetrapal_nearest_neighbour(const Tetrapal* tetrapal, const float point[3]);Tetrapal also supports nearest-neighbour queries within the triangulation. This can be faster than a standard linear search under the right circumstances. It is included for convenience.
Returns the index of the nearest neighbour to the query point defined by point[3].
int tetrapal_number_of_elements(const Tetrapal* tetrapal);Returns the number of real elements (simplices) in the triangulation, i.e. does not count symbolic 'infinite' elements.
int tetrapal_number_of_dimensions(const Tetrapal* tetrapal);Returns the number of dimensions spanned by the triangulation. Valid dimensions range from 0 to 3. A value of -1 indicates a null triangulation.
int tetrapal_element_size(const Tetrapal* tetrapal);Gets the size of an element in the triangulation, or in other words, the number of vertices that make up a simplex within the triangulation. This is the same as tetrapal_number_of_dimensions + 1.
int tetrapal_get_elements(const Tetrapal* tetrapal, int* buffer);Writes raw geometry data directly into *buffer. Each element is represented by a set of vertex indices which are packed contiguously into the buffer. Thus, the size of the buffer must be no less than the number of elements present in the triangulation multiplied by the number of indices representing each element, or tetrapal_number_of_elements × tetrapal_element_size.
Below is an example C implementation of a function that takes an input image and uses Tetrapal to perform ordered dithering with an arbitrary palette, producing an indexed image. The buffers *input and *palette are expected to contain RGB values from 0.0 to 1.0. The function sort_by_luminance simply orders the candidates and their corresponding weights in order of ascending/descending luminance.
// Normalised threshold matrix
const double bayer_matrix_4x4[4][4] =
{
{ 0.0 / 16.0, 12.0 / 16.0, 3.0 / 16.0, 15.0 / 16.0 },
{ 8.0 / 16.0, 4.0 / 16.0, 11.0 / 16.0, 7.0 / 16.0 },
{ 2.0 / 16.0, 14.0 / 16.0, 1.0 / 16.0, 13.0 / 16.0 },
{ 10.0 / 16.0, 6.0 / 16.0, 9.0 / 16.0, 5.0 / 16.0 }
};
void dither_image(const float *input, const int image_width, const int image_height,
unsigned char *output, const float *palette, const int palette_size)
{
// Candidate arrays
int candidates[4];
double weights[4];
// Triangulate the palette
Tetrapal* tetrapal = tetrapal_new(palette, palette_size);
// Iterate over pixels in the input image
for (int y = 0; y < image_height; y++)
{
for (int x = 0; x < image_width; x++)
{
// Get the current pixel from the input buffer
const int image_index = x + y * image_width;
const float* pixel = &input[image_index * 3];
// Interpolate within the triangulation to get the candidates for the current pixel and their weights
tetrapal_interpolate(tetrapal, pixel, candidates, weights);
// Sort the candidates by luminance
sort_by_luminance(candidates, weights, palette);
// Use the value in the threshold matrix to select a candidate
const double threshold = bayer_matrix_4x4[y % 4][x % 4];
double sum = 0.0;
// Accumulate the sum of weights until we pass the threshold
for (int i = 0; i < 4; i++)
{
sum += weights[i];
if (threshold < sum)
{
output[image_index] = candidates[i];
break;
}
}
}
}
// Remember to free the triangulation data
tetrapal_free(tetrapal);
}Tetrapal implements 3D Delaunay triangulation using the Bowyer-Watson incremental construction algorithm. It borrows ideas from popular computational geometry libraries such as CGAL and Geogram to ensure accuracy and correctness. This includes the use of symbolic 'infinite' vertices to guarantee convexity, as well as the ability to gracefully handle degenerate inputs in the form of lower-dimensional triangulations (2D, 1D, 0D). Coincident points are ignored during triangulation but are still stored internally to maintain consistency when indexing vertices.
A combination of extended precision integer arithmetic and static filtering via predetermined error bounds2 is used to provide robustness for geometric predicates in cases where standard arithmetic may fail to give an accurate result.
Tetrapal supports extrapolation of points outside the convex hull via projection onto the nearest surface element (which may be a triangle, line segment, or single vertex) and returning the barycentric coordinates with respect to this element. Natural-neighbour interpolation with Sibson weights is supported in 2D and 3D, but only for points that lie within the convex hull.
Point location is performed via remembering stochastic walk3, where a kd-tree approximate nearest neighbour search is used to quickly locate a starting tetrahedron close to the query point.
As its main purpose is to process colour palettes, some assumptions have been made to simplify the algorithm. However, Tetrapal can still be used as a general-purpose Delaunay triangulation library provided its limitations are observed. These include:
- That it expects the input to be normalised between 0.0 and 1.0.
- That the desired precision of the input does not exceed 1 / 65535.
- It lacks many features compared to more mature libraries.
The three tables below compare the running time of a Tetrapal-based ordered dithering algorithm against the more well-known algorithms of Thomas Knoll4 and Joel Yliluoma5. A different independent variable was chosen for each test (palette size, threshold matrix size, and input image size, respectively). The "Yliluoma's ordered dithering algorithm 2" variant of Yliluoma's algorithms was implemented. The construction of the Tetrapal data structure itself is included in the timings, and barycentric interpolaton was used for all tests. All image/palette colours were transformed to linearised sRGB space prior to processing.
| Palette Size | Tetrapal | Knoll | Yliluoma | Matrix Size | Tetrapal | Knoll | Yliluoma | Image Size | Tetrapal | Knoll | Yliluoma | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 8 | 0.091s | 0.830s | 6.923s | 2x2 | 0.127s | 0.100s | 0.364s | 128x128 | 0.010s | 0.092s | 1.089s | ||
| 16 | 0.133s | 1.457s | 14.837s | 4x4 | 0.118s | 0.392s | 2.882s | 256x256 | 0.038s | 0.365s | 4.226s | ||
| 32 | 0.146s | 2.733s | 30.529s | 8x8 | 0.122s | 1.494s | 17.032s | 512x512 | 0.146s | 1.455s | 16.587s | ||
| 64 | 0.166s | 5.240s | 61.521s | 16x16 | 0.128s | 5.676s | 92.412s | 1024x1024 | 0.523s | 5.658s | 65.638s | ||
| 128 | 0.178s | 10.194s | 126.047s | 32x32 | 0.124s | 22.421s | 470.400s | 2048x2048 | 2.039s | 22.222s | 262.450s | ||
| 256 | 0.186s | 20.072s | 257.878s | 64x64 | 0.129s | 88.046s | 2297.365s | 4096x4096 | 8.272s | 87.802s | 1056.591s |
Tetrapal is faster in almost all cases. This is because both Knoll and Yliluoma are iterative algorithms whose time complexity is a factor of both the palette size
| Algorithm | Tetrapal | Knoll | Yliluoma |
|---|---|---|---|
| Complexity |
This table shows the peak signal-to-noise ratio (PSNR) for the output of each algorithm as a rough estimate of the dither quality (higher is better). A Gaussian blur was applied to the dithered output images before measuring the PSNR. Two different 16-colour palettes were tested for each image; a custom palette by Andrew Kensler7 that remained the same for all images, and an adaptive palette generated using a variance-based colour quantisation algorithm by Xiaolin Wu8.
| Image | Tetrapal | Knoll | Yliluoma | Palette |
|---|---|---|---|---|
| Peppers (Custom) | 27.17dB | 27.03dB | 27.03dB | |
| Peppers (Adaptive) | 24.33dB | 24.37dB | 24.37dB | |
| Mandrill (Custom) | 21.14dB | 21.10dB | 21.10dB | |
| Mandrill (Adaptive) | 20.58dB | 20.61dB | 20.61dB | |
| Parrots (Custom) | 27.77dB | 27.65dB | 27.65dB | |
| Parrots (Adaptive) | 20.42dB | 20.42dB | 20.42dB |
Knoll and Yliluoma score exactly the same for each image, while the Tetrapal algorithm scores slightly better than the other two using the custom palette and slightly worse when using an adaptive palette. All algorithms scored worse using an adaptive palette, suggesting that a palette generation algorithm specifically catered towards dithering should be preferred in general. It is possible that Tetrapal is better suited to varied colour palettes with a good spread of colours, unlike the generated palettes which tended to contain many similar colours.
This table records the size in memory of the Tetrapal data structure for various palette sizes, with random colours generated uniformly inside the unit cube. The memory consumption of the Knoll and Yliluoma algorithms that were implemented were trivial (equivalent to the size of the threshold matrix in the worst case, which is usually <1KB).
| Palette Size | 8 | 16 | 32 | 64 | 128 | 256 |
|---|---|---|---|---|---|---|
| Memory | 2KB | 4KB | 8KB | 16KB | 32KB | 64KB |
Here is a visual comparison between the dithered output of Tetrapal, Knoll, and a standard implementation of ordered dithering. The 16-colour CGA palette was used.
| Algorithm | Test Image 1 | Test Image 2 |
|---|---|---|
| Tetrapal | ||
| Knoll | ||
| Standard |
These final images illustrate the visual differences between dithering via barycentric interpolation and natural neighbour interpolation. A 256x256 void-and-cluster threshold matrix was used for both images. In general, natural neighbour dithering produces perceptually smoother gradations but introduces more high-frequency noise as a result of considering a greater number of candidate colours per pixel. This can best be seen in the appearance of the sky in each of the example images.
| Barycentric | Natural Neighbour |
|---|---|
Footnotes
-
E. Gröller and W. Purgathofer, "Using tetrahedrons for dithering color pictures" (1988). ↩
-
S. Fortune and C. J. Van Wyk, "Static Analysis Yields Efficient Exact Integer Arithmetic for Computational Geometry" (1996). ↩
-
O. Devillers, S. Pion and M. Teillaud, "Walking in a Triangulation" (2006). ↩
-
T. Knoll, "Pattern Dithering" (1999). ↩
-
J. Yliluoma, "Joel Yliluoma's arbitrary-palette positional dithering algorithm" (2011). ↩
-
E. P. Mücke, I. Saias and B. Zhu, "Fast randomized point location without preprocessing in two- and three-dimensional Delaunay triangulations" (1998). ↩
-
A. Kensler, "Pixel Art Palettes for Free" (2016). ↩
-
X. Wu, "Efficient Statistical Computations For Optimal Colour Quantisation" (1995). ↩