Storing and searching hashes¶

A hash is computed once and compared many times: it is stored with the record of its image, and a new image's hash is searched for among the stored ones. This page says what has to be stored beside a hash for it to stay comparable, how each algorithm's hash is stored, what a search by comparing with every stored hash costs, measured, and how an index cuts that cost when it grows too large. The library compares two hashes; storage and indexes belong to the program, and the library has neither.

What a hash does not record¶

A uint64_t hash or a ph_digest_t holds the hash's bits and nothing about how they were computed. Two hashes are comparable only when four things are the same for both, and none of them can be read from the hash afterwards:

• The algorithm. A digest's kind refuses a comparison across kinds, but six algorithms produce bits, and the four 64-bit hashes have no kind at all: an aHash compared with a pHash returns a distance like any other. ph_algorithm_name() gives each algorithm a stable name to store, and ph_algorithm_from_name() reads it back.
• The settings that change the value, every one except the size limit (Configuring a context lists them).
• The library version. A release that changes what an algorithm computes says so in the changelog, and names the algorithms whose stored hashes are to be recomputed. ph_version_number() is one integer to store.
• The JPEG decoder. libjpeg-turbo and stb_image decode a JPEG to slightly different pixels, and hashes from the two are compared by distance, never by equality (Same hash on every machine). ph_can_use_jpeg() says which one the build has.

A collection computed in one run with one configuration stores these once, beside it, as the example below prints them on its first line; a collection that grows over years of upgrades stores them with every hash, or at least the version.

How each hash is stored¶

BMH's size is for its default 16×16 grid; every other size is fixed.

The four 64-bit hashes fit an integer column. In a column of signed 64-bit integers, a hash with its top bit set, above INT64_MAX, is stored as the negative number with the same bits: it reads back unchanged, and the XOR and bit count of a Hamming distance give the same answer on either. Copy the bits with memcpy() rather than converting the value, as C leaves the conversion of a uint64_t above INT64_MAX to int64_t to the implementation. As text, ph_hash_to_hex() writes 16 lowercase digits, the most significant first; two hashes are equal exactly when their texts are.

A digest is stored as the text of ph_digest_to_hex(), its kind followed by its bytes, and read back by ph_digest_from_hex() with the kind restored (Hashes as text). Stored as bytes instead, it takes two columns, the kind and the size bytes of data, not the whole struct: a ph_digest_t takes 136 bytes whatever its size. Reading the bytes back means filling a zeroed ph_digest_t and setting both size and kind; a digest whose kind is left at zero is compared with any other without the check.

Finding near duplicates¶

The example hashes every file given to it with pHash, on all processors, prints each hash as it would be stored, and then compares every pair, printing those within the threshold given on its command line, closest first. A file that cannot be read is reported and left out, and the others are compared.

/* Finds the pictures that appear more than once among many files: hashes every file with
 * pHash, prints each hash as it would be stored, then every pair within a threshold,
 * closest first. Build/run:
 *   cc find_duplicates.c -o find_duplicates $(pkg-config --cflags --libs libphash)
 *   ./find_duplicates 7 a.jpg b.jpg c.png ...
 * The threshold is the most bits two hashes of one picture may differ in; measure it on
 * your own images (docs/theory/comparing.md).
 */
#include <libphash.h>
#include <stdio.h>
#include <stdlib.h>

typedef struct {
    int distance;
    size_t a, b;
} pair_t;

static int by_distance(const void *x, const void *y) {
    const pair_t *p = x, *q = y;
    return p->distance - q->distance;
}

