diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 82cd9b0..d64d94a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -18,7 +18,7 @@ jobs: - uses: actions/checkout@v6 - name: Install Rust toolchain - uses: dtolnay/rust-toolchain@stable + uses: dtolnay/rust-toolchain@1.98.0 with: components: rustfmt, clippy @@ -37,6 +37,9 @@ jobs: - name: Run tests run: cargo test --workspace + - name: Run tests (sequential, parallel disabled) + run: cargo test --workspace --no-default-features --features std + - name: Run tests (with parallel feature) run: cargo test --workspace --features parallel @@ -54,7 +57,7 @@ jobs: - uses: actions/checkout@v6 - name: Install Rust toolchain - uses: dtolnay/rust-toolchain@stable + uses: dtolnay/rust-toolchain@1.98.0 with: components: clippy targets: thumbv7em-none-eabihf @@ -82,7 +85,7 @@ jobs: - uses: actions/checkout@v6 - name: Install Rust toolchain - uses: dtolnay/rust-toolchain@stable + uses: dtolnay/rust-toolchain@1.98.0 with: targets: wasm32-unknown-unknown @@ -107,7 +110,7 @@ jobs: - uses: actions/checkout@v6 - name: Install Rust toolchain - uses: dtolnay/rust-toolchain@stable + uses: dtolnay/rust-toolchain@1.98.0 - name: Rust Cache uses: Swatinem/rust-cache@v2 diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index e79e3a1..faad884 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -104,7 +104,7 @@ jobs: steps: - uses: actions/checkout@v6 - name: Install Rust - uses: dtolnay/rust-toolchain@stable + uses: dtolnay/rust-toolchain@1.98.0 - name: Publish to Crates.io run: cargo publish --token ${{ secrets.CRATES_TOKEN }} -p purecv @@ -114,7 +114,7 @@ jobs: steps: - uses: actions/checkout@v6 - name: Install Rust - uses: dtolnay/rust-toolchain@stable + uses: dtolnay/rust-toolchain@1.98.0 - name: Install wasm-pack run: curl https://rustwasm.github.io/wasm-pack/installer/init.sh -sSf | sh - name: Build and Publish to NPM diff --git a/CHANGELOG.md b/CHANGELOG.md index ce645a4..a06f71d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,48 @@ All notable changes to this project will be documented in this file. +## [0.8.0] - 2026-09-02 + +### βš™οΈ Miscellaneous Tasks + +- *(ci)* Pin stable toolchain and allow the new chunks_exact_to_as_chunks lint + +### ⚑ Performance + +- *(imgproc)* SIMD-accelerate compare_hist +- *(imgproc)* SIMD-accelerate compare_hist +- *(ci)* Run the sequential (non-parallel) test suite in CI + +### πŸ› Bug Fixes + +- *(imgproc)* Address review comments for histogram module +- *(imgproc)* Guard calc_back_project against zero-width images +- *(imgproc)* Address Qodo review comments for SIMD compare_hist +- *(wasm)* Address Qodo review findings on histogram bindings +- Address Qodo review findings on histogram examples + +### πŸ“š Documentation + +- Record the release workflow and merge policy in CLAUDE.md +- Add histogram/CLAHE examples (Rust + WASM) and update READMEs + +### πŸ•ΈοΈ WebAssembly & Emscripten + +- *(wasm)* Expose histogram module bindings + +### πŸš€ Features + +- *(imgproc)* Add histogram module +- *(imgproc)* Add parallel support to histogram module + +### 🚜 Refactor + +- *(imgproc)* Use as_chunks instead of allowing the new clippy lint + +### πŸ§ͺ Testing + +- *(imgproc)* Move histogram module tests into src/imgproc/tests.rs + ## [0.7.1] - 2026-08-10 ### βš™οΈ Miscellaneous Tasks diff --git a/CLAUDE.md b/CLAUDE.md index 3b18020..ca03057 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -51,6 +51,29 @@ Preferred scopes: `(core)` `(imgproc)` `(simd)` `(wasm)` `(parallel)` - PRs start from and target the `dev` branch (not `main`). - Keep PRs focused; one feature or fix per PR. +### Releases + +Releases go `dev` β†’ PR to `main` β†’ merge β†’ tag `vX.Y.Z` on `main`. Pushing the tag +triggers `release.yml`, which publishes to crates.io **and** npm β€” irreversible. + +**Merge release PRs with a merge commit, not "Rebase and merge" or "Squash and +merge".** Rebasing replays dev's commits as new objects on `main`, so `dev` stops +being an ancestor of `main` and the two diverge with identical content but different +SHAs. The next release PR then replays every old commit again. This happened with +v0.7.1 (#95) and had to be repaired by resetting `dev` to `main`. + +Release prep steps (see the v0.7.1 commit for a worked example): + +1. Bump the version in `Cargo.toml` (**two places** β€” `[package]` and + `[workspace.package]`) and in the root `package.json`. +2. Run `npm run build` β€” this regenerates `crates/wasm/pkg/package.json`, which is + tracked and otherwise silently drifts. +3. `npx git-cliff --config cliff.toml --tag vX.Y.Z --unreleased --prepend CHANGELOG.md`, + then add the blank line `--prepend` omits before the previous version heading. +4. Commit as `chore(release): prepare for vX.Y.Z` β€” `cliff.toml` skips this message + from the changelog, so land any other changes in their own commits *first* or they + will not appear. + ## Module Structure Convention Each top-level module (`core/`, `imgproc/`, `features2d/`, etc.) must contain its own: diff --git a/Cargo.toml b/Cargo.toml index 5a1b509..e35cce0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "purecv" -version = "0.7.1" +version = "0.8.0" authors = ["Walter Perdan "] edition = "2021" rust-version = "1.88" @@ -86,7 +86,7 @@ members = ["crates/wasm"] exclude = ["crates/no-std-smoke"] [workspace.package] -version = "0.7.1" +version = "0.8.0" authors = ["Walter Perdan "] edition = "2021" description = "A pure Rust, high-performance computer vision library focused on safety and portability." diff --git a/MAINTAINERS.md b/MAINTAINERS.md index f8ba2d9..a07207f 100644 --- a/MAINTAINERS.md +++ b/MAINTAINERS.md @@ -28,7 +28,11 @@ Publishing a new version requires a mix of manual changelog curation and automat 2. Verify that all CI checks (Formatting, Clippy, Tests for `parallel` and `simd`) are passing on the latest commit. ### Step 2: Bump the Version -Update the version number in the `Cargo.toml` file of the workspace in [package] and [workspace.package] sections. +Update the version number in the `Cargo.toml` file of the workspace in [package] and [workspace.package] sections, and in the root `package.json`. + +Also check `README.md` and `crates/wasm/README.md` for hardcoded version strings in +installation snippets (e.g. `purecv = "0.6"`) β€” these don't update automatically and +are easy to miss. Search for the old version number across both files before moving on. ### Step 3: Generate the Local Changelog We use `git-cliff` to parse the conventional commits and update the historical changelog. Run the following command in the root directory: diff --git a/README.md b/README.md index 86befb2..ddf2178 100644 --- a/README.md +++ b/README.md @@ -60,6 +60,7 @@ Unlike existing wrappers, **PureCV** is a native rewrite. It aims to provide: - **Hough Transform:** Standard (`hough_lines`) and Probabilistic (`hough_lines_p`) line detection, plus Hough Circle Transform (`hough_circles`) using internally computed Sobel gradients. Fully parallelized via the `parallel` feature. - **Resizing:** `resize` function utilizing high-performance bilinear interpolation, fully compatible with `parallel` Rayon multi-threading. - **Geometric Transformations:** `remap` (with bilinear and nearest-neighbor interpolation) and `warp_perspective` (perspective transformations) fully parallelized and SIMD-accelerated. +- **Histograms & Contrast:** `calc_hist` (multi-dimensional, uniform or non-uniform bins, optional mask and accumulation) and `calc_back_project` for histogram back-projection; `compare_hist` with all 6 OpenCV comparison methods (`Correl`, `ChiSqr`, `ChiSqrAlt`, `Intersection`, `Bhattacharyya`, `KullbackLeibler`), SIMD-accelerated via the `simd` feature; `equalize_hist` for global histogram equalization; and `Clahe` (Contrast Limited Adaptive Histogram Equalization) for `u8`/`u16` images. `calc_hist`, `calc_back_project`, `equalize_hist`, and `Clahe` are all parallelized via the `parallel` feature. ### `purecv-features2d` - **FAST Feature Detector:** Real-time corner detector (`FastFeatureDetector`) supporting Type 5_8, 7_12, and 9_16 neighborhood configurations, plus optional non-maximum suppression. @@ -86,7 +87,7 @@ Add the following to your `Cargo.toml`: ```toml [dependencies] -purecv = "0.6" +purecv = "0.8" ``` PureCV's minimum supported Rust version (MSRV) is **1.88**. @@ -104,7 +105,7 @@ PureCV's minimum supported Rust version (MSRV) is **1.88**. ### `no_std` / embedded support Build with `--no-default-features` to run on bare-metal targets such as the -ESP32 (`purecv = { version = "0.6", default-features = false }`). Only `core` +ESP32 (`purecv = { version = "0.8", default-features = false }`). Only `core` and `alloc` are required (an allocator must be provided by the target). | Module | `no_std` | Notes | @@ -120,7 +121,7 @@ features gives the scalar, single-threaded code paths. ```toml [dependencies] -purecv = { version = "0.6", default-features = false } +purecv = { version = "0.8", default-features = false } ``` ```rust @@ -149,14 +150,14 @@ To enable the `ndarray` feature: ```toml [dependencies] -purecv = { version = "0.6", features = ["ndarray"] } +purecv = { version = "0.8", features = ["ndarray"] } ``` To enable SIMD + Parallel for maximum performance: ```toml [dependencies] -purecv = { version = "0.6", features = ["parallel", "simd"] } +purecv = { version = "0.8", features = ["parallel", "simd"] } ``` ### Usage Example @@ -312,6 +313,9 @@ cargo run --example morphology # Gaussian pyramids (pyr_down, pyr_up) cargo run --example pyramids +# Histograms & contrast (calc_hist, calc_back_project, compare_hist, equalize_hist, CLAHE) +cargo run --example histogram + # Hough Transform (Lines and Circles detection) cargo run --example hough_transform @@ -344,10 +348,10 @@ cargo run --example rectification ## πŸ§ͺ Testing & Benchmarking ### Running Tests -PureCV uses a comprehensive suite of unit tests to ensure correctness and parity with OpenCV. The test suite currently includes **308 unit tests** (plus **40 doc-tests**) covering: +PureCV uses a comprehensive suite of unit tests to ensure correctness and parity with OpenCV. The test suite currently includes **342 unit tests** (plus **40 doc-tests**) covering: - **Core module:** Matrix factories, scalar arithmetic variants, bitwise scalar ops, min/max, comparison ops (`compare`, `in_range`), reduction (`reduce`, `count_non_zero`), polar/cartesian conversions, linear algebra (`determinant`, `invert`, `solve`), channel ops (`extract_channel`, `insert_channel`), `DynamicMatrix`, transforms, sorting, clustering, and RNG. -- **Imgproc module:** Filters, derivatives, edge detection, color conversions (including gray-to-RGB/BGR/RGBA/BGRA), thresholding, morphology (`erode`, `dilate`), pyramids (`pyr_down`, `pyr_up`), and kernel helpers (`get_gaussian_kernel`, `get_sobel_kernels`). +- **Imgproc module:** Filters, derivatives, edge detection, color conversions (including gray-to-RGB/BGR/RGBA/BGRA), thresholding, morphology (`erode`, `dilate`), pyramids (`pyr_down`, `pyr_up`), kernel helpers (`get_gaussian_kernel`, `get_sobel_kernels`), and histograms/CLAHE (`calc_hist`, `calc_back_project`, `compare_hist`, `equalize_hist`, `Clahe`). - **Features2d module:** Keypoint structures (`KeyPoint`), FAST corner detection (`FastFeatureDetector`), scale pyramids, and ORB feature extraction & BRIEF descriptor extraction (`Orb`). - **Video module:** Tracking and optical flow capabilities including `calc_optical_flow_pyr_lk` and `build_optical_flow_pyramid` implementations. - **Calib3d module:** SVD, homography estimation, pose estimation (`solve_pnp`), and `rodrigues`. diff --git a/benches/imgproc_bench.rs b/benches/imgproc_bench.rs index a480c1b..9afb9b4 100644 --- a/benches/imgproc_bench.rs +++ b/benches/imgproc_bench.rs @@ -41,6 +41,7 @@ use purecv::imgproc::derivatives::{laplacian, scharr, sobel}; use purecv::imgproc::edge::canny; use purecv::imgproc::feature::corner_harris; use purecv::imgproc::filter::{bilateral_filter, box_filter, gaussian_blur}; +use purecv::imgproc::histogram::{calc_back_project, calc_hist, equalize_hist, Clahe, RangeSpec}; use purecv::imgproc::hough::{hough_circles, hough_lines, hough_lines_p}; use purecv::imgproc::threshold::{threshold, ThresholdTypes}; use purecv::imgproc::{cvt_color, ColorConversionCode}; @@ -254,6 +255,100 @@ fn bench_imgproc(c: &mut Criterion) { .unwrap() }) }); + + // calc_hist benchmark setup + let mut img_hist = Matrix::::new(size, size, 1); + for (i, p) in img_hist.data.iter_mut().enumerate() { + *p = (i % 256) as u8; + } + let hist_ranges = [RangeSpec::Uniform(0.0, 256.0)]; + + c.bench_function("calc_hist_1024x1024", |b| { + b.iter(|| { + calc_hist( + black_box(&[&img_hist]), + &[0], + None, + &[256], + &hist_ranges, + false, + None, + ) + .unwrap() + }) + }); + + // calc_back_project benchmark setup + let hist_for_backproj = + calc_hist(&[&img_hist], &[0], None, &[256], &hist_ranges, false, None).unwrap(); + + c.bench_function("calc_back_project_1024x1024", |b| { + b.iter(|| { + calc_back_project( + black_box(&[&img_hist]), + &[0], + &[256], + &hist_for_backproj, + &hist_ranges, + 1.0, + ) + .unwrap() + }) + }); + + // equalize_hist benchmark setup + c.bench_function("equalize_hist_1024x1024", |b| { + b.iter(|| equalize_hist(black_box(&img_hist)).unwrap()) + }); + + // Clahe::apply_u8 benchmark setup + let clahe = Clahe::new(2.0, Size2i::new(8, 8)); + + c.bench_function("clahe_apply_u8_1024x1024", |b| { + b.iter(|| clahe.apply_u8(black_box(&img_hist)).unwrap()) + }); + // compare_hist benchmark setup + let mut hist1 = Matrix::::new(256, 1, 1); + let mut hist2 = Matrix::::new(256, 1, 1); + for (i, p) in hist1.data.iter_mut().enumerate() { + *p = (i as f32 * 0.1).sin().abs(); + } + for (i, p) in hist2.data.iter_mut().enumerate() { + *p = (i as f32 * 0.1).cos().abs(); + } + + c.bench_function("compare_hist_correl_256", |b| { + b.iter(|| { + purecv::imgproc::histogram::compare_hist( + black_box(&hist1), + black_box(&hist2), + purecv::imgproc::histogram::HistCompMethods::Correl, + ) + .unwrap() + }) + }); + + c.bench_function("compare_hist_intersection_256", |b| { + b.iter(|| { + purecv::imgproc::histogram::compare_hist( + black_box(&hist1), + black_box(&hist2), + purecv::imgproc::histogram::HistCompMethods::Intersection, + ) + .unwrap() + }) + }); + + c.bench_function("compare_hist_kullback_256", |b| { + b.iter(|| { + purecv::imgproc::histogram::compare_hist( + black_box(&hist1), + black_box(&hist2), + purecv::imgproc::histogram::HistCompMethods::KullbackLeibler, + ) + .unwrap() + }) + }); } criterion_group!(benches, bench_imgproc); diff --git a/crates/wasm/README.md b/crates/wasm/README.md index 57f945c..573f821 100644 --- a/crates/wasm/README.md +++ b/crates/wasm/README.md @@ -58,8 +58,41 @@ Because WebAssembly runs linearly in memory and holds pointers to Rust `Vec` obj Right now we have covered a large majority of operations for `core` and `imgproc`, and have started on `calib3d` and `video`: - **Core**: Arithmetic (`add`, `subtract`, `multiply`, `absdiff` etc.), Structural (`hconcat`, `vconcat`, `flip`), Geometry, constants etc. -- **ImgProc**: Filters (`blur`, `gaussian_blur`, `bilateral_filter`), Thresholding (`threshold`), Coloring (`cvt_color`), Edge Derivatives (`canny`, `sobel`, `laplacian`), Morphology (`erode`, `dilate`, `morphology_ex`, `get_structuring_element`), Pyramids (`pyr_down`, `pyr_up`, `build_pyramid`), Feature Detection (`good_features_to_track`, `corner_sub_pix`). +- **ImgProc**: Filters (`blur`, `gaussian_blur`, `bilateral_filter`), Thresholding (`threshold`), Coloring (`cvt_color`), Edge Derivatives (`canny`, `sobel`, `laplacian`), Morphology (`erode`, `dilate`, `morphology_ex`, `get_structuring_element`), Pyramids (`pyr_down`, `pyr_up`, `build_pyramid`), Feature Detection (`good_features_to_track`, `corner_sub_pix`), Histograms (`calcHistUniform`, `calcHistNonUniform`, `calcBackProjectUniform`, `calcBackProjectNonUniform`, `compareHist`, `equalizeHist`, `Clahe`). - **Video**: Optical Flow (`calc_optical_flow_pyr_lk`). - **Calib3d**: Pose Estimation (`solve_pnp`, `solve_pnp_ransac`), Homography (`find_homography`), and geometry (`rodrigues`). +### Histogram & Contrast Operations + +```javascript +import { Mat, MatVector, calcHistUniform, equalizeHist, Clahe } from '@webarkit/purecv-wasm'; + +// 1. Equalize Histogram (8-bit grayscale only) +const eqMat = equalizeHist(grayMat); + +// 2. CLAHE (Contrast Limited Adaptive Histogram Equalization) +const clahe = new Clahe(40.0, 8, 8); +const enhancedMat = clahe.apply(grayMat); + +// 3. Dense Histogram with uniform bins +const images = new MatVector(); +images.push(grayMat); +const channels = [0]; +const histSize = [256]; +const ranges = [0.0, 256.0]; // [min, max] per channel +const hist = calcHistUniform(images, channels, undefined, histSize, ranges, false, undefined); + +// To accumulate onto a previous histogram, pass accumulate=true and the +// existing histogram Mat as the last argument instead of undefined: +// calcHistUniform(images, channels, undefined, histSize, ranges, true, hist); +``` + +*Note:* `equalizeHist` and `Clahe.apply` currently support single-channel 8-bit images (`CV_8UC1`). Support for 16-bit images (`CV_16UC1`) is planned for a future release. + +*Note:* `calcBackProjectUniform`/`calcBackProjectNonUniform` take an explicit +`histSize` argument (right after `channels`) describing the shape the +histogram was built with β€” a flat histogram's bin count alone can't be +unambiguously reconstructed into a multi-dimensional shape (e.g. 8 bins +could be `[8]` or `[2, 4]`). + Note: To interface between JavaScript Typed Arrays and `purecv-wasm`, please use the available getter functions (`.data()`) which directly retrieve a Float32Array or Uint8Array view into WASM memory. diff --git a/crates/wasm/pkg/README.md b/crates/wasm/pkg/README.md index 57f945c..573f821 100644 --- a/crates/wasm/pkg/README.md +++ b/crates/wasm/pkg/README.md @@ -58,8 +58,41 @@ Because WebAssembly runs linearly in memory and holds pointers to Rust `Vec` obj Right now we have covered a large majority of operations for `core` and `imgproc`, and have started on `calib3d` and `video`: - **Core**: Arithmetic (`add`, `subtract`, `multiply`, `absdiff` etc.), Structural (`hconcat`, `vconcat`, `flip`), Geometry, constants etc. -- **ImgProc**: Filters (`blur`, `gaussian_blur`, `bilateral_filter`), Thresholding (`threshold`), Coloring (`cvt_color`), Edge Derivatives (`canny`, `sobel`, `laplacian`), Morphology (`erode`, `dilate`, `morphology_ex`, `get_structuring_element`), Pyramids (`pyr_down`, `pyr_up`, `build_pyramid`), Feature Detection (`good_features_to_track`, `corner_sub_pix`). +- **ImgProc**: Filters (`blur`, `gaussian_blur`, `bilateral_filter`), Thresholding (`threshold`), Coloring (`cvt_color`), Edge Derivatives (`canny`, `sobel`, `laplacian`), Morphology (`erode`, `dilate`, `morphology_ex`, `get_structuring_element`), Pyramids (`pyr_down`, `pyr_up`, `build_pyramid`), Feature Detection (`good_features_to_track`, `corner_sub_pix`), Histograms (`calcHistUniform`, `calcHistNonUniform`, `calcBackProjectUniform`, `calcBackProjectNonUniform`, `compareHist`, `equalizeHist`, `Clahe`). - **Video**: Optical Flow (`calc_optical_flow_pyr_lk`). - **Calib3d**: Pose Estimation (`solve_pnp`, `solve_pnp_ransac`), Homography (`find_homography`), and geometry (`rodrigues`). +### Histogram & Contrast Operations + +```javascript +import { Mat, MatVector, calcHistUniform, equalizeHist, Clahe } from '@webarkit/purecv-wasm'; + +// 1. Equalize Histogram (8-bit grayscale only) +const eqMat = equalizeHist(grayMat); + +// 2. CLAHE (Contrast Limited Adaptive Histogram Equalization) +const clahe = new Clahe(40.0, 8, 8); +const enhancedMat = clahe.apply(grayMat); + +// 3. Dense Histogram with uniform bins +const images = new MatVector(); +images.push(grayMat); +const channels = [0]; +const histSize = [256]; +const ranges = [0.0, 256.0]; // [min, max] per channel +const hist = calcHistUniform(images, channels, undefined, histSize, ranges, false, undefined); + +// To accumulate onto a previous histogram, pass accumulate=true and the +// existing histogram Mat as the last argument instead of undefined: +// calcHistUniform(images, channels, undefined, histSize, ranges, true, hist); +``` + +*Note:* `equalizeHist` and `Clahe.apply` currently support single-channel 8-bit images (`CV_8UC1`). Support for 16-bit images (`CV_16UC1`) is planned for a future release. + +*Note:* `calcBackProjectUniform`/`calcBackProjectNonUniform` take an explicit +`histSize` argument (right after `channels`) describing the shape the +histogram was built with β€” a flat histogram's bin count alone can't be +unambiguously reconstructed into a multi-dimensional shape (e.g. 8 bins +could be `[8]` or `[2, 4]`). + Note: To interface between JavaScript Typed Arrays and `purecv-wasm`, please use the available getter functions (`.data()`) which directly retrieve a Float32Array or Uint8Array view into WASM memory. diff --git a/crates/wasm/pkg/package.json b/crates/wasm/pkg/package.json index 923a255..989b496 100644 --- a/crates/wasm/pkg/package.json +++ b/crates/wasm/pkg/package.json @@ -5,7 +5,7 @@ "Walter Perdan \u003chttps://github.com/kalwalt\u003e" ], "description": "A pure Rust, high-performance computer vision library focused on safety and portability.", - "version": "0.7.1", + "version": "0.8.0", "license": "LGPL-2.1-or-later", "repository": { "type": "git", diff --git a/crates/wasm/src/lib.rs b/crates/wasm/src/lib.rs index a25f6c2..8f5a5e1 100644 --- a/crates/wasm/src/lib.rs +++ b/crates/wasm/src/lib.rs @@ -47,6 +47,7 @@ use purecv::imgproc::derivatives; use purecv::imgproc::edge; use purecv::imgproc::feature; use purecv::imgproc::filter; +use purecv::imgproc::histogram::{self, Clahe as CoreClahe, HistCompMethods, RangeSpec}; use purecv::imgproc::hough; use purecv::imgproc::morph::{self, MorphShapes, MorphTypes}; use purecv::imgproc::pyramid; @@ -2753,3 +2754,630 @@ impl ORB { Ok(obj.into()) } } +// --------------------------------------------------------------------------- +// MatVector - Collection of Mats +// --------------------------------------------------------------------------- + +/// Managed vector of Mat objects exposed to JavaScript. +/// +/// Mirrors the cv.MatVector class from OpenCV.js. +#[wasm_bindgen(js_name = "MatVector")] +pub struct MatVector { + pub(crate) inner: Vec, +} + +#[wasm_bindgen(js_class = "MatVector")] +impl MatVector { + #[wasm_bindgen(constructor)] + pub fn new() -> Self { + Self { inner: Vec::new() } + } + + /// Adds a Mat to the vector (clones the underlying matrix). + pub fn push(&mut self, mat: &Mat) { + self.inner.push(mat.inner.clone()); + } + + /// Alias for push, matching OpenCV.js conventions. + #[wasm_bindgen(js_name = "push_back")] + pub fn push_back(&mut self, mat: &Mat) { + self.push(mat); + } + + /// Returns the Mat at the given index, or an error if out of bounds. + pub fn get(&self, idx: usize) -> Result { + self.inner + .get(idx) + .cloned() + .map(|inner| Mat { inner }) + .ok_or_else(|| JsError::new("Index out of bounds")) + } + + /// Returns the number of matrices in the vector. + pub fn size(&self) -> usize { + self.inner.len() + } + + /// Returns the number of matrices in the vector. + pub fn length(&self) -> usize { + self.inner.len() + } + + /// Returns true if the vector contains no matrices. + #[wasm_bindgen(js_name = "isEmpty")] + pub fn is_empty(&self) -> bool { + self.inner.is_empty() + } +} + +impl Default for MatVector { + fn default() -> Self { + Self::new() + } +} + +// --------------------------------------------------------------------------- +// JS-side enum constants: Histogram comparison methods +// --------------------------------------------------------------------------- + +#[wasm_bindgen(js_name = "HIST_CMP_CORREL")] +pub fn hist_cmp_correl() -> i32 { + 0 +} +#[wasm_bindgen(js_name = "HIST_CMP_CHISQR")] +pub fn hist_cmp_chisqr() -> i32 { + 1 +} +#[wasm_bindgen(js_name = "HIST_CMP_CHISQR_ALT")] +pub fn hist_cmp_chisqr_alt() -> i32 { + 2 +} +#[wasm_bindgen(js_name = "HIST_CMP_INTERSECT")] +pub fn hist_cmp_intersect() -> i32 { + 3 +} +#[wasm_bindgen(js_name = "HIST_CMP_BHATTACHARYYA")] +pub fn hist_cmp_bhattacharyya() -> i32 { + 4 +} +#[wasm_bindgen(js_name = "HIST_CMP_KL_DIV")] +pub fn hist_cmp_kl_div() -> i32 { + 5 +} + +fn hist_comp_method_from_i32(m: i32) -> Result { + match m { + 0 => Ok(HistCompMethods::Correl), + 1 => Ok(HistCompMethods::ChiSqr), + 2 => Ok(HistCompMethods::ChiSqrAlt), + 3 => Ok(HistCompMethods::Intersection), + 4 => Ok(HistCompMethods::Bhattacharyya), + 5 => Ok(HistCompMethods::KullbackLeibler), + _ => Err(JsError::new(&format!("Unknown hist comp method: {m}"))), + } +} + +// --------------------------------------------------------------------------- +// Histogram operations +// --------------------------------------------------------------------------- + +/// Computes a joint dense histogram for a set of images using uniform ranges. +/// +/// * images - MatVector of input images (u8 or f32 depth, all matching). +/// * channels - List of channel indices across the images. +/// * mask - Optional single-channel u8 mask (CV_8UC1). +/// * hist_size - Array of bin counts for each histogram dimension. +/// * ranges - Flat array of [min_0, max_0, min_1, max_1, ...] per dimension. +/// * accumulate - If true, adds to `existing_hist` instead of starting from zero. +/// * existing_hist - The histogram to accumulate into when `accumulate` is true +/// (CV_32FC1, same size as the output). Ignored when `accumulate` is false. +#[wasm_bindgen(js_name = "calcHistUniform")] +pub fn calc_hist_uniform( + images: &MatVector, + channels: &[usize], + mask: Option, + hist_size: &[usize], + ranges: &[f32], + accumulate: bool, + existing_hist: Option, +) -> Result { + if images.inner.is_empty() { + return Err(JsError::new( + "calc_hist_uniform: images vector must not be empty", + )); + } + + if ranges.len() != hist_size.len() * 2 { + return Err(JsError::new( + "ranges array must have 2 elements (min, max) per dimension", + )); + } + + let mut specs = Vec::with_capacity(hist_size.len()); + for i in 0..hist_size.len() { + specs.push(RangeSpec::Uniform(ranges[i * 2], ranges[i * 2 + 1])); + } + + let mask_ref = match &mask { + Some(m) => Some(require_u8(m, "calc_hist_uniform (mask)")?), + None => None, + }; + let existing_hist_ref = match &existing_hist { + Some(h) => Some(require_f32(h, "calc_hist_uniform (existing_hist)")?), + None => None, + }; + + let first = &images.inner[0]; + let hist = match first { + DynamicMatrix { + data: DynamicData::U8(_), + } => { + let mut u8_mats = Vec::with_capacity(images.inner.len()); + for dyn_mat in &images.inner { + if let Some(m) = dyn_mat.as_matrix_u8() { + u8_mats.push(m); + } else { + return Err(JsError::new( + "calc_hist_uniform: all images must have consistent depth (u8)", + )); + } + } + histogram::calc_hist( + &u8_mats, + channels, + mask_ref, + hist_size, + &specs, + accumulate, + existing_hist_ref, + ) + .map_err(|e| JsError::new(&format!("{e}")))? + } + DynamicMatrix { + data: DynamicData::F32(_), + } => { + let mut f32_mats = Vec::with_capacity(images.inner.len()); + for dyn_mat in &images.inner { + if let Some(m) = dyn_mat.as_matrix_f32() { + f32_mats.push(m); + } else { + return Err(JsError::new( + "calc_hist_uniform: all images must have consistent depth (f32)", + )); + } + } + histogram::calc_hist( + &f32_mats, + channels, + mask_ref, + hist_size, + &specs, + accumulate, + existing_hist_ref, + ) + .map_err(|e| JsError::new(&format!("{e}")))? + } + _ => { + return Err(JsError::new( + "calc_hist_uniform: unsupported image depth (must be u8 or f32)", + )); + } + }; + + Ok(Mat { + inner: DynamicMatrix { + data: DynamicData::F32(hist), + }, + }) +} + +/// Computes a joint dense histogram for a set of images using non-uniform ranges. +/// +/// * images - MatVector of input images (u8 or f32 depth, all matching). +/// * channels - List of channel indices across the images. +/// * mask - Optional single-channel u8 mask (CV_8UC1). +/// * hist_size - Array of bin counts for each histogram dimension. +/// * ranges - Array of Float32Arrays, one per dimension, containing bin boundary coordinates. +/// * accumulate - If true, adds to `existing_hist` instead of starting from zero. +/// * existing_hist - The histogram to accumulate into when `accumulate` is true +/// (CV_32FC1, same size as the output). Ignored when `accumulate` is false. +#[wasm_bindgen(js_name = "calcHistNonUniform")] +pub fn calc_hist_non_uniform( + images: &MatVector, + channels: &[usize], + mask: Option, + hist_size: &[usize], + ranges: js_sys::Array, + accumulate: bool, + existing_hist: Option, +) -> Result { + if images.inner.is_empty() { + return Err(JsError::new( + "calc_hist_non_uniform: images vector must not be empty", + )); + } + + if ranges.length() as usize != hist_size.len() { + return Err(JsError::new( + "ranges array must have one Float32Array per dimension", + )); + } + + let mut specs = Vec::with_capacity(hist_size.len()); + for i in 0..hist_size.len() { + let js_val = ranges.get(i as u32); + let f32_array = js_sys::Float32Array::from(js_val); + let mut vec = vec![0.0f32; f32_array.length() as usize]; + f32_array.copy_to(&mut vec); + specs.push(RangeSpec::NonUniform(vec)); + } + + let mask_ref = match &mask { + Some(m) => Some(require_u8(m, "calc_hist_non_uniform (mask)")?), + None => None, + }; + let existing_hist_ref = match &existing_hist { + Some(h) => Some(require_f32(h, "calc_hist_non_uniform (existing_hist)")?), + None => None, + }; + + let first = &images.inner[0]; + let hist = match first { + DynamicMatrix { + data: DynamicData::U8(_), + } => { + let mut u8_mats = Vec::with_capacity(images.inner.len()); + for dyn_mat in &images.inner { + if let Some(m) = dyn_mat.as_matrix_u8() { + u8_mats.push(m); + } else { + return Err(JsError::new( + "calc_hist_non_uniform: all images must have consistent depth (u8)", + )); + } + } + histogram::calc_hist( + &u8_mats, + channels, + mask_ref, + hist_size, + &specs, + accumulate, + existing_hist_ref, + ) + .map_err(|e| JsError::new(&format!("{e}")))? + } + DynamicMatrix { + data: DynamicData::F32(_), + } => { + let mut f32_mats = Vec::with_capacity(images.inner.len()); + for dyn_mat in &images.inner { + if let Some(m) = dyn_mat.as_matrix_f32() { + f32_mats.push(m); + } else { + return Err(JsError::new( + "calc_hist_non_uniform: all images must have consistent depth (f32)", + )); + } + } + histogram::calc_hist( + &f32_mats, + channels, + mask_ref, + hist_size, + &specs, + accumulate, + existing_hist_ref, + ) + .map_err(|e| JsError::new(&format!("{e}")))? + } + _ => { + return Err(JsError::new( + "calc_hist_non_uniform: unsupported image depth (must be u8 or f32)", + )); + } + }; + + Ok(Mat { + inner: DynamicMatrix { + data: DynamicData::F32(hist), + }, + }) +} + +/// Computes the back projection of a histogram using uniform ranges. +/// +/// * images - MatVector of input images (u8 or f32 depth, all matching). +/// * channels - Channel indices to back-project. +/// * hist_size - Number of bins per dimension, matching how `hist` was built. +/// `hist`'s flat length alone cannot recover its shape (e.g. an 8-bin +/// histogram could be `[8]` or `[2, 4]`), so it must be supplied explicitly. +/// * hist - Input histogram Mat (CV_32FC1). +/// * ranges - Flat array of [min_0, max_0, min_1, max_1, ...] per dimension. +/// * scale - Optional scale factor for the output back projection image. +#[wasm_bindgen(js_name = "calcBackProjectUniform")] +pub fn calc_back_project_uniform( + images: &MatVector, + channels: &[usize], + hist_size: &[usize], + hist: &Mat, + ranges: &[f32], + scale: f32, +) -> Result { + if images.inner.is_empty() { + return Err(JsError::new( + "calc_back_project_uniform: images vector must not be empty", + )); + } + + let dims = channels.len(); + if ranges.len() != dims * 2 { + return Err(JsError::new( + "ranges array must have 2 elements (min, max) per channel dimension", + )); + } + + let mut specs = Vec::with_capacity(dims); + for i in 0..dims { + specs.push(RangeSpec::Uniform(ranges[i * 2], ranges[i * 2 + 1])); + } + + let h = require_f32(hist, "calc_back_project_uniform (hist)")?; + + let first = &images.inner[0]; + let bp = match first { + DynamicMatrix { + data: DynamicData::U8(_), + } => { + let mut u8_mats = Vec::with_capacity(images.inner.len()); + for dyn_mat in &images.inner { + if let Some(m) = dyn_mat.as_matrix_u8() { + u8_mats.push(m); + } else { + return Err(JsError::new( + "calc_back_project_uniform: all images must have consistent depth (u8)", + )); + } + } + histogram::calc_back_project(&u8_mats, channels, hist_size, h, &specs, scale) + .map_err(|e| JsError::new(&format!("{e}")))? + } + DynamicMatrix { + data: DynamicData::F32(_), + } => { + let mut f32_mats = Vec::with_capacity(images.inner.len()); + for dyn_mat in &images.inner { + if let Some(m) = dyn_mat.as_matrix_f32() { + f32_mats.push(m); + } else { + return Err(JsError::new( + "calc_back_project_uniform: all images must have consistent depth (f32)", + )); + } + } + histogram::calc_back_project(&f32_mats, channels, hist_size, h, &specs, scale) + .map_err(|e| JsError::new(&format!("{e}")))? + } + _ => { + return Err(JsError::new( + "calc_back_project_uniform: unsupported image depth (must be u8 or f32)", + )); + } + }; + + Ok(Mat { + inner: DynamicMatrix { + data: DynamicData::F32(bp), + }, + }) +} + +/// Computes the back projection of a histogram using non-uniform ranges. +/// +/// * images - MatVector of input images (u8 or f32 depth, all matching). +/// * channels - Channel indices to back-project. +/// * hist_size - Number of bins per dimension, matching how `hist` was built. +/// `hist`'s flat length alone cannot recover its shape (e.g. an 8-bin +/// histogram could be `[8]` or `[2, 4]`), so it must be supplied explicitly. +/// * hist - Input histogram Mat (CV_32FC1). +/// * ranges - Array of Float32Arrays containing bin boundaries. +/// * scale - Optional scale factor for the output back projection image. +#[wasm_bindgen(js_name = "calcBackProjectNonUniform")] +pub fn calc_back_project_non_uniform( + images: &MatVector, + channels: &[usize], + hist_size: &[usize], + hist: &Mat, + ranges: js_sys::Array, + scale: f32, +) -> Result { + if images.inner.is_empty() { + return Err(JsError::new( + "calc_back_project_non_uniform: images vector must not be empty", + )); + } + + let dims = channels.len(); + if ranges.length() as usize != dims { + return Err(JsError::new( + "ranges array must have one Float32Array per channel dimension", + )); + } + + let mut specs = Vec::with_capacity(dims); + for i in 0..dims { + let js_val = ranges.get(i as u32); + let f32_array = js_sys::Float32Array::from(js_val); + let mut vec = vec![0.0f32; f32_array.length() as usize]; + f32_array.copy_to(&mut vec); + specs.push(RangeSpec::NonUniform(vec)); + } + + let h = require_f32(hist, "calc_back_project_non_uniform (hist)")?; + + let first = &images.inner[0]; + let bp = match first { + DynamicMatrix { + data: DynamicData::U8(_), + } => { + let mut u8_mats = Vec::with_capacity(images.inner.len()); + for dyn_mat in &images.inner { + if let Some(m) = dyn_mat.as_matrix_u8() { + u8_mats.push(m); + } else { + return Err(JsError::new( + "calc_back_project_non_uniform: all images must have consistent depth (u8)", + )); + } + } + histogram::calc_back_project(&u8_mats, channels, hist_size, h, &specs, scale) + .map_err(|e| JsError::new(&format!("{e}")))? + } + DynamicMatrix { + data: DynamicData::F32(_), + } => { + let mut f32_mats = Vec::with_capacity(images.inner.len()); + for dyn_mat in &images.inner { + if let Some(m) = dyn_mat.as_matrix_f32() { + f32_mats.push(m); + } else { + return Err(JsError::new( + "calc_back_project_non_uniform: all images must have consistent depth (f32)", + )); + } + } + histogram::calc_back_project(&f32_mats, channels, hist_size, h, &specs, scale) + .map_err(|e| JsError::new(&format!("{e}")))? + } + _ => { + return Err(JsError::new( + "calc_back_project_non_uniform: unsupported image depth (must be u8 or f32)", + )); + } + }; + + Ok(Mat { + inner: DynamicMatrix { + data: DynamicData::F32(bp), + }, + }) +} + +/// Compares two histograms using the specified comparison method. +/// +/// * h1 - First input histogram (CV_32FC1). +/// * h2 - Second input histogram of the same size as h1 (CV_32FC1). +/// * method - Comparison method: HIST_CMP_CORREL, HIST_CMP_CHISQR, etc. +#[wasm_bindgen(js_name = "compareHist")] +pub fn compare_hist(h1: &Mat, h2: &Mat, method: i32) -> Result { + let m1 = require_f32(h1, "compareHist (h1)")?; + let m2 = require_f32(h2, "compareHist (h2)")?; + let meth = hist_comp_method_from_i32(method)?; + + histogram::compare_hist(m1, m2, meth).map_err(|e| JsError::new(&format!("{e}"))) +} + +/// Equalizes the histogram of a grayscale image. +/// +/// **Note:** Currently only 8-bit single-channel images (CV_8UC1) are supported. +/// 16-bit (CV_16UC1) images are not yet supported. +/// +/// * src - Input single-channel 8-bit image (CV_8UC1). +#[wasm_bindgen(js_name = "equalizeHist")] +pub fn equalize_hist(src: &Mat) -> Result { + let m = require_u8(src, "equalizeHist")?; + let res = histogram::equalize_hist(m).map_err(|e| JsError::new(&format!("{e}")))?; + Ok(Mat { + inner: DynamicMatrix { + data: DynamicData::U8(res), + }, + }) +} + +// --------------------------------------------------------------------------- +// Clahe wrapper +// --------------------------------------------------------------------------- + +/// Contrast Limited Adaptive Histogram Equalization. +/// +/// Wraps OpenCV-style CLAHE. +#[wasm_bindgen(js_name = "Clahe")] +pub struct WasmClahe { + inner: CoreClahe, +} + +#[wasm_bindgen(js_class = "Clahe")] +impl WasmClahe { + /// Creates a new CLAHE instance with given clip limit and grid size. + /// + /// * clip_limit - Threshold for contrast limiting (default in OpenCV is 40.0). + /// * tile_grid_width - Number of tiles horizontally (e.g. 8). + /// * tile_grid_height - Number of tiles vertically (e.g. 8). + #[wasm_bindgen(constructor)] + pub fn new(clip_limit: f64, tile_grid_width: i32, tile_grid_height: i32) -> Self { + let size = purecv::core::types::Size2i::new(tile_grid_width, tile_grid_height); + Self { + inner: histogram::create_clahe(clip_limit, size), + } + } + + /// Applies CLAHE to the input grayscale image. + /// + /// **Note:** Currently only 8-bit single-channel images (CV_8UC1) are supported. + /// 16-bit (CV_16UC1) support is not yet exposed in WebAssembly. + pub fn apply(&self, src: &Mat) -> Result { + let m = require_u8(src, "Clahe::apply")?; + let res = self + .inner + .apply_u8(m) + .map_err(|e| JsError::new(&format!("{e}")))?; + Ok(Mat { + inner: DynamicMatrix { + data: DynamicData::U8(res), + }, + }) + } + + /// Gets the clip limit threshold. + #[wasm_bindgen(getter, js_name = "clipLimit")] + pub fn clip_limit(&self) -> f64 { + self.inner.get_clip_limit() + } + + /// Sets the clip limit threshold. + #[wasm_bindgen(setter, js_name = "clipLimit")] + pub fn set_clip_limit(&mut self, limit: f64) { + self.inner.set_clip_limit(limit); + } + + /// Gets the number of tile grid columns (width). + #[wasm_bindgen(getter, js_name = "tileGridWidth")] + pub fn tile_grid_width(&self) -> i32 { + self.inner.get_tiles_grid_size().width + } + + /// Gets the number of tile grid rows (height). + #[wasm_bindgen(getter, js_name = "tileGridHeight")] + pub fn tile_grid_height(&self) -> i32 { + self.inner.get_tiles_grid_size().height + } + + /// Sets the tile grid size (columns and rows). + #[wasm_bindgen(js_name = "setTilesGridSize")] + pub fn set_tiles_grid_size(&mut self, width: i32, height: i32) { + self.inner + .set_tiles_grid_size(purecv::core::types::Size2i::new(width, height)); + } + + /// Gets the bit shift parameter. + #[wasm_bindgen(getter, js_name = "bitShift")] + pub fn bit_shift(&self) -> i32 { + self.inner.get_bit_shift() + } + + /// Sets the bit shift parameter. + #[wasm_bindgen(setter, js_name = "bitShift")] + pub fn set_bit_shift(&mut self, bit_shift: i32) { + self.inner.set_bit_shift(bit_shift); + } +} diff --git a/crates/wasm/tests/web.rs b/crates/wasm/tests/web.rs index 0c47cd8..a7f8905 100644 --- a/crates/wasm/tests/web.rs +++ b/crates/wasm/tests/web.rs @@ -1,7 +1,9 @@ #![cfg(target_arch = "wasm32")] use purecv_wasm::{ - find_homography_wasm, rodrigues_wasm, solve_pnp_wasm, Mat, Point2fVector, Point3fVector, + calc_back_project_uniform, calc_hist_uniform, compare_hist, equalize_hist, + find_homography_wasm, hist_cmp_correl, rodrigues_wasm, solve_pnp_wasm, Mat, MatVector, + Point2fVector, Point3fVector, WasmClahe, }; use wasm_bindgen_test::*; @@ -120,3 +122,91 @@ fn test_find_homography() { ) .unwrap(); } + +#[wasm_bindgen_test] +fn test_equalize_hist() { + let data = vec![0, 50, 100, 150, 200, 255]; + let src = Mat::from_u8_data(2, 3, 1, &data).unwrap(); + let eq = equalize_hist(&src).unwrap(); + assert_eq!(eq.rows(), 2); + assert_eq!(eq.cols(), 3); +} + +#[wasm_bindgen_test] +fn test_clahe() { + let data = vec![128u8; 64]; + let src = Mat::from_u8_data(8, 8, 1, &data).unwrap(); + let mut clahe = WasmClahe::new(40.0, 8, 8); + assert_eq!(clahe.clip_limit(), 40.0); + clahe.set_clip_limit(20.0); + assert_eq!(clahe.clip_limit(), 20.0); + let dst = clahe.apply(&src).unwrap(); + assert_eq!(dst.rows(), 8); + assert_eq!(dst.cols(), 8); +} + +#[wasm_bindgen_test] +fn test_calc_hist_and_compare() { + let mut mv = MatVector::new(); + let data = vec![0u8, 1, 2, 3, 4, 5, 6, 7]; + let img = Mat::from_u8_data(2, 4, 1, &data).unwrap(); + mv.push(&img); + assert_eq!(mv.length(), 1); + + let hist_size = vec![4]; + let ranges = vec![0.0, 8.0]; + let channels = vec![0]; + let h1 = calc_hist_uniform(&mv, &channels, None, &hist_size, &ranges, false, None).unwrap(); + assert_eq!(h1.rows(), 4); + assert_eq!(h1.cols(), 1); + + let score = compare_hist(&h1, &h1, hist_cmp_correl()).unwrap(); + assert!((score - 1.0).abs() < 1e-4); + + // 8 pixel values (0..8) into 4 bins over [0,8) -> 2 pixels per bin. + assert_eq!(h1.data_f32().unwrap(), vec![2.0, 2.0, 2.0, 2.0]); + + // accumulate=true should add onto existing_hist, not start from zero. + let h2 = calc_hist_uniform(&mv, &channels, None, &hist_size, &ranges, true, Some(h1)).unwrap(); + assert_eq!(h2.data_f32().unwrap(), vec![4.0, 4.0, 4.0, 4.0]); +} + +#[wasm_bindgen_test] +fn test_calc_back_project_2d_shape() { + // Regression test: calc_back_project_uniform now requires an explicit + // hist_size, so the flat 8-bin histogram below is correctly treated as + // a [2, 4] shape rather than guessed (and previously miscomputed as + // [1, 8]) from its bin count alone. + let mut mv = MatVector::new(); + let img = Mat::from_u8_data(1, 1, 2, &[1u8, 7]).unwrap(); + mv.push(&img); + + let hist_data: Vec = vec![10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0]; + let hist = Mat::from_f32_data(8, 1, 1, &hist_data).unwrap(); + + let bp = calc_back_project_uniform(&mv, &[0, 1], &[2, 4], &hist, &[0.0, 2.0, 0.0, 8.0], 1.0) + .unwrap(); + // channel 0 (value 1, range [0,2), 2 bins) -> bin 1 + // channel 1 (value 7, range [0,8), 4 bins) -> bin 3 + // flat index with strides [4, 1] -> 1*4 + 3 = 7 -> hist.data[7] = 80.0 + assert_eq!(bp.data_f32().unwrap(), vec![80.0]); +} + +#[wasm_bindgen_test] +fn test_calc_back_project_mismatched_image_size_error() { + let mut mv = MatVector::new(); + mv.push(&Mat::from_u8_data(2, 2, 1, &[0u8, 1, 2, 3]).unwrap()); + mv.push(&Mat::from_u8_data(3, 3, 1, &[0u8; 9]).unwrap()); + + let hist = Mat::from_f32_data(4, 1, 1, &[1.0, 2.0, 3.0, 4.0]).unwrap(); + assert!(calc_back_project_uniform(&mv, &[0, 1], &[4], &hist, &[0.0, 4.0], 1.0).is_err()); +} + +#[wasm_bindgen_test] +fn test_calc_hist_multichannel_mask_error() { + let mut mv = MatVector::new(); + mv.push(&Mat::from_u8_data(2, 2, 1, &[0u8, 1, 2, 3]).unwrap()); + + let mask = Mat::from_u8_data(2, 2, 3, &[1u8; 12]).unwrap(); + assert!(calc_hist_uniform(&mv, &[0], Some(mask), &[4], &[0.0, 4.0], false, None).is_err()); +} diff --git a/crates/wasm/www/example_histogram.html b/crates/wasm/www/example_histogram.html new file mode 100644 index 0000000..b82391d --- /dev/null +++ b/crates/wasm/www/example_histogram.html @@ -0,0 +1,224 @@ + + + + + + + PureCV - Histogram & CLAHE + + + +
+

