update web images.

.gitignore

@@ -15,3 +15,17 @@ thumbcache/
 legacy/
 *.pdf
 run-*.sh
+
+# demo page source assets: re-include web/img, which the broad img/ and
+# *.png patterns above would otherwise swallow
+!web/img/
+!web/img/*
+
+# demo page build artifacts, regenerated by `cargo run --release --bin
+# demoassets`, `cargo run --release --bin storekeys` and
+# `python3 pr_to_json.py` (needed when deploying the page)
+web/flickr/
+web/pr-data.json
+
+# playwright-core for verify-retrieval.mjs (npm install --no-save playwright-core)
+node_modules/

Cargo.toml

@@ -52,3 +52,18 @@ required-features = ["cli"]
 name = "bench"
 path = "src/bin/bench.rs"
 required-features = ["cli"]
+
+[[bin]]
+name = "demoassets"
+path = "src/bin/demoassets.rs"
+required-features = ["cli"]
+
+[[bin]]
+name = "storekeys"
+path = "src/bin/storekeys.rs"
+required-features = ["cli"]
+
+[[bin]]
+name = "salientlab"
+path = "src/bin/salientlab.rs"
+required-features = ["cli"]

src/bin/demoassets.rs

@@ -0,0 +1,69 @@
+//! Generates the static assets for the web demo ranking pool:
+//! a thumbnail per image plus hashes.json with the three hashes,
+//! computed from the full-size images.
+use clap::Parser;
+use image_similarity::descriptors::{Descriptor, Median, DCT, PHash};
+use image_similarity::open_image;
+use rayon::prelude::*;
+use std::fs;
+use std::io::Write;
+use std::path::PathBuf;
+
+#[derive(Parser)]
+struct Args {
+    /// Directory with source images
+    #[arg(default_value = "smolflickr")]
+    source: PathBuf,
+    /// Output directory, thumbnails go in <out>/thumbs
+    #[arg(default_value = "web/flickr")]
+    out: PathBuf,
+    /// Thumbnail bounding box in pixels
+    #[arg(long, default_value_t = 140)]
+    thumb: u32,
+}
+
+fn main() {
+    let args = Args::parse();
+    let thumbs = args.out.join("thumbs");
+    fs::create_dir_all(&thumbs).expect("creating output dir");
+
+    let mut names: Vec<String> = fs::read_dir(&args.source)
+        .expect("reading source dir")
+        .filter_map(|entry| {
+            let entry = entry.ok()?;
+            entry.path().is_file().then(|| entry.file_name().to_string_lossy().into_owned())
+        })
+        .collect();
+    names.sort();
+
+    let entries: Vec<String> = names
+        .par_iter()
+        .filter_map(|name| {
+            let img = match open_image(args.source.join(name)) {
+                Ok(img) => img,
+                Err(e) => {
+                    eprintln!("skipping {name}: {e}");
+                    return None;
+                }
+            };
+            let dct = DCT::new().describe(&img);
+            let median = Median.describe(&img);
+            let phash = PHash.describe(&img);
+
+            let thumb = img.thumbnail(args.thumb, args.thumb);
+            let mut buffer = std::io::Cursor::new(Vec::new());
+            thumb
+                .to_rgb8()
+                .write_with_encoder(image::codecs::jpeg::JpegEncoder::new_with_quality(&mut buffer, 80))
+                .expect("encoding thumbnail");
+            fs::write(thumbs.join(name), buffer.get_ref()).expect("writing thumbnail");
+
+            Some(format!("[\"{name}\",\"{dct:016x}\",\"{median:016x}\",\"{phash:016x}\"]"))
+        })
+        .collect();
+
+    let json = format!("[\n{}\n]\n", entries.join(",\n"));
+    let mut file = fs::File::create(args.out.join("hashes.json")).expect("creating hashes.json");
+    file.write_all(json.as_bytes()).expect("writing hashes.json");
+    println!("{} images -> {}", entries.len(), args.out.display());
+}

src/bin/salientlab.rs

@@ -0,0 +1,853 @@
+//! Throwaway lab for the salient-point hash variants.
+//!
+//! Two modes:
+//! - default: small-sample table of per-mutation Hamming distances plus
+//!   detector repeatability, quick iteration on a few hundred images
+//! - --pr: the full run. Calibrates on the first --calib images, hashes
+//!   every image and mutant under all variants (checkpointed to
+//!   --hashes with resume), then computes PR curves with the same
+//!   tally semantics as DescriptorStore::get_stats and plots one PNG
+//!   per variant. The dct hash rides along as a reference method.
+//!
+//! Not part of the real benchmark pipeline.
+
+use clap::Parser;
+use image_similarity::descriptors::{Descriptor, DCT};
+use image_similarity::mutators::{get_all_mutators, CropCenter, Mutator, Rotate};
+use image_similarity::open_image;
+use image_similarity::salient::{
+    binarize64, hash_strongest, pooled_vector, Feature, Salient,
+};
+use plotters::prelude::*;
+use rayon::prelude::*;
+use serde::{Deserialize, Serialize};
+use std::collections::HashSet;
+use std::path::{Path, PathBuf};
+use std::sync::atomic::{AtomicUsize, Ordering};
+use std::time::Instant;
+
+#[derive(Parser)]
+#[command(about = "Salient-point hash experiments")]
+struct Args {
+    /// Directory with base images
+    #[arg(long, default_value = "smolflickr")]
+    dir: PathBuf,
+    /// Number of images to sample (default: 200 for the table, all for --pr)
+    #[arg(long)]
+    count: Option<usize>,
+    /// Keypoints kept per image
+    #[arg(long, default_value_t = 100)]
+    features: usize,
+    /// Keypoints considered for the repeatability measurement
+    #[arg(long, default_value_t = 50)]
+    repeat_n: usize,
+    /// Full PR-curve run over the whole directory
+    #[arg(long, default_value_t = false)]
+    pr: bool,
+    /// Images used to calibrate codebook and thresholds in --pr mode
+    #[arg(long, default_value_t = 2000)]
+    calib: usize,
+    /// Checkpoint file for the hashes in --pr mode
+    #[arg(long, default_value = "salient-hashes.mp")]
+    hashes: PathBuf,
+    /// Skip hashing, compute PR curves from an existing checkpoint file
+    #[arg(long, default_value_t = false)]
+    stats_only: bool,
+}
+
+const VARIANTS: [&str; 8] = [
+    "pooled", "strongest", "consensus", "pooled-cal", "pooled-sim", "cons-cal", "words-any",
+    "words-med",
+];
+/// All hashed methods: the salient variants plus the dct reference
+const METHODS: usize = VARIANTS.len() + 1;
+
+fn method_name(v: usize) -> &'static str {
+    if v < VARIANTS.len() {
+        VARIANTS[v]
+    } else {
+        "dct-ref"
+    }
+}
+
+/// Per image-version summary statistics that the hash variants read from
+struct Stats {
+    /// Response-weighted mean descriptor, unit normalized
+    pooled: [f32; 64],
+    /// Fraction of keypoints voting 1 on each binarize64 bit position
+    votes: [f32; 64],
+    /// Keypoints assigned to each of the 64 visual words
+    counts: [u32; 64],
+}
+
+struct Calibration {
+    /// The visual-word codebook
+    centroids: Vec<[f32; 64]>,
+    /// Per-dimension corpus median of the pooled vector
+    pooled_med: [f32; 64],
+    /// Corpus median vote fraction per bit
+    vote_med: [f32; 64],
+    /// Corpus median count per visual word
+    word_med: [u32; 64],
+    /// Random +-1 hyperplanes for the simhash variant
+    planes: Vec<[f32; 64]>,
+}
+
+/// Checkpoint file contents for --pr mode. Per image: one hash per method
+/// for the base image followed by every mutant, in mutator order.
+#[derive(Serialize, Deserialize)]
+struct HashFile {
+    tags: Vec<String>,
+    results: Vec<(String, Vec<[u64; METHODS]>)>,
+}
+
+fn main() {
+    let args = Args::parse();
+    if args.pr {
+        run_pr(&args);
+    } else {
+        run_table(&args);
+    }
+}
+
+/// The standard suite (which now carries a hard rotate45) plus a mid
+/// rotation and a heavy crop, to trace the geometric falloff more finely
+fn lab_mutators() -> Vec<Box<dyn Mutator>> {
+    let mut mutators = get_all_mutators();
+    mutators.push(Box::new(Rotate { degrees: 15.0 }));
+    mutators.push(Box::new(CropCenter { keep_percent: 50 }));
+    mutators
+}
+
+fn image_paths(dir: &Path, count: usize) -> Vec<(String, PathBuf)> {
+    let mut paths: Vec<(String, PathBuf)> = std::fs::read_dir(dir)
+        .expect("cannot read image dir")
+        .filter_map(|e| e.ok())
+        .filter(|e| {
+            matches!(
+                e.path().extension().and_then(|x| x.to_str()),
+                Some("jpg") | Some("jpeg") | Some("png")
+            )
+        })
+        .filter_map(|e| e.file_name().into_string().ok().map(|n| (n, e.path())))
+        .collect();
+    paths.sort();
+    paths.truncate(count);
+    assert!(!paths.is_empty(), "no images found in {:?}", dir);
+    paths
+}
+
+// ---------------------------------------------------------------------------
+// full PR run
+
+fn run_pr(args: &Args) {
+    let paths = image_paths(&args.dir, args.count.unwrap_or(usize::MAX));
+    let salient = Salient {
+        params: image_similarity::salient::Params {
+            max_features: args.features,
+            ..Default::default()
+        },
+    };
+    let mutators = lab_mutators();
+    let tags: Vec<String> = mutators.iter().map(|m| m.tag()).collect();
+
+    let mut file = load_hashes(&args.hashes).unwrap_or(HashFile {
+        tags: tags.clone(),
+        results: Vec::new(),
+    });
+    assert_eq!(
+        file.tags, tags,
+        "{} was made with a different mutator suite, move it away first",
+        args.hashes.display()
+    );
+
+    if !args.stats_only {
+        // calibration is deterministic, so recomputing it on resume is safe
+        let calib = calibrate(&paths[..args.calib.min(paths.len())], &salient);
+
+        let done: HashSet<&String> = file.results.iter().map(|(n, _)| n).collect();
+        let todo: Vec<&(String, PathBuf)> =
+            paths.iter().filter(|(n, _)| !done.contains(n)).collect();
+        drop(done);
+        println!(
+            "{} images total, {} still to hash, {} mutations, {} methods",
+            paths.len(),
+            todo.len(),
+            mutators.len(),
+            METHODS
+        );
+
+        let t0 = Instant::now();
+        let mut last_save = Instant::now();
+        let mut hashed = 0usize;
+        for batch in todo.chunks(512) {
+            let out: Vec<(String, Vec<[u64; METHODS]>)> = batch
+                .par_iter()
+                .filter_map(|(name, path)| {
+                    let img = open_image(path).ok()?;
+                    let mut versions = Vec::with_capacity(mutators.len() + 1);
+                    versions.push(hash_methods(&img, &salient, &calib));
+                    for m in &mutators {
+                        versions.push(hash_methods(&m.mutate(&img), &salient, &calib));
+                    }
+                    Some((name.clone(), versions))
+                })
+                .collect();
+            file.results.extend(out);
+            hashed += batch.len();
+            let rate = hashed as f32 / t0.elapsed().as_secs_f32();
+            eprintln!(
+                "  {} / {} images, {:.1} img/s, ~{:.0} min left",
+                hashed,
+                todo.len(),
+                rate,
+                (todo.len() - hashed) as f32 / rate / 60.0
+            );
+            if last_save.elapsed().as_secs() >= 120 {
+                save_hashes(&args.hashes, &file);
+                last_save = Instant::now();
+            }
+        }
+        save_hashes(&args.hashes, &file);
+        println!(
+            "hashing done in {:.1} min, {} images in {}",
+            t0.elapsed().as_secs_f32() / 60.0,
+            file.results.len(),
+            args.hashes.display()
+        );
+    }
+
+    pr_stats(&file, &mutators);
+}
+
+fn calibrate(paths: &[(String, PathBuf)], salient: &Salient) -> Calibration {
+    println!("calibrating on {} images...", paths.len());
+    let t0 = Instant::now();
+    let feats: Vec<Vec<Feature>> = paths
+        .par_iter()
+        .filter_map(|(_, path)| Some(salient.features(&open_image(path).ok()?)))
+        .collect();
+    let samples: Vec<[f32; 64]> = feats
+        .iter()
+        .flat_map(|f| f.iter().take(50).map(|f| f.desc))
+        .collect();
+    let centroids = kmeans(&samples, 64, 25);
+    let stats: Vec<Stats> = feats
+        .par_iter()
+        .map(|f| stats_for(f, &centroids))
+        .collect();
+    let calib = Calibration {
+        pooled_med: column_median(&stats, |s| s.pooled),
+        vote_med: column_median(&stats, |s| s.votes),
+        word_med: column_median_u32(&stats),
+        planes: random_planes(64),
+        centroids,
+    };
+    println!(
+        "calibration took {:.1}s ({} descriptors)",
+        t0.elapsed().as_secs_f32(),
+        samples.len()
+    );
+    calib
+}
+
+fn hash_methods(img: &image::DynamicImage, salient: &Salient, calib: &Calibration) -> [u64; METHODS] {
+    let feats = salient.features(img);
+    let stats = stats_for(&feats, &calib.centroids);
+    let v = variant_hashes(&feats, &stats, calib);
+    let mut out = [0u64; METHODS];
+    out[..VARIANTS.len()].copy_from_slice(&v);
+    out[VARIANTS.len()] = DCT::new().describe(img);
+    out
+}
+
+fn load_hashes(path: &Path) -> Option<HashFile> {
+    let bytes = std::fs::read(path).ok()?;
+    match rmp_serde::from_slice(&bytes) {
+        Ok(file) => Some(file),
+        Err(e) => panic!("cannot parse {}: {e}", path.display()),
+    }
+}
+
+fn save_hashes(path: &Path, file: &HashFile) {
+    let bytes = rmp_serde::to_vec(file).expect("serialize hashes");
+    let tmp = path.with_extension("tmp");
+    std::fs::write(&tmp, bytes).expect("write hashes");
+    std::fs::rename(&tmp, path).expect("rename hashes");
+}
+
+/// Same tally as DescriptorStore::get_stats, done by brute force instead of
+/// the BK-tree: every base compares against every other hash. A hit is a
+/// true positive when it is the base's own mutant, a false positive for
+/// mutator m when it is an m-mutant of another image, and a false positive
+/// for every mutator when it is a different base image.
+fn pr_stats(file: &HashFile, mutators: &[Box<dyn Mutator>]) {
+    let max_threshold = 24usize;
+    let n = file.results.len();
+    let nm = mutators.len();
+    println!("\ncomputing PR curves over {n} images...");
+
+    for v in 0..METHODS {
+        let t0 = Instant::now();
+        let bases: Vec<u64> = file.results.iter().map(|(_, h)| h[0][v]).collect();
+        let muts: Vec<Vec<u64>> = (0..nm)
+            .map(|m| file.results.iter().map(|(_, h)| h[m + 1][v]).collect())
+            .collect();
+
+        #[derive(Clone)]
+        struct Hists {
+            // 65 slots: two hashes can differ in all 64 bits
+            base_fp: [u64; 65],
+            tp: Vec<[u64; 65]>,
+            fp: Vec<[u64; 65]>,
+        }
+        let empty = || Hists {
+            base_fp: [0; 65],
+            tp: vec![[0; 65]; nm],
+            fp: vec![[0; 65]; nm],
+        };
+        let hists = (0..n)
+            .into_par_iter()
+            .fold(empty, |mut h, i| {
+                let hb = bases[i];
+                for (j, other) in bases.iter().enumerate() {
+                    if j != i {
+                        h.base_fp[(hb ^ other).count_ones() as usize] += 1;
+                    }
+                }
+                for m in 0..nm {
+                    let mm = &muts[m];
+                    h.tp[m][(hb ^ mm[i]).count_ones() as usize] += 1;
+                    for (j, other) in mm.iter().enumerate() {
+                        if j != i {
+                            h.fp[m][(hb ^ other).count_ones() as usize] += 1;
+                        }
+                    }
+                }
+                h
+            })
+            .reduce(empty, |mut a, b| {
+                for (x, y) in a.base_fp.iter_mut().zip(&b.base_fp) {
+                    *x += y;
+                }
+                for m in 0..nm {
+                    for (x, y) in a.tp[m].iter_mut().zip(&b.tp[m]) {
+                        *x += y;
+                    }
+                    for (x, y) in a.fp[m].iter_mut().zip(&b.fp[m]) {
+                        *x += y;
+                    }
+                }
+                a
+            });
+
+        let cumulative = |at: &[u64; 65]| {
+            let mut cum = [0u64; 65];
+            let mut sum = 0;
+            for (t, count) in at.iter().enumerate() {
+                sum += count;
+                cum[t] = sum;
+            }
+            cum
+        };
+        let fp_base = cumulative(&hists.base_fp);
+
+        println!(
+            "\nStats for salient {} ({:.1}s)",
+            method_name(v),
+            t0.elapsed().as_secs_f32()
+        );
+        let curves: Vec<(String, Vec<(f64, f64)>)> = (0..nm)
+            .map(|m| {
+                let tp = cumulative(&hists.tp[m]);
+                let fp_mut = cumulative(&hists.fp[m]);
+                let curve: Vec<(f64, f64)> = (0..max_threshold)
+                    .map(|t| {
+                        let tpos = tp[t] as f64;
+                        let fpos = (fp_mut[t] + fp_base[t]) as f64;
+                        let recall = tpos / n as f64;
+                        let precision = if tpos + fpos == 0.0 {
+                            1.0
+                        } else {
+                            tpos / (tpos + fpos)
+                        };
+                        (recall, precision)
+                    })
+                    .collect();
+                (mutators[m].info(), curve)
+            })
+            .collect();
+        for (name, curve) in &curves {
+            println!("{name}:");
+            for (r, p) in curve {
+                println!("{r}, {p}");
+            }
+        }
+        plot_pr(&format!("sal-{}-pr.png", method_name(v)), method_name(v), &curves);
+    }
+    println!("\nwrote sal-<method>-pr.png for {METHODS} methods");
+}
+
+/// Same look as the PR plots from main.rs
+fn plot_pr(filename: &str, method: &str, curves: &[(String, Vec<(f64, f64)>)]) {
+    let root = BitMapBackend::new(filename, (1024, 768)).into_drawing_area();
+    root.fill(&WHITE).unwrap();
+    let mut chart = ChartBuilder::on(&root)
+        .caption(
+            format!("PR curves (salient {method})"),
+            ("sans-serif", 40).into_font(),
+        )
+        .margin(25)
+        .set_all_label_area_size(50)
+        .build_cartesian_2d(0f64..1f64, 0f64..1f64)
+        .unwrap();
+    chart
+        .configure_mesh()
+        .x_labels(10)
+        .y_labels(10)
+        .disable_mesh()
+        .x_desc("Recall")
+        .y_desc("Precision")
+        .x_label_formatter(&|x| format!("{x:.3}"))
+        .y_label_formatter(&|x| format!("{x:.3}"))
+        .draw()
+        .unwrap();
+    for (i, (name, curve)) in curves.iter().enumerate() {
+        let color = Palette99::pick(i).to_rgba();
+        chart
+            .draw_series(LineSeries::new(curve.clone(), color.filled()).point_size(2))
+            .unwrap()
+            .label(name)
+            .legend(move |(x, y)| PathElement::new(vec![(x, y), (x + 20, y)], color));
+    }
+    chart
+        .configure_series_labels()
+        .position(SeriesLabelPosition::MiddleLeft)
+        .border_style(BLACK)
+        .draw()
+        .unwrap();
+}
+
+// ---------------------------------------------------------------------------
+// small-sample distance table
+
+struct MutData {
+    tag: String,
+    w: f32,
+    h: f32,
+    feats: Vec<Feature>,
+    dct: u64,
+}
+
+struct ImageData {
+    w: f32,
+    h: f32,
+    base: Vec<Feature>,
+    dct: u64,
+    muts: Vec<MutData>,
+}
+
+fn run_table(args: &Args) {
+    let paths = image_paths(&args.dir, args.count.unwrap_or(200));
+    let salient = Salient {
+        params: image_similarity::salient::Params {
+            max_features: args.features,
+            ..Default::default()
+        },
+    };
+    let mutators = lab_mutators();
+    println!(
+        "{} images, {} mutations, {} keypoints per image",
+        paths.len(),
+        mutators.len(),
+        args.features
+    );
+
+    let t0 = Instant::now();
+    let done = AtomicUsize::new(0);
+    let data: Vec<ImageData> = paths
+        .par_iter()
+        .filter_map(|(_, path)| {
+            let img = open_image(path).ok()?;
+            let dct = DCT::new();
+            let base = salient.features(&img);
+            let muts = mutators
+                .iter()
+                .map(|m| {
+                    let mutated = m.mutate(&img);
+                    MutData {
+                        tag: m.tag(),
+                        w: mutated.width() as f32,
+                        h: mutated.height() as f32,
+                        feats: salient.features(&mutated),
+                        dct: dct.describe(&mutated),
+                    }
+                })
+                .collect();
+            let n = done.fetch_add(1, Ordering::Relaxed) + 1;
+            if n % 50 == 0 {
+                eprintln!("  {} / {} images", n, paths.len());
+            }
+            Some(ImageData {
+                w: img.width() as f32,
+                h: img.height() as f32,
+                base,
+                dct: dct.describe(&img),
+                muts,
+            })
+        })
+        .collect();
+    let avg_feats: f32 =
+        data.iter().map(|d| d.base.len() as f32).sum::<f32>() / data.len() as f32;
+    println!(
+        "extraction: {:.1}s, avg {:.0} keypoints per base image",
+        t0.elapsed().as_secs_f32(),
+        avg_feats
+    );
+
+    repeatability(&data, &mutators, args.repeat_n);
+
+    // codebook, then per-version stats, then corpus calibration
+    let samples: Vec<[f32; 64]> = data
+        .iter()
+        .flat_map(|d| d.base.iter().take(50).map(|f| f.desc))
+        .collect();
+    println!("\ntraining 64-word codebook on {} descriptors...", samples.len());
+    let centroids = kmeans(&samples, 64, 25);
+
+    let base_stats: Vec<Stats> = data
+        .par_iter()
+        .map(|d| stats_for(&d.base, &centroids))
+        .collect();
+    let mut_stats: Vec<Vec<Stats>> = data
+        .par_iter()
+        .map(|d| d.muts.iter().map(|m| stats_for(&m.feats, &centroids)).collect())
+        .collect();
+
+    let calib = Calibration {
+        pooled_med: column_median(&base_stats, |s| s.pooled),
+        vote_med: column_median(&base_stats, |s| s.votes),
+        word_med: column_median_u32(&base_stats),
+        planes: random_planes(64),
+        centroids,
+    };
+
+    let base_hashes: Vec<[u64; 8]> = data
+        .iter()
+        .zip(&base_stats)
+        .map(|(d, s)| variant_hashes(&d.base, s, &calib))
+        .collect();
+    let mut_hashes: Vec<Vec<[u64; 8]>> = data
+        .iter()
+        .zip(&mut_stats)
+        .map(|(d, stats)| {
+            d.muts
+                .iter()
+                .zip(stats)
+                .map(|(m, s)| variant_hashes(&m.feats, s, &calib))
+                .collect()
+        })
+        .collect();
+
+    println!("\nhamming distance base vs mutant, mean/p90 (dct hash as reference):");
+    print!("  {:12} {:>10}", "mutation", "dct");
+    for v in VARIANTS {
+        print!(" {:>10}", v);
+    }
+    println!();
+    for (mi, m) in mutators.iter().enumerate() {
+        print!("  {:12}", m.tag());
+        let mut dists: Vec<u32> = data
+            .iter()
+            .map(|d| (d.dct ^ d.muts[mi].dct).count_ones())
+            .collect();
+        print_cell(&mut dists, false);
+        for v in 0..VARIANTS.len() {
+            let mut dists: Vec<u32> = (0..data.len())
+                .map(|i| (base_hashes[i][v] ^ mut_hashes[i][mi][v]).count_ones())
+                .collect();
+            print_cell(&mut dists, false);
+        }
+        println!();
+    }
+
+    // impostor baseline: unrelated image pairs; p5 shows how close the
+    // nearest false positives come
+    print!("  {:12}", "impostor*");
+    let mut dists: Vec<u32> = (0..data.len())
+        .map(|i| (data[i].dct ^ data[(i + 1) % data.len()].dct).count_ones())
+        .collect();
+    print_cell(&mut dists, true);
+    for v in 0..VARIANTS.len() {
+        let mut dists: Vec<u32> = (0..data.len())
+            .map(|i| {
+                let j = (i + 1) % data.len();
+                (base_hashes[i][v] ^ base_hashes[j][v]).count_ones()
+            })
+            .collect();
+        print_cell(&mut dists, true);
+    }
+    println!("\n  *impostor row is mean/p5 over unrelated pairs: match rows want");
+    println!("   to sit far below it, especially below its p5");
+}
+
+/// One table cell: mean/p90 for match rows, mean/p5 for the impostor row
+fn print_cell(dists: &mut Vec<u32>, impostor: bool) {
+    dists.sort();
+    let mean = dists.iter().sum::<u32>() as f32 / dists.len() as f32;
+    let q = if impostor {
+        dists[dists.len() / 20]
+    } else {
+        dists[(dists.len() * 9 / 10).min(dists.len() - 1)]
+    };
+    print!(" {:>6.1}/{:<3}", mean, q);
+}
+
+fn repeatability(data: &[ImageData], mutators: &[Box<dyn Mutator>], top: usize) {
+    println!("\ndetector repeatability (top {} keypoints, 4px tolerance):", top);
+    for (mi, m) in mutators.iter().enumerate() {
+        let (mut matched, mut total) = (0usize, 0usize);
+        for img in data {
+            let md = &img.muts[mi];
+            let mapped: Vec<(f32, f32, f32)> = md
+                .feats
+                .iter()
+                .take(top)
+                .map(|f| map_to_base(&md.tag, img.w, img.h, md.w, md.h, f.kp.x, f.kp.y, f.kp.sigma))
+                .collect();
+            for f in img.base.iter().take(top) {
+                if !in_region(&md.tag, img.w, img.h, md.w, md.h, f.kp.x, f.kp.y) {
+                    continue;
+                }
+                total += 1;
+                let hit = mapped.iter().any(|(x, y, s)| {
+                    let (dx, dy) = (x - f.kp.x, y - f.kp.y);
+                    let ratio = (s / f.kp.sigma).max(f.kp.sigma / s);
+                    dx * dx + dy * dy <= 16.0 && ratio <= 1.6
+                });
+                if hit {
+                    matched += 1;
+                }
+            }
+        }
+        println!(
+            "  {:12} {:.2}",
+            m.tag(),
+            matched as f32 / total.max(1) as f32
+        );
+    }
+}
+
+// ---------------------------------------------------------------------------
+// shared variant machinery
+
+fn stats_for(feats: &[Feature], centroids: &[[f32; 64]]) -> Stats {
+    let pooled = pooled_vector(feats);
+    let mut votes = [0f32; 64];
+    for f in feats {
+        let h = binarize64(&f.desc);
+        for (bit, vote) in votes.iter_mut().enumerate() {
+            if (h >> (63 - bit)) & 1 == 1 {
+                *vote += 1.0;
+            }
+        }
+    }
+    if !feats.is_empty() {
+        for v in votes.iter_mut() {
+            *v /= feats.len() as f32;
+        }
+    }
+    let mut counts = [0u32; 64];
+    for f in feats {
+        counts[nearest(centroids, &f.desc)] += 1;
+    }
+    Stats { pooled, votes, counts }
+}
+
+fn variant_hashes(feats: &[Feature], stats: &Stats, calib: &Calibration) -> [u64; 8] {
+    let half = [0.5f32; 64];
+    let mut centered = [0f32; 64];
+    for i in 0..64 {
+        centered[i] = stats.pooled[i] - calib.pooled_med[i];
+    }
+    let mut simhash = 0u64;
+    for plane in &calib.planes {
+        simhash <<= 1;
+        let dot: f32 = plane.iter().zip(&centered).map(|(p, v)| p * v).sum();
+        if dot > 0.0 {
+            simhash |= 1;
+        }
+    }
+    let (mut any, mut med) = (0u64, 0u64);
+    for w in 0..64 {
+        any <<= 1;
+        med <<= 1;
+        if stats.counts[w] > 0 {
+            any |= 1;
+        }
+        if stats.counts[w] > calib.word_med[w] {
+            med |= 1;
+        }
+    }
+    [
+        binarize64(&stats.pooled),
+        hash_strongest(feats),
+        threshold_bits(&stats.votes, &half),
+        threshold_bits(&stats.pooled, &calib.pooled_med),
+        simhash,
+        threshold_bits(&stats.votes, &calib.vote_med),
+        any,
+        med,
+    ]
+}
+
+fn threshold_bits(v: &[f32; 64], t: &[f32; 64]) -> u64 {
+    let mut mask = 0u64;
+    for i in 0..64 {
+        mask <<= 1;
+        if v[i] > t[i] {
+            mask |= 1;
+        }
+    }
+    mask
+}
+
+fn column_median<F: Fn(&Stats) -> [f32; 64]>(stats: &[Stats], get: F) -> [f32; 64] {
+    let mut out = [0f32; 64];
+    for (i, o) in out.iter_mut().enumerate() {
+        let mut col: Vec<f32> = stats.iter().map(|s| get(s)[i]).collect();
+        col.sort_by(f32::total_cmp);
+        *o = col[col.len() / 2];
+    }
+    out
+}
+
+fn column_median_u32(stats: &[Stats]) -> [u32; 64] {
+    let mut out = [0u32; 64];
+    for (i, o) in out.iter_mut().enumerate() {
+        let mut col: Vec<u32> = stats.iter().map(|s| s.counts[i]).collect();
+        col.sort();
+        *o = col[col.len() / 2];
+    }
+    out
+}
+
+fn xorshift(state: &mut u64) -> u64 {
+    *state ^= *state << 13;
+    *state ^= *state >> 7;
+    *state ^= *state << 17;
+    *state
+}
+
+/// Fixed random +-1 hyperplanes, the classic sign-random-projection simhash
+fn random_planes(n: usize) -> Vec<[f32; 64]> {
+    let mut rng = 0x9E3779B97F4A7C15u64;
+    (0..n)
+        .map(|_| {
+            let mut plane = [0f32; 64];
+            for p in plane.iter_mut() {
+                *p = if xorshift(&mut rng) & 1 == 1 { 1.0 } else { -1.0 };
+            }
+            plane
+        })
+        .collect()
+}
+
+fn nearest(centroids: &[[f32; 64]], d: &[f32; 64]) -> usize {
+    let mut best = 0;
+    let mut best_dist = f32::MAX;
+    for (c, centroid) in centroids.iter().enumerate() {
+        let mut dist = 0.0;
+        for i in 0..64 {
+            let diff = centroid[i] - d[i];
+            dist += diff * diff;
+        }
+        if dist < best_dist {
+            best_dist = dist;
+            best = c;
+        }
+    }
+    best
+}
+
+fn kmeans(samples: &[[f32; 64]], k: usize, iters: usize) -> Vec<[f32; 64]> {
+    let mut rng = 0x243F6A8885A308D3u64;
+    let mut centroids: Vec<[f32; 64]> = (0..k)
+        .map(|_| samples[(xorshift(&mut rng) % samples.len() as u64) as usize])
+        .collect();
+    for _ in 0..iters {
+        let assign: Vec<usize> = samples
+            .par_iter()
+            .map(|s| nearest(&centroids, s))
+            .collect();
+        let mut sums = vec![[0f32; 64]; k];
+        let mut counts = vec![0usize; k];
+        for (s, &a) in samples.iter().zip(&assign) {
+            counts[a] += 1;
+            for i in 0..64 {
+                sums[a][i] += s[i];
+            }
+        }
+        for c in 0..k {
+            if counts[c] == 0 {
+                centroids[c] = samples[(xorshift(&mut rng) % samples.len() as u64) as usize];
+            } else {
+                for i in 0..64 {
+                    centroids[c][i] = sums[c][i] / counts[c] as f32;
+                }
+            }
+        }
+    }
+    centroids
+}
+
+/// Map a keypoint position in the mutant back to base image coordinates.
+/// Only the geometric mutations move points; everything else is identity.
+fn map_to_base(
+    tag: &str,
+    bw: f32,
+    bh: f32,
+    mw: f32,
+    mh: f32,
+    x: f32,
+    y: f32,
+    sigma: f32,
+) -> (f32, f32, f32) {
+    if tag == ".flip." {
+        return (bw - 1.0 - x, y, sigma);
+    }
+    if tag.starts_with(".crop") {
+        // crop_imm offsets, integer division like the mutator
+        let ox = ((bw - mw) / 2.0).floor();
+        let oy = ((bh - mh) / 2.0).floor();
+        return (x + ox, y + oy, sigma);
+    }
+    if tag == ".scale50." {
+        return (x * bw / mw, y * bh / mh, sigma * bw / mw);
+    }
+    if let Some(deg) = rot_degrees(tag) {
+        // imageproc warps output through the projection inverse, so a
+        // mutant point o sits at base position c + R(-theta) (o - c)
+        let theta = deg.to_radians();
+        let (s, c) = theta.sin_cos();
+        let (cx, cy) = (mw / 2.0, mh / 2.0);
+        let (dx, dy) = (x - cx, y - cy);
+        return (cx + c * dx + s * dy, cy - s * dx + c * dy, sigma);
+    }
+    (x, y, sigma)
+}
+
+fn rot_degrees(tag: &str) -> Option<f32> {
+    tag.strip_prefix(".rot")?.strip_suffix('.')?.parse().ok()
+}
+
+/// Whether a base keypoint can possibly survive the mutation, for a fair
+/// repeatability denominator (crops delete everything outside the window)
+fn in_region(tag: &str, bw: f32, bh: f32, mw: f32, mh: f32, x: f32, y: f32) -> bool {
+    if tag.starts_with(".crop") {
+        let ox = ((bw - mw) / 2.0).floor();
+        let oy = ((bh - mh) / 2.0).floor();
+        return x >= ox && x < ox + mw && y >= oy && y < oy + mh;
+    }
+    true
+}

src/bin/storekeys.rs

@@ -0,0 +1,38 @@
+//! Dumps the distinct hashes of a store as little-endian u64s.
+//! The web demo's at-scale BK-tree builds from this file.
+use clap::Parser;
+use image_similarity::descriptors::get_all_descriptors;
+use image_similarity::store::DescriptorStore;
+use std::fs;
+use std::path::PathBuf;
+
+#[derive(Parser)]
+struct Args {
+    /// Store to read, the file stem picks the descriptor
+    #[arg(default_value = "dct.store")]
+    store: PathBuf,
+    /// Output file of raw little-endian u64 keys
+    #[arg(default_value = "web/flickr/keys.bin")]
+    out: PathBuf,
+}
+
+fn main() {
+    let args = Args::parse();
+    let stem = args.store.file_stem().expect("store path has no file stem").to_string_lossy();
+    let descriptor = get_all_descriptors()
+        .into_iter()
+        .find(|d| d.info() == stem)
+        .unwrap_or_else(|| panic!("no descriptor named {stem}"));
+
+    let store = DescriptorStore::new(descriptor)
+        .with_file(&args.store)
+        .expect("loading store");
+
+    // sorted for a deterministic file; u64 order carries no hamming
+    // structure, so the resulting tree shape is as good as random
+    let mut keys: Vec<u64> = store.keys().collect();
+    keys.sort_unstable();
+    let bytes: Vec<u8> = keys.iter().flat_map(|k| k.to_le_bytes()).collect();
+    fs::write(&args.out, bytes).expect("writing keys");
+    println!("{} keys -> {}", keys.len(), args.out.display());
+}

src/bk.rs

@@ -44,6 +44,28 @@ impl Node {
 	fn count(&self) -> usize {
 		1 + self.children.iter().map(|(_, c)| c.count()).sum::<usize>()
 	}
+
+	/// find, but records the key of every node compared, in visit order
+	fn find_trace(&self, key: u64, radius: u64, found: &mut Vec<(u64, u64)>, visited: &mut Vec<u64>) {
+		visited.push(self.key);
+		let distance = hamming(self.key, key);
+		if distance <= radius {
+			found.push((self.key, distance));
+		}
+		for (child_distance, child) in &self.children {
+			let child_distance = *child_distance as u64;
+			if child_distance + radius >= distance && child_distance <= distance + radius {
+				child.find_trace(key, radius, found, visited);
+			}
+		}
+	}
+
+	fn walk(&self, parent: u64, distance: u8, out: &mut Vec<(u64, u64, u64)>) {
+		out.push((self.key, parent, distance as u64));
+		for (d, child) in &self.children {
+			child.walk(self.key, *d, out);
+		}
+	}
 }
 
 #[derive(Default)]
@@ -80,6 +102,27 @@ impl BkTree {
 		}
 		found
 	}
+
+	/// Like find, but also returns every key that was compared against,
+	/// in visit order. Keys not in the list were pruned away.
+	pub fn find_trace(&self, key: u64, radius: u64) -> (Vec<(u64, u64)>, Vec<u64>) {
+		let mut found = Vec::new();
+		let mut visited = Vec::new();
+		if let Some(root) = &self.root {
+			root.find_trace(key, radius, &mut found, &mut visited);
+		}
+		(found, visited)
+	}
+
+	/// Preorder walk of the tree as (key, parent key, distance to parent)
+	/// triples, for visualization. The root lists itself as its parent.
+	pub fn walk(&self) -> Vec<(u64, u64, u64)> {
+		let mut out = Vec::new();
+		if let Some(root) = &self.root {
+			root.walk(root.key, 0, &mut out);
+		}
+		out
+	}
 }
 
 #[cfg(test)]
@@ -132,4 +175,38 @@ mod tests {
 	fn empty_tree_finds_nothing() {
 		assert!(BkTree::new().find(1, 64).is_empty());
 	}
+
+	#[test]
+	fn trace_matches_find_and_visits_at_most_everything() {
+		let keys = lcg_keys(200);
+		let mut tree = BkTree::new();
+		for &key in &keys {
+			tree.insert(key);
+		}
+		for radius in [0, 4, 16] {
+			let (found, visited) = tree.find_trace(keys[3], radius);
+			assert_eq!(found, tree.find(keys[3], radius));
+			assert_eq!(visited[0], keys[0], "search starts at the root");
+			assert!(visited.len() <= keys.len());
+			// every match was necessarily compared against
+			for (key, _) in &found {
+				assert!(visited.contains(key));
+			}
+		}
+	}
+
+	#[test]
+	fn walk_lists_every_node_with_consistent_distances() {
+		let keys = lcg_keys(50);
+		let mut tree = BkTree::new();
+		for &key in &keys {
+			tree.insert(key);
+		}
+		let nodes = tree.walk();
+		assert_eq!(nodes.len(), tree.len());
+		assert_eq!(nodes[0], (keys[0], keys[0], 0), "root is its own parent");
+		for (key, parent, distance) in &nodes[1..] {
+			assert_eq!(hamming(*key, *parent), *distance);
+		}
+	}
 }

src/lib.rs

@@ -7,6 +7,7 @@ pub mod cache;
 pub mod descriptors;
 pub mod store;
 pub mod mutators;
+pub mod salient;
 pub mod wasm;
 
 /// Opens an image by sniffing its content instead of trusting the file

src/mutators/mod.rs

@@ -281,6 +281,9 @@ pub fn get_all_mutators() -> Vec<Box<dyn Mutator>> {
 		Box::new(Letterbox { bar_percent: 10 }),
 		Box::new(Logo { size_percent: 10 }),
 		Box::new(Rotate { degrees: 2.0 }),
+		// a hard rotation the frequency hashes cannot follow, where the
+		// salient-point descriptor's rotation invariance shows itself
+		Box::new(Rotate { degrees: 45.0 }),
 		Box::new(Noise { stddev: 10.0 }),
 	]
 }

src/salient.rs

@@ -0,0 +1,644 @@
+//! Salient-point detector and descriptor, built for hashing experiments.
+//!
+//! The detector is SIFT-like: extrema in a Difference-of-Gaussians scale
+//! space, filtered for contrast and edge responses. The descriptor is
+//! SURF-like: a 4x4 grid of subregions around the keypoint, each summarized
+//! by four gradient sums, 64 floats total. Every descriptor is measured in
+//! the keypoint's own scale and orientation frame, which is where the scale
+//! and rotation invariance comes from.
+//!
+//! Descriptors are canonicalized under horizontal mirroring: a flipped patch
+//! produces the same 64 floats, so hashes built from them inherit flip
+//! invariance the way the DCT hash does. See `mirror_descriptor` for the
+//! exact symmetry.
+//!
+//! Not wired into the Descriptor trait or the stores. The salientlab binary
+//! runs the experiments.
+
+use image::DynamicImage;
+use std::f32::consts::PI;
+
+const ORI_BINS: usize = 36;
+
+/// Descriptor sampling grid: DESC_SAMPLES x DESC_SAMPLES points spanning
+/// +-DESC_HALF_EXTENT keypoint sigmas, grouped into DESC_GRID x DESC_GRID
+/// subregions. These define the exact window the descriptor reads, so the
+/// browser demo draws the same box by calling Keypoint::window_corners.
+const DESC_SAMPLES: usize = 20;
+const DESC_HALF_EXTENT: f32 = 9.5;
+pub const DESC_GRID: usize = 4;
+
+/// Grayscale f32 image, values in [0, 1]
+pub struct GrayF32 {
+    w: usize,
+    h: usize,
+    data: Vec<f32>,
+}
+
+impl GrayF32 {
+    pub fn from_image(img: &DynamicImage) -> GrayF32 {
+        let luma = img.to_luma32f();
+        GrayF32 {
+            w: luma.width() as usize,
+            h: luma.height() as usize,
+            data: luma.into_raw(),
+        }
+    }
+
+    fn get(&self, x: isize, y: isize) -> f32 {
+        let x = x.clamp(0, self.w as isize - 1) as usize;
+        let y = y.clamp(0, self.h as isize - 1) as usize;
+        self.data[y * self.w + x]
+    }
+
+    fn bilinear(&self, x: f32, y: f32) -> f32 {
+        let x0 = x.floor();
+        let y0 = y.floor();
+        let fx = x - x0;
+        let fy = y - y0;
+        let (xi, yi) = (x0 as isize, y0 as isize);
+        let v00 = self.get(xi, yi);
+        let v10 = self.get(xi + 1, yi);
+        let v01 = self.get(xi, yi + 1);
+        let v11 = self.get(xi + 1, yi + 1);
+        v00 * (1.0 - fx) * (1.0 - fy)
+            + v10 * fx * (1.0 - fy)
+            + v01 * (1.0 - fx) * fy
+            + v11 * fx * fy
+    }
+
+    /// Image gradient at a fractional position, central differences
+    fn gradient(&self, x: f32, y: f32) -> (f32, f32) {
+        (
+            (self.bilinear(x + 1.0, y) - self.bilinear(x - 1.0, y)) * 0.5,
+            (self.bilinear(x, y + 1.0) - self.bilinear(x, y - 1.0)) * 0.5,
+        )
+    }
+
+    /// Separable gaussian blur, borders clamped
+    fn gauss_blur(&self, sigma: f32) -> GrayF32 {
+        let radius = (sigma * 3.0).ceil().max(1.0) as isize;
+        let mut kernel = Vec::with_capacity(2 * radius as usize + 1);
+        for i in -radius..=radius {
+            kernel.push((-(i * i) as f32 / (2.0 * sigma * sigma)).exp());
+        }
+        let sum: f32 = kernel.iter().sum();
+        for k in kernel.iter_mut() {
+            *k /= sum;
+        }
+
+        let mut tmp = vec![0.0f32; self.w * self.h];
+        for y in 0..self.h {
+            for x in 0..self.w {
+                let mut acc = 0.0;
+                for (ki, k) in kernel.iter().enumerate() {
+                    acc += k * self.get(x as isize + ki as isize - radius, y as isize);
+                }
+                tmp[y * self.w + x] = acc;
+            }
+        }
+        let tmp = GrayF32 { w: self.w, h: self.h, data: tmp };
+        let mut out = vec![0.0f32; self.w * self.h];
+        for y in 0..self.h {
+            for x in 0..self.w {
+                let mut acc = 0.0;
+                for (ki, k) in kernel.iter().enumerate() {
+                    acc += k * tmp.get(x as isize, y as isize + ki as isize - radius);
+                }
+                out[y * self.w + x] = acc;
+            }
+        }
+        GrayF32 { w: self.w, h: self.h, data: out }
+    }
+
+    /// Every other pixel, for the next octave
+    fn half(&self) -> GrayF32 {
+        let w = (self.w / 2).max(1);
+        let h = (self.h / 2).max(1);
+        let mut data = Vec::with_capacity(w * h);
+        for y in 0..h {
+            for x in 0..w {
+                data.push(self.data[(y * 2) * self.w + x * 2]);
+            }
+        }
+        GrayF32 { w, h, data }
+    }
+
+    fn sub(a: &GrayF32, b: &GrayF32) -> GrayF32 {
+        let data = a.data.iter().zip(&b.data).map(|(a, b)| a - b).collect();
+        GrayF32 { w: a.w, h: a.h, data }
+    }
+}
+
+#[derive(Clone, Copy, Debug)]
+pub struct Keypoint {
+    /// Position in original image coordinates
+    pub x: f32,
+    pub y: f32,
+    /// Scale in original image coordinates
+    pub sigma: f32,
+    /// Dominant gradient orientation, radians
+    pub angle: f32,
+    /// Absolute DoG value at the extremum
+    pub response: f32,
+}
+
+impl Keypoint {
+    /// Map a point from the keypoint's oriented frame (u along the
+    /// orientation, v perpendicular) into image coordinates. Same rotation
+    /// the descriptor uses to place its samples.
+    fn to_image(&self, u: f32, v: f32) -> (f32, f32) {
+        let (sin, cos) = self.angle.sin_cos();
+        (self.x + u * cos - v * sin, self.y + u * sin + v * cos)
+    }
+
+    /// The four corners of the oriented descriptor window in image
+    /// coordinates, from the (-u,-v) corner counter-clockwise. For drawing
+    /// the region the descriptor summarizes.
+    pub fn window_corners(&self) -> [(f32, f32); 4] {
+        let e = DESC_HALF_EXTENT * self.sigma;
+        [(-e, -e), (e, -e), (e, e), (-e, e)].map(|(u, v)| self.to_image(u, v))
+    }
+
+    /// Centres of the DESC_GRID x DESC_GRID subregions in image coordinates,
+    /// in descriptor order (row = perpendicular to orientation, column =
+    /// along it), so a caller can line them up with the 64-float descriptor.
+    pub fn subregion_centers(&self) -> [(f32, f32); DESC_GRID * DESC_GRID] {
+        let cell = DESC_SAMPLES / DESC_GRID;
+        let mut out = [(0.0, 0.0); DESC_GRID * DESC_GRID];
+        for j in 0..DESC_GRID {
+            for i in 0..DESC_GRID {
+                let su = (i * cell + cell / 2) as f32 - DESC_HALF_EXTENT;
+                let sv = (j * cell + cell / 2) as f32 - DESC_HALF_EXTENT;
+                out[j * DESC_GRID + i] = self.to_image(su * self.sigma, sv * self.sigma);
+            }
+        }
+        out
+    }
+}
+
+#[derive(Clone)]
+pub struct Feature {
+    pub kp: Keypoint,
+    pub desc: [f32; 64],
+}
+
+pub struct Params {
+    /// DoG scale samples per octave
+    pub intervals: usize,
+    /// Blur of the first level of each octave
+    pub base_sigma: f32,
+    /// Minimum absolute DoG value for a keypoint (image values are 0..1)
+    pub contrast_threshold: f32,
+    /// Maximum principal curvature ratio, rejects points on straight edges
+    pub edge_ratio: f32,
+    /// Keep only the strongest keypoints
+    pub max_features: usize,
+    /// Canonicalize descriptors under horizontal mirroring
+    pub flip_canonical: bool,
+}
+
+impl Default for Params {
+    fn default() -> Self {
+        Params {
+            intervals: 3,
+            base_sigma: 1.6,
+            contrast_threshold: 0.01,
+            edge_ratio: 10.0,
+            max_features: 100,
+            flip_canonical: true,
+        }
+    }
+}
+
+struct RawKp {
+    octave: usize,
+    level: usize,
+    x: usize,
+    y: usize,
+    bx: f32,
+    by: f32,
+    bsigma: f32,
+    response: f32,
+}
+
+pub struct Salient {
+    pub params: Params,
+}
+
+impl Default for Salient {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl Salient {
+    pub fn new() -> Salient {
+        Salient { params: Params::default() }
+    }
+
+    fn k(&self) -> f32 {
+        2f32.powf(1.0 / self.params.intervals as f32)
+    }
+
+    /// Gaussian pyramid: per octave `intervals + 3` blur levels, so the DoG
+    /// stack has `intervals` usable middle layers
+    fn pyramid(&self, img: &DynamicImage) -> Vec<Vec<GrayF32>> {
+        let s = self.params.intervals;
+        let k = self.k();
+        let sigma0 = self.params.base_sigma;
+        // assume the input carries sigma 0.5, lift it to base_sigma
+        let gray = GrayF32::from_image(img);
+        let mut base = gray.gauss_blur((sigma0 * sigma0 - 0.25).max(0.01).sqrt());
+
+        let mut octaves = Vec::new();
+        while base.w.min(base.h) >= 24 && octaves.len() < 6 {
+            let mut levels = vec![base];
+            for i in 1..s + 3 {
+                let sig_prev = sigma0 * k.powi(i as i32 - 1);
+                let sig_total = sigma0 * k.powi(i as i32);
+                let sig_diff = (sig_total * sig_total - sig_prev * sig_prev).sqrt();
+                levels.push(levels[i - 1].gauss_blur(sig_diff));
+            }
+            // levels[s] carries 2 * base_sigma: halved it seeds the next octave
+            base = levels[s].half();
+            octaves.push(levels);
+        }
+        octaves
+    }
+
+    fn detect(&self, pyramid: &[Vec<GrayF32>]) -> Vec<RawKp> {
+        let k = self.k();
+        let mut raw = Vec::new();
+        for (o, levels) in pyramid.iter().enumerate() {
+            let dogs: Vec<GrayF32> = levels
+                .windows(2)
+                .map(|w| GrayF32::sub(&w[1], &w[0]))
+                .collect();
+            let scale_up = (1usize << o) as f32;
+            for l in 1..dogs.len() - 1 {
+                let d = &dogs[l];
+                for y in 1..d.h - 1 {
+                    for x in 1..d.w - 1 {
+                        let v = d.data[y * d.w + x];
+                        if v.abs() < self.params.contrast_threshold {
+                            continue;
+                        }
+                        if !is_extremum(&dogs, l, x, y, v) {
+                            continue;
+                        }
+                        if is_edge(d, x, y, self.params.edge_ratio) {
+                            continue;
+                        }
+                        let osigma = self.params.base_sigma * k.powi(l as i32);
+                        raw.push(RawKp {
+                            octave: o,
+                            level: l,
+                            x,
+                            y,
+                            bx: x as f32 * scale_up,
+                            by: y as f32 * scale_up,
+                            bsigma: osigma * scale_up,
+                            response: v.abs(),
+                        });
+                    }
+                }
+            }
+        }
+        raw
+    }
+
+    /// Detect, orient and describe the strongest keypoints of an image
+    pub fn features(&self, img: &DynamicImage) -> Vec<Feature> {
+        let pyramid = self.pyramid(img);
+        let mut raw = self.detect(&pyramid);
+        raw.sort_by(|a, b| b.response.total_cmp(&a.response));
+
+        // greedy dedupe: extrema often fire in adjacent scale layers
+        let mut picked: Vec<RawKp> = Vec::new();
+        for kp in raw {
+            if picked.len() >= self.params.max_features {
+                break;
+            }
+            let dup = picked.iter().any(|p| {
+                let (dx, dy) = (p.bx - kp.bx, p.by - kp.by);
+                let ratio = (p.bsigma / kp.bsigma).max(kp.bsigma / p.bsigma);
+                dx * dx + dy * dy < 4.0 && ratio < 1.6
+            });
+            if !dup {
+                picked.push(kp);
+            }
+        }
+
+        picked
+            .iter()
+            .map(|kp| {
+                let img = &pyramid[kp.octave][kp.level];
+                let osigma = self.params.base_sigma * self.k().powi(kp.level as i32);
+                let angle = orientation(img, kp.x, kp.y, osigma);
+                let mut desc = describe(img, kp.x as f32, kp.y as f32, osigma, angle);
+                if self.params.flip_canonical {
+                    desc = canonicalize(desc);
+                }
+                Feature {
+                    kp: Keypoint {
+                        x: kp.bx,
+                        y: kp.by,
+                        sigma: kp.bsigma,
+                        angle,
+                        response: kp.response,
+                    },
+                    desc,
+                }
+            })
+            .collect()
+    }
+}
+
+/// Strict local extremum over the 26 neighbors in space and scale
+fn is_extremum(dogs: &[GrayF32], l: usize, x: usize, y: usize, v: f32) -> bool {
+    for layer in &dogs[l - 1..=l + 1] {
+        for dy in -1isize..=1 {
+            for dx in -1isize..=1 {
+                let n = layer.get(x as isize + dx, y as isize + dy);
+                if std::ptr::eq(layer, &dogs[l]) && dx == 0 && dy == 0 {
+                    continue;
+                }
+                if (v > 0.0 && n >= v) || (v < 0.0 && n <= v) {
+                    return false;
+                }
+            }
+        }
+    }
+    true
+}
+
+/// Principal curvature ratio test on the 2x2 spatial Hessian of the DoG,
+/// rejects points that sit on an edge and slide along it
+fn is_edge(d: &GrayF32, x: usize, y: usize, r: f32) -> bool {
+    let (x, y) = (x as isize, y as isize);
+    let v = d.get(x, y);
+    let dxx = d.get(x + 1, y) - 2.0 * v + d.get(x - 1, y);
+    let dyy = d.get(x, y + 1) - 2.0 * v + d.get(x, y - 1);
+    let dxy = (d.get(x + 1, y + 1) - d.get(x + 1, y - 1) - d.get(x - 1, y + 1)
+        + d.get(x - 1, y - 1))
+        / 4.0;
+    let tr = dxx + dyy;
+    let det = dxx * dyy - dxy * dxy;
+    det <= 0.0 || tr * tr * r >= (r + 1.0) * (r + 1.0) * det
+}
+
+/// Dominant gradient orientation: 36-bin histogram of gradient angles in a
+/// gaussian-weighted neighborhood, smoothed, peak refined by a parabola
+fn orientation(img: &GrayF32, x: usize, y: usize, sigma: f32) -> f32 {
+    let sig_w = 1.5 * sigma;
+    let radius = (3.0 * sig_w).round() as isize;
+    let mut hist = [0f32; ORI_BINS];
+    for dy in -radius..=radius {
+        for dx in -radius..=radius {
+            let (px, py) = (x as isize + dx, y as isize + dy);
+            let gx = (img.get(px + 1, py) - img.get(px - 1, py)) * 0.5;
+            let gy = (img.get(px, py + 1) - img.get(px, py - 1)) * 0.5;
+            let mag = (gx * gx + gy * gy).sqrt();
+            if mag == 0.0 {
+                continue;
+            }
+            let mut phi = gy.atan2(gx);
+            if phi < 0.0 {
+                phi += 2.0 * PI;
+            }
+            let w = (-((dx * dx + dy * dy) as f32) / (2.0 * sig_w * sig_w)).exp();
+            let bin = ((phi / (2.0 * PI) * ORI_BINS as f32) as usize).min(ORI_BINS - 1);
+            hist[bin] += w * mag;
+        }
+    }
+    // two passes of circular [1 2 1] smoothing
+    for _ in 0..2 {
+        let orig = hist;
+        for b in 0..ORI_BINS {
+            let l = orig[(b + ORI_BINS - 1) % ORI_BINS];
+            let r = orig[(b + 1) % ORI_BINS];
+            hist[b] = (l + 2.0 * orig[b] + r) / 4.0;
+        }
+    }
+    let peak = (0..ORI_BINS)
+        .max_by(|a, b| hist[*a].total_cmp(&hist[*b]))
+        .unwrap();
+    let l = hist[(peak + ORI_BINS - 1) % ORI_BINS];
+    let c = hist[peak];
+    let r = hist[(peak + 1) % ORI_BINS];
+    let denom = l - 2.0 * c + r;
+    let delta = if denom.abs() > 1e-12 { 0.5 * (l - r) / denom } else { 0.0 };
+    (peak as f32 + 0.5 + delta) * 2.0 * PI / ORI_BINS as f32
+}
+
+/// SURF-style descriptor: 20x20 samples spaced `sigma` apart, rotated to the
+/// keypoint orientation, grouped into 4x4 subregions. Each subregion gets
+/// [sum du, sum dv, sum |du|, sum |dv|] where (du, dv) is the gradient
+/// rotated into the keypoint frame. Unit normalized.
+///
+/// Layout: index = (j * 4 + i) * 4 + c, with i along the orientation axis,
+/// j perpendicular to it, c the component.
+fn describe(img: &GrayF32, x: f32, y: f32, sigma: f32, angle: f32) -> [f32; 64] {
+    let (sin, cos) = angle.sin_cos();
+    let sig_w = 3.3 * sigma;
+    let cell = DESC_SAMPLES / DESC_GRID;
+    let mut desc = [0f32; 64];
+    for sj in 0..DESC_SAMPLES {
+        for si in 0..DESC_SAMPLES {
+            let u = (si as f32 - DESC_HALF_EXTENT) * sigma;
+            let v = (sj as f32 - DESC_HALF_EXTENT) * sigma;
+            let px = x + u * cos - v * sin;
+            let py = y + u * sin + v * cos;
+            let (gx, gy) = img.gradient(px, py);
+            let du = cos * gx + sin * gy;
+            let dv = -sin * gx + cos * gy;
+            let w = (-(u * u + v * v) / (2.0 * sig_w * sig_w)).exp();
+            let base = ((sj / cell) * DESC_GRID + si / cell) * 4;
+            desc[base] += w * du;
+            desc[base + 1] += w * dv;
+            desc[base + 2] += w * du.abs();
+            desc[base + 3] += w * dv.abs();
+        }
+    }
+    let norm: f32 = desc.iter().map(|v| v * v).sum::<f32>().sqrt();
+    if norm > 0.0 {
+        for v in desc.iter_mut() {
+            *v /= norm;
+        }
+    }
+    desc
+}
+
+/// The descriptor a horizontally mirrored copy of the patch would produce.
+///
+/// Under a mirror the dominant orientation maps to pi - theta. Working out
+/// the sample grid in that frame: the u axis is preserved, the v axis
+/// negates. So subregion (i, j) swaps with (i, 3 - j), the dv sum changes
+/// sign, and the other three components are untouched.
+pub fn mirror_descriptor(d: &[f32; 64]) -> [f32; 64] {
+    let mut out = [0f32; 64];
+    for j in 0..4 {
+        for i in 0..4 {
+            let src = (j * 4 + i) * 4;
+            let dst = ((3 - j) * 4 + i) * 4;
+            out[dst] = d[src];
+            out[dst + 1] = -d[src + 1];
+            out[dst + 2] = d[src + 2];
+            out[dst + 3] = d[src + 3];
+        }
+    }
+    out
+}
+
+/// Pick the mirror representative deterministically, so a patch and its
+/// mirror image canonicalize to the same descriptor
+fn canonicalize(d: [f32; 64]) -> [f32; 64] {
+    let dv_total: f32 = (0..16).map(|s| d[s * 4 + 1]).sum();
+    let mirror = if dv_total.abs() > 1e-6 {
+        dv_total < 0.0
+    } else {
+        // fallback: skew of the |dv| mass along the v axis
+        let skew: f32 = (0..16).map(|s| ((s / 4) as f32 - 1.5) * d[s * 4 + 3]).sum();
+        skew < 0.0
+    };
+    if mirror {
+        mirror_descriptor(&d)
+    } else {
+        d
+    }
+}
+
+/// 64 bits from one descriptor: the 32 signed sums contribute their sign,
+/// the 32 magnitude sums contribute a comparison against their own median.
+/// Same recipe as the DCT hash: signs plus ordinal facts, nothing absolute.
+pub fn binarize64(desc: &[f32; 64]) -> u64 {
+    let mut mags: Vec<f32> = (0..16)
+        .flat_map(|s| [desc[s * 4 + 2], desc[s * 4 + 3]])
+        .collect();
+    mags.sort_by(f32::total_cmp);
+    let median = mags[16];
+    let mut mask = 0u64;
+    for (i, &v) in desc.iter().enumerate() {
+        mask <<= 1;
+        let bit = if i % 4 < 2 { v > 0.0 } else { v > median };
+        if bit {
+            mask |= 1;
+        }
+    }
+    mask
+}
+
+/// Response-weighted mean of all descriptors, unit normalized.
+/// Symmetric pooling, so it inherits every per-keypoint invariance.
+pub fn pooled_vector(feats: &[Feature]) -> [f32; 64] {
+    let mut acc = [0f32; 64];
+    for f in feats {
+        for (a, d) in acc.iter_mut().zip(&f.desc) {
+            *a += f.kp.response * d;
+        }
+    }
+    let norm: f32 = acc.iter().map(|v| v * v).sum::<f32>().sqrt();
+    if norm > 0.0 {
+        for v in acc.iter_mut() {
+            *v /= norm;
+        }
+    }
+    acc
+}
+
+/// Variant 1: binarized pooled descriptor
+pub fn hash_pooled(feats: &[Feature]) -> u64 {
+    binarize64(&pooled_vector(feats))
+}
+
+/// Variant 2: binarized descriptor of the single strongest keypoint
+pub fn hash_strongest(feats: &[Feature]) -> u64 {
+    match feats.first() {
+        Some(f) => binarize64(&f.desc),
+        None => 0,
+    }
+}
+
+/// Variant 3: per-bit majority vote over all binarized descriptors
+pub fn hash_consensus(feats: &[Feature]) -> u64 {
+    if feats.is_empty() {
+        return 0;
+    }
+    let hashes: Vec<u64> = feats.iter().map(|f| binarize64(&f.desc)).collect();
+    let mut mask = 0u64;
+    for bit in (0..64).rev() {
+        let votes = hashes.iter().filter(|h| (*h >> bit) & 1 == 1).count();
+        mask <<= 1;
+        if votes * 2 > hashes.len() {
+            mask |= 1;
+        }
+    }
+    mask
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::descriptors::testimg;
+
+    fn dist(a: u64, b: u64) -> u32 {
+        (a ^ b).count_ones()
+    }
+
+    #[test]
+    fn mirror_is_an_involution() {
+        let mut d = [0f32; 64];
+        for (i, v) in d.iter_mut().enumerate() {
+            *v = (i as f32 * 0.37).sin();
+        }
+        assert_eq!(mirror_descriptor(&mirror_descriptor(&d)), d);
+    }
+
+    #[test]
+    fn canonical_ignores_mirroring() {
+        let mut d = [0f32; 64];
+        for (i, v) in d.iter_mut().enumerate() {
+            *v = (i as f32 * 0.37).sin();
+        }
+        assert_eq!(canonicalize(d), canonicalize(mirror_descriptor(&d)));
+    }
+
+    #[test]
+    fn binarize_packs_all_bits() {
+        // strongest positive signal on the first and last positions
+        let mut d = [0f32; 64];
+        d[0] = 1.0; // first sign bit, lands on bit 63
+        d[63] = 1.0; // a magnitude above the (zero) median, lands on bit 0
+        let h = binarize64(&d);
+        assert_eq!(h, (1 << 63) | 1);
+    }
+
+    #[test]
+    fn features_deterministic() {
+        // synthetic gradients are too smooth to produce DoG extrema
+        let img = image::open("img/meowl.jpg").unwrap();
+        let s = Salient::new();
+        let a = s.features(&img);
+        let b = s.features(&img);
+        assert_eq!(a.len(), b.len());
+        assert!(!a.is_empty());
+        for (fa, fb) in a.iter().zip(&b) {
+            assert_eq!(fa.desc, fb.desc);
+        }
+    }
+
+    #[test]
+    fn flip_invariant_on_real_photo() {
+        let img = image::open("img/meowl.jpg").unwrap();
+        let s = Salient::new();
+        let a = s.features(&img);
+        let b = s.features(&img.fliph());
+        assert!(!a.is_empty());
+        // detection mirrors exactly, orientation binning wobbles a little
+        let d = dist(hash_pooled(&a), hash_pooled(&b));
+        assert!(d <= 6, "pooled flip distance {d}");
+        let d = dist(hash_consensus(&a), hash_consensus(&b));
+        assert!(d <= 6, "consensus flip distance {d}");
+    }
+}

src/store.rs

@@ -262,6 +262,11 @@ impl DescriptorStore {
 		self.names.len()
 	}
 
+	/// All distinct hashes in the store
+	pub fn keys(&self) -> impl Iterator<Item = u64> + '_ {
+		self.map.keys().copied()
+	}
+
 	pub fn is_empty(&self) -> bool {
 		self.names.is_empty()
 	}

src/wasm/mod.rs

@@ -1,7 +1,9 @@
 //! Wasm bindings for the browser demos.
 //! Hashes come back as BigInt, images go in and out as encoded bytes.
+use crate::bk::BkTree;
 use crate::descriptors::{Descriptor, Median, DCT, PHash};
 use crate::mutators::{self, Mutator};
+use crate::salient::{self, Params, Salient};
 use wasm_bindgen::prelude::*;
 use std::io::Cursor;
 
@@ -71,6 +73,227 @@ pub fn phash_lowfreq(bytes: &[u8]) -> Result<Vec<f64>, JsError> {
 	Ok(PHash.lowfreq(&img).to_vec())
 }
 
+/// DCT hash of every mutant in the standard 16-mutator suite.
+/// The image is shrunk first so the heavy mutators stay quick, which can
+/// wobble the hashes a bit or two compared to a full-size run.
+#[wasm_bindgen]
+pub fn suite_hashes(bytes: &[u8]) -> Result<Vec<u64>, JsError> {
+	let img = load(bytes)?;
+	let img = if img.width().max(img.height()) > 512 {
+		img.thumbnail(512, 512)
+	} else {
+		img
+	};
+	let dct = DCT::new();
+	Ok(mutators::get_all_mutators()
+		.iter()
+		.map(|m| dct.describe(&m.mutate(&img)))
+		.collect())
+}
+
+/// Labels of the suite mutators, newline separated, same order as suite_hashes
+#[wasm_bindgen]
+pub fn suite_labels() -> String {
+	mutators::get_all_mutators()
+		.iter()
+		.map(|m| m.info())
+		.collect::<Vec<_>>()
+		.join("\n")
+}
+
+/// A queryable set of hashes: the BK-tree plus a flat copy of the keys
+/// as a linear scan baseline.
+#[wasm_bindgen]
+pub struct BkIndex {
+	tree: BkTree,
+	keys: Vec<u64>,
+}
+
+#[wasm_bindgen]
+impl BkIndex {
+	#[wasm_bindgen(constructor)]
+	pub fn new() -> BkIndex {
+		BkIndex { tree: BkTree::new(), keys: Vec::new() }
+	}
+
+	/// Inserts one key, ignoring duplicates. Fine for demo-sized trees,
+	/// use insert_bytes for bulk loads.
+	pub fn insert(&mut self, key: u64) {
+		if !self.keys.contains(&key) {
+			self.keys.push(key);
+			self.tree.insert(key);
+		}
+	}
+
+	/// Bulk load of little-endian u64 keys (the keys.bin format),
+	/// assumed to be distinct.
+	pub fn insert_bytes(&mut self, bytes: &[u8]) {
+		for chunk in bytes.chunks_exact(8) {
+			let key = u64::from_le_bytes(chunk.try_into().unwrap());
+			self.keys.push(key);
+			self.tree.insert(key);
+		}
+	}
+
+	pub fn len(&self) -> u32 {
+		self.keys.len() as u32
+	}
+
+	/// All keys within the radius, as flattened (key, distance) pairs
+	pub fn find(&self, key: u64, radius: u32) -> Vec<u64> {
+		self.tree
+			.find(key, radius as u64)
+			.into_iter()
+			.flat_map(|(k, d)| [k, d])
+			.collect()
+	}
+
+	/// Number of hashes the tree compared against during a find
+	pub fn find_compared(&self, key: u64, radius: u32) -> u32 {
+		self.tree.find_trace(key, radius as u64).1.len() as u32
+	}
+
+	/// The keys compared against during a find, in visit order
+	pub fn trace(&self, key: u64, radius: u32) -> Vec<u64> {
+		self.tree.find_trace(key, radius as u64).1
+	}
+
+	/// Linear scan baseline, same result as find
+	pub fn scan(&self, key: u64, radius: u32) -> Vec<u64> {
+		let radius = radius as u64;
+		self.keys
+			.iter()
+			.map(|&k| (k, ((k ^ key).count_ones()) as u64))
+			.filter(|&(_, d)| d <= radius)
+			.flat_map(|(k, d)| [k, d])
+			.collect()
+	}
+
+	/// The tree as flattened (key, parent key, distance to parent) triples
+	/// in preorder; the root lists itself as parent.
+	pub fn structure(&self) -> Vec<u64> {
+		self.tree
+			.walk()
+			.into_iter()
+			.flat_map(|(k, p, d)| [k, p, d])
+			.collect()
+	}
+}
+
+impl Default for BkIndex {
+	fn default() -> Self {
+		Self::new()
+	}
+}
+
+/// Detected salient points of an image plus their descriptors, for the
+/// part-two demos. Built once, then queried for the overlay, the ranking
+/// and the per-point neighbourhood description. Coordinates are in the
+/// decoded image's own pixels (see width/height), scale them to the canvas.
+#[wasm_bindgen]
+pub struct SalientImage {
+    feats: Vec<salient::Feature>,
+    width: u32,
+    height: u32,
+}
+
+#[wasm_bindgen]
+impl SalientImage {
+    /// Detect and describe up to `max_features` points, strongest first.
+    #[wasm_bindgen(constructor)]
+    pub fn new(bytes: &[u8], max_features: u32) -> Result<SalientImage, JsError> {
+        let img = load(bytes)?;
+        let detector = Salient {
+            params: Params {
+                max_features: max_features.max(1) as usize,
+                ..Default::default()
+            },
+        };
+        let feats = detector.features(&img);
+        Ok(SalientImage { feats, width: img.width(), height: img.height() })
+    }
+
+    /// Number of keypoints kept.
+    pub fn len(&self) -> u32 {
+        self.feats.len() as u32
+    }
+
+    pub fn is_empty(&self) -> bool {
+        self.feats.is_empty()
+    }
+
+    /// Width of the coordinate space the keypoints live in.
+    pub fn width(&self) -> u32 {
+        self.width
+    }
+
+    /// Height of the coordinate space the keypoints live in.
+    pub fn height(&self) -> u32 {
+        self.height
+    }
+
+    /// Every keypoint as [x, y, sigma, angle, response], strongest first.
+    /// The response is the ranking score (how blob-like and high-contrast
+    /// the point is); it is what "importance" means here.
+    pub fn keypoints(&self) -> Vec<f32> {
+        self.feats
+            .iter()
+            .flat_map(|f| [f.kp.x, f.kp.y, f.kp.sigma, f.kp.angle, f.kp.response])
+            .collect()
+    }
+
+    /// The 64-float descriptor of keypoint `i`: 16 subregions, each
+    /// [sum du, sum dv, sum |du|, sum |dv|] in the point's oriented frame.
+    pub fn descriptor(&self, i: u32) -> Vec<f32> {
+        self.feats
+            .get(i as usize)
+            .map(|f| f.desc.to_vec())
+            .unwrap_or_default()
+    }
+
+    /// The four corners of keypoint `i`'s oriented sampling window in image
+    /// coordinates, as [x, y] x 4. Draw this to show the described region.
+    pub fn window(&self, i: u32) -> Vec<f32> {
+        self.feats
+            .get(i as usize)
+            .map(|f| {
+                f.kp.window_corners()
+                    .iter()
+                    .flat_map(|&(x, y)| [x, y])
+                    .collect()
+            })
+            .unwrap_or_default()
+    }
+
+    /// The 16 subregion centres of keypoint `i` in image coordinates, as
+    /// [x, y] x 16, in the same order as the descriptor's 16 blocks.
+    pub fn subregions(&self, i: u32) -> Vec<f32> {
+        self.feats
+            .get(i as usize)
+            .map(|f| {
+                f.kp.subregion_centers()
+                    .iter()
+                    .flat_map(|&(x, y)| [x, y])
+                    .collect()
+            })
+            .unwrap_or_default()
+    }
+
+    /// The 64-bit strongest-keypoint hash (the variant that survived at 25k).
+    pub fn hash(&self) -> u64 {
+        salient::hash_strongest(&self.feats)
+    }
+
+    /// The 64-bit hash of an arbitrary keypoint `i`, so a demo can show a
+    /// specific point's fingerprint rather than only the strongest one.
+    pub fn descriptor_hash(&self, i: u32) -> u64 {
+        self.feats
+            .get(i as usize)
+            .map(|f| salient::binarize64(&f.desc))
+            .unwrap_or(0)
+    }
+}
+
 /// Apply a single mutator to an image, returns PNG bytes.
 /// The meaning of `amount` depends on the kind.
 #[wasm_bindgen]

web/style.css

@@ -2,19 +2,33 @@ body {
 	margin: 0;
 	padding: 0;
 	color: var(--foreground);
-	background: var(--background) url('bg.jpeg');
-	background-size: cover;
+	background: var(--background);
 	font-family: sans-serif;
-	height: 100%;
-	overflow: hidden;
 	font-size: 14px;
 }
 
-#wrapper {
-	overflow-x: hidden;
-	overflow-y: auto;
+/* Fixed background layer, same construction as the landing page */
+
+#bg {
+	position: fixed;
+	width: 100vw;
+	height: 100vh;
+	z-index: -1;
+	background: var(--background);
+}
+
+#bg img {
+	display: block;
 	width: 100%;
 	height: 100%;