int main(int argc, char **argv) {
    if (argc < 3) {
        fprintf(stderr, "usage: %s <max-bits> <image>...\n", argv[0]);
        return 1;
    }
    char *end;
    long max_bits = strtol(argv[1], &end, 10);
    if (*end != '\0' || max_bits < 0 || max_bits > 64) {
        fprintf(stderr, "max-bits must be a number from 0 to 64\n");
        return 1;
    }
    size_t n = (size_t)(argc - 2);
    ph_batch_item_t *items = calloc(n, sizeof(*items));
    if (!items) {
        return 1;
    }
    for (size_t i = 0; i < n; i++) {
        items[i].path = argv[i + 2];
    }

    /* Every file hashed with pHash on all CPUs; a file that cannot be read fails alone. */
    ph_error_t err = ph_hash_files(items, n, PH_HASH_PHASH, 0);
    if (err != PH_SUCCESS) {
        fprintf(stderr, "batch failed: %s\n", ph_get_error_string(err));
        free(items);
        return 1;
    }

    /* A stored hash carries what decides its value: the algorithm, the library version,
     * and whether JPEG went through libjpeg-turbo or stb_image. The settings are the
     * defaults here; a program that changes them stores them too. */
    printf("# algorithm=%s version=%d native-jpeg=%d\n", ph_algorithm_name(PH_ALGO_PHASH),
           ph_version_number(), ph_can_use_jpeg());
    for (size_t i = 0; i < n; i++) {
        char hex[17];
        if (items[i].status != PH_SUCCESS) {
            fprintf(stderr, "%s: %s\n", items[i].path, ph_get_error_string(items[i].status));
        } else if (ph_hash_to_hex(items[i].hashes[0], hex, sizeof(hex)) == PH_SUCCESS) {
            printf("%s  %s\n", hex, items[i].path);
        }
    }

    /* Every pair compared once: n * (n - 1) / 2 Hamming distances, a linear scan per file. */
    size_t found = 0, room = 16;
    pair_t *pairs = malloc(room * sizeof(*pairs));
    for (size_t i = 0; i < n && pairs; i++) {
        for (size_t j = i + 1; j < n; j++) {
            if (items[i].status != PH_SUCCESS || items[j].status != PH_SUCCESS) {
                continue;
            }
            int d = ph_hamming_distance(items[i].hashes[0], items[j].hashes[0]);
            if (d > max_bits) {
                continue;
            }
            if (found == room) {
                pair_t *more = realloc(pairs, 2 * room * sizeof(*pairs));
                if (!more) {
                    free(pairs);
                    pairs = NULL;
                    break;
                }
                pairs = more;
                room *= 2;
            }
            pairs[found++] = (pair_t){d, i, j};
        }
    }
    if (!pairs) {
        free(items);
        return 1;
    }
    qsort(pairs, found, sizeof(*pairs), by_distance);
    for (size_t k = 0; k < found; k++) {
        printf("%2d bits  %s  %s\n", pairs[k].distance, items[pairs[k].a].path,
               items[pairs[k].b].path);
    }
    printf("pairs within %ld bits: %zu\n", max_bits, found);

    free(pairs);
    free(items);
    return 0;
}

Comparing every pair of a collection of \(n\) hashes is \(n(n-1)/2\) comparisons, and looking up one new image is \(n\): a linear scan, the simplest search there is, and exact, since it compares with every stored hash.

What a linear scan costs¶

Each row is one comparison function, timed over a large stored set of its algorithm's digests, with one query compared with every one of them.¹ The two columns on the right multiply the time of one comparison out, for a query against a million stored hashes and for every pair of a hundred thousand.

The minimum of the runs, on the machine named at the foot of the page:

• A 64-bit hash is the cheapest to search. For ph_hamming_distance(), the call costs more than the comparison it makes: the same XOR and bit count written into the search loop, where the compiler can inline them, runs several times faster over the same array. A program that scans a large collection keeps its 64-bit hashes in a uint64_t array and counts the bits in its own loop.
• A digest costs more than its bytes. ph_hamming_distance_digest() and the other digest functions check both digests' kind and size before comparing, and every ph_digest_t in an array takes its 136 bytes. An 8-byte digest is several times slower to scan than the same hash as a uint64_t, and BMH's and mHash's digests, at 32 and 72 bytes, cost about the same as each other.
• Radial and ColorHash are the exceptions. Radial's peak correlation is a correlation at every cyclic shift of one digest against the other, a product for every pair of their coefficients, 1600 of them for 40 bytes, where a Hamming distance is a few word operations (Radial: peak correlation); ColorHash's intersection passes over its 108 bins twice and divides. Either comparison costs as much as a few dozen Hamming distances, and a search by it takes as many times longer.
• Every pair grows with the square of the collection. A query against a million stored hashes is a million comparisons; every pair of a hundred thousand is five billion, and the last column is longer than the one before it by the same factor of five thousand. Past that, a program either narrows the candidates before comparing them (an index, below) or compares only each new image with the collection, not the collection with itself.

/* site_stages scan: what one comparison costs in a search that compares a query with every
 * stored hash, a linear scan. Each case compares one query with a stored set of one
 * algorithm's shape: SCAN_HASHES 64-bit hashes, or SCAN_DIGESTS digests in an array of
 * ph_digest_t. The values are pseudo-random from a fixed seed, as these comparisons cost
 * the same whatever the values (Radial's apart, which refuses a digest without variance
 * before comparing, and random bytes always have some). Each case counts the stored values
 * within a fixed distance and adds it to scan_sink, so that the compiler cannot leave a
 * comparison out. */
enum {
    SCAN_HASHES = 1 << 20,
    SCAN_DIGESTS = 1 << 16,
    SCAN_WITHIN_BITS = 10,
};

static volatile size_t scan_sink;

typedef struct {
    uint64_t query;
    const uint64_t *hashes;
    ph_digest_t query_digest;
    const ph_digest_t *digests;
} scan_t;