Initializing WASM...

+
+ +
+
+

Histogram & CLAHE

+

calc_hist, compare_hist, equalize_hist & Contrast Limited Adaptive Histogram Equalization

+
+ +
+
+

Click or drag image here

+ +
+ +
+
+ + +
+
+ + +
+
+
+ +
+
+ +
+
+ +
+
256-bin uniform histogram of the grayscale source (calc_hist)
+ +
+ +
+
compare_hist: source histogram vs. CLAHE-equalized histogram
+
+ +
+
+
+ + + + diff --git a/crates/wasm/www/example_histogram.js b/crates/wasm/www/example_histogram.js new file mode 100644 index 0000000..956984a --- /dev/null +++ b/crates/wasm/www/example_histogram.js @@ -0,0 +1,221 @@ +/* + * example_histogram.js + * purecv + * + * This file is part of purecv - WebARKit. + * + * purecv is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * purecv is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with purecv. If not, see . + * + * As a special exception, the copyright holders of this library give you + * permission to link this library with independent modules to produce an + * executable, regardless of the license terms of these independent modules, and to + * copy and distribute the resulting executable under terms of your choice, + * provided that you also meet, for each linked independent module, the terms and + * conditions of the license of that module. An independent module is a module + * which is neither derived from nor based on this library. If you modify this + * library, you may extend this exception to your version of the library, but you + * are not obligated to do so. If you do not wish to do so, delete this exception + * statement from your version. + * + * Copyright 2026 WebARKit. + * + * Author(s): Walter Perdan @kalwalt https://github.com/kalwalt + * + */ + +import { initWasm, loadImage, getScaledDimensions, canvasToMat, matToCanvas } from './cv_demo_utils.js'; + +let sourceImage = null; +let cv = null; +const compareRow = document.getElementById('compare-row'); +const histCanvas = document.getElementById('histogram-canvas'); +const scoresEl = document.getElementById('scores'); +const clipSlider = document.getElementById('clip-limit'); +const clipDisplay = document.getElementById('val-clip'); +const tileSlider = document.getElementById('tile-grid'); +const tilesDisplay = document.getElementById('val-tiles'); +const tilesDisplay2 = document.getElementById('val-tiles-2'); + +const COMPARE_METHODS = [ + { name: 'Correl', ctor: (cv) => cv.HIST_CMP_CORREL() }, + { name: 'ChiSqr', ctor: (cv) => cv.HIST_CMP_CHISQR() }, + { name: 'ChiSqrAlt', ctor: (cv) => cv.HIST_CMP_CHISQR_ALT() }, + { name: 'Intersect', ctor: (cv) => cv.HIST_CMP_INTERSECT() }, + { name: 'Bhattacharyya', ctor: (cv) => cv.HIST_CMP_BHATTACHARYYA() }, + { name: 'KL Divergence', ctor: (cv) => cv.HIST_CMP_KL_DIV() }, +]; + +async function start() { + try { + cv = await initWasm(); + document.getElementById('loader').classList.add('hidden'); + + sourceImage = await loadImage('https://raw.githubusercontent.com/opencv/opencv/master/samples/data/butterfly.jpg'); + processImage(); + } catch (err) { + console.error("WASM Initialization failed:", err); + document.getElementById('loader').innerHTML = `

Error loading WASM: ${err.message}

`; + } +} + +function addCanvasBox(label) { + const box = document.createElement('div'); + box.className = 'level-box'; + + const canvas = document.createElement('canvas'); + const text = document.createElement('span'); + text.className = 'level-label'; + text.innerText = label; + + box.appendChild(canvas); + box.appendChild(text); + compareRow.appendChild(box); + return canvas; +} + +function drawHistogram(canvas, histData, color) { + const w = canvas.clientWidth || 512; + const h = 160; + canvas.width = w; + canvas.height = h; + const ctx = canvas.getContext('2d'); + ctx.clearRect(0, 0, w, h); + + const max = Math.max(...histData, 1); + const binWidth = w / histData.length; + + ctx.fillStyle = color; + for (let i = 0; i < histData.length; i++) { + const barHeight = (histData[i] / max) * (h - 4); + ctx.fillRect(i * binWidth, h - barHeight, Math.max(binWidth, 1), barHeight); + } +} + +function processImage() { + if (!sourceImage || !cv) return; + + compareRow.innerHTML = ''; + scoresEl.innerHTML = ''; + + const { width, height } = getScaledDimensions(sourceImage, 512); + const tempCanvas = document.createElement('canvas'); + tempCanvas.width = width; + tempCanvas.height = height; + const tempCtx = tempCanvas.getContext('2d', "willReadFrequently: true"); + tempCtx.drawImage(sourceImage, 0, 0, width, height); + + const rgbaMat = canvasToMat(cv, tempCanvas, tempCtx); + const gray = cv.cvtColor(rgbaMat, cv.COLOR_RGBA2GRAY()); + rgbaMat.free(); + + const clipLimit = parseFloat(clipSlider.value); + const tileGrid = parseInt(tileSlider.value); + + // Declared outside the try so a mid-way failure can still free whatever + // was already allocated (WASM objects are not garbage-collected). + let images = null; + let hist = null; + let equalized = null; + let clahe = null; + let claheOut = null; + let claheImages = null; + let claheHist = null; + + try { + // --- calc_hist: 256-bin uniform histogram of the grayscale source --- + images = new cv.MatVector(); + images.push(gray); + const histSize = [256]; + const ranges = [0.0, 256.0]; + hist = cv.calcHistUniform(images, [0], undefined, histSize, ranges, false, undefined); + drawHistogram(histCanvas, hist.dataF32(), '#43e97b'); + + // --- equalize_hist: global histogram equalization --- + equalized = cv.equalizeHist(gray); + + // --- CLAHE: contrast-limited adaptive histogram equalization --- + clahe = new cv.Clahe(clipLimit, tileGrid, tileGrid); + claheOut = clahe.apply(gray); + + // --- compare_hist: source vs. CLAHE-equalized histograms --- + claheImages = new cv.MatVector(); + claheImages.push(claheOut); + claheHist = cv.calcHistUniform(claheImages, [0], undefined, histSize, ranges, false, undefined); + + for (const method of COMPARE_METHODS) { + const score = cv.compareHist(hist, claheHist, method.ctor(cv)); + const item = document.createElement('div'); + item.className = 'score-item'; + item.innerHTML = `${method.name}${score.toFixed(4)}`; + scoresEl.appendChild(item); + } + + // --- render the three grayscale variants side by side --- + matToCanvas(gray, addCanvasBox(`Grayscale: ${gray.cols}x${gray.rows}`)); + matToCanvas(equalized, addCanvasBox('equalize_hist')); + matToCanvas(claheOut, addCanvasBox(`CLAHE (clip=${clipLimit}, ${tileGrid}x${tileGrid})`)); + } catch (e) { + console.error("Histogram/CLAHE error:", e); + } finally { + gray.free(); + images?.free(); + hist?.free(); + equalized?.free(); + clahe?.free(); + claheOut?.free(); + claheImages?.free(); + claheHist?.free(); + } +} + +clipSlider.oninput = () => { + clipDisplay.innerText = clipSlider.value; + processImage(); +}; + +tileSlider.oninput = () => { + tilesDisplay.innerText = tileSlider.value; + tilesDisplay2.innerText = tileSlider.value; + processImage(); +}; + +const fileInput = document.getElementById('file-input'); +const dropZone = document.getElementById('drop-zone'); + +async function loadFile(file) { + if (!file) return; + const url = URL.createObjectURL(file); + try { + sourceImage = await loadImage(url); + processImage(); + } finally { + URL.revokeObjectURL(url); + } +} + +dropZone.onclick = () => fileInput.click(); +fileInput.onchange = (e) => loadFile(e.target.files[0]); + +dropZone.ondragover = (e) => { + e.preventDefault(); + dropZone.classList.add('drag-over'); +}; +dropZone.ondragleave = () => dropZone.classList.remove('drag-over'); +dropZone.ondrop = (e) => { + e.preventDefault(); + dropZone.classList.remove('drag-over'); + loadFile(e.dataTransfer.files[0]); +}; + +start(); diff --git a/crates/wasm/www/index.html b/crates/wasm/www/index.html index b52b2c8..db72be8 100644 --- a/crates/wasm/www/index.html +++ b/crates/wasm/www/index.html @@ -124,6 +124,11 @@

ORB Keypoints & BRIEF

Image Pyramids

Construct and visualize a Gaussian pyramid stack at multiple scales.

+ + Imgproc +

Histogram & CLAHE

+

Compute histograms, compare distributions, and enhance contrast with global and adaptive (CLAHE) equalization.

+ Hough

Line Detection

diff --git a/examples/histogram.rs b/examples/histogram.rs new file mode 100644 index 0000000..fdc9f28 --- /dev/null +++ b/examples/histogram.rs @@ -0,0 +1,205 @@ +/* + * histogram.rs + * purecv + * + * This file is part of purecv - WebARKit. + * + * purecv is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * purecv is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with purecv. If not, see . + * + * As a special exception, the copyright holders of this library give you + * permission to link this library with independent modules to produce an + * executable, regardless of the license terms of these independent modules, and to + * copy and distribute the resulting executable under terms of your choice, + * provided that you also meet, for each linked independent module, the terms and + * conditions of the license of that module. An independent module is a module + * which is neither derived from nor based on this library. If you modify this + * library, you may extend this exception to your version of the library, but you + * are not obligated to do so. If you do not wish to do so, delete this exception + * statement from your version. + * + * Copyright 2026 WebARKit. + * + * Author(s): Walter Perdan @kalwalt https://github.com/kalwalt + * + */ + +use image::{DynamicImage, GenericImageView, ImageBuffer, Luma}; +use purecv::core::logging::tags; +use purecv::core::{Matrix, Size2i}; +use purecv::imgproc::{ + calc_back_project, calc_hist, compare_hist, create_clahe, cvt_color_rgb_to_gray, equalize_hist, + HistCompMethods, RangeSpec, +}; +use purecv::version; +use purecv::{cv_log_error, cv_log_info}; +use std::path::Path; + +fn main() -> Result<(), Box> { + purecv::core::logging::init_basic_logger()?; + + cv_log_info!(tags::PURECV, "--- Histogram & CLAHE Example ---"); + version::print_version(); + + // 1. Load the image + let img_path = "examples/data/butterfly.jpg"; + if !Path::new(img_path).exists() { + cv_log_error!( + tags::IMGPROC, + "{} not found. Run from the project root.", + img_path + ); + return Ok(()); + } + + let img = image::open(img_path)?; + let (width, height) = img.dimensions(); + cv_log_info!( + tags::IMGPROC, + "loaded image: {} ({}x{})", + img_path, + width, + height + ); + + let rgb_img = img.to_rgb8(); + let mat_rgb = Matrix::from_vec(height as usize, width as usize, 3, rgb_img.into_raw()); + let mat_gray = cvt_color_rgb_to_gray(&mat_rgb)?; + + // Create output directory if it doesn't exist. + std::fs::create_dir_all("examples/data/out")?; + + // --- calc_hist: uniform bins --- + + cv_log_info!(tags::IMGPROC, "computing a 256-bin uniform histogram..."); + let hist_size = [256usize]; + let ranges = [RangeSpec::Uniform(0.0, 256.0)]; + let hist = calc_hist(&[&mat_gray], &[0], None, &hist_size, &ranges, false, None)?; + print_histogram_summary(&hist); + + // --- calc_hist: non-uniform bins --- + // + // Four unevenly-spaced bins: a wide shadow bucket, two mid-tone buckets, + // and a narrow highlight bucket. + cv_log_info!(tags::IMGPROC, "computing a 4-bin non-uniform histogram..."); + let non_uniform_size = [4usize]; + let boundaries = vec![0.0, 96.0, 160.0, 224.0, 256.0]; + let non_uniform_ranges = [RangeSpec::NonUniform(boundaries)]; + let non_uniform_hist = calc_hist( + &[&mat_gray], + &[0], + None, + &non_uniform_size, + &non_uniform_ranges, + false, + None, + )?; + println!( + " non-uniform bins [0,96) [96,160) [160,224) [224,256): {:?}", + non_uniform_hist.data + ); + + // --- calc_back_project --- + // + // Projects the histogram back onto the source image: each pixel is + // replaced by its own bin's count, scaled for visibility. Bright + // regions in the output correspond to common intensities in the image. + // The scale is derived from the histogram's own peak so the output uses + // the full 0-255 display range, matching OpenCV's typical demo pattern. + cv_log_info!(tags::IMGPROC, "back-projecting the histogram..."); + let max_bin = hist.data.iter().cloned().fold(0.0f32, f32::max); + let bp_scale = if max_bin > 0.0 { 255.0 / max_bin } else { 1.0 }; + let back_projected = + calc_back_project(&[&mat_gray], &[0], &hist_size, &hist, &ranges, bp_scale)?; + let mut bp_u8 = Matrix::::new(back_projected.rows, back_projected.cols, 1); + for (dst, &src) in bp_u8.data.iter_mut().zip(back_projected.data.iter()) { + *dst = src.clamp(0.0, 255.0) as u8; + } + save_matrix_gray(&bp_u8, "examples/data/out/output_back_project.png")?; + + // --- compare_hist --- + // + // Compare the source histogram against the histogram of a brightened + // copy of the same image, across every HistCompMethods variant. + cv_log_info!(tags::IMGPROC, "comparing histograms..."); + let mut brightened = Matrix::::new(mat_gray.rows, mat_gray.cols, 1); + for (dst, &src) in brightened.data.iter_mut().zip(mat_gray.data.iter()) { + *dst = src.saturating_add(40); + } + let brightened_hist = calc_hist(&[&brightened], &[0], None, &hist_size, &ranges, false, None)?; + + for method in [ + HistCompMethods::Correl, + HistCompMethods::ChiSqr, + HistCompMethods::ChiSqrAlt, + HistCompMethods::Intersection, + HistCompMethods::Bhattacharyya, + HistCompMethods::KullbackLeibler, + ] { + let score = compare_hist(&hist, &brightened_hist, method)?; + println!(" {:?}: {:.4}", method, score); + } + + // --- equalize_hist --- + + cv_log_info!(tags::IMGPROC, "equalizing histogram..."); + let equalized = equalize_hist(&mat_gray)?; + save_matrix_gray(&equalized, "examples/data/out/output_equalize_hist.png")?; + + // --- CLAHE --- + // + // Contrast Limited Adaptive Histogram Equalization avoids over-amplifying + // noise in flat regions, unlike global equalize_hist above. Two tile + // grid sizes are compared: a coarser 4x4 grid and a finer 8x8 grid. + cv_log_info!(tags::IMGPROC, "applying CLAHE (4x4 tiles)..."); + let clahe_coarse = create_clahe(40.0, Size2i::new(4, 4)); + let clahe_coarse_out = clahe_coarse.apply_u8(&mat_gray)?; + save_matrix_gray(&clahe_coarse_out, "examples/data/out/output_clahe_4x4.png")?; + + cv_log_info!(tags::IMGPROC, "applying CLAHE (8x8 tiles)..."); + let clahe_fine = create_clahe(40.0, Size2i::new(8, 8)); + let clahe_fine_out = clahe_fine.apply_u8(&mat_gray)?; + save_matrix_gray(&clahe_fine_out, "examples/data/out/output_clahe_8x8.png")?; + + cv_log_info!( + tags::IMGPROC, + "done! Check the output_*.png files under examples/data/out/." + ); + Ok(()) +} + +/// Prints the min/max/mean bin count and the 3 most populated bins of a +/// flattened `f32` histogram (as returned by `calc_hist`). +fn print_histogram_summary(hist: &Matrix) { + let data = &hist.data; + let total_bins = data.len(); + let sum: f32 = data.iter().sum(); + let mean = sum / total_bins as f32; + let max_bin = data + .iter() + .enumerate() + .max_by(|a, b| a.1.partial_cmp(b.1).unwrap()) + .map(|(idx, &v)| (idx, v)) + .unwrap_or((0, 0.0)); + println!( + " {total_bins} bins, total count {sum:.0}, mean {mean:.2}, most populated bin {} (count {:.0})", + max_bin.0, max_bin.1 + ); +} + +fn save_matrix_gray(mat: &Matrix, filename: &str) -> image::ImageResult<()> { + let img: ImageBuffer, Vec> = + ImageBuffer::from_raw(mat.cols as u32, mat.rows as u32, mat.data.clone()) + .expect("Failed to create image buffer"); + DynamicImage::ImageLuma8(img).save(filename) +} diff --git a/package.json b/package.json index dada5ed..a9f0e00 100644 --- a/package.json +++ b/package.json @@ -1,6 +1,6 @@ { "name": "purecv", - "version": "0.7.1", + "version": "0.8.0", "description": "A pure Rust, high-performance computer vision library focused on safety and portability.", "private": true, "scripts": { diff --git a/src/imgproc.rs b/src/imgproc.rs index 6475188..e664829 100644 --- a/src/imgproc.rs +++ b/src/imgproc.rs @@ -57,6 +57,7 @@ pub mod edge; pub mod feature; pub mod filter; pub mod geometric; +pub mod histogram; pub mod hough; pub mod morph; pub mod pyramid; @@ -79,6 +80,10 @@ pub use edge::*; pub use feature::*; pub use filter::*; pub use geometric::{remap, warp_perspective, InterpolationFlags}; +pub use histogram::{ + calc_back_project, calc_hist, compare_hist, create_clahe, equalize_hist, Clahe, + HistCompMethods, RangeSpec, +}; pub use hough::*; pub use morph::{dilate, erode, get_structuring_element, morphology_ex, MorphShapes, MorphTypes}; pub use pyramid::{build_pyramid, pyr_down, pyr_up}; diff --git a/src/imgproc/color.rs b/src/imgproc/color.rs index 650bd23..0d60ac0 100644 --- a/src/imgproc/color.rs +++ b/src/imgproc/color.rs @@ -61,7 +61,7 @@ fn rgb_to_gray_row(out_row: &mut [u8], in_row: &[u8]) { } #[cfg(not(feature = "simd"))] { - for (out_pixel, in_val) in out_row.iter_mut().zip(in_row.chunks_exact(3)) { + for (out_pixel, in_val) in out_row.iter_mut().zip(in_row.as_chunks::<3>().0.iter()) { let r = in_val[0] as f32; let g = in_val[1] as f32; let b = in_val[2] as f32; @@ -79,7 +79,7 @@ fn bgr_to_gray_row(out_row: &mut [u8], in_row: &[u8]) { } #[cfg(not(feature = "simd"))] { - for (out_pixel, in_val) in out_row.iter_mut().zip(in_row.chunks_exact(3)) { + for (out_pixel, in_val) in out_row.iter_mut().zip(in_row.as_chunks::<3>().0.iter()) { let b = in_val[0] as f32; let g = in_val[1] as f32; let r = in_val[2] as f32; @@ -97,7 +97,7 @@ fn rgba_to_gray_row(out_row: &mut [u8], in_row: &[u8]) { } #[cfg(not(feature = "simd"))] { - for (out_pixel, in_val) in out_row.iter_mut().zip(in_row.chunks_exact(4)) { + for (out_pixel, in_val) in out_row.iter_mut().zip(in_row.as_chunks::<4>().0.iter()) { let r = in_val[0] as f32; let g = in_val[1] as f32; let b = in_val[2] as f32; @@ -115,7 +115,7 @@ fn bgra_to_gray_row(out_row: &mut [u8], in_row: &[u8]) { } #[cfg(not(feature = "simd"))] { - for (out_pixel, in_val) in out_row.iter_mut().zip(in_row.chunks_exact(4)) { + for (out_pixel, in_val) in out_row.iter_mut().zip(in_row.as_chunks::<4>().0.iter()) { let b = in_val[0] as f32; let g = in_val[1] as f32; let r = in_val[2] as f32; @@ -347,7 +347,9 @@ pub fn cvt_color_gray_to_rgb(input: &Matrix) -> Result, &'static .par_chunks_exact_mut(out_row_len) .zip(input.data.par_chunks_exact(in_row_len)) .for_each(|(out_row, in_row)| { - for (out_val, in_pixel) in out_row.chunks_exact_mut(3).zip(in_row.iter()) { + for (out_val, in_pixel) in + out_row.as_chunks_mut::<3>().0.iter_mut().zip(in_row.iter()) + { let v = *in_pixel; out_val[0] = v; out_val[1] = v; @@ -363,7 +365,9 @@ pub fn cvt_color_gray_to_rgb(input: &Matrix) -> Result, &'static .chunks_exact_mut(out_row_len) .zip(input.data.chunks_exact(in_row_len)) .for_each(|(out_row, in_row)| { - for (out_val, in_pixel) in out_row.chunks_exact_mut(3).zip(in_row.iter()) { + for (out_val, in_pixel) in + out_row.as_chunks_mut::<3>().0.iter_mut().zip(in_row.iter()) + { let v = *in_pixel; out_val[0] = v; out_val[1] = v; @@ -392,7 +396,9 @@ pub fn cvt_color_gray_to_bgr(input: &Matrix) -> Result, &'static .par_chunks_exact_mut(out_row_len) .zip(input.data.par_chunks_exact(in_row_len)) .for_each(|(out_row, in_row)| { - for (out_val, in_pixel) in out_row.chunks_exact_mut(3).zip(in_row.iter()) { + for (out_val, in_pixel) in + out_row.as_chunks_mut::<3>().0.iter_mut().zip(in_row.iter()) + { let v = *in_pixel; out_val[0] = v; out_val[1] = v; @@ -408,7 +414,9 @@ pub fn cvt_color_gray_to_bgr(input: &Matrix) -> Result, &'static .chunks_exact_mut(out_row_len) .zip(input.data.chunks_exact(in_row_len)) .for_each(|(out_row, in_row)| { - for (out_val, in_pixel) in out_row.chunks_exact_mut(3).zip(in_row.iter()) { + for (out_val, in_pixel) in + out_row.as_chunks_mut::<3>().0.iter_mut().zip(in_row.iter()) + { let v = *in_pixel; out_val[0] = v; out_val[1] = v; @@ -437,7 +445,9 @@ pub fn cvt_color_gray_to_rgba(input: &Matrix) -> Result, &'static .par_chunks_exact_mut(out_row_len) .zip(input.data.par_chunks_exact(in_row_len)) .for_each(|(out_row, in_row)| { - for (out_val, in_pixel) in out_row.chunks_exact_mut(4).zip(in_row.iter()) { + for (out_val, in_pixel) in + out_row.as_chunks_mut::<4>().0.iter_mut().zip(in_row.iter()) + { let v = *in_pixel; out_val[0] = v; out_val[1] = v; @@ -454,7 +464,9 @@ pub fn cvt_color_gray_to_rgba(input: &Matrix) -> Result, &'static .chunks_exact_mut(out_row_len) .zip(input.data.chunks_exact(in_row_len)) .for_each(|(out_row, in_row)| { - for (out_val, in_pixel) in out_row.chunks_exact_mut(4).zip(in_row.iter()) { + for (out_val, in_pixel) in + out_row.as_chunks_mut::<4>().0.iter_mut().zip(in_row.iter()) + { let v = *in_pixel; out_val[0] = v; out_val[1] = v; @@ -484,7 +496,9 @@ pub fn cvt_color_gray_to_bgra(input: &Matrix) -> Result, &'static .par_chunks_exact_mut(out_row_len) .zip(input.data.par_chunks_exact(in_row_len)) .for_each(|(out_row, in_row)| { - for (out_val, in_pixel) in out_row.chunks_exact_mut(4).zip(in_row.iter()) { + for (out_val, in_pixel) in + out_row.as_chunks_mut::<4>().0.iter_mut().zip(in_row.iter()) + { let v = *in_pixel; out_val[0] = v; out_val[1] = v; @@ -501,7 +515,9 @@ pub fn cvt_color_gray_to_bgra(input: &Matrix) -> Result, &'static .chunks_exact_mut(out_row_len) .zip(input.data.chunks_exact(in_row_len)) .for_each(|(out_row, in_row)| { - for (out_val, in_pixel) in out_row.chunks_exact_mut(4).zip(in_row.iter()) { + for (out_val, in_pixel) in + out_row.as_chunks_mut::<4>().0.iter_mut().zip(in_row.iter()) + { let v = *in_pixel; out_val[0] = v; out_val[1] = v; diff --git a/src/imgproc/histogram.rs b/src/imgproc/histogram.rs new file mode 100644 index 0000000..0fc2e7a --- /dev/null +++ b/src/imgproc/histogram.rs @@ -0,0 +1,1393 @@ +/* + * histogram.rs + * purecv + * + * This file is part of purecv - WebARKit. + * + * purecv is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * purecv is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with purecv. If not, see . + * + * As a special exception, the copyright holders of this library give you + * permission to link this library with independent modules to produce an + * executable, regardless of the license terms of these independent modules, and to + * copy and distribute the resulting executable under terms of your choice, + * provided that you also meet, for each linked independent module, the terms and + * conditions of the license of that module. An independent module is a module + * which is neither derived from nor based on this library. If you modify this + * library, you may extend this exception to your version of the library, but you + * are not obligated to do so. If you do not wish to do so, delete this exception + * statement from your version. + * + * Copyright 2026 WebARKit. + * + * Author(s): XiaoPengYouCode https://github.com/XiaoPengYouCode + * + */ + +use core::cmp::Ordering; + +use alloc::{vec, vec::Vec}; +#[allow(unused_imports)] +use num_traits::Float; +use num_traits::ToPrimitive; + +use crate::core::error::Result; +use crate::core::logging::tags; +use crate::core::types::{BorderTypes, Size2i}; +use crate::core::utils::border_interpolate; +use crate::core::Matrix; +use crate::cv_bail; +use crate::cv_err; + +#[cfg(feature = "parallel")] +use rayon::prelude::*; + +// --------------------------------------------------------------------------- +// Enums +// --------------------------------------------------------------------------- + +/// Comparison method for `compare_hist`. +/// +/// Mirrors OpenCV's `HistCompMethods`. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum HistCompMethods { + Correl = 0, + ChiSqr = 1, + ChiSqrAlt = 2, + Intersection = 3, + Bhattacharyya = 4, + KullbackLeibler = 5, +} + +/// The range specification for a single histogram dimension. +#[derive(Debug, Clone)] +pub enum RangeSpec { + Uniform(f32, f32), + NonUniform(Vec), +} + +// --------------------------------------------------------------------------- +// Multi-image channel resolution +// --------------------------------------------------------------------------- + +/// Resolves a global channel index across multiple images into +/// `(image_index, local_channel_index)`. +/// +/// OpenCV semantics: channels are numbered sequentially across images. +fn resolve_channel(global_ch: usize, images: &[&Matrix]) -> Result<(usize, usize)> { + let mut remaining = global_ch; + for (img_idx, img) in images.iter().enumerate() { + if remaining < img.channels { + return Ok((img_idx, remaining)); + } + remaining -= img.channels; + } + cv_bail!( + tags::IMGPROC, + InvalidInput, + "resolve_channel: global channel {} exceeds total channels {}", + global_ch, + images.iter().map(|i| i.channels).sum::() + ); +} + +/// Reads a single pixel value from multi-image channel indexing and converts to `f32`. +#[inline(always)] +fn read_pixel_f32( + images: &[&Matrix], + global_ch: usize, + y: usize, + x: usize, +) -> Option { + let (img_idx, local_ch) = resolve_channel(global_ch, images).ok()?; + images[img_idx].get(y, x, local_ch)?.to_f32() +} + +// --------------------------------------------------------------------------- +// Bin mapping helpers +// --------------------------------------------------------------------------- + +#[inline(always)] +fn uniform_bin(val: f32, range_lo: f32, range_hi: f32, hist_size: usize) -> Option { + if hist_size == 0 { + return None; + } + if val < range_lo || val >= range_hi { + return None; + } + let idx = ((val - range_lo) * (hist_size as f32) / (range_hi - range_lo)) as i32; + Some(idx.clamp(0, hist_size as i32 - 1) as usize) +} + +#[inline(always)] +fn nonuniform_bin(val: f32, boundaries: &[f32], hist_size: usize) -> Option { + if hist_size == 0 || boundaries.len() <= hist_size { + return None; + } + if val < boundaries[0] || val >= boundaries[hist_size] { + return None; + } + // Find the first boundary that is strictly greater than val, + // then subtract 1 to get the bin index. + match boundaries[..=hist_size] + .binary_search_by(|b| b.partial_cmp(&val).unwrap_or(Ordering::Less)) + { + Ok(idx) => { + // val == boundaries[idx]: bin is idx (left-inclusive) + // bin i = [boundaries[i], boundaries[i+1]) + // If val == b1, it belongs to bin 1: [b1, b2) + // binary_search finds idx where b[idx] == val + // So bin = idx. But if idx == 0 and val == b0, bin = 0. + if idx <= hist_size { + Some(idx) + } else { + None + } + } + Err(idx) => { + // boundaries[idx-1] <= val < boundaries[idx] + if idx > 0 && idx <= hist_size { + Some(idx - 1) + } else { + None + } + } + } +} + +fn map_bin(val: f32, range: &RangeSpec, hist_size: usize) -> Option { + match range { + RangeSpec::Uniform(lo, hi) => uniform_bin(val, *lo, *hi, hist_size), + RangeSpec::NonUniform(boundaries) => nonuniform_bin(val, boundaries, hist_size), + } +} + +/// Validates that `hist_size` and `ranges` describe a well-formed set of +/// histogram dimensions, shared by `calc_hist` and `calc_back_project` so +/// the two can't drift: `hist_size[d] > 0`, `Uniform(lo, hi)` satisfies +/// `lo < hi` (also rejects NaN, since `partial_cmp` returns `None` for it), +/// and `NonUniform(boundaries)` has exactly `hist_size[d] + 1` strictly +/// increasing entries. +/// +/// `fn_name` is used as the error-message prefix so callers keep their own +/// identity in the message (e.g. `"calc_hist: ..."` vs +/// `"calc_back_project: ..."`). +fn validate_hist_ranges(fn_name: &str, hist_size: &[usize], ranges: &[RangeSpec]) -> Result<()> { + for (d, &sz) in hist_size.iter().enumerate() { + if sz == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "{}: hist_size[{}] must be > 0 (got 0)", + fn_name, + d + ); + } + match &ranges[d] { + RangeSpec::Uniform(lo, hi) => { + if lo.partial_cmp(hi) != Some(Ordering::Less) { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "{}: Uniform range [{}] must satisfy lo < hi (got {} >= {})", + fn_name, + d, + lo, + hi + ); + } + } + RangeSpec::NonUniform(boundaries) => { + let expected_len = sz.checked_add(1).ok_or_else(|| { + cv_err!( + tags::IMGPROC, + InvalidInput, + "{}: hist_size[{}] ({}) is too large (hist_size[{}]+1 overflows)", + fn_name, + d, + sz, + d + ) + })?; + if boundaries.len() != expected_len { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "{}: NonUniform boundaries[{}] length {} must be hist_size[{}]+1 ({})", + fn_name, + d, + boundaries.len(), + d, + expected_len + ); + } + for k in 0..boundaries.len() - 1 { + if boundaries[k].partial_cmp(&boundaries[k + 1]) != Some(Ordering::Less) { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "{}: NonUniform boundaries[{}][{}] ({}) must be < boundaries[{}][{}] ({})", + fn_name, + d, + k, + boundaries[k], + d, + k + 1, + boundaries[k + 1] + ); + } + } + } + } + } + Ok(()) +} + +// --------------------------------------------------------------------------- +// calc_hist +// --------------------------------------------------------------------------- + +/// Calculates a multi-dimensional histogram of a set of images. +/// +/// * `images` - Slice of input images (may have multiple channels). +/// * `channels` - Global channel indices used for histogram computation. +/// Channel numbering spans across images: if `images[0]` has 3 channels +/// and `images[1]` has 1 channel, then channel 3 refers to `images[1]`'s channel 0. +/// * `mask` - Optional mask. Must have the same size as `images[0]`. +/// * `hist_size` - Number of bins per dimension. +/// * `ranges` - One `RangeSpec` per dimension. `Uniform(lo, hi)` for uniform bins, +/// `NonUniform(boundaries)` for explicit boundaries. +/// * `accumulate` - If true, adds to the returned histogram rather than clearing. +/// On first call, pass `None` for `hist` or create a zero histogram. +/// +/// Returns a flattened `Matrix` histogram of size `(product(hist_size), 1, 1)`. +pub fn calc_hist( + images: &[&Matrix], + channels: &[usize], + mask: Option<&Matrix>, + hist_size: &[usize], + ranges: &[RangeSpec], + accumulate: bool, + hist: Option<&Matrix>, +) -> Result> { + let dims = hist_size.len(); + if dims == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_hist: hist_size must not be empty" + ); + } + if images.is_empty() { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_hist: images must not be empty" + ); + } + if channels.len() != dims { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_hist: channels length ({}) must match hist_size length ({})", + channels.len(), + dims + ); + } + if ranges.len() != dims { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_hist: ranges length ({}) must match hist_size length ({})", + ranges.len(), + dims + ); + } + + let rows = images[0].rows; + let cols = images[0].cols; + + for img in images.iter() { + if img.rows != rows || img.cols != cols { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_hist: all images must have the same size" + ); + } + } + + if let Some(m) = mask { + if m.rows != rows || m.cols != cols { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_hist: mask must have the same size as images" + ); + } + if m.channels != 1 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_hist: mask must be single-channel (got {})", + m.channels + ); + } + } + + // Validate channel indices + for &ch in channels { + let _ = resolve_channel(ch, images)?; + } + + // Validate hist_size and ranges (prevents panics in bin mapping) + validate_hist_ranges("calc_hist", hist_size, ranges)?; + + let total_bins: usize = hist_size.iter().product(); + + // Validate accumulate histogram size to prevent OOB + if accumulate { + if let Some(h) = hist { + let hist_len = h.data.len(); + if hist_len != total_bins { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_hist: accumulate hist length {} does not match product(hist_size) {}", + hist_len, + total_bins + ); + } + } + } + + let mut hist_data = if accumulate { + hist.map(|h| h.data.clone()) + .unwrap_or_else(|| vec![0.0f32; total_bins]) + } else { + vec![0.0f32; total_bins] + }; + + // Strides (row-major, last dim varies fastest) + let mut strides = vec![1usize; dims]; + for i in (0..dims - 1).rev() { + strides[i] = strides[i + 1] * hist_size[i + 1]; + } + + let accumulate_row = |local: &mut [f32], y: usize| { + for x in 0..cols { + if let Some(m) = mask { + if let Some(&v) = m.get(y, x, 0) { + if v == 0 { + continue; + } + } + } + + let mut bin_idx = 0usize; + let mut out_of_range = false; + + for d in 0..dims { + let val = read_pixel_f32(images, channels[d], y, x).unwrap_or(0.0); + match map_bin(val, &ranges[d], hist_size[d]) { + Some(b) => bin_idx += b * strides[d], + None => { + out_of_range = true; + break; + } + } + } + + if !out_of_range { + local[bin_idx] += 1.0; + } + } + }; + + #[cfg(feature = "parallel")] + { + let partial = (0..rows) + .into_par_iter() + .fold( + || vec![0.0f32; total_bins], + |mut local, y| { + accumulate_row(&mut local, y); + local + }, + ) + .reduce( + || vec![0.0f32; total_bins], + |mut a, b| { + for (av, bv) in a.iter_mut().zip(b.iter()) { + *av += bv; + } + a + }, + ); + for (hv, pv) in hist_data.iter_mut().zip(partial.iter()) { + *hv += pv; + } + } + + #[cfg(not(feature = "parallel"))] + { + for y in 0..rows { + accumulate_row(&mut hist_data, y); + } + } + + Ok(Matrix::from_vec(total_bins, 1, 1, hist_data)) +} + +// --------------------------------------------------------------------------- +// calc_back_project +// --------------------------------------------------------------------------- + +/// Calculates the back projection of a histogram. +/// +/// * `images` - Slice of input images. +/// * `channels` - Global channel indices (same semantics as `calc_hist`). +/// * `hist_size` - Number of bins per dimension, matching how `hist` was built +/// by `calc_hist`. `hist`'s flattened bin count alone cannot recover this +/// shape unambiguously (e.g. 8 bins could be `[8]`, `[2, 4]`, `[4, 2]`, ...), +/// so it must be supplied explicitly rather than guessed. +/// * `hist` - Input histogram (`f32`, flattened multi-dimensional, size `product(hist_size)`). +/// * `ranges` - One `RangeSpec` per dimension. +/// * `scale` - Scale factor for the output values. +/// +/// Returns a single-channel `Matrix` of the same size as `images[0]`. +/// Values are scaled and clamped to `[0.0, 255.0]` matching OpenCV's u8 output range. +pub fn calc_back_project( + images: &[&Matrix], + channels: &[usize], + hist_size: &[usize], + hist: &Matrix, + ranges: &[RangeSpec], + scale: f32, +) -> Result> { + let dims = channels.len(); + if dims == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_back_project: channels must not be empty" + ); + } + if images.is_empty() { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_back_project: images must not be empty" + ); + } + if ranges.len() != dims { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_back_project: ranges length ({}) must match channels length ({})", + ranges.len(), + dims + ); + } + if hist_size.len() != dims { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_back_project: hist_size length ({}) must match channels length ({})", + hist_size.len(), + dims + ); + } + validate_hist_ranges("calc_back_project", hist_size, ranges)?; + + let rows = images[0].rows; + let cols = images[0].cols; + + for img in images.iter() { + if img.rows != rows || img.cols != cols { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_back_project: all images must have the same size" + ); + } + } + + // Validate channel indices + for &ch in channels { + let _ = resolve_channel(ch, images)?; + } + + let total_bins: usize = hist_size.iter().product(); + if total_bins == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_back_project: hist_size must not contain a zero dimension" + ); + } + let hist_len = hist.rows * hist.cols * hist.channels; + if hist_len != total_bins { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "calc_back_project: hist length {} does not match product(hist_size) {}", + hist_len, + total_bins + ); + } + + let mut strides = vec![1usize; dims]; + for i in (0..dims - 1).rev() { + strides[i] = strides[i + 1] * hist_size[i + 1]; + } + + let mut dst = Matrix::::new(rows, cols, 1); + + // `chunks_mut`/`par_chunks_mut` panic on a zero chunk size regardless of + // slice length, so a zero-width input must bail out before reaching them. + if cols == 0 { + return Ok(dst); + } + + let process_row = |y: usize, dst_row: &mut [f32]| { + for (x, out_pixel) in dst_row.iter_mut().enumerate() { + let mut bin_idx = 0usize; + let mut out_of_range = false; + + for d in 0..dims { + let val = read_pixel_f32(images, channels[d], y, x).unwrap_or(0.0); + match map_bin(val, &ranges[d], hist_size[d]) { + Some(b) => bin_idx += b * strides[d], + None => { + out_of_range = true; + break; + } + } + } + + *out_pixel = if out_of_range { + 0.0f32 + } else { + (hist.data[bin_idx] * scale).clamp(0.0, 255.0) + }; + } + }; + + #[cfg(feature = "parallel")] + { + dst.data + .par_chunks_mut(cols) + .enumerate() + .for_each(|(y, dst_row)| { + process_row(y, dst_row); + }); + } + + #[cfg(not(feature = "parallel"))] + { + for (y, dst_row) in dst.data.chunks_mut(cols).enumerate() { + process_row(y, dst_row); + } + } + + Ok(dst) +} + +// --------------------------------------------------------------------------- +// compare_hist +// --------------------------------------------------------------------------- + +/// Compares two dense histograms using the specified method. +/// +/// Both histograms must be single-channel `f32` with the same size. +pub fn compare_hist(h1: &Matrix, h2: &Matrix, method: HistCompMethods) -> Result { + if h1.channels != 1 || h2.channels != 1 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "compare_hist: histograms must be single-channel" + ); + } + + let len = h1.data.len(); + if len != h2.data.len() { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "compare_hist: histograms must have the same size ({} vs {})", + len, + h2.data.len() + ); + } + if len == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "compare_hist: histograms must not be empty" + ); + } + + #[cfg(feature = "simd")] + { + if let Some(res) = + crate::imgproc::simd::simd_compare_hist_f32(&h1.data, &h2.data, method as u8) + { + return Ok(res); + } + } + + let n = len as f64; + + match method { + HistCompMethods::Correl => { + let mut s1 = 0.0f64; + let mut s2 = 0.0f64; + let mut s11 = 0.0f64; + let mut s12 = 0.0f64; + let mut s22 = 0.0f64; + + for i in 0..len { + let a = h1.data[i] as f64; + let b = h2.data[i] as f64; + s1 += a; + s2 += b; + s11 += a * a; + s22 += b * b; + s12 += a * b; + } + + let scale = 1.0 / n; + let num = s12 - s1 * s2 * scale; + let denom2 = (s11 - s1 * s1 * scale) * (s22 - s2 * s2 * scale); + Ok(if denom2.abs() > f64::EPSILON { + num / denom2.sqrt() + } else { + 1.0 + }) + } + HistCompMethods::ChiSqr => { + let mut result = 0.0f64; + for i in 0..len { + let a = h1.data[i] as f64; + let b = h2.data[i] as f64; + if a.abs() > f64::EPSILON { + let diff = a - b; + result += diff * diff / a; + } + } + Ok(result) + } + HistCompMethods::ChiSqrAlt => { + let mut result = 0.0f64; + for i in 0..len { + let a = h1.data[i] as f64; + let b = h2.data[i] as f64; + let sum = a + b; + if sum.abs() > f64::EPSILON { + let diff = a - b; + result += diff * diff / sum; + } + } + Ok(result * 2.0) + } + HistCompMethods::Intersection => { + let mut result = 0.0f64; + for i in 0..len { + let a = h1.data[i] as f64; + let b = h2.data[i] as f64; + result += a.min(b); + } + Ok(result) + } + HistCompMethods::Bhattacharyya => { + let mut s1 = 0.0f64; + let mut s2 = 0.0f64; + let mut bc = 0.0f64; + + for i in 0..len { + let a = h1.data[i] as f64; + let b = h2.data[i] as f64; + s1 += a; + s2 += b; + bc += (a * b).sqrt(); + } + + let norm = s1 * s2; + let norm_factor = if norm.abs() > f64::EPSILON { + 1.0 / norm.sqrt() + } else { + 1.0 + }; + + Ok(((1.0 - bc * norm_factor).max(0.0)).sqrt()) + } + HistCompMethods::KullbackLeibler => { + let mut result = 0.0f64; + for i in 0..len { + let p = h1.data[i] as f64; + let q = h2.data[i] as f64; + if p.abs() > f64::EPSILON { + let q_adj = if q.abs() <= f64::EPSILON { 1e-10 } else { q }; + result += p * (p / q_adj).ln(); + } + } + Ok(result) + } + } +} + +// --------------------------------------------------------------------------- +// equalize_hist +// --------------------------------------------------------------------------- + +/// Equalizes the histogram of a grayscale image. +/// +/// * `src` - Source 8-bit single-channel image. +pub fn equalize_hist(src: &Matrix) -> Result> { + if src.channels != 1 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "equalize_hist: source must be single-channel (got {})", + src.channels + ); + } + + const HIST_SZ: usize = 256; + let mut hist = [0u32; HIST_SZ]; + + for &val in src.data.iter() { + hist[val as usize] += 1; + } + + let mut i = 0; + while i < HIST_SZ && hist[i] == 0 { + i += 1; + } + + if i == HIST_SZ { + return Ok(Matrix::new(src.rows, src.cols, 1)); + } + + let total = src.rows * src.cols; + if hist[i] == total as u32 { + let mut dst = Matrix::::new(src.rows, src.cols, 1); + for pixel in dst.data.iter_mut() { + *pixel = i as u8; + } + return Ok(dst); + } + + let scale = (HIST_SZ as f64 - 1.0) / (total as f64 - hist[i] as f64); + let mut lut = [0u8; HIST_SZ]; + let mut sum = 0u32; + + lut[i] = 0; + i += 1; + for j in i..HIST_SZ { + sum += hist[j]; + lut[j] = (sum as f64 * scale).round().clamp(0.0, 255.0) as u8; + } + + let mut dst = Matrix::::new(src.rows, src.cols, 1); + + #[cfg(feature = "parallel")] + { + dst.data + .par_iter_mut() + .zip(src.data.par_iter()) + .for_each(|(d, &s)| { + *d = lut[s as usize]; + }); + } + + #[cfg(not(feature = "parallel"))] + { + for (d, &s) in dst.data.iter_mut().zip(src.data.iter()) { + *d = lut[s as usize]; + } + } + + Ok(dst) +} + +// --------------------------------------------------------------------------- +// CLAHE +// --------------------------------------------------------------------------- + +/// Contrast Limited Adaptive Histogram Equalization. +#[derive(Debug, Clone)] +pub struct Clahe { + clip_limit: f64, + tiles_x: usize, + tiles_y: usize, + bit_shift: i32, +} + +impl Clahe { + pub fn new(clip_limit: f64, tile_grid_size: Size2i) -> Self { + let tiles_x = tile_grid_size.width.max(0) as usize; + let tiles_y = tile_grid_size.height.max(0) as usize; + Self { + clip_limit, + tiles_x, + tiles_y, + bit_shift: 0, + } + } + + pub fn set_clip_limit(&mut self, clip_limit: f64) { + self.clip_limit = clip_limit; + } + + pub fn get_clip_limit(&self) -> f64 { + self.clip_limit + } + + pub fn set_tiles_grid_size(&mut self, tile_grid_size: Size2i) { + self.tiles_x = tile_grid_size.width.max(0) as usize; + self.tiles_y = tile_grid_size.height.max(0) as usize; + } + + pub fn get_tiles_grid_size(&self) -> Size2i { + Size2i::new(self.tiles_x as i32, self.tiles_y as i32) + } + + pub fn set_bit_shift(&mut self, bit_shift: i32) { + self.bit_shift = bit_shift; + } + + pub fn get_bit_shift(&self) -> i32 { + self.bit_shift + } + + /// Applies CLAHE to a single-channel `u8` or `u16` image. + /// + /// The `src` must be `Matrix` (CV_8UC1) or `Matrix` (CV_16UC1). + /// For u16, the function uses bit_shift to reduce the histogram size. + pub fn apply_u8(&self, src: &Matrix) -> Result> { + if src.channels != 1 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u8: source must be single-channel (got {})", + src.channels + ); + } + self.apply_impl_u8(src) + } + + pub fn apply_u16(&self, src: &Matrix) -> Result> { + if src.channels != 1 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u16: source must be single-channel (got {})", + src.channels + ); + } + self.apply_impl_u16(src) + } + + fn apply_impl_u8(&self, src: &Matrix) -> Result> { + if self.tiles_x == 0 || self.tiles_y == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u8: tiles_x and tiles_y must be > 0 (got {}x{})", + self.tiles_x, + self.tiles_y + ); + } + if self.bit_shift < 0 || self.bit_shift > 7 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u8: bit_shift must be in 0..=7 for u8 (got {})", + self.bit_shift + ); + } + if src.rows == 0 || src.cols == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u8: source must not be empty" + ); + } + let hist_size = 256usize >> self.bit_shift; + // Defensive: unreachable since bit_shift is validated to 0..=7 above. + if hist_size == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u8: invalid hist_size 0 (bit_shift {})", + self.bit_shift + ); + } + + // OpenCV parity: pad bottom/right with BORDER_REFLECT_101 so dimensions become divisible + let pad_bottom = (self.tiles_y - src.rows % self.tiles_y) % self.tiles_y; + let pad_right = (self.tiles_x - src.cols % self.tiles_x) % self.tiles_x; + let need_pad = pad_bottom > 0 || pad_right > 0; + + let padded = if need_pad { + pad_reflect101(src, pad_bottom, pad_right) + } else { + src.clone() + }; + + let tile_rows = padded.rows / self.tiles_y; + let tile_cols = padded.cols / self.tiles_x; + // Defensive: unreachable given a non-empty src and validated tiles, + // since padding above makes the padded dims divisible by the grid. + if tile_rows == 0 || tile_cols == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u8: tile size must be > 0 (padded {}x{} / tiles {}x{} -> tile {}x{})", + padded.rows, + padded.cols, + self.tiles_x, + self.tiles_y, + tile_cols, + tile_rows + ); + } + let tile_size_total = tile_rows * tile_cols; + + let lut_scale = (hist_size as f64 - 1.0) / tile_size_total as f64; + + let mut clip_limit = 0i32; + if self.clip_limit > 0.0 { + clip_limit = (self.clip_limit * tile_size_total as f64 / hist_size as f64) as i32; + clip_limit = clip_limit.max(1); + } + + let num_tiles = match self.tiles_x.checked_mul(self.tiles_y) { + Some(v) => v, + None => cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u8: tiles_x * tiles_y overflow" + ), + }; + let lut_len = match num_tiles.checked_mul(hist_size) { + Some(v) => v, + None => cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u8: lut size overflow" + ), + }; + let mut lut = vec![0u8; lut_len]; + + let build_tile_lut = |tile_idx: usize, tile_lut: &mut [u8]| { + let ty = tile_idx / self.tiles_x; + let tx = tile_idx % self.tiles_x; + let y0 = ty * tile_rows; + let x0 = tx * tile_cols; + + let mut tile_hist = vec![0i32; hist_size]; + for dy in 0..tile_rows { + for dx in 0..tile_cols { + let val = *padded.get(y0 + dy, x0 + dx, 0).unwrap_or(&0) as usize; + let bin = val >> self.bit_shift; + if bin < hist_size { + tile_hist[bin] += 1; + } + } + } + + clip_and_redistribute(&mut tile_hist, clip_limit, hist_size); + + let mut sum = 0i32; + for bin in 0..hist_size { + sum += tile_hist[bin]; + tile_lut[bin] = (sum as f64 * lut_scale).round().clamp(0.0, 255.0) as u8; + } + }; + + #[cfg(feature = "parallel")] + { + lut.par_chunks_mut(hist_size) + .enumerate() + .for_each(|(tile_idx, tile_lut)| { + build_tile_lut(tile_idx, tile_lut); + }); + } + + #[cfg(not(feature = "parallel"))] + { + for (tile_idx, tile_lut) in lut.chunks_mut(hist_size).enumerate() { + build_tile_lut(tile_idx, tile_lut); + } + } + + let mut dst = Matrix::::new(src.rows, src.cols, 1); + interpolate_tiles_u8( + src, + &mut dst, + &lut, + tile_rows, + tile_cols, + self.bit_shift, + self.tiles_x, + self.tiles_y, + hist_size, + ); + Ok(dst) + } + + fn apply_impl_u16(&self, src: &Matrix) -> Result> { + if self.tiles_x == 0 || self.tiles_y == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u16: tiles_x and tiles_y must be > 0 (got {}x{})", + self.tiles_x, + self.tiles_y + ); + } + if self.bit_shift < 0 || self.bit_shift > 15 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u16: bit_shift must be in 0..=15 for u16 (got {})", + self.bit_shift + ); + } + if src.rows == 0 || src.cols == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u16: source must not be empty" + ); + } + let hist_size = 65536usize >> self.bit_shift; + // Defensive: unreachable since bit_shift is validated to 0..=15 above. + if hist_size == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u16: invalid hist_size 0 (bit_shift {})", + self.bit_shift + ); + } + + let pad_bottom = (self.tiles_y - src.rows % self.tiles_y) % self.tiles_y; + let pad_right = (self.tiles_x - src.cols % self.tiles_x) % self.tiles_x; + let need_pad = pad_bottom > 0 || pad_right > 0; + + let padded = if need_pad { + pad_reflect101(src, pad_bottom, pad_right) + } else { + src.clone() + }; + + let tile_rows = padded.rows / self.tiles_y; + let tile_cols = padded.cols / self.tiles_x; + // Defensive: unreachable given a non-empty src and validated tiles, + // since padding above makes the padded dims divisible by the grid. + if tile_rows == 0 || tile_cols == 0 { + cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u16: tile size must be > 0 (padded {}x{} / tiles {}x{} -> tile {}x{})", + padded.rows, + padded.cols, + self.tiles_x, + self.tiles_y, + tile_cols, + tile_rows + ); + } + let tile_size_total = tile_rows * tile_cols; + + let lut_scale = (hist_size as f64 - 1.0) / tile_size_total as f64; + + let mut clip_limit = 0i32; + if self.clip_limit > 0.0 { + clip_limit = (self.clip_limit * tile_size_total as f64 / hist_size as f64) as i32; + clip_limit = clip_limit.max(1); + } + + let num_tiles = match self.tiles_x.checked_mul(self.tiles_y) { + Some(v) => v, + None => cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u16: tiles_x * tiles_y overflow" + ), + }; + let lut_len = match num_tiles.checked_mul(hist_size) { + Some(v) => v, + None => cv_bail!( + tags::IMGPROC, + InvalidInput, + "Clahe::apply_u16: lut size overflow" + ), + }; + let mut lut = vec![0u16; lut_len]; + + let build_tile_lut = |tile_idx: usize, tile_lut: &mut [u16]| { + let ty = tile_idx / self.tiles_x; + let tx = tile_idx % self.tiles_x; + let y0 = ty * tile_rows; + let x0 = tx * tile_cols; + + let mut tile_hist = vec![0i32; hist_size]; + for dy in 0..tile_rows { + for dx in 0..tile_cols { + let val = *padded.get(y0 + dy, x0 + dx, 0).unwrap_or(&0) as usize; + let bin = val >> self.bit_shift; + if bin < hist_size { + tile_hist[bin] += 1; + } + } + } + + clip_and_redistribute(&mut tile_hist, clip_limit, hist_size); + + let mut sum = 0i32; + for bin in 0..hist_size { + sum += tile_hist[bin]; + tile_lut[bin] = (sum as f64 * lut_scale).round().clamp(0.0, 65535.0) as u16; + } + }; + + #[cfg(feature = "parallel")] + { + lut.par_chunks_mut(hist_size) + .enumerate() + .for_each(|(tile_idx, tile_lut)| { + build_tile_lut(tile_idx, tile_lut); + }); + } + + #[cfg(not(feature = "parallel"))] + { + for (tile_idx, tile_lut) in lut.chunks_mut(hist_size).enumerate() { + build_tile_lut(tile_idx, tile_lut); + } + } + + let mut dst = Matrix::::new(src.rows, src.cols, 1); + interpolate_tiles_u16( + src, + &mut dst, + &lut, + tile_rows, + tile_cols, + self.bit_shift, + self.tiles_x, + self.tiles_y, + hist_size, + ); + Ok(dst) + } +} + +fn clip_and_redistribute(tile_hist: &mut [i32], clip_limit: i32, hist_size: usize) { + if clip_limit <= 0 { + return; + } + // Exact port of OpenCV's CLAHE_CalcLut_Body redistribution (clahe.cpp): + // uniform batch plus fixed-step residual. + let mut clipped = 0i32; + for bin in tile_hist.iter_mut() { + if *bin > clip_limit { + clipped += *bin - clip_limit; + *bin = clip_limit; + } + } + + let redist_batch = clipped / hist_size as i32; + let residual = clipped - redist_batch * hist_size as i32; + + for bin in tile_hist.iter_mut() { + *bin += redist_batch; + } + + if residual > 0 { + let step = (hist_size as i32 / residual).max(1); + let mut i = 0i32; + let mut rem = residual; + while i < hist_size as i32 && rem > 0 { + tile_hist[i as usize] += 1; + i += step; + rem -= 1; + } + } +} + +/// Computes the pair of tile indices and blend weight for CLAHE's bilinear +/// tile interpolation along one axis. +/// +/// `coord_f` is the pixel's tile-space coordinate, already offset by -0.5 so +/// tile centers land on integers. Returns `(idx1, idx2, weight1, weight2)` +/// where `weight1 + weight2 == 1.0` and both weights are in `[0, 1]`. +/// +/// The fraction is derived from the *unclamped* floor of `coord_f` so it +/// always lands in `[0, 1)` - a genuine convex blend. Only the indices used +/// for LUT lookup are clamped to the valid tile range; clamping the index +/// *before* computing the fraction (as an earlier version of this code did) +/// produces fractions outside `[0, 1]` at the first/last tile, extrapolating +/// beyond that tile's LUT instead of blending within it. +#[inline(always)] +pub(crate) fn tile_interp_weights(coord_f: f64, num_tiles: usize) -> (usize, usize, f64, f64) { + let idx1_raw = coord_f.floor() as i32; + let weight2 = coord_f - idx1_raw as f64; + let weight1 = 1.0 - weight2; + let max_idx = num_tiles as i32 - 1; + let idx1 = idx1_raw.clamp(0, max_idx) as usize; + let idx2 = (idx1_raw + 1).clamp(0, max_idx) as usize; + (idx1, idx2, weight1, weight2) +} + +#[allow(clippy::too_many_arguments)] +fn interpolate_tiles_u8( + src: &Matrix, + dst: &mut Matrix, + lut: &[u8], + tile_rows: usize, + tile_cols: usize, + bit_shift: i32, + tiles_x: usize, + tiles_y: usize, + hist_size: usize, +) { + let inv_tw = 1.0f64 / tile_cols as f64; + let inv_th = 1.0f64 / tile_rows as f64; + let lut_idx = + |ty: usize, tx: usize, bin: usize| ty * tiles_x * hist_size + tx * hist_size + bin; + let cols = src.cols; + + let process_row = |y: usize, dst_row: &mut [u8]| { + let tyf = y as f64 * inv_th - 0.5; + let (ty1, ty2, ya1, ya) = tile_interp_weights(tyf, tiles_y); + + for (x, out_pixel) in dst_row.iter_mut().enumerate() { + let txf = x as f64 * inv_tw - 0.5; + let (tx1, tx2, xa1, xa) = tile_interp_weights(txf, tiles_x); + + let src_val = *src.get(y, x, 0).unwrap_or(&0) as usize; + let bin = (src_val >> bit_shift).min(hist_size - 1); + + let v00 = lut[lut_idx(ty1, tx1, bin)] as f64; + let v01 = lut[lut_idx(ty1, tx2, bin)] as f64; + let v10 = lut[lut_idx(ty2, tx1, bin)] as f64; + let v11 = lut[lut_idx(ty2, tx2, bin)] as f64; + + let val = (v00 * xa1 + v01 * xa) * ya1 + (v10 * xa1 + v11 * xa) * ya; + *out_pixel = (val.round() as u8) << bit_shift; + } + }; + + #[cfg(feature = "parallel")] + { + dst.data + .par_chunks_mut(cols) + .enumerate() + .for_each(|(y, dst_row)| { + process_row(y, dst_row); + }); + } + + #[cfg(not(feature = "parallel"))] + { + for (y, dst_row) in dst.data.chunks_mut(cols).enumerate() { + process_row(y, dst_row); + } + } +} + +#[allow(clippy::too_many_arguments)] +fn interpolate_tiles_u16( + src: &Matrix, + dst: &mut Matrix, + lut: &[u16], + tile_rows: usize, + tile_cols: usize, + bit_shift: i32, + tiles_x: usize, + tiles_y: usize, + hist_size: usize, +) { + let inv_tw = 1.0f64 / tile_cols as f64; + let inv_th = 1.0f64 / tile_rows as f64; + let lut_idx = + |ty: usize, tx: usize, bin: usize| ty * tiles_x * hist_size + tx * hist_size + bin; + let cols = src.cols; + + let process_row = |y: usize, dst_row: &mut [u16]| { + let tyf = y as f64 * inv_th - 0.5; + let (ty1, ty2, ya1, ya) = tile_interp_weights(tyf, tiles_y); + + for (x, out_pixel) in dst_row.iter_mut().enumerate() { + let txf = x as f64 * inv_tw - 0.5; + let (tx1, tx2, xa1, xa) = tile_interp_weights(txf, tiles_x); + + let src_val = *src.get(y, x, 0).unwrap_or(&0) as usize; + let bin = (src_val >> bit_shift).min(hist_size - 1); + + let v00 = lut[lut_idx(ty1, tx1, bin)] as f64; + let v01 = lut[lut_idx(ty1, tx2, bin)] as f64; + let v10 = lut[lut_idx(ty2, tx1, bin)] as f64; + let v11 = lut[lut_idx(ty2, tx2, bin)] as f64; + + let val = (v00 * xa1 + v01 * xa) * ya1 + (v10 * xa1 + v11 * xa) * ya; + *out_pixel = (val.round() as u16) << bit_shift; + } + }; + + #[cfg(feature = "parallel")] + { + dst.data + .par_chunks_mut(cols) + .enumerate() + .for_each(|(y, dst_row)| { + process_row(y, dst_row); + }); + } + + #[cfg(not(feature = "parallel"))] + { + for (y, dst_row) in dst.data.chunks_mut(cols).enumerate() { + process_row(y, dst_row); + } + } +} + +fn pad_reflect101( + src: &Matrix, + pad_bottom: usize, + pad_right: usize, +) -> Matrix { + let new_rows = src.rows + pad_bottom; + let new_cols = src.cols + pad_right; + let mut dst = Matrix::::new(new_rows, new_cols, 1); + for y in 0..new_rows { + let sy = border_interpolate(y as i32, src.rows as i32, BorderTypes::Reflect101) as usize; + for x in 0..new_cols { + let sx = + border_interpolate(x as i32, src.cols as i32, BorderTypes::Reflect101) as usize; + dst.set(y, x, 0, *src.get(sy, sx, 0).unwrap_or(&T::default())); + } + } + dst +} + +pub fn create_clahe(clip_limit: f64, tile_grid_size: Size2i) -> Clahe { + Clahe::new(clip_limit, tile_grid_size) +} diff --git a/src/imgproc/simd.rs b/src/imgproc/simd.rs index 5b731db..c790377 100644 --- a/src/imgproc/simd.rs +++ b/src/imgproc/simd.rs @@ -265,3 +265,109 @@ mod tests { assert_eq!(dst[3], 4.0); } } +/// SIMD-friendly (LLVM auto-vectorized) compare_hist methods. +#[cfg(feature = "simd")] +pub(crate) fn simd_compare_hist_f32(h1: &[f32], h2: &[f32], method: u8) -> Option { + match method { + 0 => { + // Correl + let mut s1 = 0.0f64; + let mut s2 = 0.0f64; + let mut s11 = 0.0f64; + let mut s12 = 0.0f64; + let mut s22 = 0.0f64; + for (a, b) in h1.iter().copied().zip(h2.iter().copied()) { + let a = a as f64; + let b = b as f64; + s1 += a; + s2 += b; + s11 += a * a; + s22 += b * b; + s12 += a * b; + } + let n = h1.len() as f64; + let scale = 1.0 / n; + let num = s12 - s1 * s2 * scale; + let denom2 = (s11 - s1 * s1 * scale) * (s22 - s2 * s2 * scale); + Some(if denom2.abs() > f64::EPSILON { + num / denom2.sqrt() + } else { + 1.0 + }) + } + 1 => { + // ChiSqr + let mut result = 0.0f64; + for (a, b) in h1.iter().copied().zip(h2.iter().copied()) { + let a = a as f64; + let b = b as f64; + let diff = a - b; + let a_adj = if a.abs() <= f64::EPSILON { 1.0 } else { a }; + let val = diff * diff / a_adj; + let add_val = if a.abs() > f64::EPSILON { val } else { 0.0 }; + result += add_val; + } + Some(result) + } + 2 => { + // ChiSqrAlt + let mut result = 0.0f64; + for (a, b) in h1.iter().copied().zip(h2.iter().copied()) { + let a = a as f64; + let b = b as f64; + let sum = a + b; + let diff = a - b; + let sum_adj = if sum.abs() <= f64::EPSILON { 1.0 } else { sum }; + let val = diff * diff / sum_adj; + let add_val = if sum.abs() > f64::EPSILON { val } else { 0.0 }; + result += add_val; + } + Some(result * 2.0) + } + 3 => { + // Intersection + let mut result = 0.0f64; + for (a, b) in h1.iter().copied().zip(h2.iter().copied()) { + let a = a as f64; + let b = b as f64; + result += a.min(b); + } + Some(result) + } + 4 => { + // Bhattacharyya + let mut s1 = 0.0f64; + let mut s2 = 0.0f64; + let mut bc = 0.0f64; + for (a, b) in h1.iter().copied().zip(h2.iter().copied()) { + let a = a as f64; + let b = b as f64; + s1 += a; + s2 += b; + bc += (a * b).sqrt(); + } + let norm = s1 * s2; + let norm_factor = if norm.abs() > f64::EPSILON { + 1.0 / norm.sqrt() + } else { + 1.0 + }; + Some(((1.0 - bc * norm_factor).max(0.0)).sqrt()) + } + 5 => { + // KullbackLeibler + let mut result = 0.0f64; + for (a, b) in h1.iter().copied().zip(h2.iter().copied()) { + let a = a as f64; + let b = b as f64; + let q_adj = if b.abs() <= f64::EPSILON { 1e-10 } else { b }; + let p_adj = if a.abs() <= f64::EPSILON { 1.0 } else { a }; + let val = a * (p_adj / q_adj).ln(); + let add_val = if a.abs() > f64::EPSILON { val } else { 0.0 }; + result += add_val; + } + Some(result) + } + _ => None, + } +} diff --git a/src/imgproc/tests.rs b/src/imgproc/tests.rs index b758584..58fc9da 100644 --- a/src/imgproc/tests.rs +++ b/src/imgproc/tests.rs @@ -1316,4 +1316,773 @@ mod imgproc_tests { assert_eq!(res.data, vec![0.0, 10.0, 0.0, 30.0]); } + + #[test] + fn test_compare_hist_simd_coverage() { + let mut h1 = Matrix::::new(1, 10, 1); + let mut h2 = Matrix::::new(1, 10, 1); + for i in 0..10 { + *h1.at_mut(0, i, 0).unwrap() = i as f32; + *h2.at_mut(0, i, 0).unwrap() = (9 - i) as f32; + } + + let c = compare_hist(&h1, &h2, HistCompMethods::Correl).unwrap(); + assert!(c > -1.1 && c < 1.1); + + let c = compare_hist(&h1, &h2, HistCompMethods::ChiSqr).unwrap(); + assert!(c >= 0.0); + + let c = compare_hist(&h1, &h2, HistCompMethods::Intersection).unwrap(); + assert!(c >= 0.0); + + let c = compare_hist(&h1, &h2, HistCompMethods::Bhattacharyya).unwrap(); + assert!(c >= 0.0); + } + + #[test] + fn test_calc_hist_uniform_1d() { + let data: Vec = (0..16).collect(); + let img = Matrix::from_vec(4, 4, 1, data); + let hist = calc_hist( + &[&img], + &[0], + None, + &[16], + &[RangeSpec::Uniform(0.0, 16.0)], + false, + None, + ) + .unwrap(); + assert_eq!(hist.data.len(), 16); + for &v in hist.data.iter() { + assert_eq!(v, 1.0); + } + } + + #[test] + fn test_calc_hist_uniform_1d_fewer_bins() { + let data: Vec = (0..16).collect(); + let img = Matrix::from_vec(4, 4, 1, data); + let hist = calc_hist( + &[&img], + &[0], + None, + &[4], + &[RangeSpec::Uniform(0.0, 16.0)], + false, + None, + ) + .unwrap(); + assert_eq!(hist.data.len(), 4); + for &v in hist.data.iter() { + assert_eq!(v, 4.0); + } + } + + #[test] + fn test_calc_hist_mask() { + let img = Matrix::from_vec(3, 4, 1, vec![0u8, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]); + // Partial mask: checkerboard + let mask = Matrix::from_vec(3, 4, 1, vec![0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1]); + let h = calc_hist( + &[&img], + &[0], + Some(&mask), + &[4], + &[RangeSpec::Uniform(0.0, 12.0)], + false, + None, + ) + .unwrap(); + assert_eq!(h.data, vec![1.0, 2.0, 1.0, 2.0]); + + // All-zero mask: nothing counted + let img2 = Matrix::from_vec(2, 2, 1, vec![0u8, 1, 2, 3]); + let mask_zero = Matrix::from_vec(2, 2, 1, vec![0u8; 4]); + let h = calc_hist( + &[&img2], + &[0], + Some(&mask_zero), + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .unwrap(); + assert_eq!(h.data, vec![0.0, 0.0, 0.0, 0.0]); + + // All-one mask: all counted + let mask_one = Matrix::from_vec(2, 2, 1, vec![1u8; 4]); + let h = calc_hist( + &[&img2], + &[0], + Some(&mask_one), + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .unwrap(); + assert_eq!(h.data, vec![1.0, 1.0, 1.0, 1.0]); + } + + #[test] + fn test_calc_hist_multichannel_mask_error() { + // Regression test: a mask must be single-channel. Previously only + // rows/cols were checked, so a multichannel mask was silently + // accepted and only its channel 0 was ever read. + let img = Matrix::from_vec(2, 2, 1, vec![0u8, 1, 2, 3]); + let mask = Matrix::from_vec(2, 2, 3, vec![1u8; 12]); + assert!(calc_hist( + &[&img], + &[0], + Some(&mask), + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .is_err()); + } + + #[test] + fn test_calc_hist_multi_image_channel_indexing() { + // images[0] has 2 channels, images[1] has 1 channel + let img0 = Matrix::from_vec(2, 1, 2, vec![10, 20, 30, 40]); + let img1 = Matrix::from_vec(2, 1, 1, vec![100, 200]); + + let hist = calc_hist( + &[&img0, &img1], + &[0, 2], + None, + &[2, 2], + &[ + RangeSpec::Uniform(0.0, 50.0), + RangeSpec::Uniform(0.0, 250.0), + ], + false, + None, + ) + .unwrap(); + + // pixel(0,0): ch0=10->bin0, ch2=100->bin0 -> idx=0 + // pixel(1,0): ch0=30->bin1, ch2=200->bin1 -> idx=3 + assert_eq!(hist.data[0], 1.0); + assert_eq!(hist.data[1], 0.0); + assert_eq!(hist.data[2], 0.0); + assert_eq!(hist.data[3], 1.0); + } + + #[test] + fn test_calc_hist_nonuniform() { + let data: Vec = (0..20).collect(); + let img = Matrix::from_vec(4, 5, 1, data); + // Non-uniform: boundaries [0, 5, 20] -> 2 bins: [0,5) and [5,20) + let hist = calc_hist( + &[&img], + &[0], + None, + &[2], + &[RangeSpec::NonUniform(vec![0.0, 5.0, 20.0])], + false, + None, + ) + .unwrap(); + assert_eq!(hist.data[0], 5.0); // values 0..4 + assert_eq!(hist.data[1], 15.0); // values 5..19 + } + + #[test] + fn test_calc_back_project() { + let data: Vec = vec![0, 1, 2, 3, 4, 5, 6, 7]; + let img = Matrix::from_vec(2, 4, 1, data); + let hist = Matrix::from_vec(4, 1, 1, vec![10.0, 20.0, 30.0, 40.0]); + let bp = calc_back_project( + &[&img], + &[0], + &[4], + &hist, + &[RangeSpec::Uniform(0.0, 8.0)], + 1.0, + ) + .unwrap(); + assert_eq!( + bp.data, + vec![10.0, 10.0, 20.0, 20.0, 30.0, 30.0, 40.0, 40.0] + ); + } + + #[test] + fn test_calc_back_project_2d_shape() { + // Regression test: without an explicit hist_size, a flat 8-bin + // histogram used to have its shape *guessed* from its length alone + // (infer_hist_size), which silently produced the wrong shape for any + // non-perfect-power dims (e.g. [2, 4] was inferred as [1, 8]) and + // corrupted the bin math. hist_size is now required explicitly. + let img = Matrix::from_vec(1, 1, 2, vec![1u8, 7]); + let hist = Matrix::from_vec( + 8, + 1, + 1, + vec![10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0], + ); + let bp = calc_back_project( + &[&img], + &[0, 1], + &[2, 4], + &hist, + &[RangeSpec::Uniform(0.0, 2.0), RangeSpec::Uniform(0.0, 8.0)], + 1.0, + ) + .unwrap(); + // channel 0 (value 1, range [0,2), 2 bins) -> bin 1 + // channel 1 (value 7, range [0,8), 4 bins) -> bin 3 + // flat index with strides [4, 1] -> 1*4 + 3 = 7 -> hist.data[7] = 80.0 + assert_eq!(bp.data, vec![80.0]); + } + + #[test] + fn test_calc_back_project_hist_size_mismatch_error() { + let img = Matrix::from_vec(1, 4, 1, vec![0u8, 1, 2, 3]); + let hist = Matrix::from_vec(4, 1, 1, vec![10.0, 20.0, 30.0, 40.0]); + assert!(calc_back_project( + &[&img], + &[0], + &[3], // doesn't match hist's 4 bins + &hist, + &[RangeSpec::Uniform(0.0, 4.0)], + 1.0, + ) + .is_err()); + } + + #[test] + fn test_calc_back_project_invalid_uniform_range_error() { + // Regression test: calc_back_project used to skip the range + // validation calc_hist already had (lo < hi, NaN rejection, + // non-uniform boundary shape/ordering), silently producing a + // plausible-but-wrong projection instead of an error. + let img = Matrix::from_vec(1, 4, 1, vec![0u8, 1, 2, 3]); + let hist = Matrix::from_vec(4, 1, 1, vec![10.0, 20.0, 30.0, 40.0]); + + // NaN bounds: lo.partial_cmp(hi) is None, must be rejected. + assert!(calc_back_project( + &[&img], + &[0], + &[4], + &hist, + &[RangeSpec::Uniform(f32::NAN, f32::NAN)], + 1.0, + ) + .is_err()); + + // lo >= hi. + assert!(calc_back_project( + &[&img], + &[0], + &[4], + &hist, + &[RangeSpec::Uniform(4.0, 0.0)], + 1.0, + ) + .is_err()); + } + + #[test] + fn test_calc_back_project_invalid_nonuniform_boundaries_error() { + let img = Matrix::from_vec(1, 4, 1, vec![0u8, 1, 2, 3]); + let hist = Matrix::from_vec(4, 1, 1, vec![10.0, 20.0, 30.0, 40.0]); + + // Wrong length: needs hist_size[0] + 1 = 5 boundaries, only 3 given. + assert!(calc_back_project( + &[&img], + &[0], + &[4], + &hist, + &[RangeSpec::NonUniform(vec![0.0, 2.0, 4.0])], + 1.0, + ) + .is_err()); + + // Not strictly increasing. + assert!(calc_back_project( + &[&img], + &[0], + &[4], + &hist, + &[RangeSpec::NonUniform(vec![0.0, 2.0, 1.0, 3.0, 4.0])], + 1.0, + ) + .is_err()); + } + + #[test] + fn test_calc_back_project_hist_size_overflow_error() { + // Regression test: validate_hist_ranges computed hist_size[d] + 1 + // unchecked, so hist_size[d] == usize::MAX panicked with "attempt to + // add with overflow" instead of returning an InvalidInput error. + // Confirmed the pre-fix code panics here via std::panic::catch_unwind. + let img = Matrix::from_vec(1, 1, 1, vec![0u8]); + let hist = Matrix::from_vec(1, 1, 1, vec![1.0f32]); + let result = calc_back_project( + &[&img], + &[0], + &[usize::MAX], + &hist, + &[RangeSpec::NonUniform(vec![])], + 1.0, + ); + assert!(result.is_err()); + } + + #[test] + fn test_compare_hist_correl_identical() { + let h1 = Matrix::from_vec(4, 1, 1, vec![1.0, 2.0, 3.0, 4.0]); + let h2 = Matrix::from_vec(4, 1, 1, vec![1.0, 2.0, 3.0, 4.0]); + let corr = compare_hist(&h1, &h2, HistCompMethods::Correl).unwrap(); + assert!((corr - 1.0).abs() < 1e-10); + } + + #[test] + fn test_compare_hist_correl_opposite() { + let h1 = Matrix::from_vec(4, 1, 1, vec![1.0, 2.0, 3.0, 4.0]); + let h2 = Matrix::from_vec(4, 1, 1, vec![4.0, 3.0, 2.0, 1.0]); + let corr = compare_hist(&h1, &h2, HistCompMethods::Correl).unwrap(); + assert!((corr - (-1.0)).abs() < 1e-10); + } + + #[test] + fn test_compare_hist_intersect() { + let h1 = Matrix::from_vec(4, 1, 1, vec![1.0, 2.0, 3.0, 4.0]); + let h2 = Matrix::from_vec(4, 1, 1, vec![4.0, 3.0, 2.0, 1.0]); + let inter = compare_hist(&h1, &h2, HistCompMethods::Intersection).unwrap(); + assert!((inter - 6.0).abs() < 1e-10); + } + + #[test] + fn test_compare_hist_chi_sqr() { + let h1 = Matrix::from_vec(3, 1, 1, vec![1.0, 2.0, 3.0]); + let h2 = Matrix::from_vec(3, 1, 1, vec![1.0, 2.0, 3.0]); + let chi = compare_hist(&h1, &h2, HistCompMethods::ChiSqr).unwrap(); + assert!((chi - 0.0).abs() < 1e-10); + } + + #[test] + fn test_compare_hist_bhattacharyya_identical() { + let h1 = Matrix::from_vec(4, 1, 1, vec![0.25, 0.25, 0.25, 0.25]); + let h2 = Matrix::from_vec(4, 1, 1, vec![0.25, 0.25, 0.25, 0.25]); + let bc = compare_hist(&h1, &h2, HistCompMethods::Bhattacharyya).unwrap(); + assert!(bc.abs() < 1e-10); + } + + #[test] + fn test_compare_hist_kl_divergence() { + let h1 = Matrix::from_vec(3, 1, 1, vec![0.5, 0.3, 0.2]); + let h2 = Matrix::from_vec(3, 1, 1, vec![0.5, 0.3, 0.2]); + let kl = compare_hist(&h1, &h2, HistCompMethods::KullbackLeibler).unwrap(); + assert!(kl.abs() < 1e-10); + } + + #[test] + fn test_equalize_hist_uniform() { + let img = Matrix::from_vec(4, 4, 1, vec![128u8; 16]); + let dst = equalize_hist(&img).unwrap(); + for &v in dst.data.iter() { + assert_eq!(v, 128); + } + } + + #[test] + fn test_equalize_hist_gradient() { + let data: Vec = (0..=255).collect(); + let img = Matrix::from_vec(16, 16, 1, data); + let dst = equalize_hist(&img).unwrap(); + assert_eq!(dst.rows, 16); + assert_eq!(dst.cols, 16); + assert_eq!(dst.channels, 1); + let min_val = *dst.data.iter().min().unwrap(); + let max_val = *dst.data.iter().max().unwrap(); + assert_eq!(min_val, 0); + assert_eq!(max_val, 255); + } + + #[test] + fn test_clahe_basic() { + let data: Vec = (0..=255).collect(); + let img = Matrix::from_vec(16, 16, 1, data); + let clahe = create_clahe(40.0, Size2i::new(4, 4)); + let dst = clahe.apply_u8(&img).unwrap(); + assert_eq!(dst.rows, 16); + assert_eq!(dst.cols, 16); + assert_eq!(dst.channels, 1); + assert_eq!(dst.data.len(), 256); + } + + #[test] + fn test_clahe_u16() { + let data: Vec = (0..=1023).collect(); + let img = Matrix::from_vec(32, 32, 1, data); + let clahe = create_clahe(40.0, Size2i::new(4, 4)); + let dst = clahe.apply_u16(&img).unwrap(); + assert_eq!(dst.rows, 32); + assert_eq!(dst.cols, 32); + assert_eq!(dst.channels, 1); + } + + #[test] + fn test_calc_hist_2d() { + let img = Matrix::from_vec(3, 3, 1, vec![0, 1, 2, 3, 4, 5, 6, 7, 8]); + let hist = calc_hist( + &[&img, &img], + &[0, 0], + None, + &[3, 3], + &[RangeSpec::Uniform(0.0, 9.0), RangeSpec::Uniform(0.0, 9.0)], + false, + None, + ) + .unwrap(); + assert_eq!(hist.data.len(), 9); + assert_eq!(hist.data[0], 3.0); + assert_eq!(hist.data[1], 0.0); + assert_eq!(hist.data[2], 0.0); + assert_eq!(hist.data[3], 0.0); + assert_eq!(hist.data[4], 3.0); + assert_eq!(hist.data[5], 0.0); + assert_eq!(hist.data[6], 0.0); + assert_eq!(hist.data[7], 0.0); + assert_eq!(hist.data[8], 3.0); + } + + #[test] + fn test_calc_hist_empty_images_error() { + assert!(calc_hist::( + &[], + &[0], + None, + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .is_err()); + } + + #[test] + fn test_calc_hist_mismatched_channels_error() { + let img = Matrix::from_vec(2, 2, 1, vec![0u8; 4]); + assert!(calc_hist( + &[&img], + &[0, 1], + None, + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .is_err()); + } + + #[test] + fn test_calc_back_project_empty_channels_error() { + let img = Matrix::from_vec(2, 2, 1, vec![0u8; 4]); + let hist = Matrix::from_vec(4, 1, 1, vec![0.0; 4]); + assert!(calc_back_project::( + &[&img], + &[], + &[], + &hist, + &[RangeSpec::Uniform(0.0, 4.0)], + 1.0, + ) + .is_err()); + } + + #[test] + fn test_calc_back_project_zero_width() { + // Regression test: chunks_mut/par_chunks_mut panic on a zero chunk + // size regardless of slice length, so a zero-width image must not + // reach the row-chunking dispatch. + let img = Matrix::::from_vec(3, 0, 1, vec![]); + let hist = Matrix::from_vec(4, 1, 1, vec![0.0; 4]); + let dst = calc_back_project( + &[&img], + &[0], + &[4], + &hist, + &[RangeSpec::Uniform(0.0, 4.0)], + 1.0, + ) + .unwrap(); + assert_eq!(dst.rows, 3); + assert_eq!(dst.cols, 0); + assert_eq!(dst.data.len(), 0); + } + + #[test] + fn test_calc_hist_u16_input() { + let data: Vec = vec![0, 100, 200, 300, 400, 500]; + let img = Matrix::from_vec(2, 3, 1, data); + let hist = calc_hist( + &[&img], + &[0], + None, + &[3], + &[RangeSpec::Uniform(0.0, 600.0)], + false, + None, + ) + .unwrap(); + assert_eq!(hist.data, vec![2.0, 2.0, 2.0]); + } + + #[test] + fn test_calc_hist_f32_input() { + let data: Vec = vec![0.5, 1.5, 2.5, 3.5]; + let img = Matrix::from_vec(2, 2, 1, data); + let hist = calc_hist( + &[&img], + &[0], + None, + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .unwrap(); + assert_eq!(hist.data, vec![1.0, 1.0, 1.0, 1.0]); + } + + #[test] + fn test_calc_hist_boundary_exclusion() { + let img = Matrix::from_vec(1, 4, 1, vec![0u8, 4, 8, 12]); + let hist = calc_hist( + &[&img], + &[0], + None, + &[4], + &[RangeSpec::Uniform(0.0, 16.0)], + false, + None, + ) + .unwrap(); + assert_eq!(hist.data, vec![1.0, 1.0, 1.0, 1.0]); + + // Value at exact hi boundary should be excluded + let img2 = Matrix::from_vec(1, 2, 1, vec![4u8, 4]); + let hist2 = calc_hist( + &[&img2], + &[0], + None, + &[2], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .unwrap(); + assert_eq!(hist2.data, vec![0.0, 0.0]); + + // Value just below hi should be included + let img3 = Matrix::from_vec(1, 1, 1, vec![3u8]); + let hist3 = calc_hist( + &[&img3], + &[0], + None, + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .unwrap(); + assert_eq!(hist3.data[3], 1.0); + } + + #[test] + fn test_calc_hist_accumulate() { + let img = Matrix::from_vec(1, 4, 1, vec![0u8, 1, 2, 3]); + let h1 = calc_hist( + &[&img], + &[0], + None, + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .unwrap(); + assert_eq!(h1.data, vec![1.0, 1.0, 1.0, 1.0]); + + let h2 = calc_hist( + &[&img], + &[0], + None, + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + true, + Some(&h1), + ) + .unwrap(); + assert_eq!(h2.data, vec![2.0, 2.0, 2.0, 2.0]); + + let h3 = calc_hist( + &[&img], + &[0], + None, + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + true, + Some(&h2), + ) + .unwrap(); + assert_eq!(h3.data, vec![3.0, 3.0, 3.0, 3.0]); + } + + #[test] + fn test_compare_hist_edge_cases() { + let z = Matrix::from_vec(3, 1, 1, vec![0.0, 0.0, 0.0]); + assert!((compare_hist(&z, &z, HistCompMethods::Correl).unwrap() - 1.0).abs() < 1e-10); + assert!((compare_hist(&z, &z, HistCompMethods::Intersection).unwrap()).abs() < 1e-10); + assert!( + (compare_hist(&z, &z, HistCompMethods::Bhattacharyya).unwrap() - 1.0).abs() < 1e-10 + ); + + let p1 = Matrix::from_vec(4, 1, 1, vec![0.5, 0.25, 0.125, 0.125]); + let p2 = Matrix::from_vec(4, 1, 1, vec![0.5, 0.25, 0.125, 0.125]); + assert!((compare_hist(&p1, &p2, HistCompMethods::Correl).unwrap() - 1.0).abs() < 1e-10); + assert!((compare_hist(&p1, &p2, HistCompMethods::ChiSqr).unwrap()).abs() < 1e-10); + } + + #[test] + fn test_calc_back_project_scale() { + let img = Matrix::from_vec(1, 4, 1, vec![0u8, 1, 2, 3]); + let hist = Matrix::from_vec(4, 1, 1, vec![10.0, 20.0, 30.0, 40.0]); + let bp = calc_back_project( + &[&img], + &[0], + &[4], + &hist, + &[RangeSpec::Uniform(0.0, 4.0)], + 2.0, + ) + .unwrap(); + assert_eq!(bp.data, vec![20.0, 40.0, 60.0, 80.0]); + } + + #[test] + fn test_calc_hist_errors() { + let img = Matrix::from_vec(2, 2, 1, vec![0u8; 4]); + assert!(calc_hist::( + &[], + &[0], + None, + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .is_err()); + assert!(calc_hist( + &[&img], + &[0, 1], + None, + &[4], + &[RangeSpec::Uniform(0.0, 4.0)], + false, + None, + ) + .is_err()); + let img2 = Matrix::from_vec(3, 3, 1, vec![0u8; 9]); + assert!(calc_hist( + &[&img, &img2], + &[0, 0], + None, + &[2, 2], + &[RangeSpec::Uniform(0.0, 2.0), RangeSpec::Uniform(0.0, 2.0)], + false, + None, + ) + .is_err()); + let hist = Matrix::from_vec(4, 1, 1, vec![0.0; 4]); + assert!(calc_back_project::( + &[&img], + &[], + &[], + &hist, + &[RangeSpec::Uniform(0.0, 4.0)], + 1.0, + ) + .is_err()); + } + + #[test] + fn test_calc_hist_multichannel_select() { + let img = Matrix::from_vec( + 1, + 3, + 3, + vec![ + 10, 100, 200, // pixel 0 + 20, 150, 250, // pixel 1 + 30, 50, 100, // pixel 2 + ], + ); + let hist = calc_hist( + &[&img], + &[1], + None, + &[3], + &[RangeSpec::Uniform(0.0, 200.0)], + false, + None, + ) + .unwrap(); + assert_eq!(hist.data, vec![1.0, 1.0, 1.0]); + } + + #[test] + fn test_clahe_no_clip() { + let data: Vec = (0..=255).collect(); + let img = Matrix::from_vec(16, 16, 1, data); + let clahe = create_clahe(0.0, Size2i::new(4, 4)); + let dst = clahe.apply_u8(&img).unwrap(); + assert_eq!(dst.data.len(), 256); + } + + #[test] + fn test_clahe_tile_interp_weights_are_bounded() { + // Regression test: the interpolation fraction used to be derived + // *after* clamping the tile index, which produced weights outside + // [0, 1] (extrapolation) at the top/left edges - e.g. at y=0 with + // tile_rows=4, tyf=-0.5 used to yield ya=-0.5, ya1=1.5 instead of a + // convex blend. Check the invariant holds across the coordinate + // range that actually occurs (coord_f = pixel*inv_tile_extent - 0.5, + // for pixel in 0..src_extent), for several tile counts. + use crate::imgproc::histogram::tile_interp_weights; + + for num_tiles in [1usize, 2, 3, 5] { + let tile_extent = 4usize; // pixels per tile + let inv = 1.0 / tile_extent as f64; + let src_extent = num_tiles * tile_extent; + + for pixel in 0..src_extent { + let coord_f = pixel as f64 * inv - 0.5; + let (idx1, idx2, w1, w2) = tile_interp_weights(coord_f, num_tiles); + + assert!( + (0.0..=1.0).contains(&w1) && (0.0..=1.0).contains(&w2), + "weights out of [0,1] at pixel={pixel}, num_tiles={num_tiles}: w1={w1}, w2={w2}" + ); + assert!( + (w1 + w2 - 1.0).abs() < 1e-12, + "weights don't sum to 1 at pixel={pixel}: w1={w1}, w2={w2}" + ); + assert!(idx1 < num_tiles && idx2 < num_tiles); + } + } + + // The exact case from the bug report: y=0, tile_rows=4 -> tyf=-0.5. + let (idx1, idx2, w1, w2) = tile_interp_weights(-0.5, 2); + assert_eq!((idx1, idx2), (0, 0)); // both clamp to the first tile + assert!((w1 - 0.5).abs() < 1e-12 && (w2 - 0.5).abs() < 1e-12); + } }