+	object-fit: cover;
+	object-position: 50% 50%;
+	opacity: 0;
+	transition: opacity 3s linear;
+}
+
+#bg img.loaded {
+	opacity: 1;
 }
 
 h1, h2, h3, h4, h5, h6 {
@@ -25,7 +39,7 @@ h1, h2, h3, h4, h5, h6 {
 
 h1 {
 	margin-top: 0;
-	font-size: 24px;
+	font-size: 28px;
 }
 
 h2 {
@@ -34,45 +48,98 @@ h2 {
 
 p {
 	margin-top: 2px;
+	line-height: 1.5;
 }
 
 a {
 	color: var(--color2);
+	text-decoration: none;
+	font-weight: bold;
 }
 
-strong {
-	color: var(--color4);
+a:hover {
+	color: var(--color3);
 }
 
 em {
 	color: var(--color3);
 }
 
+sup.fn {
+	font-size: 10px;
+}
+
+sup.fn a {
+	font-weight: normal;
+}
+
 pre {
-	margin: 0;
+	margin: 8px 0;
+	max-width: 640px;
 }
 
 pre > code {
-	border: 1px solid;
+	border: 1px solid color-mix(in srgb, var(--foreground) 25%, transparent);
+	display: block;
+	padding: 8px 12px;
+	overflow-x: auto;
+	font-size: 13px;
 }
 
 article {
 	max-width: 1024px;
-	min-height: 128px;
 	margin: 25px auto;
 	padding: 20px;
 	background: color-mix(in srgb, var(--background) 90%, transparent);
 	box-shadow: 5px 5px 10px rgba(0, 0, 0, 0.9);
 }
 
-code.block {
-	display: block;
+@media (max-width: 1064px) {
+	article {
+		margin: 15px 12px;
+	}
 }
 
 img {
 	max-width: 100%;
 }
 
+math[display="block"] {
+	margin: 12px 0 12px 24px;
+	font-size: 17px;
+	color: var(--color4);
+}
+
+aside.explainer {
+	border-left: 3px solid var(--color4);
+	background: color-mix(in srgb, var(--foreground) 5%, transparent);
+	padding: 8px 12px;
+	margin: 12px 0;
+	font-size: 13px;
+}
+
+.rewrite-target {
+	text-decoration: underline dashed var(--color3);
+	text-underline-offset: 3px;
+}
+
+.rewrite-target::after {
+	content: ' ✎';
+	color: var(--color3);
+	font-size: 11px;
+}
+
+.attribution {
+	font-size: 12px;
+	opacity: 0.7;
+}
+
+.big-number {
+	font-size: 26px;
+	color: var(--color1);
+	margin: 8px 0 4px 0;
+}
+
 /* Demo */
 
 .needs-image {
@@ -88,135 +155,107 @@ img {
 	min-width: 256px;
 }
 
-/* Resize demo */
+.sample {
+	height: 48px;
+	margin: 0 4px 4px 0;
+	cursor: pointer;
+	vertical-align: middle;
+}
 
-#resize {
+.sample:hover {
+	outline: 2px solid var(--color1);
+}
+
+/* Comparison slider */
+
+.compare {
 	position: relative;
-	height: 512px;
+	max-width: 512px;
+	overflow: hidden;
+	touch-action: none;
+	cursor: ew-resize;
+	user-select: none;
+}
+
+.compare img {
+	display: block;
 	width: 100%;
-	transition: height 2s linear;
+	pointer-events: none;
 }
 
-#resize.resize {
-	height: 256px;
-}
-
-#resize .image-original {
+/* The top image stretches to the base image's box, so differently
+   sized results (8x8 thumbs, crops) stay comparable under the divider */
+.compare .compare-b {
 	position: absolute;
-	width: 100%;
-	height: 512px;
+	top: 0;
+	left: 0;
+	height: 100%;
+	object-fit: fill;
+	clip-path: inset(0 0 0 50%);
 }
 
-#resize .image-original img {
-	max-height: 50vh;
-	min-height: 257px;
-	max-width: 100%;
-	min-width: 257px;
-	/* height: 512px; */
-}
-
-#resize-original {
-	position: relative;
-	height: 512px;
-}
-
-#resize-original img {
+.compare .compare-divider {
 	position: absolute;
-	max-width: 90vw;
-	max-height: 512px;
-	min-width: 257px;
-	min-height: 257px;
-	opacity: 1;
-	transition:
-		width 2s linear,
-		height 2s linear,
-		filter 3s linear 2s,
-		opacity 3s ease-in 5s;
+	top: 0;
+	bottom: 0;
+	left: 50%;
+	width: 2px;
+	margin-left: -1px;
+	background: var(--color1);
 }
 
-#resize.resize #resize-original img {
-	width: 256px !important;
-	height: 256px !important;
-	filter: grayscale(1) blur(10px);
-	opacity: 0;
-}
-
-#resize-original img.grayscale {
-	filter: grayscale(1) blur(10px);
-}
-
-#resize-original img.fade {
-	opacity: 0;
-}
-
-#image-resize img {
-	width: 256px;
-	height: 256px;
-	image-rendering: pixelated;
-	image-rendering: -moz-crisp-edges;
-}
-
-#image-dct {
-	margin-left: 32px;
-	display:grid;
-	grid-template-columns: repeat(8, 32px);
-}
-
-#image-dct > div {
-	height: 32px;
-	line-height:32px;
-	font-size: 11px;
+.compare .compare-handle {
+	position: absolute;
+	top: 50%;
+	left: 50%;
+	transform: translate(-50%, -50%);
+	width: 28px;
+	height: 28px;
+	line-height: 28px;
 	text-align: center;
-	border: 1px solid var(--color1);
-	box-sizing: border-box;
-}
-
-.similarity-example {
-	display: flex;
-	flex-flow: row nowrap;
-	align-items: flex-end;
-}
-
-figure {
-	margin: 10px;
-	width: 50%;
-}
-
-/* Loader */ 
-
-#loader.show {
 	border-radius: 50%;