/* The library's comparison of two 64-bit hashes, called once per stored hash. */
static size_t scan_hamming(const scan_t *s) {
    size_t within = 0;
    for (size_t i = 0; i < SCAN_HASHES; i++) {
        within += ph_hamming_distance(s->query, s->hashes[i]) <= SCAN_WITHIN_BITS;
    }
    return within;
}

#if defined(__GNUC__) || defined(__clang__)
/* The same distance computed in the caller's loop, where the compiler can inline the
 * popcount and vectorize the loop: what a search loop of its own costs. */
static size_t scan_popcount(const scan_t *s) {
    size_t within = 0;
    for (size_t i = 0; i < SCAN_HASHES; i++) {
        within += __builtin_popcountll(s->query ^ s->hashes[i]) <= SCAN_WITHIN_BITS;
    }
    return within;
}
#endif

/* The function a digest's kind calls for, once per stored digest. */
static size_t scan_digests(const scan_t *s) {
    size_t within = 0;
    const ph_digest_t *q = &s->query_digest;
    for (size_t i = 0; i < SCAN_DIGESTS; i++) {
        const ph_digest_t *d = &s->digests[i];
        double v = 0.0;
        switch (q->kind) {
            case PH_DIGEST_KIND_BITS:
                within += ph_hamming_distance_digest(q, d) <= SCAN_WITHIN_BITS;
                break;
            case PH_DIGEST_KIND_COEFFICIENTS:
                within += ph_radial_similarity(q, d, &v) == PH_SUCCESS && v >= 0.9;
                break;
            case PH_DIGEST_KIND_HISTOGRAM:
                within += ph_histogram_intersection(q, d, &v) == PH_SUCCESS && v >= 0.9;
                break;
            default:
                within += ph_l2_distance(q, d) <= 10.0;
                break;
        }
    }
    return within;
}

/* The cases: a name, the scan, how many values it compares, and for a digest the
 * algorithm whose shape (ph_digest_info()) the stored digests take. */
static const struct {
    const char *name;
    size_t (*scan)(const scan_t *s);
    size_t count;
    ph_algorithm_t algo;
} scan_cases[] = {
    {"hamming", scan_hamming, SCAN_HASHES, PH_ALGO_PHASH},
#if defined(__GNUC__) || defined(__clang__)
    {"popcount", scan_popcount, SCAN_HASHES, PH_ALGO_PHASH},
#endif
    {"phash", scan_digests, SCAN_DIGESTS, PH_ALGO_PHASH},
    {"bmh", scan_digests, SCAN_DIGESTS, PH_ALGO_BMH},
    {"mhash", scan_digests, SCAN_DIGESTS, PH_ALGO_MHASH},
    {"radial", scan_digests, SCAN_DIGESTS, PH_ALGO_RADIAL},
    {"color_hash", scan_digests, SCAN_DIGESTS, PH_ALGO_COLOR_HASH},
    {"color_moments", scan_digests, SCAN_DIGESTS, PH_ALGO_COLOR_MOMENTS},
};

/* xorshift64: a fixed sequence of pseudo-random values. */
static uint64_t scan_random(void) {
    static uint64_t x = 0x9E3779B97F4A7C15u;
    x ^= x << 13;
    x ^= x >> 7;
    x ^= x << 17;
    return x;
}

/* A digest of the shape `algo` gives with the default settings, its bytes random. */
static int random_digest(ph_algorithm_t algo, ph_digest_t *d) {
    size_t size;
    ph_digest_kind_t kind;
    if (ph_digest_info(NULL, algo, &size, &kind) != PH_SUCCESS) {
        return 1;
    }
    memset(d, 0, sizeof(*d));
    d->size = (uint8_t)size;
    d->kind = (uint8_t)kind;
    for (size_t i = 0; i < size; i++) {
        d->data[i] = (uint8_t)scan_random();
    }
    return 0;
}

/* site_stages scan: every case of scan_cases[], timed as the cases of `time` are; the
 * minimum and the median of its runs, in nanoseconds per comparison, as one JSON object. */
