From 33a88f885b728a94e1c45e87ffbbe4d4022f0f46 Mon Sep 17 00:00:00 2001 From: cospectrum Date: Sun, 19 Jan 2025 22:25:55 +0300 Subject: [PATCH] add `pHash` (#709) * add phash * update msrv --- .github/workflows/check.yml | 2 +- Cargo.toml | 14 +- examples/template_matching.rs | 2 +- src/image_hash/bits.rs | 106 ++++++++++++ src/image_hash/mod.rs | 11 ++ src/image_hash/phash.rs | 146 ++++++++++++++++ src/image_hash/signals.rs | 313 ++++++++++++++++++++++++++++++++++ src/lib.rs | 1 + src/utils.rs | 18 +- 9 files changed, 601 insertions(+), 12 deletions(-) create mode 100644 src/image_hash/bits.rs create mode 100644 src/image_hash/mod.rs create mode 100644 src/image_hash/phash.rs create mode 100644 src/image_hash/signals.rs diff --git a/.github/workflows/check.yml b/.github/workflows/check.yml index 4d85cc2b..c8523e17 100644 --- a/.github/workflows/check.yml +++ b/.github/workflows/check.yml @@ -98,7 +98,7 @@ jobs: # https://docs.github.com/en/actions/learn-github-actions/contexts#context-availability strategy: matrix: - msrv: ["1.79.0"] # Don't forget to update the `rust-version` in Cargo.toml as well + msrv: ["1.80.1"] # Don't forget to update the `rust-version` in Cargo.toml as well name: ubuntu / ${{ matrix.msrv }} steps: - uses: actions/checkout@v4 diff --git a/Cargo.toml b/Cargo.toml index f0f90c56..3520d054 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,7 +3,7 @@ name = "imageproc" version = "0.26.0" authors = ["theotherphil"] # note: when changed, also update `msrv` in `.github/workflows/check.yml` -rust-version = "1.79.0" +rust-version = "1.80.1" edition = "2021" license = "MIT" description = "Image processing operations" @@ -22,20 +22,21 @@ ab_glyph = { version = "0.2.23", default-features = false, features = ["std"] } approx = { version = "0.5", default-features = false } image = { version = "0.25.0", default-features = false } itertools = { version = "0.13.0", default-features = false, features = [ - "use_std", + "use_std", ] } nalgebra = { version = "0.32", default-features = false, features = ["std"] } num = { version = "0.4.1", default-features = false } rand = { version = "0.8.5", default-features = false, features = [ - "std", - "std_rng", + "std", + "std_rng", ] } rand_distr = { version = "0.4.3", default-features = false } rayon = { version = "1.8.0", optional = true, default-features = false } sdl2 = { version = "0.36", optional = true, default-features = false, features = [ - "bundled", + "bundled", ] } katexit = { version = "0.1.4", optional = true, default-features = false } +rustdct = "0.7.1" [target.'cfg(target_arch = "wasm32")'.dependencies] getrandom = { version = "0.2", default-features = false, features = ["js"] } @@ -54,6 +55,9 @@ features = ["property-testing", "katexit"] opt-level = 3 debug = true +[profile.dev] +opt-level = 1 + [profile.bench] opt-level = 3 debug = true diff --git a/examples/template_matching.rs b/examples/template_matching.rs index 320743a4..f0a25474 100644 --- a/examples/template_matching.rs +++ b/examples/template_matching.rs @@ -45,7 +45,7 @@ If the optional boolean argument parallel is given, match_template will be calle let template_y = args[4].parse().unwrap(); let template_w = args[5].parse().unwrap(); let template_h = args[6].parse().unwrap(); - let parallel = args.get(7).map_or(false, |s| s.parse().unwrap()); + let parallel = args.get(7).is_some_and(|s| s.parse().unwrap()); TemplateMatchingArgs { input_path, diff --git a/src/image_hash/bits.rs b/src/image_hash/bits.rs new file mode 100644 index 00000000..bee8bf4f --- /dev/null +++ b/src/image_hash/bits.rs @@ -0,0 +1,106 @@ +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +pub(super) struct Bits64(u64); + +impl Bits64 { + const N: usize = 64; + + pub fn new(v: impl IntoIterator) -> Self { + let mut bits = Self::zeros(); + let mut n = 0; + for bit in v { + if bit { + bits.set_bit_at(n); + } else { + bits.unset_bit_at(n); + }; + n += 1; + } + assert_eq!(n, Self::N); + bits + } + pub fn zeros() -> Self { + Self(0) + } + #[allow(dead_code)] + pub fn to_bitarray(self) -> [bool; Self::N] { + let mut bits = [false; Self::N]; + for (n, bit) in bits.iter_mut().enumerate() { + *bit = self.bit_at(n) + } + bits + } + pub fn hamming_distance(self, other: Bits64) -> u32 { + self.xor(other).0.count_ones() + } + fn xor(self, other: Self) -> Self { + Self(self.0 ^ other.0) + } + fn bit_at(self, n: usize) -> bool { + assert!(n < Self::N); + let bit = self.0 & (1 << n); + bit != 0 + } + fn set_bit_at(&mut self, n: usize) { + assert!(n < Self::N); + self.0 |= 1 << n; + } + fn unset_bit_at(&mut self, n: usize) { + assert!(n < Self::N); + self.0 &= !(1 << n); + } +} + +#[cfg(test)] +mod tests { + use std::collections::HashMap; + + use super::*; + + #[test] + fn test_bits64_ops() { + let mut bits = Bits64::zeros(); + bits.set_bit_at(0); + assert_eq!(bits, Bits64(1)); + bits.set_bit_at(1); + assert_eq!(bits, Bits64(1 + 2)); + bits.unset_bit_at(0); + assert_eq!(bits, Bits64(2)); + bits.unset_bit_at(1); + assert_eq!(bits, Bits64::zeros()); + + bits.set_bit_at(2); + assert_eq!(bits, Bits64(4)); + } + #[test] + fn test_bitarray() { + let mut v = [false; Bits64::N]; + v[3] = true; + v[6] = true; + let bits = Bits64::new(v); + assert_eq!(bits.to_bitarray(), v); + } + #[test] + fn test_bits64_new() { + const N: usize = 64; + + let mut v = [false; N]; + v[0] = true; + assert_eq!(Bits64::new(v), Bits64(1)); + v[1] = true; + assert_eq!(Bits64::new(v), Bits64(1 + 2)); + } + #[test] + #[should_panic] + fn test_bits64_new_fail() { + const N: usize = 64; + let it = (1..N).map(|x| x % 2 == 0); + let _bits = Bits64::new(it); + } + + #[test] + fn test_hash() { + let one = Bits64(1); + let map = HashMap::from([(one, "1")]); + assert_eq!(map[&one], "1"); + } +} diff --git a/src/image_hash/mod.rs b/src/image_hash/mod.rs new file mode 100644 index 00000000..ccc18122 --- /dev/null +++ b/src/image_hash/mod.rs @@ -0,0 +1,11 @@ +//! [Perceptual hashing] algorithms for images. +//! +//! [Perceptual hashing]: https://en.wikipedia.org/wiki/Perceptual_hashing + +mod bits; +mod phash; +mod signals; + +use bits::Bits64; + +pub use phash::{phash, PHash}; diff --git a/src/image_hash/phash.rs b/src/image_hash/phash.rs new file mode 100644 index 00000000..1adb0031 --- /dev/null +++ b/src/image_hash/phash.rs @@ -0,0 +1,146 @@ +use super::{signals, Bits64}; +use crate::definitions::Image; +use image::{imageops, math::Rect, Luma}; +use std::borrow::Cow; + +/// Stores the result of the [`phash`]. +/// Implements [`Hash`] trait. +/// +/// # Note +/// The hash value may vary between versions. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +pub struct PHash(Bits64); + +impl PHash { + /// Compute the [hamming distance] between hashes. + /// + /// # Example + /// ```no_run + /// use imageproc::image_hash; + /// + /// # fn main() { + /// # let img1 = image::open("first.png").unwrap().to_luma32f(); + /// # let img2 = image::open("second.png").unwrap().to_luma32f(); + /// let hash1 = image_hash::phash(&img1); + /// let hash2 = image_hash::phash(&img2); + /// dbg!(hash1.hamming_distance(hash2)); + /// # } + /// ``` + /// + /// [hamming distance]: https://en.wikipedia.org/wiki/Hamming_distance + pub fn hamming_distance(self, PHash(other): PHash) -> u32 { + self.0.hamming_distance(other) + } +} + +/// Compute the [pHash] using [DCT]. +/// +/// # Example +/// +/// ```no_run +/// use imageproc::image_hash; +/// +/// # fn main() { +/// let img1 = image::open("first.png").unwrap().to_luma32f(); +/// let img2 = image::open("second.png").unwrap().to_luma32f(); +/// let hash1 = image_hash::phash(&img1); +/// let hash2 = image_hash::phash(&img2); +/// dbg!(hash1.hamming_distance(hash2)); +/// # } +/// ``` +/// +/// [pHash]: https://phash.org/docs/pubs/thesis_zauner.pdf +/// [DCT]: https://en.wikipedia.org/wiki/Discrete_cosine_transform +pub fn phash(img: &Image>) -> PHash { + const N: u32 = 8; + const HASH_FACTOR: u32 = 4; + let img = imageops::resize( + img, + HASH_FACTOR * N, + HASH_FACTOR * N, + imageops::FilterType::Lanczos3, + ); + let dct = signals::dct2d(Cow::Owned(img)); + let topleft = Rect { + x: 1, + y: 1, + width: N, + height: N, + }; + let topleft_dct = crate::compose::crop(&dct, topleft); + debug_assert_eq!(topleft_dct.dimensions(), (N, N)); + assert_eq!(topleft_dct.len(), (N * N) as usize); + let mean = + topleft_dct.iter().copied().reduce(|a, b| a + b).unwrap() / (topleft_dct.len() as f32); + let bits = topleft_dct.iter().map(|&x| x > mean); + PHash(Bits64::new(bits)) +} + +#[cfg(test)] +mod tests { + use super::*; + #[test] + fn test_phash() { + let img1 = gray_image!(type: f32, + 1., 2., 3.; + 4., 5., 6. + ); + let mut img2 = img1.clone(); + *img2.get_pixel_mut(0, 0) = Luma([0f32]); + let mut img3 = img2.clone(); + *img3.get_pixel_mut(0, 1) = Luma([0f32]); + + let hash1 = phash(&img1); + let hash2 = phash(&img2); + let hash3 = phash(&img3); + + assert_eq!(0, hash1.hamming_distance(hash1)); + assert_eq!(0, hash2.hamming_distance(hash2)); + assert_eq!(0, hash3.hamming_distance(hash3)); + + assert_eq!(hash1.hamming_distance(hash2), hash2.hamming_distance(hash1)); + + assert!(hash1.hamming_distance(hash2) > 0); + assert!(hash1.hamming_distance(hash3) > 0); + assert!(hash2.hamming_distance(hash3) > 0); + + assert!(hash1.hamming_distance(hash2) < hash1.hamming_distance(hash3)); + } +} + +#[cfg(not(miri))] +#[cfg(test)] +mod proptests { + use super::*; + use crate::proptest_utils::arbitrary_image; + use proptest::prelude::*; + + const N: usize = 100; + + proptest! { + #[test] + fn proptest_phash(img in arbitrary_image(0..N, 0..N)) { + let hash = phash(&img); + assert_eq!(0, hash.hamming_distance(hash)); + } + } +} + +#[cfg(not(miri))] +#[cfg(test)] +mod benches { + use super::*; + use crate::utils::luma32f_bench_image; + use test::{black_box, Bencher}; + + const N: u32 = 600; + + #[bench] + fn bench_phash(b: &mut Bencher) { + let img = luma32f_bench_image(N, N); + b.iter(|| { + let img = black_box(&img); + black_box(phash(img)); + }); + } +} diff --git a/src/image_hash/signals.rs b/src/image_hash/signals.rs new file mode 100644 index 00000000..b4f4e5b6 --- /dev/null +++ b/src/image_hash/signals.rs @@ -0,0 +1,313 @@ +use crate::definitions::Image; +use image::Luma; +use std::borrow::Cow; + +/// Computes the 2 dimensional [DCT]. +/// +/// [DCT]: https://en.wikipedia.org/wiki/Discrete_cosine_transform +pub(super) fn dct2d(img: Cow>>) -> Image> { + #[allow(non_snake_case)] + let T = |img: Cow>| -> Image<_> { + let (w, h) = img.as_ref().dimensions(); + if w == h { + let mut img = img.into_owned(); + transpose_inplace(&mut img); + return img; + } + transpose(img.as_ref()) + }; + + let mut planner = rustdct::DctPlanner::new(); + let rows_ctx = planner.plan_dct2(img.width() as usize); + let cols_ctx = planner.plan_dct2(img.height() as usize); + let cap = { + let rows_cap = rows_ctx.len() + rows_ctx.get_scratch_len(); + let cols_cap = cols_ctx.len() + cols_ctx.get_scratch_len(); + rows_cap.max(cols_cap) + }; + let mut arena = vec![0f32; cap]; + + let dct = dct1d(&T(img), cols_ctx.as_ref(), &mut arena); + dct1d(&T(Cow::Owned(dct)), rows_ctx.as_ref(), &mut arena) +} + +/// Computes the 1 dimensional [DCT] for each row. +/// The required `arena.len()` is equal to `ctx.len() + ctx.get_scratch_len()`. +/// +/// [DCT]: https://en.wikipedia.org/wiki/Discrete_cosine_transform +// TODO: compute inplace +fn dct1d( + img: &Image>, + ctx: &dyn rustdct::TransformType2And3, + arena: &mut [f32], +) -> Image> { + let width = usize::try_from(img.width()).unwrap(); + let height = usize::try_from(img.height()).unwrap(); + assert_eq!(width, ctx.len()); + + let arena_len = ctx.len() + ctx.get_scratch_len(); + assert!(arena_len <= arena.len()); + let (dct_buf, scratch_space) = arena[..arena_len].split_at_mut(ctx.len()); + + let mut dct = |inout: &mut [f32]| { + debug_assert_eq!(inout.len(), ctx.len()); + debug_assert_eq!(scratch_space.len(), ctx.get_scratch_len()); + ctx.process_dct2_with_scratch(inout, scratch_space); + }; + + let mut img_buf = Vec::with_capacity(width * height); + for row in img.rows() { + let row = row.into_iter().map(|p| p.0[0]); + debug_assert_eq!(row.len(), dct_buf.len()); + for (dst, src) in dct_buf.iter_mut().zip(row) { + *dst = src; + } + dct(dct_buf); + img_buf.extend_from_slice(dct_buf); + } + debug_assert_eq!(img_buf.len(), width * height); + Image::from_vec(img.width(), img.height(), img_buf).unwrap() +} + +fn transpose_inplace(img: &mut Image>) { + assert_eq!( + img.width(), + img.height(), + "inplace transposition supported only for square images." + ); + let n = usize::try_from(img.width()).unwrap(); + let at = |row, col| { + debug_assert!(row < n); + debug_assert!(col < n); + row + n * col + }; + let buf = img.as_mut(); + assert_eq!(buf.len(), n * n); + for row in 0..n { + for col in row..n { + let a = buf[at(row, col)]; + let b = buf[at(col, row)]; + buf[at(row, col)] = b; + buf[at(col, row)] = a; + } + } +} + +fn transpose(img: &Image>) -> Image> { + let nwidth = img.height(); + let nheight = img.width(); + Image::from_fn(nwidth, nheight, |x, y| *img.get_pixel(y, x)) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_transpose_inplace() { + #[allow(non_snake_case)] + let T = |img: &Image<_>| { + let mut img = img.clone(); + transpose_inplace(&mut img); + img + }; + let img = gray_image!(type: f32, 1.); + assert_eq!(T(&img), img); + + let img = gray_image!(type: f32, + 1., 2., 3.; + 4., 5., 6.; + 7., 8., 9. + ); + let img_t = gray_image!(type: f32, + 1., 4., 7.; + 2., 5., 8.; + 3., 6., 9. + ); + assert_pixels_eq!(T(&img), img_t); + assert_pixels_eq!(T(&T(&img)), img); + } + #[test] + fn test_transpose() { + #[allow(non_snake_case)] + let T = |img: &Image<_>| transpose(img); + + let img = gray_image!(type: f32, + 1., 2., 3.; + 4., 5., 6. + ); + let img_t = gray_image!(type: f32, + 1., 4.; + 2., 5.; + 3., 6. + ); + assert_pixels_eq!(T(&img), img_t); + assert_pixels_eq!(T(&T(&img)), img); + } + #[test] + fn test_dct1d() { + let mut arena = vec![0f32; 1_000]; + let mut planner = rustdct::DctPlanner::new(); + #[allow(non_snake_case)] + let mut DCT = |img: &Image<_>| -> Image<_> { + let ctx = planner.plan_dct2(img.width() as usize); + dct1d(img, ctx.as_ref(), arena.as_mut()) + }; + let img = gray_image!(type: f32, + 1., 2., 3. + ); + let expected = gray_image!(type: f32, + 6.0, -1.7320508, 0.0 + ); + assert_pixels_eq!(DCT(&img), expected); + + let img = gray_image!(type: f32, + 1., 2., 3.; + 4., 5., 6. + ); + let expected = gray_image!(type: f32, + 6.0, -1.7320508, 0.0; + 15.0, -1.7320508, 0.0 + ); + assert_pixels_eq!(DCT(&img), expected); + + let img = gray_image!(type: f32, + 1., 2., 3.; + 4., 5., 6.; + 7., 8., 9. + ); + let expected = gray_image!(type: f32, + 6.0, -1.7320508, 0.0; + 15.0, -1.7320508, 0.0; + 24.0, -1.7320508, 0.0 + ); + assert_pixels_eq!(DCT(&img), expected); + } + #[test] + fn test_dct2d() { + #[allow(non_snake_case)] + let DCT = |img: &Image<_>| dct2d(Cow::Borrowed(img)); + + let img = gray_image!(type: f32, + 1., 2., 3. + ); + let expected = gray_image!(type: f32, + 6.0, -1.7320508, 0.0 + ); + assert_pixels_eq!(DCT(&img), expected); + + let img = gray_image!(type: f32, + 1., 2., 3.; + 4., 5., 6. + ); + let expected = gray_image!(type: f32, + 21.0, -3.464_101_6, 0.0; + -6.3639607, 0.0, 0.0 + ); + assert_pixels_eq!(DCT(&img), expected); + + let img = gray_image!(type: f32, + 1., 2., 3.; + 4., 5., 6.; + 7., 8., 9. + ); + let expected = gray_image!(type: f32, + 45.0, -5.196_152, 0.0; + -15.588_457, 0.0, 0.0; + 0.0, 0.0, 0.0 + ); + assert_pixels_eq!(DCT(&img), expected); + } +} + +#[cfg(not(miri))] +#[cfg(test)] +mod proptests { + use super::*; + use crate::proptest_utils::arbitrary_image; + use proptest::prelude::*; + + const N: usize = 100; + + proptest! { + #[test] + fn proptest_transpose(img in arbitrary_image(0..N, 0..N)) { + let img_t = transpose(&img); + assert_eq!(img.width(), img_t.height()); + assert_eq!(img.height(), img_t.width()); + assert_pixels_eq!(img, transpose(&img_t)); + } + #[test] + fn proptest_transpose_inplace(img in arbitrary_image(0..N, 0..N)) { + #[allow(non_snake_case)] + let T = |img: &Image<_>| { + let mut img = img.clone(); + transpose_inplace(&mut img); + img + }; + if img.width() != img.height() { + return Ok(()); + } + assert_pixels_eq!(T(&img), transpose(&img)); + assert_pixels_eq!(T(&T(&img)), img); + } + #[test] + fn proptest_dct1d(img in arbitrary_image(0..N, 0..N)) { + let ctx = rustdct::DctPlanner::new().plan_dct2(img.width() as usize); + let mut arena = vec![0f32; ctx.len() + ctx.get_scratch_len()]; + let dct = dct1d(&img, ctx.as_ref(), arena.as_mut()); + assert_eq!(dct.dimensions(), img.dimensions()); + } + #[test] + fn proptest_dct2d(img in arbitrary_image(0..N, 0..N)) { + let dct = dct2d(Cow::Borrowed(&img)); + assert_eq!(dct.dimensions(), img.dimensions()); + } + } +} + +#[cfg(not(miri))] +#[cfg(test)] +mod benches { + use super::*; + use crate::utils::luma32f_bench_image; + use test::{black_box, Bencher}; + + const N: u32 = 600; + + #[bench] + fn bench_transpose(b: &mut Bencher) { + let img = luma32f_bench_image(N, N); + b.iter(|| { + black_box(transpose(&img)); + }); + } + #[bench] + fn bench_transpose_inplace(b: &mut Bencher) { + let mut img = luma32f_bench_image(N, N); + b.iter(|| { + black_box(&mut img); + transpose_inplace(&mut img); + }); + } + + #[bench] + fn bench_dct1d(b: &mut Bencher) { + let img = luma32f_bench_image(N, N); + let ctx = rustdct::DctPlanner::new().plan_dct2(img.width() as usize); + let mut arena = vec![0f32; ctx.len() + ctx.get_scratch_len()]; + b.iter(|| { + let dct = dct1d(black_box(&img), ctx.as_ref(), &mut arena); + black_box(dct); + }); + } + + #[bench] + fn bench_dct2d(b: &mut Bencher) { + let img = luma32f_bench_image(N, N); + b.iter(|| { + let dct = dct2d(black_box(Cow::Borrowed(&img))); + black_box(dct); + }); + } +} diff --git a/src/lib.rs b/src/lib.rs index 1a0cdc89..e595bd78 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -37,6 +37,7 @@ pub mod gradients; pub mod haar; pub mod hog; pub mod hough; +pub mod image_hash; pub mod integral_image; pub mod kernel; pub mod local_binary_patterns; diff --git a/src/utils.rs b/src/utils.rs index 033fa286..efd116bc 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -1,15 +1,16 @@ //! Utils for testing and debugging. +use crate::definitions::Image; use image::{ open, DynamicImage, GenericImage, GenericImageView, GrayImage, Luma, Pixel, Rgb, RgbImage, }; - use itertools::Itertools; + use std::cmp::{max, min}; use std::collections::HashSet; -use std::fmt; -use std::fmt::Write; +use std::hint::black_box; use std::path::Path; +use std::{fmt, fmt::Write}; /// Helper for defining greyscale images. /// @@ -675,7 +676,14 @@ pub fn gray_bench_image(width: u32, height: u32) -> GrayImage { image.put_pixel(x, y, Luma([intensity])); } } - image + black_box(image) +} + +/// `Luma` image to use in benchmarks. See comment on `gray_bench_image`. +pub fn luma32f_bench_image(width: u32, height: u32) -> Image> { + use image::DynamicImage; + let img = gray_bench_image(width, height); + DynamicImage::ImageLuma8(img).to_luma32f() } /// RGB image to use in benchmarks. See comment on `gray_bench_image`. @@ -690,7 +698,7 @@ pub fn rgb_bench_image(width: u32, height: u32) -> RgbImage { image.put_pixel(x, y, Rgb([r, g, b])); } } - image + black_box(image) } #[cfg(test)]