-  width: 256px;
-  height: 256px;
-  font-size: 10px;
-  position: fixed;
-  z-index: 99;
-  top: calc(50vh - 128px);
-  left: calc(50vw - 128px); 
-  border: 16px solid transparent;
-  border-top: 16px solid darkred;
-  transform: translateZ(0);
-  animation: loader 2s infinite ease-in-out;
+	background: var(--color1);
+	color: var(--background);
+	font-size: 15px;
 }
 
-@keyframes loader {
-  0% {
-    transform: rotate(0deg);
-  }
-  50% {
-    transform: rotate(360deg);
-  }
-  100% {
-    transform: rotate(360deg);
-  }
+.compare .compare-label {
+	position: absolute;
+	top: 6px;
+	padding: 2px 6px;
+	font-size: 11px;
+	background: color-mix(in srgb, var(--background) 75%, transparent);
+	color: var(--foreground);
 }
-/* Skeleton */
 
-.todo {
-	opacity: 0.55;
-	font-style: italic;
-	border-left: 3px solid var(--color3);
-	padding-left: 8px;
+.compare .compare-label.a {
+	left: 6px;
+}
+
+.compare .compare-label.b {
+	right: 6px;
+}
+
+/* Busy indicator while the worker renders a mutation */
+
+.compare.busy::after {
+	content: 'mutating…';
+	position: absolute;
+	right: 8px;
+	bottom: 8px;
+	padding: 3px 9px;
+	background: color-mix(in srgb, var(--background) 80%, transparent);
+	color: var(--color1);
+	font-size: 12px;
+	animation: busy-pulse 1s ease-in-out infinite;
+}
+
+.distances.busy td {
+	opacity: 0.35;
+}
+
+@keyframes busy-pulse {
+	50% { opacity: 0.35; }
 }
 
 /* Hash visualisation */
@@ -233,6 +272,14 @@ figure {
 	letter-spacing: 2px;
 }
 
+.hash-value.small {
+	display: block;
+	margin-top: 4px;
+	font-size: 12px;
+	letter-spacing: 1px;
+	opacity: 0.8;
+}
+
 .bit-grid {
 	display: grid;
 	grid-template-columns: repeat(8, 20px);
@@ -248,6 +295,14 @@ figure {
 	border: 1px solid var(--background);
 }
 
+.bit-grid.small {
+	grid-template-columns: repeat(8, 13px);
+}
+
+.bit-grid.small .bit {
+	height: 13px;
+}
+
 .bit-grid .bit.on {
 	background: var(--foreground);
 }
@@ -256,6 +311,20 @@ figure {
 	background: color-mix(in srgb, var(--foreground) 12%, transparent);
 }
 
+/* Changed bits keep their on/off brightness and get a red tint + border,
+   so light-changed and dark-changed cells stay distinguishable */
+.bit-grid .bit.diff {
+	border-color: hsl(4 70% 55%);
+}
+
+.bit-grid .bit.on.diff {
+	background: color-mix(in srgb, var(--foreground) 55%, hsl(4 70% 50%));
+}
+
+.bit-grid .bit.off.diff {
+	background: hsl(4 60% 50% / 0.25);
+}
+
 .heat-grid {
 	display: grid;
 	grid-template-columns: repeat(8, 32px);
@@ -266,7 +335,15 @@ figure {
 	box-sizing: border-box;
 }
 
-.pixelated, #image-resize img, #phash-resize img {
+.heat-grid.mini {
+	grid-template-columns: repeat(8, 22px);
+}
+
+.heat-grid.mini > div {
+	height: 22px;
+}
+
+.pixelated, #phash-resize img, .input-thumb img {
 	image-rendering: pixelated;
 	image-rendering: -moz-crisp-edges;
 }
@@ -276,71 +353,520 @@ figure {
 	height: 128px;
 }
 
-.sample {
-	height: 48px;
-	margin: 0 4px;
-	cursor: pointer;
-	vertical-align: middle;
+.input-thumb img {
+	width: 176px;
+	height: 176px;
 }
 
-/* Mutation compare */
+/* DCT extras */
 
-.compare {
+.grid-overlay {
 	position: relative;
-	max-width: 640px;
+	width: 256px;
 }
 
-.compare img {
-	display: block;
-	width: 100%;
+.grid-overlay.mini {
+	width: 176px;
 }
 
-.compare #compare-b {
+#dct-basis-wrap {
+	position: relative;
+	display: inline-block;
+}
+
+.basis-hover {
+	display: none;
+	position: absolute;
+	outline: 2px solid var(--color1);
+	outline-offset: -1px;
+	pointer-events: none;
+}
+
+#dct-basis {
+	cursor: crosshair;
+}
+
+#dct-formula .formula-result {
+	display: none;
+}
+
+#dct-formula.live .formula-result {
+	display: inline;
+}
+
+#dct-formula.live .var-u,
+#dct-formula.live .var-v,
+#dct-formula.live .formula-value {
+	color: var(--color1);
+}
+
+#image-dct .hover-cell {
+	outline: 2px solid var(--color1);
+	outline-offset: -2px;
+}
+
+#dct-zigzag {
 	position: absolute;
 	top: 0;
 	left: 0;
