Skip to content

mHash — Marr–Hildreth Hash

Bring the image to 512×512, spread its gray levels evenly, and pass it through the Mexican-hat filter of Marr and Hildreth, which answers where the image has edges and fine detail at the filter's scale. Sum the answer over a 31×31 grid of blocks, and record, for each block of 64 small windows of that grid, whether it answers more strongly than its neighbors in the window. The result is 576 bits, a fingerprint of where the image's detail lies, much finer than the 64 bits of the other hashes and much more sensitive to where exactly it lies.

Call ph_compute_mhash(), or PH_ALGO_MHASH in ph_compute_digest()
Output digest of 72 bytes, 576 bits, compared with ph_hamming_distance_digest()
Source the operator: Marr and Hildreth, 1980; the hash built on it: pHash's ph_mh_imagehash(), which has no paper (provenance)
Cost 1.8 ms on a 400×400 image, 91 ms on 20 Mpx, after decoding1

The steps

The steps below are those of the default parameters, a 512×512 image and a kernel reaching 8 pixels from its center; ph_context_set_mhash_params() sets others, below. Every picture on this page is computed by the library from the same photograph, and the digest at the end is what ph_compute_mhash() returns for it.

The stages of mHash on the example photograph: the decoded image, the blurred and resized grayscale image, the equalized image, its response to the kernel, the block sums and the 576 bits The stages of mHash on the example photograph: the decoded image, the blurred and resized grayscale image, the equalized image, its response to the kernel, the block sums and the 576 bits

1. Grayscale. Each pixel becomes one luminance value,

\[ Y = \left\lfloor \frac{38R + 75G + 15B}{128} \right\rfloor , \]

the same integer approximation of the ITU-R BT.601 weights as for every grayscale hash. ph_context_set_gray_weights() changes them; image preparation has why these weights.

2. Blur. The grayscale image is blurred at full resolution by a Gaussian of standard deviation 1 pixel, cut off at 3 pixels on either side and scaled so its seven weights add up to 1, first along the rows and then along the columns; the edges repeat their last pixel, and each result is rounded to a whole gray level.

3. Normalize to 512×512. The blurred image is resampled to 512×512 with the Mitchell filter, regardless of its aspect ratio: a wide picture is squeezed, a tall one stretched. Every image then reaches the filter at the same size, so the filter's scale is the same fraction of every picture.

4. Equalize. Each gray level \(v\) is replaced by its rank among the image's \(N = 512^2\) pixels,

\[ v \;\mapsto\; \left\lfloor 255 \cdot \frac{C(v) - C_{\min}}{N - C_{\min}} + \frac12 \right\rfloor , \]

where \(C(v)\) counts the pixels at level \(v\) or below and \(C_{\min}\) the pixels at the darkest level present. The darkest level becomes 0, the brightest 255, and the levels between spread so that each takes its share of the range.

The image before and after equalizing, with the histogram of its gray levels and their cumulative share The image before and after equalizing, with the histogram of its gray levels and their cumulative share

The cumulative share becomes a straight line: whatever the exposure, the picture uses the whole range evenly. A change that keeps the order of the gray levels, such as a brighter exposure or a tone curve, changes the equalized image only through rounding.

5. The kernel. The filter is the Laplacian of Gaussian, sampled as the Mexican hat

\[ K(x, y) = (2 - A)\, e^{-A/2}, \qquad A = \frac{x^2 + y^2}{\alpha^{2\ell}}, \qquad x, y = -s \dots s, \quad s = \lfloor 4 \alpha^{\ell} \rfloor , \]

with the defaults \(\alpha = 2\) and level \(\ell = 1\): \(s = 8\), a 17×17 kernel. Its center is positive, a ring around it negative, and its weights nearly cancel:

The 17×17 kernel at the defaults sums to 0.017, against a weight of 2 at its center.

The 17×17 kernel, positive center in blue and negative ring in orange, and its middle row against the continuous curve The 17×17 kernel, positive center in blue and negative ring in orange, and its middle row against the continuous curve

The kernel crosses zero at \(A = 2\), a distance of \(\sqrt 2\, \alpha^{\ell} \approx 2.8\) pixels from its center. On a smooth area it answers almost nothing; on a line or a spot about that wide it answers strongly, positive on one side of a light-dark change and negative on the other.

6. The response. Each pixel of the equalized image \(I\) is replaced by the kernel's weighted sum around it, the edges repeating their last pixel:

\[ R(y, x) = \sum_{j=-s}^{s} \sum_{i=-s}^{s} K(i, j)\, I(y + j,\, x + i) . \]

The response at every pixel of the 512×512 image, with the grid of 31×31 blocks drawn over it, and the strips at the right and bottom that no block covers shaded The response at every pixel of the 512×512 image, with the grid of 31×31 blocks drawn over it, and the strips at the right and bottom that no block covers shaded

The response is a map of detail: the spines of the cactus and the rim of the pot light up, the wall and the table stay pale.

7. Block sums. The response is summed over 16×16 blocks, 31 of them across and 31 down:

\[ B(b_y, b_x) = \sum_{y = 16 b_y}^{16 b_y + 15} \; \sum_{x = 16 b_x}^{16 b_x + 15} R(y, x), \qquad b_y, b_x = 0 \dots 30 . \]

31 blocks of 16 cover 496 of the 512 pixels; the last 16 rows and columns enter only through the kernel's reach of the blocks next to them. The library does not compute \(R\) at every pixel; it computes these sums directly, faster and with the same result.

8. Windows and bits. The grid is read through 3×3 windows placed every 4 blocks, 8 across and 8 down: 64 windows, nine block sums each. The windows do not touch: every fourth row and column of the grid, seven of each, belongs to no window, so the digest reads 576 of the 961 block sums.

The 31×31 block sums with the 64 windows outlined and numbered, the rows and columns between them belonging to none The 31×31 block sums with the 64 windows outlined and numbered, the rows and columns between them belonging to none

Each block sum in a window is compared with the mean of the nine, and sets its bit when it is above:

\[ b = \begin{cases} 1 & B > \bar B_{\text{window}} \\ 0 & B \le \bar B_{\text{window}} \end{cases} \]

One window of the example: nine block sums, their mean, and the bits of the sums above it One window of the example: nine block sums, their mean, and the bits of the sums above it

Each bit therefore says where the detail lies within its window, which sets the bits independently of the response's overall strength: the source scales the response to [0, 1] before this step, and the library leaves that out, since scaling every sum by the same positive factor and shifting it by the same amount cannot change which side of their mean any of them falls.

Folding the block sum into the kernel

Computing \(R\) at every pixel of the 496×496 area and then summing it costs \(496^2 \times 289\), some 71 million multiply-adds. The kernel weight does not depend on the pixel, so it can come out of the block sum:

\[ B(b_y, b_x) = \sum_{j,i} K(i, j) \sum_{(y,x) \in \text{block}} I(y + j,\, x + i) , \]

and what is left inside is the sum of a 16×16 square of the image, shifted by \((j, i)\). An integral image, the table of sums of every top-left rectangle, gives any such square in four lookups. That makes \(31^2 \times 289\), some 280 thousand multiply-adds, 256 times fewer.

flowchart LR
    I[Equalized image] --> R["Response at every pixel<br/>496² × 289"] --> S[Sum over each block] --> B[Block sums]
    I --> T[Integral image] --> Q["289 shifted 16×16 squares<br/>per block, 31² × 289"] --> B

The square sums are exact integers, and only 289 products are added per block, in double precision. Evaluated the long way, the response sums 289 large products that nearly cancel, and in single precision the rounding error survives into the block sum. Most of the time it does not reach a bit; where it does, it is decisive. Here is the long way against the library's digest, with the response computed in double and in single precision:

Images Response computed in Against the folded block sums
the example photograph double the same digest
the example photograph float the same digest
synthetic corpus, 24 images; not photographs double the same digest on every image
synthetic corpus, 24 images; not photographs float another digest on 1 of 24 images, 108 bits apart
200 public-domain photographs from Wikimedia Commons double the same digest on every image
200 public-domain photographs from Wikimedia Commons float another digest on 1 of 200 images, 1 bit apart