int mode_scan(int argc, char **argv) {
    (void)argc;
    (void)argv;
    uint64_t *hashes = malloc(SCAN_HASHES * sizeof(*hashes));
    ph_digest_t *digests = malloc(SCAN_DIGESTS * sizeof(*digests));
    if (!hashes || !digests) {
        free(hashes);
        free(digests);
        return fail("out of memory", NULL);
    }
    scan_t s = {scan_random(), hashes, {{0}, 0, 0, {0}}, digests};
    for (size_t i = 0; i < SCAN_HASHES; i++) {
        hashes[i] = scan_random();
    }
    json_t j = json_begin(stdout);
    json_int(&j, "hashes", SCAN_HASHES);
    json_int(&j, "digests", SCAN_DIGESTS);
    json_t cases = json_object(&j, "cases");
    static double ns[TIME_MAX_RUNS];
    int status = 0;
    for (size_t c = 0; c < COUNT(scan_cases) && !status; c++) {
        for (size_t i = 0; i < SCAN_DIGESTS && !status; i++) {
            status = random_digest(scan_cases[c].algo, &digests[i]);
        }
        if (status || random_digest(scan_cases[c].algo, &s.query_digest)) {
            status = fail("no digest shape for", scan_cases[c].name);
            break;
        }
        int runs = 0;
        double total = 0.0;
        for (int k = -1; k < TIME_MAX_RUNS; k++) { /* k = -1: the warm-up */
            double t0 = seconds();
            scan_sink += scan_cases[c].scan(&s);
            double t = seconds() - t0;
            if (k >= 0) {
                ns[runs++] = t * 1e9 / (double)scan_cases[c].count;
                total += t;
                if (runs >= TIME_MIN_RUNS && total >= TIME_MIN_SECONDS) {
                    break;
                }
            }
        }
        qsort(ns, (size_t)runs, sizeof(ns[0]), cmp_double);
        json_t o = json_object(&cases, scan_cases[c].name);
        json_double(&o, "min_ns", ns[0]);
        json_double(&o, "median_ns", ns[runs / 2]);
        json_int(&o, "runs", runs);
        json_int(&o, "size", scan_cases[c].scan == scan_digests ? s.query_digest.size : 8);
        json_close_object(&o);
    }
    json_close_object(&cases);
    json_end(&j);
    free(hashes);
    free(digests);
    return status;
}

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")

cmake --preset release && cmake --build --preset release --target site_stages
build/release/site_stages scan

An index¶

An index finds the stored hashes within a threshold of a query without comparing it with all of them. The library has none, and the two ideas below are for a program to build over its own storage. Both are for the bit hashes, whose Hamming distance is a true distance; they find exactly the hashes a linear scan would.

Multi-index hashing cuts every 64-bit hash into blocks, for example four of 16 bits, and keeps one table per block that maps the block's value to the hashes that have it. If two hashes differ in at most \(t\) bits, the differences fall into four blocks, so at least one block differs in at most \(\lfloor t/4 \rfloor\) bits. A query therefore looks up, in each of the four tables, every block value within \(\lfloor t/4 \rfloor\) bits of its own; whatever it finds is a candidate, and a full Hamming distance decides.

flowchart LR
    Q[Query hash] --> B1[Bits 63–48]
    Q --> B2[Bits 47–32]
    Q --> B3[Bits 31–16]
    Q --> B4[Bits 15–0]
    B1 & B2 & B3 & B4 --> L["Look up values within ⌊t/4⌋ bits<br>in that block's table"]
    L --> C[Candidates]
    C --> V["Hamming distance ≤ t?"]

How far a lookup reaches depends on the threshold. Here are the thresholds that accept 95 % of the copies of each corpus (Choosing a threshold), and what each lookup has to read:

At the strict thresholds of the photographs, a lookup reads a few values in each of the four tables. At the looser thresholds of the synthetic images, it reads thousands in each, and an index of four blocks saves less. More, shorter blocks lower the radius, down to exact lookups once there are more blocks than the threshold has bits; but a shorter block is shared by more stored hashes, and each of them is a candidate to verify.

A BK-tree uses the triangle inequality of the Hamming distance. Every node holds a hash, and its children are filed under their distance from it. A query at distance \(d\) from a node can only have answers among the children filed under \(d - t\) to \(d + t\), so the others, with everything under them, are skipped. How much is skipped depends on the threshold and on how the stored hashes spread: a threshold near half the bits skips almost nothing. It is measured on the collection, against the linear scan it replaces.

flowchart TD
    R["Root hash"] -->|"distance 3"| A[Hash]
    R -->|"distance 9"| B[Hash]
    R -->|"distance 21"| C[Hash]
    A -->|"distance 2"| D[Hash]
    B -->|"distance 5"| E[Hash]

Choosing the threshold¶

The threshold belongs to the collection, and a threshold measured on one set of images is a guess on another (Choosing a threshold). The example helps measure one: run it with a loose threshold over images whose duplicates are known, and read the pairs from the closest out. The distance at which pairs stop being the same picture is where the threshold goes, and Measuring a threshold on your own collection turns that into a share of copies accepted and a share of different pairs let through.

When the library or the settings change¶

A stored hash stays comparable as long as the four things it does not record stay the same. When one of them changes:

• A setting or the algorithm: every stored hash of that algorithm is recomputed with the new configuration. Hashes from two configurations are never compared with each other; the number such a comparison returns measures nothing.
• The library version: only the algorithms the changelog names change their values, and only their hashes are recomputed. Migrating from 1.x describes three ways to do it, from recomputing everything at once to recomputing each hash the next time it is read.
• The JPEG decoder: the hashes of some JPEG files move by small amounts (Same hash on every machine), which a search by threshold tolerates and a search by equality does not.