Hilbert-Curve Dithering

DEFINITION

visit the cells in Hilbert order, not in scan order
b = r^(1/(q−1)),  wᵢ = bⁱ ⁄ Σⱼ bʲ,  i = 0 … q−1
aₖ = sₖ + Σᵢ wᵢ·eₖ₋q₊ᵢ,  oₖ = nearest level to aₖ,  eₖ = aₖ − oₖ

NOTES

Plate 52 pushes each rounding error onto the four neighbours a raster scan has not reached yet, and the scan is the flaw: a mistake made at the left of a row is still travelling when the row ends, and the eye reads what comes back as worms leaning in the direction of the scan. Riemersma changed the order of the visits instead of the arithmetic. Walk the image along the Hilbert curve of plate 18 — the line that reaches every cell of a square and keeps cells that are close on the line close in the square — and carry a short list of the errors just made, each weighted by how recently it was made, the weights falling geometrically from newest to oldest. Every cell is then judged against its own tone plus everything the walk has just got wrong nearby, and nearby now means nearby in both directions at once rather than along one row. Measured here at the published settings, the sixteen cells such a list covers lie at most 7.28 cells apart and 5.57 apart on average, which is under twice the square root of sixteen and is the square-root growth the curve promises; the same sixteen consecutive cells of a boustrophedon scan lie 15.00 and 14.51 apart, and of a plain raster scan up to 127, the whole width of the lattice. The bias the scan direction leaves behind falls with it, from 0.43 to 0.04 on a flat field. The bright clump crawling the plate is that list, drawn at one instant of the walk with every cell shaded by the weight it carries, and the faint thread behind it is the curve it came along. Turn the memory down and the plate forgets an error before it has any chance to cancel it; turn the decay ratio up and only the newest error counts at all, which is rounding with the diffusion switched off. Bayer at plate 51 needs no memory whatsoever, and the trade between the three of them is what this corner of the atlas is for.

PROVENANCE

Origin

Thiadmer Riemersma, "A Balanced Dither Algorithm", C/C++ Users Journal 16(12), December 1998. The author has maintained an expanded version of the article at compuphase.com/riemer.htm ever since, and that page is the src link above and the text this implementation was written from. Both checked 17 August 2026.

What the record says, and where it disagrees with itself

The author names the article "A Balanced Dither Algorithm". The ACM Digital Library and the Dr Dobbs archive that absorbed the journal both index the same piece as "A Balanced Dithering Technique". Neither title is wrong to cite and the discrepancy is recorded here rather than quietly resolved. There is no DOI: the ACM identifier 10.5555/306490.306494 carries the placeholder prefix ACM assigns to items it holds without one, and queried against the Crossref API on 17 August 2026 it answers 404. So it is not offered here as a doi:, and src points instead at the page the author keeps, which is the only stable resolvable address this work has.

Standing

Public domain in the sense this atlas uses the phrase: a traversal order and a weighting rule, set out in a published article. No patent is claimed by the article or the page and none was found. The compuphase page is itself licensed CC BY-SA 3.0, which covers the text on that page and the C listing it carries, neither of which is in this implementation — see below.

Implemented from the description, not from the listing

The compuphase page ships a complete C program. It was not read and nothing here is adapted from it. What was taken is the published prose: that the image is traversed along a Hilbert curve, that a list of the most recent quantisation errors is carried, that the list is short — the author calls sixteen entries a practical minimum and sets the oldest entry at one sixteenth the weight of the newest — and the decay rule itself, that for a ratio r over a list of length q the geometric base is b = r^(1/(q−1)) and item i carries weight proportional to bⁱ. The Hilbert walk below is the standard d to (x, y) bit recursion, written out again here rather than shared with plate 18, on the same reasoning that gives three separate plates their own copy of the colour matching fit: a specimen has to be readable on its own in this drawer.

Why the weights are divided by their sum

The published rule fixes the shape of the list but not its total, and the prose says only to add the weighted errors to the pixel. Left at the total the formula produces, the sixteen entries Riemersma names at a ratio of sixteen sum to 5.618, so every error is fed back more than five times over. Measured through this exact loop on a flat grey, that puts the mean output 0.25 of full scale away from the tone it exists to reproduce — a quarter of the range, on the one input where a dither has nothing to do but be right. Dividing by the sum leaves every ratio in the list untouched, and therefore the whole published shape untouched, and brings the worst tone error over greys from 0.05 to 0.95 down to 0.0001. Floyd and Steinberg publish the same conservation as part of their weights: 7 + 3 + 5 + 1 over 16 is exactly 1. It is what makes an error-diffusion filter reproduce tone at all, and it is normalised here on that ground rather than on any reading of a listing that was never opened.

The square is cropped, not skipped, and that was measured both ways