The synthetic image where single precision decides 108 bits is four flat squares of color. In a window inside a flat area the nine block sums are equal, every comparison with their mean is a tie, and rounding alone breaks it. The folded sums break such ties the same way on every run; a direct evaluation in single precision breaks them by the order of its additions.

Digest layout

Bits run window by window, in the order of their numbers above, and within a window row by row: bit \(9w + 3r + c\) is row \(r\), column \(c\) of window \(w\). They are packed most significant bit first, eight to a byte, so the first window takes all of byte 0 and the top bit of byte 1.

The first 18 bits: the nine of window 0 and the nine of window 1, with the bytes they fall into The first 18 bits: the nine of window 0 and the nine of window 1, with the bytes they fall into

The order is pHash's; neither paper specifies one. Two digests are compared bit by bit, so only the count of differing bits matters, not where they are.

What changes the hash

The hash records where fine detail lies, block by block, to within a few pixels. Anything that moves the detail across the grid changes it; anything that keeps the detail in place changes it little.

Bits that differ from the original under nine transforms, on the example photograph Bits that differ from the original under nine transforms, on the example photograph

  • Rotation and cropping move the hash furthest: a turn of a degree or two, or a crop of a few percent, moves the detail by a few pixels at the image's edges, as much as the kernel's center is wide, and the block sums near those edges change. A few degrees more, or a larger crop, and the digest is as far from the original as an unrelated image's, about half of its 576 bits.
  • Brightness, contrast and gamma keep the order of the gray levels, which equalizing turns into the same image up to rounding, until the brightest pixels clip. They move a few percent of the bits.
  • Recompression, noise, blur and downscaling change the fine detail itself, which is what the kernel measures; the hash moves with their strength, a few percent of the bits at moderate settings.

These are measurements on one photograph; the same edits over two whole corpora follow below.

The numbers behind the chart
Transform Strength Bits that differ
JPEG quality 95 1
JPEG quality 75 0
JPEG quality 50 35
JPEG quality 30 58
JPEG quality 15 69
JPEG quality 5 106
Downscale 0.75 17
Downscale 0.5 28
Downscale 0.25 90
Downscale 0.125 151
Rotation 1 106
Rotation 2 170
Rotation 5 268
Rotation 10 291
Rotation 20 280
Rotation 45 290
Brightness 0.5 26
Brightness 0.7 30
Brightness 0.85 23
Brightness 1.15 17
Brightness 1.3 56
Brightness 1.5 71
Contrast 0.5 26
Contrast 0.7 19
Contrast 0.85 17
Contrast 1.15 15
Contrast 1.3 24
Contrast 1.5 39
Gamma 0.5 30
Gamma 0.7 24
Gamma 0.85 22
Gamma 1.2 19
Gamma 1.5 17
Gamma 2 20
Gaussian blur 0.5 4
Gaussian blur 1 12
Gaussian blur 2 32
Gaussian blur 4 83
Gaussian blur 8 154
Noise 2 18
Noise 5 35
Noise 10 49
Noise 20 70
Noise 40 93
Crop 5 200
Crop 10 295
Crop 20 270
Crop 30 281
How this was measured

Each transform is applied alone to the photograph above, and the edited copy is saved losslessly, as PPM (as a JPEG of the given quality for that panel). The library hashes the original and every copy, and compares each copy with the original by ph_similarity_digest(). The chart shows \(576 \times (1 - \text{similarity})\), which is the Hamming distance: how many of the 576 bits the edit flipped. Below is the code that ran, not a copy of it.

The edits, in Python with Pillow:

tools/site/measure/transforms.py
def transforms():
    """(name, x-axis label, [(strength, PIL image -> (image, (file extension, quality)))]).

    Each edit is applied alone to the original. A variant is saved as PPM, which is
    lossless (and quick to write for a corpus of thousands of variants), so the edit is
    the only difference; JPEG quality is the one edit that is the saving itself.
    """
    def jpeg(q):
        return lambda im: (im, ("jpg", q))

    def scale(s):
        return lambda im: (im.resize((max(1, round(im.width * s)), max(1, round(im.height * s))),
                                     Image.LANCZOS), ("ppm", None))

    def rotate(deg):
        # The corners the turn uncovers are filled with the image's mean color, so they
        # add no edge of their own.
        def f(im):
            fill = tuple(int(v) for v in np.asarray(im).reshape(-1, 3).mean(axis=0))
            return im.rotate(deg, resample=Image.BICUBIC, fillcolor=fill), ("ppm", None)
        return f

    def brightness(k):
        return lambda im: (ImageEnhance.Brightness(im).enhance(k), ("ppm", None))

    def contrast(k):
        return lambda im: (ImageEnhance.Contrast(im).enhance(k), ("ppm", None))

    def gamma(g):
        def f(im):
            a = (np.asarray(im).astype(np.float64) / 255.0) ** g
            return Image.fromarray((a * 255.0).round().astype(np.uint8)), ("ppm", None)
        return f

    def blur(r):
        return lambda im: (im.filter(ImageFilter.GaussianBlur(r)), ("ppm", None))

    def noise(sigma):
        # A fixed seed: the same noise on every build, so the chart does not move.
        def f(im):
            rng = np.random.default_rng(1)
            a = np.asarray(im).astype(np.float64) + rng.normal(0, sigma, (im.height, im.width, 3))
            return Image.fromarray(np.clip(a, 0, 255).round().astype(np.uint8)), ("ppm", None)
        return f

    def crop(p):
        def f(im):
            dx, dy = round(im.width * p / 2), round(im.height * p / 2)
            return im.crop((dx, dy, im.width - dx, im.height - dy)), ("ppm", None)
        return f

    return [
        ("JPEG quality", "quality", [(q, jpeg(q)) for q in (95, 75, 50, 30, 15, 5)]),
        ("Downscale", "scale", [(s, scale(s)) for s in (0.75, 0.5, 0.25, 0.125)]),
        ("Rotation", "degrees", [(d, rotate(d)) for d in (1, 2, 5, 10, 20, 45)]),
        ("Brightness", "factor", [(k, brightness(k)) for k in (0.5, 0.7, 0.85, 1.15, 1.3, 1.5)]),
        ("Contrast", "factor", [(k, contrast(k)) for k in (0.5, 0.7, 0.85, 1.15, 1.3, 1.5)]),
        ("Gamma", "exponent", [(g, gamma(g)) for g in (0.5, 0.7, 0.85, 1.2, 1.5, 2.0)]),
        ("Gaussian blur", "radius, px", [(r, blur(r)) for r in (0.5, 1, 2, 4, 8)]),
        ("Noise", "σ, gray levels", [(s, noise(s)) for s in (2, 5, 10, 20, 40)]),
        ("Crop", "border removed, %", [(p, crop(p / 100)) for p in (5, 10, 20, 30)]),
    ]

Writing the copies and handing them to the measuring tool:

tools/site/measure/robustness.py
def measure_variants(tool, base, ref):
    """{algorithm: {transform: [(strength, value)]}}: every edit of `base` (a Pillow
    image), of both groups of transforms.py, compared with `ref`, the same pixels saved
    losslessly.

    Every variant is written to a file, and `site_stages measure` hashes them with the
    library and compares each variant with the reference, by each algorithm's own
    metric; `value` is None where the comparison does not apply.
    """
    with tempfile.TemporaryDirectory() as tmp:
        files, keys = [], []
        for name, _, steps in transforms() + content_edits():
            for strength, op in steps:
                files.append(write_edit(op, base, os.path.join(tmp, str(len(files)))))
                keys.append((name, strength))
        rows = run_lines(tool, "measure", ref, *files)
    result = {}
    for (name, strength), row in zip(keys, rows):
        for algo, value in row.items():
            if algo != "file":
                result.setdefault(algo, {}).setdefault(name, []).append((strength, value))
    return result


def measure_robustness(tool, image):
    """measure_variants() for one image file."""
    base = Image.open(image).convert("RGB")
    with tempfile.TemporaryDirectory() as tmp:
        ref = os.path.join(tmp, "reference.ppm")
        base.save(ref)
        return measure_variants(tool, base, ref)