-	clip-path: inset(0 0 0 50%);
+	pointer-events: none;
 }
 
-#compare-slider, #mutator-amount {
-	width: 320px;
+#dct-zigzag polyline {
+	stroke: var(--color1);
+	stroke-width: 3;
+	stroke-linejoin: round;
+}
+
+#dct-zigzag polyline.play {
+	transition: stroke-dashoffset 2.5s ease-in-out;
+}
+
+#dct-basis {
+	background: var(--background);
+}
+
+/* Reconstruction demos */
+
+.recon-row {
+	display: flex;
+	flex-flow: row wrap;
+	gap: 24px;
+	align-items: flex-start;
+	margin: 10px 0;
+}
+
+.recon-row figure {
+	margin: 0;
+	text-align: center;
+}
+
+.recon-row figcaption {
+	font-size: 12px;
+	opacity: 0.8;
+	margin-top: 4px;
+}
+
+.recon-row canvas {
+	width: 176px;
+	height: 176px;
+	image-rendering: pixelated;
+	image-rendering: -moz-crisp-edges;
+}
+
+#recon-k {
+	width: 256px;
+	vertical-align: middle;
+}
+
+.recon-note {
+	font-size: 12px;
+	opacity: 0.7;
+	margin-left: 8px;
+}
+
+#recon-zigzag {
+	position: absolute;
+	top: 0;
+	left: 0;
+	pointer-events: none;
+}
+
+#recon-zigzag polyline {
+	stroke-width: 2.5;
+	stroke-linejoin: round;
+	transition: stroke-dashoffset 0.12s linear;
+}
+
+#recon-zigzag .zig-head {
+	stroke: var(--color1);
+}
+
+#recon-zigzag .zig-tail {
+	stroke: color-mix(in srgb, var(--foreground) 45%, transparent);
+}
+
+/* Tables */
+
+table.facts {
+	border-collapse: collapse;
+	margin: 10px 0;
+}
+
+table.facts th, table.facts td {
+	text-align: left;
+	padding: 4px 14px 4px 0;
+	border-bottom: 1px solid color-mix(in srgb, var(--foreground) 20%, transparent);
+	font-size: 13px;
+}
+
+table.facts th {
+	color: var(--color1);
+	font-weight: normal;
 }
 
 .distances td {
 	text-align: center;
 	font-size: 18px;
 	min-width: 64px;
+}
+
+.distances td.match {
+	color: hsl(135 50% 45%);
+}
+
+.distances td.nomatch {
+	color: hsl(4 70% 55%);
+}
+
+/* LSH demo */
+
+.lsh-thumb img {
+	height: 110px;
+	margin-bottom: 6px;
+	display: block;
+}
+
+.lsh-distance {
+	font-size: 16px;
+}
+
+.lsh-distance .match {
+	font-size: 22px;
+	color: hsl(135 50% 45%);
+}
+
+.lsh-distance .nomatch {
+	font-size: 22px;
+	color: hsl(4 70% 55%);
+}
+
+/* Flip demo */
+
+#flip-stage img {
+	max-width: 256px;
+	max-height: 256px;
+	transition: transform 0.7s ease-in-out;
+}
+
+#flip-stage img.flipped {
+	transform: scaleX(-1);
+}
+
+#ghost-canvas {
+	transition: transform 0.4s ease-in-out;
+}
+
+#ghost-canvas.flipped {
+	transform: scaleX(-1);
+}
+
+.flip-table th {
 	color: var(--color1);
+	font-weight: normal;
+	text-align: left;
+	padding-right: 12px;
+}
+
+.flip-table td {
+	vertical-align: top;
+	padding: 6px 12px 6px 0;
+}
+
+.flip-distance {
+	font-size: 24px;
+}
+
+.flip-distance.match {
+	color: hsl(135 50% 45%);
+}
+
+.flip-distance.nomatch {
+	color: hsl(4 70% 55%);
 }
 
 /* Ranking */
 