The curve is defined on a 2ⁿ square and this frame is not one. The usual treatment is to run the walk over the covering square and skip the cells that fall outside the picture, and it was written that way first. Measured, it breaks the one property the curve is here for: at the aspect the drawer uses, 14 of 10,879 consecutive pairs of kept cells are not neighbours at all, the worst of them 57 cells apart, and the worst span of sixteen consecutive kept cells goes from 7.28 on the square to 63.13 — a raster scan number, on the plate whose argument is that it does not have raster scan numbers. So the square is scaled to cover the frame instead and the frame crops it. Every cell the walk visits is a cell of the complete curve, in order, with no jumps anywhere; the cells that fall off the top and bottom edges are computed and not shown. It costs about half as many cells again at the drawer aspect and nothing at all on a square canvas, and it keeps the traversal the exact object the plate is about.

Checked, not just plotted

Three measurements, all on the default 128-cell square. First, the walk is a walk: every one of the 16,384 cells is visited exactly once and every consecutive pair is a lattice neighbour, asserted rather than assumed. Second, locality: over every starting position, the greatest distance between any two of sixteen consecutive cells is 7.28 cells and the mean is 5.57, under 2√16 = 8; eight consecutive cells span 5.00 at worst and thirty-two span 11.18, against 5.66 and 11.31 for the same bound. A boustrophedon scan spans 7.00, 15.00 and 31.00 for the same three counts — length, not area — and a plain raster scan spans the full 127 at every one of them. Third, directional emphasis: dithering a flat grey, which is isotropic, and asking whether the output prefers an axis, the horizontal-against-vertical correlation gap over greys from 0.10 to 0.90 reaches 0.4262 for Floyd-Steinberg and 0.0437 walking the curve; on the rotating field this plate actually dithers, swept over twelve bearings, 0.1314 against 0.0324. What that does not say: no deterministic error diffusion escapes locking onto a periodic pattern on a perfectly flat field and this one does not either — Floyd-Steinberg collapses to an exact checkerboard at one half and reaches a diagonal gap of 1.12 at one third, where walking the curve reaches 0.54 at one quarter. The claim measurement earns is that the scan direction stops being visible, not that the lock-in goes away.

Constants

mem is the length of the error list and ratio is how much more the newest entry weighs than the oldest — the two numbers the published implementation fixes at 16 and 16:1, which are the defaults here. mem runs under the stated minimum of sixteen deliberately: at four entries the plate visibly forgets, which is the exhibit, and it is still a dither rather than a blank. lv is the number of tone levels, errordiff carrying the same constant under the same name. order is the lattice, 2ⁿ cells on a side, locked from regenerate for the reason in taxonomy.js; it defaults to 7 and reaches 8, which is the gray-scott arrangement — a constant that costs time defaults low for the grid frame budget and still ranges high for someone looking closely, because the drawer repaints outside both tiers. Liveness, measured the way regenerate actually jitters, all constants at once, 40 tuples across 40 seeds: 40 of 40 live, every one of them putting at least 46.5% of the lattice on a level other than the majority one — the right measure here, because a dead dither is one level everywhere and a pixel-contrast probe reads a healthy number off a single stray cell. Swept exhaustively rather than randomly, all 3,072 grid points of the declared box at a 4-step in mem and a 2-step in ratio, the worst minority share anywhere is 46.3% and the reproduced tone stays within 0.0015 of the source. There is no dead ground in this box, which is what a threshold rule with a conserved feedback loop ought to give and is worth having checked rather than assumed.

Source

https://www.compuphase.com/riemer.htm

HOUDINI · VEX

The same published mathematics as a Detail Wrangle body. Paste it into a Wrangle with Run Over set to Detail; every constant is the published value plus a tweak channel, so Create Spare Parameters gives a slider that starts where the paper does.

// FORMA — PL. 131 · HILBERT-CURVE DITHERING — Thiadmer Riemersma, 1998
//   visit the cells in Hilbert order, not in scan order
//   b = r^(1/(q−1)),  wᵢ = bⁱ ⁄ Σⱼ bʲ,  i = 0 … q−1
//   aₖ = sₖ + Σᵢ wᵢ·eₖ₋q₊ᵢ,  oₖ = nearest level to aₖ,  eₖ = aₖ − oₖ
// Paste into a Detail Wrangle (Run Over: Detail), no inputs needed.
// Written from the published mathematics, not adapted from any code.
// Constants arrive at their published values. Press the node's Create
// Spare Parameters button and every tweak becomes a slider — starting
// at 0, the published figure, and moving in the constant's own units.
// https://forma-gen.com/#plate=riemersma

float p_mem   = 16 + chf('mem_tweak');        // error list length · live 4 .. 32
float p_ratio = 16 + chf('ratio_tweak');      // weight ratio, newest : oldest · live 2 .. 64
float p_lv    = 2 + chf('lv_tweak');          // tone levels · live 2 .. 5
float p_order = 7 + chf('order_tweak');       // lattice order — 2ⁿ cells square · live 6 .. 8