The measuring tool, linked against the library: it computes every algorithm's digest of the original, then of each copy, and compares the two by that algorithm's own metric.

tools/site/stages/measure.c
/* site_stages measure <reference> <variant>...: every algorithm's digest of the reference,
 * then one JSON line per variant with its comparison to the reference, per algorithm (the
 * similarity in [0, 1] for a bit hash, the peak correlation for Radial, the histogram
 * intersection for ColorHash, the L2 distance for ColorMoments). A `--load=<settings>`
 * between the images loads the ones after it so (take_load_settings()). */
int mode_measure(int argc, char **argv) {
    ph_digest_t ref[PH_ALGORITHM_COUNT], d[PH_ALGORITHM_COUNT];
    int have_ref[PH_ALGORITHM_COUNT], have[PH_ALGORITHM_COUNT];
    int status = 0, images = 0;
    for (int v = 0; v < argc && !status; v++) {
        if (take_load_settings(argv[v], &status)) {
            continue;
        }
        ph_context_t *ctx = NULL;
        if (load_image(&ctx, argv[v])) {
            return 1;
        }
        digest_all(ctx, images ? d : ref, images ? have : have_ref);
        ph_free(ctx);
        if (images++ == 0) {
            continue;
        }
        json_t j = json_begin(stdout);
        json_string(&j, "file", argv[v]);
        compare_all(&j, ref, have_ref, d, have);
        json_end(&j);
    }
    return status;
}
tools/site/stages/measure.c
/* One algorithm's comparison of two digests by its own metric; 0 when it does not apply. */
static int compare(const ph_digest_t *a, const ph_digest_t *b, double *out) {
    switch (a->kind) {
        case PH_DIGEST_KIND_BITS:
            *out = ph_similarity_digest(a, b);
            return *out >= 0.0;
        case PH_DIGEST_KIND_COEFFICIENTS:
            return ph_radial_similarity(a, b, out) == PH_SUCCESS;
        case PH_DIGEST_KIND_HISTOGRAM:
            return ph_histogram_intersection(a, b, out) == PH_SUCCESS;
        case PH_DIGEST_KIND_VECTOR16:
            *out = ph_l2_distance(a, b);
            return *out >= 0.0;
        default:
            return 0;
    }
}

The number on the chart:

tools/site/measure/metric.py
def bits_that_differ(similarity, bits_total):
    """The number on a bit hash's chart. ph_similarity_digest() is 1 - distance / bits,
    so this is the Hamming distance itself: how many of the hash's bits the edit flipped."""
    return (1.0 - similarity) * bits_total

To reproduce, make site builds the tool and redraws every figure. The tool can also be run by hand on any pair of images; it prints one JSON line per variant, with each algorithm's comparison:

cmake --preset release && cmake --build --preset release --target site_stages
build/release/site_stages measure reference.png variant.png …

Over two corpora

The same nine edits, applied to every image of two corpora: 200 public-domain photographs from Wikimedia Commons (the list), and the 24 generated images the library's tests measure (methodology), which are not photographs. Each line is the median over the corpus and the band around it holds the middle half of its images.

Bits that differ from the original under nine transforms, median and middle half over each corpus Bits that differ from the original under nine transforms, median and middle half over each corpus

The photographs repeat the single image: rotation and cropping reach half of the bits at a few degrees and a few percent, the other edits stay within a small share of them. The synthetic images also react to noise. Most of them are flat areas of a few colors, whose narrow histograms equalizing stretches over the whole range, and with them any noise on the flat areas: the response then answers to the noise.

Copies and different images

Robustness alone proves little: a hash that never changes is perfectly robust. What makes a hash useful is that copies of an image land closer together than different images do. Here a copy is an original after one moderate edit (one strength of each of the nine, listed under "How this was measured"), and different images are every pair of distinct originals in the corpus.

Distribution of bits that differ, for copies and for pairs of different images, in both corpora Distribution of bits that differ, for copies and for pairs of different images, in both corpora

Different images gather tightly around 288, half of the 576 bits: with that many bits, two unrelated digests differ in very nearly half of them. Copies split in two. Most land within a few dozen bits, but the rotated and cropped copies land near half, among the different images. The threshold that accepts 95 % of the copies has to reach up to them, and d′, the gap between the two means in units of their spread, is the lowest of the bit hashes on this site. That the threshold still lets few different pairs through is the tight cluster's doing, not the copies'.

The numbers behind the charts
Corpus Copies Different pairs d′ Threshold Different pairs within it
synthetic corpus, 24 images; not photographs 216 276 2.94 255 6.9%
200 public-domain photographs from Wikimedia Commons 1800 19900 3.46 251 0.7%

Median bits that differ over each corpus, with the 25th and 75th percentiles:

Transform Strength synthetic corpus, 24 images; not photographs 200 public-domain photographs from Wikimedia Commons
JPEG quality 95 12.5 (2.5–48.5) 2 (1–5)
JPEG quality 75 39 (12.3–65) 8 (4–20)
JPEG quality 50 54.5 (15.2–91.5) 12 (7–31.2)
JPEG quality 30 70 (22–123.3) 18 (9–42)
JPEG quality 15 78 (1–169.5) 31.5 (18.7–63)
JPEG quality 5 108 (37.7–229.5) 71.5 (47–108)
Downscale 0.75 33 (20.5–54) 7 (4–10.2)
Downscale 0.5 68 (25.8–116) 16 (11–23)
Downscale 0.25 151 (59.5–224.5) 48 (38.5–62)
Downscale 0.125 208 (144.7–272.2) 101 (86–120)
Rotation 1 72.5 (30.5–124.3) 144 (125–159.3)
Rotation 2 111.5 (74.2–192) 214 (191.5–228)
Rotation 5 228.5 (124–258) 273 (259–286)
Rotation 10 257.5 (129.2–288.2) 283 (272.7–293)
Rotation 20 270.5 (186.7–287.3) 286 (279–295)
Rotation 45 278 (207.2–290.2) 287 (276–293)
Brightness 0.5 11 (0.7–40.7) 7 (3–22)
Brightness 0.7 8 (0–34.5) 6 (2–16.3)
Brightness 0.85 5.5 (0–31.5) 5 (2–13.2)
Brightness 1.15 6 (0–65) 8 (3.7–20.2)
Brightness 1.3 21 (0–134.7) 14 (7–34)
Brightness 1.5 35.5 (5.8–167.8) 25.5 (13–48.5)
Contrast 0.5 11 (0.7–49.3) 7 (3–20.2)
Contrast 0.7 10 (0–38.5) 6 (2–17)
Contrast 0.85 5.5 (0–32.5) 5.5 (2–13.2)
Contrast 1.15 7.5 (0–23.5) 8 (3–18)
Contrast 1.3 12 (0–57) 12 (5.8–24.3)
Contrast 1.5 27 (0.7–115.5) 18.5 (8–31)
Gamma 0.5 13 (0–135.8) 15 (9–25)
Gamma 0.7 6.5 (0–97.5) 10 (5–18)
Gamma 0.85 3.5 (0–69.3) 6 (3–16)
Gamma 1.2 3.5 (0–75.5) 7 (4–14.3)
Gamma 1.5 9 (0–153.2) 12.5 (8–20)
Gamma 2 22 (0–180.2) 19 (12.7–30.2)
Gaussian blur 0.5 4.5 (0–15.8) 3 (1–5)
Gaussian blur 1 25.5 (10–43.8) 7 (4–12)
Gaussian blur 2 83.5 (29.2–125) 18.5 (13–24.3)
Gaussian blur 4 169 (65.3–242) 46 (37.5–59)
Gaussian blur 8 208.5 (152.2–273.5) 99 (82–114.3)
Noise 2 91.5 (17.2–205.5) 7 (3–24.3)
Noise 5 132.5 (34.5–235.5) 11 (5–36)
Noise 10 170 (60.5–254) 17 (9–52.3)
Noise 20 196.5 (109–259.3) 31 (17.8–74.7)
Noise 40 237 (148.5–265.5) 53 (34–102)
Crop 5 173.5 (94.3–230.3) 247 (229.3–259)
Crop 10 240.5 (136.5–277.5) 286 (269.5–298)
Crop 20 275 (157.3–305.3) 284.5 (273–293)
Crop 30 265 (178.2–288) 286 (276–295)
How this was measured

