Radial — Radial Variance Hash¶

/* 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;
}

/* 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;
    }
}

/* 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;
}

/* site_stages radial-profiles <reference> <image>...: for every image, the default
 * Radial's variance profile and digest, recomputed from its stages and checked against
 * ph_compute_radial_hash(), and the digest compared with the reference's by
 * ph_radial_similarity() (null where the library refuses, an image with no angular
 * structure): one JSON line per image, the reference first. */
int mode_radial_profiles(int argc, char **argv) {
    ph_digest_t ref = {0};
    for (int i = 0; i < argc; i++) {
        ph_context_t *ctx = NULL;
        if (load_image(&ctx, argv[i])) {
            return 1;
        }
        radial_stages_t s;
        ph_digest_t d;
        int bad = radial_run(ctx, &radial_variants[RADIAL_DEFAULT], &s) ||
                  ph_compute_radial_hash(ctx, &d) != PH_SUCCESS;
        ph_free(ctx);
        if (bad) {
            radial_free(&s);
            return 1;
        }
        if (i == 0) {
            ref = d;
        }
        json_t j = json_begin(stdout);
        json_string(&j, "file", argv[i]);
        radial_json(&j, "radial", &s, 1);
        double pcc = 0.0;
        if (ph_radial_similarity(&ref, &d, &pcc) == PH_SUCCESS) {
            json_double(&j, "similarity", pcc);
        } else {
            json_null(&j, "similarity");
        }
        json_end(&j);
        radial_free(&s);
    }
    return 0;
}

/* One Radial computation with the given settings, recomputed from its stages and checked
 * against ph_compute_radial_hash() with the same settings: the grayscale image blurred at
 * `sigma` and put through `gamma`; the variance along each of `projections` lines through
 * the center, `samples` points on each; the profile's mean and spread; the profile
 * standardized; its first PH_RADIAL_COEFFS DCT coefficients; and the digest, the
 * coefficients mapped onto 0..255 by their own minimum and maximum (all zero when the
 * profile has no spread worth the name). The context is left at the default settings. */
typedef struct {
    radial_variant_t settings;
    double center_x, center_y, radius;
    double *variance, *standardized; /* projections each */
    double mean, spread_sq;
    int structure;
    double coefficients[PH_RADIAL_COEFFS];
    uint8_t digest[PH_RADIAL_COEFFS];
    uint8_t *blurred; /* width x height */
} radial_stages_t;

static void radial_free(radial_stages_t *s) {
    free(s->variance);
    free(s->standardized);
    free(s->blurred);
    s->variance = s->standardized = NULL;
    s->blurred = NULL;
}

static int radial_run(ph_context_t *ctx, const radial_variant_t *v, radial_stages_t *s) {
    const int w = ctx->image.width, h = ctx->image.height, projections = v->projections;
    const size_t npix = (size_t)w * (size_t)h;
    memset(s, 0, sizeof(*s));
    s->settings = *v;
    s->variance = malloc((size_t)projections * sizeof(double));
    s->standardized = calloc((size_t)projections, sizeof(double));
    s->blurred = malloc(npix);
    float *scratch = malloc(npix * sizeof(float));
    const uint8_t *gray = ph_get_gray(ctx);
    if (!s->variance || !s->standardized || !s->blurred || !scratch || !gray ||
        !radial_settings(ctx, v)) {
        free(scratch);
        return fail("cannot set up Radial", NULL);
    }
    ph_gaussian_blur_sigma(gray, w, h, v->sigma, scratch, s->blurred);
    free(scratch);
    ph_apply_gamma(ctx, s->blurred, w, h);

    s->center_x = w / 2.0;
    s->center_y = h / 2.0;
    s->radius = (w < h ? w : h) / 2.0;
    double sum = 0.0, sum_sq = 0.0;
    for (int i = 0; i < projections; i++) {
        double theta = i * M_PI / projections;
        s->variance[i] =
            ph_projection_variance(s->blurred, w, h, s->center_x, s->center_y, s->radius,
                                   (float)cos(theta), (float)sin(theta), v->samples);
        sum += s->variance[i];
        sum_sq += s->variance[i] * s->variance[i];
    }
    s->mean = sum / projections;
    s->spread_sq = sum_sq / projections - s->mean * s->mean;
    s->structure = s->mean > PH_RADIAL_MIN_MEAN_VARIANCE &&
                   s->spread_sq > PH_RADIAL_MIN_RELATIVE_SPREAD * s->mean * s->mean;
    if (s->structure) {
        for (int i = 0; i < projections; i++) {
            s->standardized[i] = (s->variance[i] - s->mean) / sqrt(s->spread_sq);
        }
        if (ph_dct1d_partial(s->standardized, projections, PH_RADIAL_COEFFS, s->coefficients) !=
            PH_SUCCESS) {
            return fail("ph_dct1d_partial failed", NULL);
        }
        double lo = s->coefficients[0], hi = s->coefficients[0];
        for (int k = 1; k < PH_RADIAL_COEFFS; k++) {
            lo = s->coefficients[k] < lo ? s->coefficients[k] : lo;
            hi = s->coefficients[k] > hi ? s->coefficients[k] : hi;
        }
        for (int k = 0; k < PH_RADIAL_COEFFS; k++) {
            s->digest[k] = (uint8_t)(255.0 * (s->coefficients[k] - lo) / (hi - lo));
        }
    }

    ph_digest_t lib;
    int bad = ph_compute_radial_hash(ctx, &lib) != PH_SUCCESS;
    if (!radial_settings(ctx, &radial_variants[RADIAL_DEFAULT]) || bad) {
        return fail("ph_compute_radial_hash failed", NULL);
    }
    char what[64];
    snprintf(what, sizeof(what), "Radial (%d x %d, sigma %g, gamma %g)", projections, v->samples,
             (double)v->sigma, (double)v->gamma);
    return check_digest(what, s->digest, PH_RADIAL_COEFFS, &lib);
}

/* 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;
}