// The plate's own colour: FORMA's COLOUR accent as a cosine ramp,
// brightest near t = 0 and t = 1, near-black around t = 0.5.
vector forma_ramp(float t){
  return set(
    0.46 + 0.5 * cos(6.28318530718 * (t + 0)),
    0.2002 + 0.2176 * cos(6.28318530718 * (t + 0.05)),
    0.3896 + 0.4235 * cos(6.28318530718 * (t + 0.1)));
}

// Riemersma dither: the same quantise-and-carry as the Floyd-Steinberg port
// two plates over, walked along a Hilbert curve instead of along the rows.
// A short list of the most recent errors travels with the walk, weighted
// geometrically so the newest entry outweighs the oldest by p_ratio, and each
// cell is judged against its own tone plus everything on that list.
//
// The weights are divided by their sum. That is not in the published prose,
// which fixes only their shape, and it is not optional: the series as written
// totals 5.618 at the published 16 entries and 16:1, feeding every error back
// more than five times, which puts a flat grey a quarter of full scale off its
// own tone. Dividing changes no ratio inside the list. Floyd and Steinberg
// publish the same conservation inside theirs — 7 + 3 + 5 + 1 over 16 is 1.
//
// The source is a sine grating that on the plate rotates on the clock through
// its own phase; a cook has neither, so it is frozen at 30 degrees — the same
// angle errordiff cooks at, deliberately, so the two ports are one test card
// under two traversals and can be laid side by side. A staging decision on the
// record: the mathematics is the traversal and the carry, not the card.
//
// The lattice is the full 2^order square. On the page the frame crops that
// square; a cook has no frame, so nothing is cropped and nothing is skipped.
// Points come out in the order the curve visits them, so @ptnum IS the
// position along the walk — an Add SOP set to polygon-by-point-order draws the
// traversal itself. One point per non-black cell, the plate's own tone-scaled
// ramp as Cd and the tone as Alpha. The crawling error-list cursor the plate
// draws over the dither is the animation's, not the construction's, and is not
// ported, on turmite's precedent. Deterministic by construction: there is no
// randomness anywhere in this algorithm.

// d -> (x, y) on the 2^k square by the standard bit recursion: at each scale
// two bits of the index choose the quadrant and the lower bits are reflected
// into it. Written with /2 and %2 rather than shifts — the low bit of a xor is
// just parity, and VEX has no << to lean on anyway. Every input arrives as a
// parameter, because a VEX function cannot see the snippet scope.
vector forma_riem_cell(int d; int n){
    int x = 0, y = 0, tt = d;
    for (int s = 1; s < n; s *= 2){
        int rx = (tt / 2) % 2;
        int ry = (tt + rx) % 2;     // low bit of tt ^ rx — xor is parity there
        if (ry == 0){
            if (rx == 1){ x = s - 1 - x; y = s - 1 - y; }
            int tmp = x;  x = y;  y = tmp;
        }
        x += s * rx;  y += s * ry;
        tt /= 4;
    }
    return set(float(x), float(y), 0.0);
}

int order = int(rint(p_order));
int n = 1;
for (int g = 0; g < order; g++) n *= 2;
int L = max(2, int(rint(p_lv)));
int q = max(2, int(rint(p_mem)));
float r = max(1.0001, p_ratio);

// the geometric list, then normalised
float b = exp(log(r) / float(q - 1));
float wt[];
resize(wt, q);
float wsum = 0.0;
for (int i = 0; i < q; i++){
    wt[i] = pow(b, float(i));
    wsum += wt[i];
}
for (int i = 0; i < q; i++) wt[i] /= wsum;

float err[];
resize(err, q);
for (int i = 0; i < q; i++) err[i] = 0.0;

float ang = M_PI / 6.0;
float co = cos(ang), si = sin(ang);

// the one ordered pass that is the whole algorithm, in curve order
for (int d = 0; d < n * n; d++){
    vector c = forma_riem_cell(d, n);
    int x = int(c.x), y = int(c.y);
    float u = float(x) / float(n), v = float(y) / float(n);
    float tone = 0.5 + 0.5 * sin(6.283 * 1.5 * (u * co + v * si));

    float carried = 0.0;
    for (int k = 0; k < q; k++) carried += err[k] * wt[k];
    float a = tone + carried;

    float lit = min(float(L - 1), max(0.0, rint(a * float(L - 1))));
    float tq = lit / float(L - 1);

    // the error just made joins the list and the oldest falls off the end
    for (int k = 0; k < q - 1; k++) err[k] = err[k + 1];
    err[q - 1] = a - tq;

    if (tq <= 0.0) continue;    // the ground is the ink, not a point
    // canvas y runs down; negated so the grating leans as the plate shows it
    int pt = addpoint(0, set(float(x - n / 2), float(n / 2 - y), 0.0));
    setpointattrib(0, "Cd", pt, forma_ramp(0.86 + 0.2 * tq) * tq);
    setpointattrib(0, "Alpha", pt, tq);
}