For each corpus, every image goes through the nine edits above, and each copy is compared with its original exactly as on the single photograph. Every pair of distinct originals is compared too, by site_stages pairs. Below is the code that ran.

The photographs are downloaded once and checked against the SHA-256 the manifest records for each:

tools/site/fetch_corpus.py
def download(row, unreachable):
    """The photograph's bytes as the manifest names them, from the copy of the corpus or
    else from Commons; (None, why) when neither has them. A source that cannot be reached
    is added to `unreachable` and not tried again."""
    why = "no source"
    for source in (MIRROR + row["sha256"] + extension(row), row["url"]):
        host = urllib.parse.urlsplit(source).netloc
        if host in unreachable:
            continue
        try:
            data = get(source)
        except urllib.error.HTTPError as e:
            why = f"{host}: {e}"
            continue
        except OSError as e:
            unreachable[host] = str(e)
            why = f"cannot reach {host} ({e})"
            continue
        finally:
            if source == row["url"]:
                time.sleep(0.2)  # Commons asks clients to pace their requests
        if hashlib.sha256(data).hexdigest() == row["sha256"]:
            return data, None
        why = f"{host}: the download differs from the manifest"
    return None, why


def unpack_archive(rows, d):
    """Writes the rows' photographs that the copy's archive holds, verified, into d.
    An archive that cannot be downloaded or read writes nothing; the files are then
    downloaded one by one."""
    wanted = {row["sha256"] + extension(row): row for row in rows}
    try:
        req = urllib.request.Request(MIRROR_ARCHIVE, headers={"User-Agent": USER_AGENT})
        with urllib.request.urlopen(req, timeout=300) as r, \
                tarfile.open(fileobj=r, mode="r|gz") as tar:
            for member in tar:
                row = wanted.get(os.path.basename(member.name))
                if row is None or not member.isfile():
                    continue
                data = tar.extractfile(member).read()
                if hashlib.sha256(data).hexdigest() != row["sha256"]:
                    continue
                path = os.path.join(d, local_name(row))
                with open(path + ".part", "wb") as f:
                    f.write(data)
                os.replace(path + ".part", path)
    except (OSError, tarfile.TarError) as e:
        print(f"warning: photo corpus: the archive of the copy: {e}", file=sys.stderr)


_fetched = {}


def fetch(cache=None):
    """The verified local paths of the manifest's photographs, downloading what is missing.
    Returns (paths, problems); a photograph that is unavailable is a problem, not an error.
    One process downloads once: a second call returns the first one's answer."""
    d = cache_dir(cache)
    if d in _fetched:
        return _fetched[d]
    os.makedirs(d, exist_ok=True)

    def cached(row):
        path = os.path.join(d, local_name(row))
        if not os.path.exists(path):
            return False
        with open(path, "rb") as f:
            return hashlib.sha256(f.read()).hexdigest() == row["sha256"]

    rows = read_manifest()
    missing = [row for row in rows if not cached(row)]
    if len(missing) > 10:
        unpack_archive(missing, d)
    paths, problems, unreachable = [], [], {}
    for row in rows:
        path = os.path.join(d, local_name(row))
        if row not in missing or cached(row):
            paths.append(path)
            continue
        data, why = download(row, unreachable)
        if data is None:
            problems.append(f"{row['file']}: {why}")
            continue
        with open(path + ".part", "wb") as f:
            f.write(data)
        os.replace(path + ".part", path)
        paths.append(path)
    _fetched[d] = (paths, problems)
    return paths, problems

Measuring a corpus, or taking the result from the cache when nothing that decides it has changed:

tools/site/measure/corpus.py
def measure_corpus(tool, corpus):
    """{"n", "label", "robust": {algo: {transform: [[strength, [value per image]]]}},
    "edits": {the same, for the edits that change the picture}, "different": {algo:
    [value per pair]}} for one corpus, from the cache when nothing that decides it has
    changed. Values are each algorithm's own metric, as `site_stages measure` and
    `site_stages pairs` print it; None where it does not apply."""
    paths = images(tool, corpus)

    def measure():
        with tempfile.TemporaryDirectory() as tmp:
            # The originals are decoded once, by Pillow, and written losslessly, so that
            # the reference an edit is compared with and the image that is edited are the
            # same pixels.
            refs = [os.path.join(tmp, f"{i:04d}.ppm") for i in range(len(paths))]
            with concurrent.futures.ProcessPoolExecutor() as pool:
                per_image = list(pool.map(_measure_one, [(tool, p, r) for p, r in
                                                         zip(paths, refs)]))
            rows = run_lines(tool, "pairs", *refs) if len(refs) >= 2 else []

        robust, edits = {}, {}
        content = {name for name, _, _ in content_edits()}
        for result in per_image:
            for algo, by_transform in result.items():
                for name, points in by_transform.items():
                    group = edits if name in content else robust
                    slot = group.setdefault(algo, {}).setdefault(
                        name, [[s, []] for s, _ in points])
                    for (strength, value), cell in zip(points, slot):
                        cell[1].append(value)
        different = {}
        for row in rows:
            for algo, value in row.items():
                if algo not in ("a", "b"):
                    different.setdefault(algo, []).append(value)
        return {"n": len(paths), "label": label(corpus, len(paths)), "robust": robust,
                "edits": edits, "different": different}

    return cached(corpus, cache_key(paths), measure, f"{corpus} corpus, {len(paths)} images")

Comparing every pair of originals:

tools/site/stages/measure.c
/* site_stages pairs <image>...: every image's digests, then one JSON line per pair of
 * distinct images (i < j, in the order given, as "a" and "b") with their comparison,
 * per algorithm. A `--load=<settings>` before the images loads them so. */
int mode_pairs(int argc, char **argv) {
    ph_digest_t(*dig)[PH_ALGORITHM_COUNT] = calloc((size_t)argc, sizeof(*dig));
    int (*have)[PH_ALGORITHM_COUNT] = calloc((size_t)argc, sizeof(*have));
    int bad = (!dig || !have) && fail("out of memory", NULL);
    int n = 0;
    for (int i = 0; i < argc && !bad; i++) {
        if (take_load_settings(argv[i], &bad)) {
            continue;
        }
        ph_context_t *ctx = NULL;
        bad = load_image(&ctx, argv[i]);
        if (!bad) {
            digest_all(ctx, dig[n], have[n]);
            n++;
            ph_free(ctx);
        }
    }
    for (int i = 0; i < n && !bad; i++) {
        for (int k = i + 1; k < n; k++) {
            json_t j = json_begin(stdout);
            json_int(&j, "a", i);
            json_int(&j, "b", k);
            compare_all(&j, dig[i], have[i], dig[k], have[k]);
            json_end(&j);
        }
    }
    free(dig);
    free(have);
    return bad;
}

The edits that make a copy, and the numbers on the second chart:

tools/site/measure/separability.py
# A "copy" is an original after one moderate edit: one strength of each transform, the
# kind of change a picture goes through between two places it is stored.
COPY_STRENGTHS = {
    "JPEG quality": 50,
    "Downscale": 0.5,
    "Rotation": 2,
    "Brightness": 1.15,
    "Contrast": 1.15,
    "Gamma": 1.2,
    "Gaussian blur": 1,
    "Noise": 5,
    "Crop": 5,
}
tools/site/measure/separability.py
def separability(copies, different, lower_is_closer, recall=RECALL):
    """d', the threshold that keeps `recall` of the copies (RECALL unless a page asks for
    another share), and the share of different pairs that threshold also accepts.

    d' is the gap between the two means in units of their pooled spread, as
    tests/src/test_hash_properties.c defines it (separability()), signed so that a
    positive d' means copies are closer than different images.
    """
    c, d = np.asarray(copies, float), np.asarray(different, float)
    sign = 1.0 if lower_is_closer else -1.0
    pooled = np.sqrt((c.std() ** 2 + d.std() ** 2) / 2.0)
    dprime = sign * (d.mean() - c.mean()) / pooled if pooled > 0 else float("inf")
    # The threshold: the closest value that still accepts `recall` of the copies.
    t = sign * np.quantile(sign * c, recall, method="inverted_cdf")
    accepted = (d <= t) if lower_is_closer else (d >= t)
    return dprime, t, accepted.mean()

