Oregami
Repositories/oxedyne/fe2o3

oxedyne/fe2o3/fe2o3_hash/src/phash.rs

14.6 KiB, 30 runs

created by r1870400018:17752, which is this file's identity for as long as the history lasts, whatever it is later renamed to

download · who wrote it · its history

1//! Perceptual image hashing for near-duplicate detection.
2//!
3//! Two 64 bit hashes are offered. [`PerceptualHash::dhash`] is a difference hash: the image is
4//! reduced to a nine by eight grid and each pixel compared with its right hand neighbour. It
5//! costs almost nothing and is the right first pass over a large collection.
6//! [`PerceptualHash::phash`] takes a discrete cosine transform of a thirty-two square reduction
7//! and thresholds the low frequency coefficients about their median. It survives compression,
8//! rescaling and mild tonal shifts better than the difference hash, and is the right instrument
9//! for confirming a candidate the cheap pass has thrown up.
10//!
11//! The recommended use is therefore two stages: hash everything with [`PerceptualHash::dhash`],
12//! shortlist by small [`PerceptualHash::distance`], then confirm each shortlisted pair with
13//! [`PerceptualHash::phash`].
14//!
15//! # Choosing a threshold
16//!
17//! Thresholds belong to the collection, so measure rather than assume. As a starting point,
18//! over a spread of photographs put through a half-size reduction, a heavy re-encode, a ten per
19//! cent brightening and a lossless to lossy conversion, the distance between a photograph and
20//! its own variants never exceeded four for the difference hash or two for the cosine transform
21//! hash, while the closest unrelated pair was nineteen and twenty-four respectively. A first
22//! pass at ten and a confirmation at eight therefore sit in a wide empty gap.
23//!
24//! Cropping is the case that defeats both. A ten per cent centre crop moved the distance to a
25//! median of eleven and a worst case of twenty-eight, which overlaps the unrelated population.
26//! Neither hash is a crop detector, and a collection full of crops needs a different instrument.
27//!
28//! Following this library's practice for primitives, the caller owns the decode. These functions
29//! accept a greyscale luma grid, never a file path and never compressed bytes, so the choice of
30//! image decoder stays with the application. [`luma_from_rgb`] and [`luma_from_rgba`] convert an
31//! interleaved buffer for callers whose decoder hands back colour.
32//!
33//! # Example
34//! ```
35//! use oxedyne_fe2o3_core::prelude::*;
36//! use oxedyne_fe2o3_hash::phash::{LumaGrid, PerceptualHash};
37//!
38//! fn near_duplicate(px: &[u8], w: usize, h: usize) -> Outcome<(u64, u64)> {
39//! let grid = res!(LumaGrid::new(px, w, h));
40//! let d = res!(PerceptualHash::dhash(&grid));
41//! let p = res!(PerceptualHash::phash(&grid));
42//! assert_eq!(res!(d.distance(&d)), 0);
43//! assert!(d.distance(&p).is_err()); // Unlike kinds cannot be compared.
44//! Ok((d.bits(), p.bits()))
45//! }
46//!
47//! let px = vec![0u8; 64 * 64];
48//! assert!(near_duplicate(&px, 64, 64).is_ok());
49//! ```
50//!
51//! [Written with AI entirely](https://need2know.ai/entirely-ai/code)\
52//! Anthropic Claude
53
54use oxedyne_fe2o3_core::prelude::*;
55
56use std::{
57 fmt,
58 f64::consts::PI,
59};
60
61
62// The reduction each hash works on: nine by eight for the difference hash, a thirty two square
63// for the transform hash, of which the top left eight square is thresholded.
64pub const DHASH_W: usize = 9;
65pub const DHASH_H: usize = 8;
66pub const PHASH_N: usize = 32;
67pub const PHASH_K: usize = 8;
68
69const MAX_DIM: usize = 1 << 20; // guards against an absurd resample
70
71/// A borrowed greyscale grid: eight bits per pixel, row major, no row padding.
72#[derive(Clone, Copy, Debug)]
73pub struct LumaGrid<'a> {
74 dat: &'a [u8],
75 w: usize,
76 h: usize,
77}
78
79impl<'a> LumaGrid<'a> {
80
81 /// Wraps a luma buffer, checking that its length matches the stated dimensions.
82 pub fn new(dat: &'a [u8], w: usize, h: usize) -> Outcome<Self> {
83 if w == 0 || h == 0 {
84 return Err(err!(
85 "A luma grid must have a positive width and height, found {} by {}.", w, h;
86 Input, Invalid, TooSmall));
87 }
88 if w > MAX_DIM || h > MAX_DIM {
89 return Err(err!(
90 "A luma grid dimension of {} by {} exceeds the {} pixel limit.", w, h, MAX_DIM;
91 Input, Invalid, TooBig));
92 }
93 let need = w * h;
94 if dat.len() < need {
95 return Err(err!(
96 "A {} by {} luma grid needs {} bytes, {} were supplied.", w, h, need, dat.len();
97 Input, Invalid, TooSmall, Size));
98 }
99 Ok(Self { dat, w, h })
100 }
101
102 pub fn width(&self) -> usize {
103 self.w
104 }
105
106 pub fn height(&self) -> usize {
107 self.h
108 }
109
110 pub fn data(&self) -> &'a [u8] {
111 self.dat
112 }
113
114 /// Reduces the grid to `tw` by `th` samples by averaging over the source area of each.
115 ///
116 /// Every source pixel contributes in proportion to its overlap with the target cell, so the
117 /// result is stable against small changes in the source dimensions. Enlargement is
118 /// permitted and simply replicates, though hashing an image smaller than the reduction
119 /// carries little information.
120 pub fn resample(&self, tw: usize, th: usize) -> Outcome<Vec<f64>> {
121 if tw == 0 || th == 0 {
122 return Err(err!(
123 "A resample target must have a positive width and height, found {} by {}.",
124 tw, th;
125 Input, Invalid, TooSmall));
126 }
127 let sx = self.w as f64 / tw as f64; // Source pixels per target cell, horizontally
128 let sy = self.h as f64 / th as f64; // Source pixels per target cell, vertically
129 let mut out = vec![0.0f64; tw * th];
130 for ty in 0..th {
131 let y0 = ty as f64 * sy;
132 let y1 = (ty + 1) as f64 * sy;
133 let iy0 = y0.floor() as usize;
134 let iy1 = (y1.ceil() as usize).min(self.h);
135 for tx in 0..tw {
136 let x0 = tx as f64 * sx;
137 let x1 = (tx + 1) as f64 * sx;
138 let ix0 = x0.floor() as usize;
139 let ix1 = (x1.ceil() as usize).min(self.w);
140 let mut acc = 0.0f64;
141 let mut wt = 0.0f64;
142 for y in iy0..iy1 {
143 let hy = (y1.min((y + 1) as f64) - y0.max(y as f64)).max(0.0);
144 if hy == 0.0 {
145 continue;
146 }
147 let row = y * self.w;
148 for x in ix0..ix1 {
149 let hx = (x1.min((x + 1) as f64) - x0.max(x as f64)).max(0.0);
150 if hx == 0.0 {
151 continue;
152 }
153 let a = hx * hy;
154 acc += a * self.dat[row + x] as f64;
155 wt += a;
156 }
157 }
158 out[ty * tw + tx] = if wt > 0.0 {
159 acc / wt
160 } else {
161 // The cell fell between sample centres, which only happens when the target
162 // is larger than the source; take the nearest source pixel.
163 let y = iy0.min(self.h - 1);
164 let x = ix0.min(self.w - 1);
165 self.dat[y * self.w + x] as f64
166 };
167 }
168 }
169 Ok(out)
170 }
171}
172
173/// A 64 bit perceptual hash, tagged by the algorithm that produced it.
174///
175/// Two hashes are comparable only when they were produced the same way, which the tag enforces.
176#[derive(Clone, Copy, Debug, Eq, Hash, Ord, PartialEq, PartialOrd)]
177pub enum PerceptualHash {
178 DHash(u64), // cheap, the right first pass
179 PHash(u64), // slower, the right confirmation
180}
181
182impl fmt::Display for PerceptualHash {
183 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
184 match self {
185 Self::DHash(b) => write!(f, "d:{:016x}", b),
186 Self::PHash(b) => write!(f, "p:{:016x}", b),
187 }
188 }
189}
190
191impl PerceptualHash {
192
193 /// Computes the difference hash of a luma grid.
194 ///
195 /// The grid is reduced to nine by eight samples and each sample compared with the one to its
196 /// right, giving sixty-four bits written most significant first, row by row.
197 pub fn dhash(grid: &LumaGrid) -> Outcome<Self> {
198 let px = res!(grid.resample(DHASH_W, DHASH_H));
199 let mut bits = 0u64;
200 for y in 0..DHASH_H {
201 for x in 0..DHASH_H {
202 let left = px[y * DHASH_W + x];
203 let right = px[y * DHASH_W + x + 1];
204 bits = (bits << 1) | u64::from(left < right);
205 }
206 }
207 Ok(Self::DHash(bits))
208 }
209
210 /// Computes the discrete cosine transform hash of a luma grid.
211 ///
212 /// The grid is reduced to thirty-two square, transformed, and the eight by eight low
213 /// frequency block thresholded about the median of its alternating-current terms. The
214 /// constant term is excluded from that median, since overall brightness would otherwise
215 /// shift every bit, but it still contributes its own bit.
216 pub fn phash(grid: &LumaGrid) -> Outcome<Self> {
217 let px = res!(grid.resample(PHASH_N, PHASH_N));
218 let co = res!(dct2d(&px, PHASH_N));
219
220 // Gather the low frequency block, keeping the alternating-current terms apart so the
221 // median is not dragged by the constant term.
222 let mut block = [0.0f64; PHASH_K * PHASH_K];
223 let mut ac = Vec::with_capacity(PHASH_K * PHASH_K - 1);
224 for v in 0..PHASH_K {
225 for u in 0..PHASH_K {
226 let c = co[v * PHASH_N + u];
227 block[v * PHASH_K + u] = c;
228 if !(u == 0 && v == 0) {
229 ac.push(c);
230 }
231 }
232 }
233 let med = median(&mut ac);
234
235 let mut bits = 0u64;
236 for c in block.iter() {
237 bits = (bits << 1) | u64::from(*c > med);
238 }
239 Ok(Self::PHash(bits))
240 }
241
242 pub fn bits(&self) -> u64 {
243 match self {
244 Self::DHash(b) => *b,
245 Self::PHash(b) => *b,
246 }
247 }
248
249 /// Were both hashes produced by the same algorithm?
250 pub fn same_kind(&self, other: &Self) -> bool {
251 matches!(
252 (self, other),
253 (Self::DHash(..), Self::DHash(..)) | (Self::PHash(..), Self::PHash(..))
254 )
255 }
256
257 /// The result runs from zero, for identical hashes, to sixty-four. A distance near thirty
258 /// two means no relationship at all, since that is what a random pair gives.
259 pub fn distance(&self, other: &Self) -> Outcome<u32> {
260 if !self.same_kind(other) {
261 return Err(err!(
262 "A {} hash cannot be compared with a {} hash; the bit positions mean different \
263 things.", self.label(), other.label();
264 Input, Invalid, Mismatch));
265 }
266 Ok(hamming(self.bits(), other.bits()))
267 }
268
269 /// Returns the algorithm name, for messages.
270 pub fn label(&self) -> &'static str {
271 match self {
272 Self::DHash(..) => "difference",
273 Self::PHash(..) => "cosine transform",
274 }
275 }
276}
277
278/// Returns the number of differing bits between two hashes.
279pub fn hamming(a: u64, b: u64) -> u32 {
280 (a ^ b).count_ones()
281}
282
283/// Converts an interleaved buffer to luma, taking the first three channels of each pixel.
284///
285/// The weights are those of Rec. 601, the convention every common image tool applies when it
286/// desaturates. `stride` is the number of bytes per pixel and must be at least three.
287pub fn luma_from_interleaved(
288 dat: &[u8],
289 w: usize,
290 h: usize,
291 stride: usize,
292)
293 -> Outcome<Vec<u8>>
294{
295 if stride < 3 {
296 return Err(err!(
297 "An interleaved buffer needs at least three bytes per pixel, {} were declared.",
298 stride;
299 Input, Invalid, TooSmall));
300 }
301 if w == 0 || h == 0 {
302 return Err(err!(
303 "An interleaved buffer must have a positive width and height, found {} by {}.", w, h;
304 Input, Invalid, TooSmall));
305 }
306 let need = w * h * stride;
307 if dat.len() < need {
308 return Err(err!(
309 "A {} by {} buffer at {} bytes per pixel needs {} bytes, {} were supplied.",
310 w, h, stride, need, dat.len();
311 Input, Invalid, TooSmall, Size));
312 }
313 let mut out = Vec::with_capacity(w * h);
314 for i in 0..(w * h) {
315 let p = i * stride;
316 let y = 0.299 * dat[p] as f64
317 + 0.587 * dat[p + 1] as f64
318 + 0.114 * dat[p + 2] as f64;
319 out.push(y.round().clamp(0.0, 255.0) as u8);
320 }
321 Ok(out)
322}
323
324/// Converts a packed red, green, blue buffer to luma.
325pub fn luma_from_rgb(dat: &[u8], w: usize, h: usize) -> Outcome<Vec<u8>> {
326 luma_from_interleaved(dat, w, h, 3)
327}
328
329/// Converts a packed red, green, blue, alpha buffer to luma, ignoring the alpha channel.
330pub fn luma_from_rgba(dat: &[u8], w: usize, h: usize) -> Outcome<Vec<u8>> {
331 luma_from_interleaved(dat, w, h, 4)
332}
333
334/// Computes the orthonormal two dimensional type two discrete cosine transform of a square grid.
335///
336/// The transform is separable, so it is applied along the rows and then along the columns, which
337/// costs `n` cubed multiplications rather than `n` to the fourth.
338fn dct2d(px: &[f64], n: usize) -> Outcome<Vec<f64>> {
339 if n == 0 {
340 return Err(err!(
341 "A discrete cosine transform needs a positive side length, {} was given.", n;
342 Input, Invalid, TooSmall));
343 }
344 if px.len() < n * n {
345 return Err(err!(
346 "A {} square discrete cosine transform needs {} samples, {} were supplied.",
347 n, n * n, px.len();
348 Input, Invalid, TooSmall, Size));
349 }
350 // Basis table: cos((2i + 1) k pi / 2n), indexed as [i * n + k].
351 let mut cos = vec![0.0f64; n * n];
352 for i in 0..n {
353 for k in 0..n {
354 cos[i * n + k] = (((2 * i + 1) as f64) * (k as f64) * PI / (2.0 * n as f64)).cos();
355 }
356 }
357 let s0 = (1.0 / n as f64).sqrt(); // Scale of the constant term
358 let sk = (2.0 / n as f64).sqrt(); // Scale of every other term
359
360 // Rows first.
361 let mut tmp = vec![0.0f64; n * n];
362 for y in 0..n {
363 for u in 0..n {
364 let mut acc = 0.0f64;
365 for x in 0..n {
366 acc += px[y * n + x] * cos[x * n + u];
367 }
368 tmp[y * n + u] = acc * if u == 0 { s0 } else { sk };
369 }
370 }
371 // Then columns.
372 let mut out = vec![0.0f64; n * n];
373 for u in 0..n {
374 for v in 0..n {
375 let mut acc = 0.0f64;
376 for y in 0..n {
377 acc += tmp[y * n + u] * cos[y * n + v];
378 }
379 out[v * n + u] = acc * if v == 0 { s0 } else { sk };
380 }
381 }
382 Ok(out)
383}
384
385/// Returns the median of a slice, sorting it in place; an empty slice gives zero.
386fn median(v: &mut [f64]) -> f64 {
387 if v.is_empty() {
388 return 0.0;
389 }
390 v.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
391 let n = v.len();
392 if n % 2 == 1 {
393 v[n / 2]
394 } else {
395 (v[n / 2 - 1] + v[n / 2]) / 2.0
396 }
397}