-.thumb-grid {
+.pool {
+	display: grid;
+	grid-template-columns: repeat(auto-fill, minmax(72px, 1fr));
+	gap: 4px;
+	max-height: 330px;
+	overflow-y: auto;
+	margin: 10px 0;
+}
+
+.pool img {
+	width: 100%;
+	aspect-ratio: 1;
+	object-fit: cover;
+	cursor: pointer;
+	display: block;
+}
+
+.pool img:hover {
+	outline: 2px solid var(--color1);
+}
+
+.pool img.query {
+	outline: 3px solid var(--color1);
+}
+
+.ranking-results h3 {
+	margin-top: 14px;
+}
+
+.thumb-row {
 	display: flex;
 	flex-flow: row wrap;
 	gap: 8px;
 }
 
-.thumb-grid figure {
-	width: 96px;
+.thumb-row figure {
+	width: 84px;
 	margin: 0;
-	cursor: pointer;
 	text-align: center;
 }
 
-.thumb-grid figure.query {
-	outline: 2px solid var(--color1);
+.thumb-row figure.hit {
+	outline: 2px solid hsl(135 50% 45%);
 }
 
-.thumb-grid img {
-	width: 96px;
-	height: 96px;
+.thumb-row img {
+	width: 84px;
+	height: 84px;
 	object-fit: cover;
+	display: block;
 }
 
-.thumb-grid figcaption {
+.thumb-row figcaption {
 	font-size: 11px;
 	overflow: hidden;
 	text-overflow: ellipsis;
 	white-space: nowrap;
 }
+
+/* Precision-recall charts */
+
+.pr-controls {
+	margin: 10px 0;
+}
+
+#pr-legend {
+	display: flex;
+	flex-flow: row wrap;
+	gap: 4px 14px;
+	margin-bottom: 8px;
+}
+
+.pr-group {
+	display: flex;
+	flex-flow: row wrap;
+	gap: 4px;
+}
+
+.pr-chip {
+	display: inline-flex;
+	align-items: center;
+	gap: 5px;
+	border: 1px solid color-mix(in srgb, var(--foreground) 20%, transparent);
+	background: none;
+	color: var(--foreground);
+	font-size: 11px;
+	padding: 2px 7px;
+	cursor: pointer;
+	transition: opacity 0.15s;
+}
+
+.pr-chip.dim {
+	opacity: 0.3;
+}
+
+.pr-chip.pinned {
+	border-color: var(--color1);
+}
+
+.pr-slider {
+	display: inline-block;
+	margin: 4px 0;
+	color: var(--color4);
+}
+
+.pr-slider input {
+	width: 220px;
+	vertical-align: middle;
+	margin-left: 8px;
+}
+
+.pr-charts {
+	display: flex;
+	flex-flow: row wrap;
+	gap: 12px;
+}
+
+.pr-charts figure {
+	margin: 0;
+	flex: 1 1 300px;
+	max-width: 340px;
+}
+
+.pr-charts figcaption {
+	color: var(--color1);
+	font-size: 14px;
+	margin-bottom: 2px;
+}
+
+.pr-svg {
+	width: 100%;
+	height: auto;
+	background: color-mix(in srgb, var(--foreground) 3%, transparent);
+}
+
+.pr-grid {
+	stroke: color-mix(in srgb, var(--foreground) 12%, transparent);
+	stroke-width: 1;
+}
+
+.pr-tick {
+	fill: color-mix(in srgb, var(--foreground) 55%, transparent);
+	font-size: 9px;
+}
+
+.pr-axis {
+	fill: color-mix(in srgb, var(--foreground) 70%, transparent);
+	font-size: 11px;
+}
+
+.pr-line {
+	fill: none;
+	stroke-width: 1.8;
+	transition: opacity 0.15s;
+}
+
+.pr-dot {
+	stroke: var(--background);
+	stroke-width: 2;
+	transition: opacity 0.15s;
+}
+
+.pr-line.faded, .pr-dot.faded {
+	opacity: 0.12;
+}
+
+#pr-tooltip {
+	position: fixed;
+	z-index: 10;
+	background: color-mix(in srgb, var(--background) 92%, transparent);
+	border: 1px solid var(--color1);
+	color: var(--foreground);
+	font-size: 12px;
+	padding: 4px 8px;
+	pointer-events: none;
+	white-space: nowrap;
+}
+
+/* Retrieval demos */
+
+.dist-figure, .toy-figure {
+	margin: 10px 0;
+}
+
+.dist-figure figcaption, .toy-figure figcaption {
+	font-size: 12px;
+	opacity: 0.8;
+	margin-top: 4px;
+}
+
+.dist-svg {
+	max-width: 720px;
+}
+
+.dist-bar {
+	fill: color-mix(in srgb, var(--foreground) 45%, transparent);
+}
+
+.dist-copy {
+	fill: hsl(135 50% 45%);
+	stroke: var(--background);
+	stroke-width: 1;
+}
+
+.dist-region {
+	fill: color-mix(in srgb, var(--color1) 10%, transparent);
+}
+
+.dist-thresh {
+	stroke: var(--color1);
+	stroke-width: 2;
+}
+
+.dist-readout, .toy-stats {
+	margin-left: 12px;
+	color: var(--color4);
+}
+
+.toy-node {
+	fill: color-mix(in srgb, var(--foreground) 30%, transparent);
+	transition: opacity 0.3s, fill 0.3s;
+}
+
+.toy-node.seen {
+	fill: color-mix(in srgb, var(--foreground) 10%, var(--background));
+	stroke: var(--foreground);
+	stroke-width: 2;
+}
+
+.toy-node.match {
+	fill: hsl(135 50% 45%);
+}
+
+.toy-node.pruned {
+	opacity: 0.15;
+}
+
+.toy-edge {
+	stroke: color-mix(in srgb, var(--foreground) 45%, transparent);
+	stroke-width: 1.5;
+	transition: opacity 0.3s;
+}
+
+.toy-edge-label {
+	fill: color-mix(in srgb, var(--foreground) 70%, transparent);
+	font-size: 9px;
+	transition: opacity 0.3s;
+}
+
+.toy-edge.pruned, .toy-edge-label.pruned {
+	opacity: 0.15;
+}
+
+.big-stats .stat-note {
+	font-size: 12px;
+	opacity: 0.7;
+	margin: 0;
+	max-width: 220px;
+}
+
+.big-curve-svg {
+	width: 250px;
+}
+
+.big-curve-line {
+	stroke: var(--color2);
+}
+
+.big-curve-dot {
+	fill: var(--color1);
+}
+
+/* Lists */
+
+ul.loose-ends li, ol.references li {
+	margin-bottom: 6px;
+	line-height: 1.5;
+}
+
+ol.references li:target {
+	outline: 1px dashed var(--color1);
+	outline-offset: 4px;
+}

web/worker.js

@@ -1,15 +1,32 @@
-import init, { hash_all, resize_preview, dct_coefficients, dct_masks, phash_lowfreq, mutate } from './pkg/image_similarity.js';
+import init, { hash_all, resize_preview, dct_coefficients, dct_masks, phash_lowfreq, mutate, suite_hashes, suite_labels, BkIndex } from './pkg/image_similarity.js';
 
 const ready = init();
 
+// The retrieval demos keep their trees in wasm memory between calls:
+// a small one that gets visualized and the full experiment index.
+let toyIndex = null;
+let bigIndex = null;
+
+// Average runtime of fn, sampled until enough wall time passed to
+// outrun the clamped performance.now resolution
+function time(fn) {
+	const start = performance.now();
+	let reps = 0;
+	do {
+		fn();
+		reps++;
+	} while (performance.now() - start < 20);
+	return (performance.now() - start) / reps;
+}
+
 onmessage = async (e) => {
 	const { id, op, args } = e.data;
 	await ready;
 	try {
-		const bytes = new Uint8Array(args.buffer);
 		let result;
 		switch (op) {
-			case 'pipeline':
+			case 'pipeline': {
+				const bytes = new Uint8Array(args.buffer);
 				result = {
 					resize8: resize_preview(bytes, 8),
 					resize32: resize_preview(bytes, 32),
@@ -19,14 +36,63 @@ onmessage = async (e) => {
 					hashes: hash_all(bytes),
 				};
 				break;
+			}
 			case 'hashes':
-				result = { hashes: hash_all(bytes) };
+				result = { hashes: hash_all(new Uint8Array(args.buffer)) };
 				break;
 			case 'mutate': {
-				const png = mutate(bytes, args.kind, args.amount);
+				const png = mutate(new Uint8Array(args.buffer), args.kind, args.amount);
 				result = { png, hashes: hash_all(png) };
 				break;
 			}
+			case 'suite':
+				result = {
+					hashes: suite_hashes(new Uint8Array(args.buffer)),
+					labels: suite_labels().split('\n'),
+				};
+				break;
+			case 'toy_build': {
+				toyIndex?.free();
+				toyIndex = new BkIndex();
+				for (const key of args.keys) toyIndex.insert(key);
+				result = { structure: toyIndex.structure(), size: toyIndex.len() };
+				break;
+			}
+			case 'toy_query':
+				result = {
+					found: toyIndex.find(args.key, args.radius),
+					visited: toyIndex.trace(args.key, args.radius),
+				};
+				break;
+			case 'big_build': {
+				bigIndex?.free();
+				bigIndex = new BkIndex();
+				bigIndex.insert_bytes(new Uint8Array(args.buffer));
+				result = { size: bigIndex.len() };
+				break;
+			}
+			case 'big_query': {
+				const found = bigIndex.find(args.key, args.radius);
+				result = {
+					found,
+					compared: bigIndex.find_compared(args.key, args.radius),
+					treeMs: time(() => bigIndex.find(args.key, args.radius)),
+					scanMs: time(() => bigIndex.scan(args.key, args.radius)),
+				};
+				break;
+			}
+			// compared counts and found counts for every radius up to max,
+			// so the radius slider works without a round trip per step
+			case 'big_profile': {
+				const compared = [];
+				const found = [];
+				for (let t = 0; t <= args.max; t++) {
+					compared.push(bigIndex.find_compared(args.key, t));
+					found.push(bigIndex.find(args.key, t).length / 2);
+				}
+				result = { compared, found };
+				break;
+			}
 			default:
 				throw new Error(`unknown op ${op}`);
 		}