The synthetic corpus is generated by make_base() in tests/src/synthetic_corpus.h, the same code the tests use. make site measures both corpora and draws every figure.

Edits that change the picture

The edits above are those a copy goes through. The ones here change what the picture shows: a patch of the image's own top-left corner pasted over its center, covering 1 to 16 % of the frame; the hue of every pixel turned, which recolors the picture; and a quarter turn, a half turn and a mirror image. Whether a hash should notice them depends on what it is for: a search for copies wants to ignore a recoloring, a search for retouched pictures wants to catch the patch.

The example photograph with a patch over 4 and 16 % of the frame, its hue turned by 90 and 180 degrees, turned by 90 and 180 degrees, and mirrored The example photograph with a patch over 4 and 16 % of the frame, its hue turned by 90 and 180 degrees, turned by 90 and 180 degrees, and mirrored

Each panel shows how far the edited images land from their originals over both corpora. The dashed line is the threshold of the chart above, the one that accepts 95 % of the copies, in each corpus's color: an edited image below it would be taken for a copy.

Bits that differ from the original for a patch, a hue rotation, and turns and a mirror, median and middle half over each corpus, against the threshold that accepts 95 % of copies Bits that differ from the original for a patch, a hue rotation, and turns and a mirror, median and middle half over each corpus, against the threshold that accepts 95 % of copies

  • Every patch passes for a copy. The threshold sits so high, to take in the rotated and cropped copies, that even a patch over 16 % of the frame lands far below it.
  • Turning the hue moves the hash little on photographs, whose muted colors keep their luminance under a hue rotation, and far on the flat, saturated colors of the synthetic images. Equalizing does not protect against it: a recoloring changes the order of the gray levels, not only their spacing.
  • A turn or a mirror image moves the hash as far as an unrelated image.
The numbers behind the chart

Median bits that differ, and the share of edited images within the threshold that accepts 95% of the copies; synthetic: synthetic corpus, 24 images; not photographs; photographs: 200 public-domain photographs from Wikimedia Commons.

Edit Synthetic: median (25th–75th) accepted as a copy Photographs: median (25th–75th) accepted as a copy
patch over 1% of the frame 9.5 (2.5–22.2) 100% 5 (3–8) 100%
patch over 4% of the frame 21.5 (17.5–29.2) 100% 19 (16.8–23) 100%
patch over 9% of the frame 47 (20.7–60.8) 100% 37 (33–44) 100%
patch over 16% of the frame 59.5 (49.3–68.5) 100% 59 (53.8–64) 100%
hue turned 30° 67 (4–225.5) 88% 31 (18–48) 100%
hue turned 90° 268.5 (87.7–429.8) 46% 25 (13–37) 100%
hue turned 180° 298 (106–402.8) 42% 49 (28.7–75.2) 100%
turned 90° 287.5 (119.5–297.2) 38% 284 (273–293) 8%
turned 180° 218.5 (82.2–328.5) 54% 288 (274.8–298.2) 10%
mirrored 273.5 (125–360.2) 46% 278 (257–294) 22%
How this was measured

The same measurement as for the copies, over the same images: each edit is applied alone to the original, the copy is saved losslessly, and the library compares it with the original. The threshold is the one of the chart above, computed from the copies. Below is the code of the edits.

tools/site/measure/transforms.py
def content_edits():
    """The edits that change the picture, in the shape of transforms(); a strength is a
    number, or a word where the edit has no scale.

    A patch replaces a square at the center, covering the given share of the frame, with
    the square of the same size at the top-left corner of the same image: a local edit
    with content the image itself supplies. A hue rotation turns every pixel's hue by the
    given angle and keeps its saturation and value, a recoloring. The turns and the mirror
    change only the orientation.
    """
    def patch(area):
        def f(im):
            side = area ** 0.5
            w, h = round(im.width * side), round(im.height * side)
            x0, y0 = (im.width - w) // 2, (im.height - h) // 2
            out = im.copy()
            out.paste(im.crop((0, 0, w, h)), (x0, y0))
            return out, ("ppm", None)
        return f

    def hue(degrees):
        def f(im):
            hsv = np.asarray(im.convert("HSV")).copy()
            hsv[..., 0] = (hsv[..., 0].astype(int) + round(degrees * 256 / 360)) % 256
            return Image.fromarray(hsv, "HSV").convert("RGB"), ("ppm", None)
        return f

    def turn(method):
        return lambda im: (im.transpose(method), ("ppm", None))

    return [
        ("Local patch", "share of the frame, %", [(a, patch(a / 100)) for a in (1, 4, 9, 16)]),
        ("Hue rotation", "degrees", [(d, hue(d)) for d in (30, 90, 180)]),
        ("Turn and mirror", "", [("90°", turn(Image.Transpose.ROTATE_90)),
                                 ("180°", turn(Image.Transpose.ROTATE_180)),
                                 ("mirror", turn(Image.Transpose.FLIP_LEFT_RIGHT))]),
    ]

A small edit against a rescale

A hash that should find edited copies of a picture has to move further for an edit than for the changes a copy goes through. Here are mHash's median distances for both, the edits of a copy at the strength that makes one, and the patches:

Median bits that differ for each edit of a copy and for patches over 1 to 16 % of the frame, sorted, over both corpora Median bits that differ for each edit of a copy and for patches over 1 to 16 % of the frame, sorted, over both corpora

On the photographs, a patch over 4 % of the frame moves mHash about as far as halving the image's size, and a 2-degree turn or a 5 % crop moves it several times further than a patch over 16 %. On the synthetic images, recompression, noise and downscaling move it further than small patches. The order is the wrong way round for finding edits: mHash describes where detail lies, and any small shift of the detail outweighs a local change of it. For aHash, dHash, pHash and wHash the order is the right one: halving the size moves none of their bits, a patch moves some.

How the nine algorithms compare is on choosing an algorithm.

Cost

mHash is the most expensive of the hashes on a small image. The grayscale conversion, the blur and the resampling read every pixel of the image and grow with it; the equalizing and the integral image work at 512×512 whatever the image; the block sums take a fixed 280 thousand multiply-adds. On a large image the blur, at full resolution, is most of the time.

How this was measured

Measured on Intel(R) Xeon(R) Platinum 8370C CPU @ 2.80GHz, 4 cores, Linux x86_64; compiled by 13.3.0, optimized, with libjpeg-turbo decoding the JPEGs; the minimum of the runs. Decoding, every hash at its default settings, and in bold this page's.

Case 400×400 20 Mpx
decode 0.53 ms 105 ms
ahash 0.21 ms 23 ms
dhash 0.23 ms 27 ms
phash 0.21 ms 23 ms
whash 0.21 ms 23 ms
bmh 0.21 ms 23 ms
mhash 1.8 ms 91 ms
radial 1.4 ms 120 ms
color_hash 0.24 ms 31 ms
color_moments 0.52 ms 66 ms
mhash_size_62 0.85 ms 92 ms
mhash_size_128 0.91 ms 94 ms
mhash_size_256 1.1 ms 93 ms
mhash_size_512 1.8 ms 92 ms
mhash_size_1024 4.7 ms 97 ms
mhash_size_2048 17 ms 109 ms
mhash_size_4096 66 ms 159 ms

Each case runs once to warm up, then until it has run at least five times and for at least a second; the time is the minimum. The grayscale image and the area grid the library caches are dropped before every run, so each run is the first hash on a loaded image. The cases, and the timing loop:

tools/site/stages/timing.c
/* Seconds on a monotonic clock, as tests/src/bench_hash.c reads it. The calendar clock
 * (timespec_get() with TIME_UTC) would step back when the system corrects its time, and a
 * run spanning the step would come out negative and become the minimum. On Darwin,
 * CLOCK_MONOTONIC counts whole microseconds, 2 % of the cheapest hash; mach_absolute_time()
 * counts nanoseconds. */
static double seconds(void) {
#ifdef __APPLE__
    static mach_timebase_info_data_t tb;
    if (tb.denom == 0) {
        mach_timebase_info(&tb);
    }
    return (double)mach_absolute_time() * tb.numer / tb.denom * 1e-9;
#else
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return (double)ts.tv_sec + (double)ts.tv_nsec * 1e-9;
#endif
}

/* One timed case: `run` does the work once on a loaded image and returns nonzero on
 * failure; `arg` selects a variant of it. A case named after an algorithm, or starting
 * with its name and "_", is shown on that algorithm's page. */
typedef struct {
    const char *name;
    int (*run)(ph_context_t *ctx, const char *path, int arg);
    int arg;
} time_case_t;

static int time_decode(ph_context_t *ctx, const char *path, int arg) {
    (void)arg;
    return ph_load_from_file(ctx, path) != PH_SUCCESS;
}

/* The first hash computed on a loaded image, which is what a caller who loads an image
 * and hashes it pays: the grayscale conversion and the shared area grid are dropped
 * before every run, so each run computes them (the dropping is untimed). */
static int time_hash(ph_context_t *ctx, const char *path, int algo) {
    (void)path;
    ph_digest_t d;
    return ph_compute_digest(ctx, (ph_algorithm_t)algo, &d) != PH_SUCCESS;
}

/* The four 64-bit hashes in one ph_compute_multi() call, which computes the grayscale
 * image and the area grid once for all of them. */
static int time_multi(ph_context_t *ctx, const char *path, int arg) {
    (void)path;
    (void)arg;
    uint64_t h[PH_HASH_FLAGS_COUNT];
    return ph_compute_multi(ctx, PH_HASH_AHASH | PH_HASH_DHASH | PH_HASH_PHASH | PH_HASH_WHASH,
                            h) != PH_SUCCESS;
}

static int time_whash_full(ph_context_t *ctx, const char *path, int arg) {
    (void)path;
    (void)arg;
    uint64_t h;
    int bad = ph_context_set_whash_mode(ctx, PH_WHASH_FULL) != PH_SUCCESS ||
              ph_compute_whash(ctx, &h) != PH_SUCCESS;
    return ph_context_set_whash_mode(ctx, PH_WHASH_FAST) != PH_SUCCESS || bad;
}

static int time_mhash_size(ph_context_t *ctx, const char *path, int size) {
    (void)path;
    ph_digest_t d;
    int bad = ph_context_set_mhash_params(ctx, PH_MH_ALPHA, PH_MH_LEVEL, size) != PH_SUCCESS ||
              ph_compute_mhash(ctx, &d) != PH_SUCCESS;
    return ph_context_set_mhash_params(ctx, PH_MH_ALPHA, PH_MH_LEVEL, PH_MH_IMAGE_SIZE) !=
               PH_SUCCESS ||
           bad;
}

/* pHash with a dct_size other than the default, and the default block. */
static int time_phash_dct(ph_context_t *ctx, const char *path, int size) {
    (void)path;
    uint64_t h;
    int bad = ph_context_set_phash_params(ctx, size, PH_DCT_REDUCTION_SIZE) != PH_SUCCESS ||
              ph_compute_phash(ctx, &h) != PH_SUCCESS;
    return ph_context_set_phash_params(ctx, PH_DCT_SIZE, PH_DCT_REDUCTION_SIZE) != PH_SUCCESS ||
           bad;
}

/* Radial with the settings of radial_variants[variant]. */
static int time_radial(ph_context_t *ctx, const char *path, int variant) {
    (void)path;
    ph_digest_t d;
    int bad = !radial_settings(ctx, &radial_variants[variant]) ||
              ph_compute_radial_hash(ctx, &d) != PH_SUCCESS;
    return !radial_settings(ctx, &radial_variants[RADIAL_DEFAULT]) || bad;
}

/* Loading the image, then, unless `algo` is LOAD_ONLY, computing that algorithm: what a
 * caller pays for a load and a hash. `arg` packs the load -- a ph_decode_scale_t, plus
 * LOAD_GRAY for ph_context_set_load_grayscale() and LOAD_FRESH for a context created and
 * freed around the load, as a caller who creates one per image does -- and the algorithm
 * above it. Without LOAD_FRESH, a context of its own per kind of load, kept across runs,
 * so the cases after it hash the image as it was first loaded. */
enum {
    LOAD_GRAY = 4,
    LOAD_FRESH = 8,
    LOAD_ALGO = 16,
    LOAD_ONLY = -1,
};

static int time_load(ph_context_t *ctx, const char *path, int arg) {
    (void)ctx;
    static ph_context_t *kept[LOAD_FRESH];
    int kind = arg % LOAD_FRESH, fresh = arg & LOAD_FRESH, algo = arg / LOAD_ALGO - 1;
    ph_context_t *c = fresh ? NULL : kept[kind];
    if (!c && ph_create(&c) != PH_SUCCESS) {
        return 1;
    }
    ph_digest_t d;
    int bad = ph_context_set_decode_scale(c, (ph_decode_scale_t)(kind & 3)) != PH_SUCCESS ||
              ph_context_set_load_grayscale(c, (kind & LOAD_GRAY) != 0) != PH_SUCCESS ||
              ph_load_from_file(c, path) != PH_SUCCESS ||
              (algo != LOAD_ONLY && ph_compute_digest(c, (ph_algorithm_t)algo, &d) != PH_SUCCESS);
    if (fresh) {
        ph_free(c);
    } else {
        kept[kind] = c;
    }
    return bad;
}

#define LOADED(name, flags, algo) {name, time_load, (flags) + ((algo) + 1) * LOAD_ALGO}

static const time_case_t time_cases[] = {
    {"decode", time_decode, 0},
    {"ahash", time_hash, PH_ALGO_AHASH},
    {"dhash", time_hash, PH_ALGO_DHASH},
    {"phash", time_hash, PH_ALGO_PHASH},
    {"whash", time_hash, PH_ALGO_WHASH},
    {"whash_full", time_whash_full, 0},
    {"bmh", time_hash, PH_ALGO_BMH},
    {"mhash", time_hash, PH_ALGO_MHASH},
    {"radial", time_hash, PH_ALGO_RADIAL},
    {"color_hash", time_hash, PH_ALGO_COLOR_HASH},
    {"color_moments", time_hash, PH_ALGO_COLOR_MOMENTS},
    {"multi", time_multi, 0},
    {"mhash_size_62", time_mhash_size, 62},
    {"mhash_size_128", time_mhash_size, 128},
    {"mhash_size_256", time_mhash_size, 256},
    {"mhash_size_512", time_mhash_size, 512},
    {"mhash_size_1024", time_mhash_size, 1024},
    {"mhash_size_2048", time_mhash_size, 2048},
    {"mhash_size_4096", time_mhash_size, 4096},
    {"phash_dct_8", time_phash_dct, 8},
    {"phash_dct_16", time_phash_dct, 16},
    {"phash_dct_24", time_phash_dct, 24},
    {"radial_sigma_1", time_radial, RADIAL_SIGMA_1},
    {"radial_sigma_8", time_radial, RADIAL_SIGMA_8},
    {"radial_gamma_2", time_radial, RADIAL_GAMMA_2},
    {"radial_grid_40x32", time_radial, RADIAL_GRID_40X32},
    {"radial_grid_90x64", time_radial, RADIAL_GRID_90X64},
    {"radial_grid_360x256", time_radial, RADIAL_GRID_360X256},
    {"radial_grid_1440x1024", time_radial, RADIAL_GRID_1440X1024},
    {"radial_grid_4096x4096", time_radial, RADIAL_GRID_4096X4096},
    LOADED("scale_half_decode", PH_DECODE_SCALE_HALF, LOAD_ONLY),
    LOADED("scale_quarter_decode", PH_DECODE_SCALE_QUARTER, LOAD_ONLY),
    LOADED("scale_eighth_decode", PH_DECODE_SCALE_EIGHTH, LOAD_ONLY),
    LOADED("scale_full_phash", PH_DECODE_SCALE_FULL, PH_ALGO_PHASH),
    LOADED("scale_half_phash", PH_DECODE_SCALE_HALF, PH_ALGO_PHASH),
    LOADED("scale_quarter_phash", PH_DECODE_SCALE_QUARTER, PH_ALGO_PHASH),
    LOADED("scale_eighth_phash", PH_DECODE_SCALE_EIGHTH, PH_ALGO_PHASH),
    LOADED("scale_full_mhash", PH_DECODE_SCALE_FULL, PH_ALGO_MHASH),
    LOADED("scale_half_mhash", PH_DECODE_SCALE_HALF, PH_ALGO_MHASH),
    LOADED("scale_quarter_mhash", PH_DECODE_SCALE_QUARTER, PH_ALGO_MHASH),
    LOADED("scale_eighth_mhash", PH_DECODE_SCALE_EIGHTH, PH_ALGO_MHASH),
    LOADED("scale_full_radial", PH_DECODE_SCALE_FULL, PH_ALGO_RADIAL),
    LOADED("scale_half_radial", PH_DECODE_SCALE_HALF, PH_ALGO_RADIAL),
    LOADED("scale_quarter_radial", PH_DECODE_SCALE_QUARTER, PH_ALGO_RADIAL),
    LOADED("scale_eighth_radial", PH_DECODE_SCALE_EIGHTH, PH_ALGO_RADIAL),
    LOADED("gray_decode", LOAD_GRAY, LOAD_ONLY),
    LOADED("gray_phash", LOAD_GRAY, PH_ALGO_PHASH),
    LOADED("gray_mhash", LOAD_GRAY, PH_ALGO_MHASH),
    LOADED("gray_radial", LOAD_GRAY, PH_ALGO_RADIAL),
    LOADED("fresh_decode", LOAD_FRESH, LOAD_ONLY),
    LOADED("fresh_phash", LOAD_FRESH, PH_ALGO_PHASH),
};

static int cmp_double(const void *a, const void *b) {
    double x = *(const double *)a, y = *(const double *)b;
    return (x > y) - (x < y);
}

/* One warm-up run, then at least TIME_MIN_RUNS runs and TIME_MIN_SECONDS, at most
 * TIME_MAX_RUNS; the minimum is what the pages quote. */
enum {
    TIME_MIN_RUNS = 5,
    TIME_MAX_RUNS = 300,
};

#define TIME_MIN_SECONDS 1.0

/* Whether case `name` is among the `n` names of `only`; every case is when `n` is 0. */
static int time_chosen(const char *name, int n, char **only) {
    for (int i = 0; i < n; i++) {
        if (strcmp(name, only[i]) == 0) {
            return 1;
        }
    }
    return n == 0;
}

/* site_stages time <image> [case...]: every case of time_cases[] on the image, or the cases
 * named; the minimum and the median of its runs, in milliseconds, as one JSON object, with
 * the build that ran them. */
int mode_time(int argc, char **argv) {
    for (int i = 1; i < argc; i++) {
        int known = 0;
        for (size_t c = 0; c < COUNT(time_cases); c++) {
            known |= strcmp(argv[i], time_cases[c].name) == 0;
        }
        if (!known) {
            return fail("no such timed case", argv[i]);
        }
    }
    ph_context_t *ctx = NULL;
    if (load_image(&ctx, argv[0])) {
        return 1;
    }
    json_t j = json_begin(stdout);
    json_int(&j, "width", ctx->image.width);
    json_int(&j, "height", ctx->image.height);
    json_string(&j, "build_info", ph_get_build_info());
#ifdef __VERSION__
    json_string(&j, "compiler", __VERSION__);
#endif
    json_t cases = json_object(&j, "cases");
    static double ms[TIME_MAX_RUNS];
    for (size_t c = 0; c < COUNT(time_cases); c++) {
        const time_case_t *tc = &time_cases[c];
        if (!time_chosen(tc->name, argc - 1, argv + 1)) {
            continue;
        }
        int runs = 0;
        double total = 0.0;
        for (int k = -1; k < TIME_MAX_RUNS; k++) { /* k = -1: the warm-up */
            ph_drop_gray_cache(ctx);
            double t0 = seconds();
            if (tc->run(ctx, argv[0], tc->arg)) {
                ph_free(ctx);
                return fail("timed case failed", tc->name);
            }
            double t = seconds() - t0;
            if (k >= 0) {
                ms[runs++] = t * 1e3;
                total += t;
                if (runs >= TIME_MIN_RUNS && total >= TIME_MIN_SECONDS) {
                    break;
                }
            }
        }
        qsort(ms, (size_t)runs, sizeof(ms[0]), cmp_double);
        json_t o = json_object(&cases, tc->name);
        json_double(&o, "min_ms", ms[0]);
        json_double(&o, "median_ms", ms[runs / 2]);
        json_int(&o, "runs", runs);
        json_close_object(&o);
    }
    json_close_object(&cases);
    json_end(&j);
    ph_free(ctx);
    return 0;
}

The times are measured again only when the library, the tool or the machine changes:

tools/site/measure/timing.py
def measure_timing(tool):
    """{"machine": {...}, "small": {...}, "large": {...}, "scan": {...}}, each image's
    entry as `site_stages time` prints it (width, height, build_info, compiler, and per
    case the minimum and median in milliseconds) and the comparisons of a linear search as
    `site_stages scan` prints them (per case, nanoseconds per comparison), from the cache
    when nothing that decides them has changed, the machine included."""
    host = machine()
    paths = [os.path.join(ROOT, p) for p in IMAGES.values()]
    key = fingerprint(KEY_FILES, [tool] + paths, json.dumps(host, sort_keys=True))

    def measure():
        data = {"machine": host}
        for name, path in zip(IMAGES, paths):
            out = subprocess.run([tool, "time", path], check=True, capture_output=True,
                                 text=True).stdout
            data[name] = json.loads(out)
        data["scan"] = json.loads(subprocess.run([tool, "scan"], check=True,
                                                 capture_output=True, text=True).stdout)
        return data

    return cached("timing", key, measure, "timing")

Parameters

ph_context_set_mhash_params() sets three values.

Parameter Range Default Meaning
alpha above 1 2 the base of the kernel's scale
level 0 or above 1 its exponent: the kernel reaches \(\lfloor 4\alpha^{\ell} \rfloor\) pixels from its center, at most 32 (a 65×65 kernel)
size 62 to 4096 512 the side of step 3; the grid stays 31×31, so a block is \(\lfloor \text{size}/31 \rfloor\) pixels

The digest is 576 bits whatever the values, and digests are comparable only when computed with the same ones. What the hash sees is decided by one ratio: the kernel's scale against the block's size, both in pixels of the normalized image.

The kernel's scale

Raising alpha or level widens the kernel, which then answers to coarser detail. Here are levels 0, 1 and 2 at alpha 2, their kernels drawn to the same scale, and the block sums each gives:

Kernels of 9×9, 17×17 and 33×33 at levels 0, 1 and 2, drawn to scale, and the block sums each gives, with the bits each digest differs from the default's Kernels of 9×9, 17×17 and 33×33 at levels 0, 1 and 2, drawn to scale, and the block sums each gives, with the bits each digest differs from the default's

The block sums keep their overall pattern, the cactus against the plain wall, and change in detail; a few dozen bits move each way.

The size

size changes the block, not the kernel: at 256 a block is 8 pixels, at 1024 it is 33, against a kernel that stays 17 pixels wide. The coarser the kernel against the block, the more the block sums follow the picture's broad layout; the finer, the more they follow texture. At 1024 the digest of the example is nearly as far from the default's as an unrelated image's.

Block sums at sizes 256, 512 and 1024, with the bits each digest differs from the default's Block sums at sizes 256, 512 and 1024, with the bits each digest differs from the default's

The work that depends on size grows with its square; the blur, which does not, is the larger part on a large image. Here is the time against size on the two timed images:

Time of mHash against size on a 400×400 and a 20-megapixel image, both axes logarithmic Time of mHash against size on a 400×400 and a 20-megapixel image, both axes logarithmic

On the small image the time grows slowly up to the default and roughly with the square of size above it; on the large one the blur's cost hides it up to a size several times the default.

The numbers behind the chart
size 400×400 5472×3648
62 0.85 ms 92 ms
128 0.91 ms 94 ms
256 1.1 ms 93 ms
512 1.8 ms 92 ms
1024 4.7 ms 97 ms
2048 17 ms 109 ms
4096 66 ms 159 ms

The defaults are pHash's, which fixes the size at 512 and does not expose it. Across 24 settings of level and size the property tests' separability moves within a narrow range with no trend in either, and no setting beats the source's (provenance has the sweep).

Settings that affect it

Five context settings change the pixels mHash reads.

Setting Effect on mHash
ph_context_set_gray_weights() the grayscale formula of step 1
ph_context_set_load_grayscale() a JPEG decoder converts to grayscale itself (measured), which moves gray levels by one here and there
ph_context_set_auto_orient() on by default: the hash describes the image as displayed, after its EXIF rotation
ph_context_set_alpha_mode() how transparent pixels are composited before grayscale (which background)
ph_context_set_decode_scale() a JPEG read by libjpeg-turbo is decoded at ½, ¼ or ⅛ of its size (measured); the blur's σ is fixed in pixels, so a smaller decode changes what mHash reads, and the digest moves at every scale, toward its threshold for copies

A gray level here and there reaches the response through the kernel, and the block sums of fine detail are sensitive to it:

On the example photograph the two digests are 14 of 576 bits apart.

In code

mHash returns a digest, like BMH, Radial and the two color hashes. This example computes all five for two images and compares each pair with the function its kind calls for; mHash's goes through ph_hamming_distance_digest():

examples/digest_and_metrics.c
/* Hashes two images with every digest algorithm, compares each pair with the function
 * that digest's kind calls for, shows what a mismatched comparison returns, and stores a
 * digest as text and reads it back. Build/run:
 *   cc digest_and_metrics.c -o digest_and_metrics $(pkg-config --cflags --libs libphash)
 *   ./digest_and_metrics a.jpg b.jpg
 */
#include <libphash.h>
#include <stdio.h>
#include <string.h>

/* Loads `path` into `ctx` in color, so the color hashes can run on it too. */
static int load(ph_context_t *ctx, const char *path) {
    ph_error_t err = ph_load_from_file(ctx, path);
    if (err != PH_SUCCESS) {
        fprintf(stderr, "failed to load %s: %s (%s)\n", path, ph_get_error_string(err),
                ph_get_last_error_message(ctx));
        return 0;
    }
    return 1;
}

/* Compares two digests of one algorithm with the function its kind calls for. */
static void compare(const char *name, const ph_digest_t *a, const ph_digest_t *b) {
    double score = 0.0;
    switch ((ph_digest_kind_t)a->kind) {
        case PH_DIGEST_KIND_BITS:
            printf("%-13s %3d of %d bits differ, similarity %.3f\n", name,
                   ph_hamming_distance_digest(a, b), a->size * 8, ph_similarity_digest(a, b));
            break;
        case PH_DIGEST_KIND_COEFFICIENTS:
            if (ph_radial_similarity(a, b, &score) == PH_SUCCESS) {
                printf("%-13s peak correlation %.3f (same image at %.1f and above)\n", name, score,
                       PH_RADIAL_PCC_THRESHOLD);
            } else {
                printf("%-13s no angular structure to compare\n", name);
            }
            break;
        case PH_DIGEST_KIND_HISTOGRAM:
            if (ph_histogram_intersection(a, b, &score) == PH_SUCCESS) {
                printf("%-13s histogram intersection %.3f\n", name, score);
            }
            break;
        case PH_DIGEST_KIND_VECTOR16:
        case PH_DIGEST_KIND_VECTOR:
            printf("%-13s L2 distance %.2f\n", name, ph_l2_distance(a, b));
            break;
        default:
            printf("%-13s kind %d: no comparison\n", name, a->kind);
            break;
    }
}


int main(int argc, char **argv) {
    if (argc != 3) {
        fprintf(stderr, "usage: %s <image-a> <image-b>\n", argv[0]);
        return 1;
    }

    static const ph_algorithm_t algos[] = {PH_ALGO_MHASH, PH_ALGO_BMH, PH_ALGO_RADIAL,
                                           PH_ALGO_COLOR_HASH, PH_ALGO_COLOR_MOMENTS};

    enum {
        N_ALGOS = sizeof(algos) / sizeof(algos[0]),
    };

    ph_digest_t da[N_ALGOS], db[N_ALGOS];

    ph_context_t *ctx = NULL;
    if (ph_create(&ctx) != PH_SUCCESS) {
        fprintf(stderr, "ph_create failed\n");
        return 1;
    }

    /* One context, reused: each load replaces the previous image and its cached work. */
    for (int img = 0; img < 2; img++) {
        if (!load(ctx, argv[1 + img])) {
            ph_free(ctx);
            return 1;
        }
        for (int i = 0; i < N_ALGOS; i++) {
            ph_error_t err = ph_compute_digest(ctx, algos[i], img == 0 ? &da[i] : &db[i]);
            if (err != PH_SUCCESS) {
                fprintf(stderr, "%s failed: %s\n", ph_algorithm_name(algos[i]),
                        ph_get_error_string(err));
                ph_free(ctx);
                return 1;
            }
        }
    }

    for (int i = 0; i < N_ALGOS; i++) {
        compare(ph_algorithm_name(algos[i]), &da[i], &db[i]);
    }

    /* The wrong metric is refused, not answered with a meaningless number: a radial
     * digest is quantized coefficients, not a bit vector. */
    printf("hamming distance of two radial digests: %d (refused)\n",
           ph_hamming_distance_digest(&da[2], &db[2]));

    /* Store a digest as text and read it back; the text carries the kind. */
    char hex[PH_DIGEST_HEX_BUFFER_SIZE];
    ph_digest_t back;
    if (ph_digest_to_hex(&da[1], hex, sizeof(hex)) != PH_SUCCESS ||
        ph_digest_from_hex(hex, &back) != PH_SUCCESS) {
        fprintf(stderr, "hex round trip failed\n");
        ph_free(ctx);
        return 1;
    }
    printf("BMH as text: %s\n", hex);
    printf("read back: %s\n", (back.size == da[1].size && back.kind == da[1].kind &&
                               memcmp(back.data, da[1].data, back.size) == 0)
                                  ? "identical"
                                  : "DIFFERENT");

    ph_free(ctx);
    return 0;
}

Where it comes from

Two sources cover two halves of mHash. The operator is the Laplacian of Gaussian of D. Marr and E. Hildreth, "Theory of edge detection" (1980), sampled as the Mexican hat. The hash built on it has no paper: it is ph_mh_imagehash() of the pHash library, by Evan Klinger and David Starkweather, and Zauner's thesis (2010), which describes it, says the construction "has not been proposed previously". pHash's code is therefore the primary source for every step and constant on this page: the blur at 1, the size 512, the 256 levels, the scale \(\alpha = 2\), \(\ell = 1\), the 16-pixel blocks, the 3×3 windows at stride 4.

Despite the name, neither pHash nor this library detects edges as Marr and Hildreth do, by the zero-crossings of the response: the response is summed over blocks and compared within windows. The operator is theirs; the edge detector is not implemented.

The digest is not bit for bit pHash's. pHash blurs with Deriche's recursive filter and resamples with a quintic filter, where this library uses a direct Gaussian and the Mitchell filter; it also computes the block sums the long way. The construction and every parameter are reproduced; the arithmetic is not (provenance § 5).


  1. Measured when this site was built, on Intel(R) Xeon(R) Platinum 8370C CPU @ 2.80GHz, 4 cores, Linux x86_64; compiled by 13.3.0, optimized, with libjpeg-turbo decoding the JPEGs: the minimum of at least five runs after a warm-up. What each time covers, and the code, are under “How this was measured” beside it. ↩