diff --git a/README.md b/README.md index 20a9eac..7b9af18 100644 --- a/README.md +++ b/README.md @@ -44,7 +44,7 @@ Unlike existing wrappers, **PureCV** is a native rewrite. It aims to provide: - **Utilities:** `add_weighted`, `check_range`, `absdiff`, `get_tick_count`, `get_tick_frequency`. - **Mathematical Constants:** OpenCV-compatible constants — `CV_PI`, `CV_PI_2`, `CV_2PI`, `CV_PI_4`, `CV_LOG2`, `CV_LN2`, `CV_E`, `CV_LN10`, `CV_SQRT2` — backed by `std::f64::consts` for maximum precision. - **ndarray Interop:** Optional, zero-cost conversions to/from `ndarray::Array3` via the `ndarray` feature flag. -- **SIMD Acceleration** (`simd` feature): Trait-based dispatch via `pulp` for `f32`, `f64`, and `u8` types. Accelerated operations include `add`, `sub`, `mul`, `div`, `min`, `max`, `sqrt`, `dot`, `sum`, `add_weighted`, `convert_scale_abs`, `magnitude`, `simd_row_min_max`, `simd_min_max_col`, and `simd_gaussian_5tap_h/v`. Falls back to scalar loops at zero cost when disabled. +- **SIMD Acceleration** (`simd` feature): Trait-based dispatch via `pulp` for `f32`, `f64`, and `u8` types. Accelerated operations include `add`, `sub`, `mul`, `div`, `min`, `max`, `sqrt`, `dot`, `sum`, `add_weighted`, `convert_scale_abs`, `magnitude`, `simd_row_min_max`, `simd_min_max_col`, `simd_gaussian_5tap_h/v`, and `simd_remap_bilinear_row`/`simd_remap_nearest_row`. Falls back to scalar loops at zero cost when disabled. ### `purecv-imgproc` - **Color Conversions:** High-performance `cvt_color` supporting RGB, BGR, Gray, RGBA, BGRA and more. Up to **6.6× speedup** with Parallel + SIMD. SIMD-accelerated paths (`simd` feature) use fixed-point integer arithmetic (coefficients 77/150/29 ≈ 0.299/0.587/0.114 × 256) for all `*_to_gray` conversions — portable to x86 SSE/AVX, ARM NEON, and WASM `simd128` via `pulp`. @@ -56,6 +56,7 @@ Unlike existing wrappers, **PureCV** is a native rewrite. It aims to provide: - **Feature Detection:** `corner_harris`, `corner_min_eigen_val` (Shi-Tomasi), `good_features_to_track`, `corner_sub_pix` refinement, and structure tensor computation via `corner_eigen_vals_and_vecs`. Supports both Harris and Shi-Tomasi responses with non-maximum suppression. - **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. ### `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. @@ -67,6 +68,8 @@ Unlike existing wrappers, **PureCV** is a native rewrite. It aims to provide: - **Optical Flow:** Pyramidal Lucas-Kanade optical flow implementation with `calc_optical_flow_pyr_lk` and `build_optical_flow_pyramid`. Includes robust window-based tracking, sub-pixel accuracy, spatial gradient optimization, and iterative refinement. ### `purecv-calib3d` +- **Camera Undistortion:** `init_undistort_rectify_map` to compute lens undistortion and rectification maps. +- **Epipolar Geometry:** `find_fundamental_mat` supporting the normalized 8-Point algorithm and robust RANSAC estimation with Sampson distance. - **Pose Estimation:** Camera pose estimation using `solve_pnp` (Iterative) and `solve_pnp_ransac`. - **Homography:** Direct Linear Transformation (DLT) and RANSAC-based `find_homography` for robust planar perspective mapping. - **Geometry:** `rodrigues` for converting between rotation vectors and 3x3 rotation matrices. @@ -250,6 +253,9 @@ cargo run --example optical_flow_video # Pose Estimation (solve_pnp and rodrigues) cargo run --example pose_estimation + +# Camera Undistortion and Perspective Warping +cargo run --example rectification ``` ## 🧪 Testing & Benchmarking diff --git a/crates/wasm/Cargo.toml b/crates/wasm/Cargo.toml index d268e52..959ca8b 100644 --- a/crates/wasm/Cargo.toml +++ b/crates/wasm/Cargo.toml @@ -21,6 +21,7 @@ js-sys = "0.3.69" web-sys = { version = "0.3.69", features = ["console"] } serde = { version = "1.0", features = ["derive"] } serde-wasm-bindgen = "0.6" +num-traits = "0.2" purecv = { path = "../../", default-features = false } [features] diff --git a/crates/wasm/scripts/build-dual.ps1 b/crates/wasm/scripts/build-dual.ps1 index 1745a45..9ba90a5 100644 --- a/crates/wasm/scripts/build-dual.ps1 +++ b/crates/wasm/scripts/build-dual.ps1 @@ -3,6 +3,7 @@ $ErrorActionPreference = "Stop" + # Get the directory where the script is located $scriptDir = $PSScriptRoot $wasmDir = Join-Path $scriptDir ".." | Resolve-Path diff --git a/crates/wasm/src/lib.rs b/crates/wasm/src/lib.rs index 8d49f32..a25f6c2 100644 --- a/crates/wasm/src/lib.rs +++ b/crates/wasm/src/lib.rs @@ -1086,6 +1086,136 @@ pub fn laplacian( }) } +/// Helper to convert a JS integer into an `InterpolationFlags`. +fn interpolation_flags_from_i32(val: i32) -> Result { + match val { + 0 => Ok(purecv::imgproc::InterpolationFlags::Nearest), + 1 => Ok(purecv::imgproc::InterpolationFlags::Linear), + _ => Err(JsError::new(&format!("Unknown interpolation type: {val}"))), + } +} + +/// Applies a generic geometrical transformation to an image. +#[wasm_bindgen(js_name = "remap")] +pub fn remap_wasm( + src: &Mat, + map1: &Mat, + map2: &Mat, + interpolation: i32, + border_mode: i32, + border_value: &Scalar, +) -> Result { + use purecv::core::types::Scalar as CoreScalar; + use purecv::imgproc::remap; + + let m1 = require_f32(map1, "remap (map1)")?; + let m2 = require_f32(map2, "remap (map2)")?; + let interp = interpolation_flags_from_i32(interpolation)?; + let border = border_type_from_i32(border_mode)?; + + let data = match &src.inner.data { + DynamicData::U8(m) => { + let bv = CoreScalar::new( + num_traits::FromPrimitive::from_f64(border_value.v0).unwrap_or(0), + num_traits::FromPrimitive::from_f64(border_value.v1).unwrap_or(0), + num_traits::FromPrimitive::from_f64(border_value.v2).unwrap_or(0), + num_traits::FromPrimitive::from_f64(border_value.v3).unwrap_or(0), + ); + let result = + remap(m, m1, m2, interp, border, bv).map_err(|e| JsError::new(&format!("{e}")))?; + DynamicData::U8(result) + } + DynamicData::F32(m) => { + let bv = CoreScalar::new( + num_traits::FromPrimitive::from_f64(border_value.v0).unwrap_or(0.0), + num_traits::FromPrimitive::from_f64(border_value.v1).unwrap_or(0.0), + num_traits::FromPrimitive::from_f64(border_value.v2).unwrap_or(0.0), + num_traits::FromPrimitive::from_f64(border_value.v3).unwrap_or(0.0), + ); + let result = + remap(m, m1, m2, interp, border, bv).map_err(|e| JsError::new(&format!("{e}")))?; + DynamicData::F32(result) + } + DynamicData::F64(m) => { + let bv = CoreScalar::new( + border_value.v0, + border_value.v1, + border_value.v2, + border_value.v3, + ); + let result = + remap(m, m1, m2, interp, border, bv).map_err(|e| JsError::new(&format!("{e}")))?; + DynamicData::F64(result) + } + _ => return Err(JsError::new("remap: unsupported depth type")), + }; + + Ok(Mat { + inner: DynamicMatrix { data }, + }) +} + +/// Applies a perspective transformation to an image. +#[wasm_bindgen(js_name = "warpPerspective")] +pub fn warp_perspective_wasm( + src: &Mat, + m: &Mat, + dsize_width: i32, + dsize_height: i32, + flags: i32, + border_mode: i32, + border_value: &Scalar, +) -> Result { + use purecv::core::types::{Scalar as CoreScalar, Size2i}; + use purecv::imgproc::warp_perspective; + + let m_mat = require_f64(m, "warpPerspective (M)")?; + let interp = interpolation_flags_from_i32(flags)?; + let border = border_type_from_i32(border_mode)?; + let dsize = Size2i::new(dsize_width, dsize_height); + + let data = match &src.inner.data { + DynamicData::U8(m) => { + let bv = CoreScalar::new( + num_traits::FromPrimitive::from_f64(border_value.v0).unwrap_or(0), + num_traits::FromPrimitive::from_f64(border_value.v1).unwrap_or(0), + num_traits::FromPrimitive::from_f64(border_value.v2).unwrap_or(0), + num_traits::FromPrimitive::from_f64(border_value.v3).unwrap_or(0), + ); + let result = warp_perspective(m, m_mat, dsize, interp, border, bv) + .map_err(|e| JsError::new(&format!("{e}")))?; + DynamicData::U8(result) + } + DynamicData::F32(m) => { + let bv = CoreScalar::new( + num_traits::FromPrimitive::from_f64(border_value.v0).unwrap_or(0.0), + num_traits::FromPrimitive::from_f64(border_value.v1).unwrap_or(0.0), + num_traits::FromPrimitive::from_f64(border_value.v2).unwrap_or(0.0), + num_traits::FromPrimitive::from_f64(border_value.v3).unwrap_or(0.0), + ); + let result = warp_perspective(m, m_mat, dsize, interp, border, bv) + .map_err(|e| JsError::new(&format!("{e}")))?; + DynamicData::F32(result) + } + DynamicData::F64(m) => { + let bv = CoreScalar::new( + border_value.v0, + border_value.v1, + border_value.v2, + border_value.v3, + ); + let result = warp_perspective(m, m_mat, dsize, interp, border, bv) + .map_err(|e| JsError::new(&format!("{e}")))?; + DynamicData::F64(result) + } + _ => return Err(JsError::new("warpPerspective: unsupported depth type")), + }; + + Ok(Mat { + inner: DynamicMatrix { data }, + }) +} + // --------------------------------------------------------------------------- // Blur / filter operations (require f32 depth) // --------------------------------------------------------------------------- @@ -2220,6 +2350,113 @@ pub fn rodrigues_wasm(src: &Mat, dst: &mut Mat) -> Result<(), JsError> { rodrigues(src_mat, dst_mat).map_err(|e| JsError::new(&format!("{e}"))) } +#[wasm_bindgen(js_name = "initUndistortRectifyMap")] +pub fn init_undistort_rectify_map_wasm( + camera_matrix: &Mat, + dist_coeffs: &Mat, + r: Option, + new_camera_matrix: &Mat, + size_width: i32, + size_height: i32, +) -> Result { + use purecv::calib3d::undistort::init_undistort_rectify_map; + use purecv::core::types::Size2i; + + let cam_mat = require_f64(camera_matrix, "initUndistortRectifyMap (camera_matrix)")?; + let dist = require_f64(dist_coeffs, "initUndistortRectifyMap (dist_coeffs)")?; + + let r_mat = match &r { + Some(mat) => Some(require_f64(mat, "initUndistortRectifyMap (R)")?), + None => None, + }; + + let new_cam_mat = require_f64( + new_camera_matrix, + "initUndistortRectifyMap (new_camera_matrix)", + )?; + + let (map1, map2) = init_undistort_rectify_map( + cam_mat, + dist, + r_mat, + new_cam_mat, + Size2i::new(size_width, size_height), + ) + .map_err(|e| JsError::new(&format!("{e}")))?; + + let obj = js_sys::Object::new(); + let js_map1 = Mat { + inner: purecv::core::dynamic::DynamicMatrix { + data: purecv::core::dynamic::DynamicData::F32(map1), + }, + }; + let js_map2 = Mat { + inner: purecv::core::dynamic::DynamicMatrix { + data: purecv::core::dynamic::DynamicData::F32(map2), + }, + }; + + js_sys::Reflect::set(&obj, &JsValue::from_str("map1"), &JsValue::from(js_map1)) + .map_err(|_| JsError::new("Failed to set map1"))?; + js_sys::Reflect::set(&obj, &JsValue::from_str("map2"), &JsValue::from(js_map2)) + .map_err(|_| JsError::new("Failed to set map2"))?; + + Ok(obj.into()) +} + +#[wasm_bindgen(js_name = "findFundamentalMat")] +pub fn find_fundamental_mat_wasm( + points1: &Point2fVector, + points2: &Point2fVector, + method: i32, + ransac_reproj_threshold: f64, + confidence: f64, + max_iters: usize, +) -> Result { + use purecv::calib3d::fundamental::{find_fundamental_mat, FundamentalMatMethod}; + + let method_enum = match method { + 2 => FundamentalMatMethod::FM_8POINT, + 8 => FundamentalMatMethod::FM_RANSAC, + _ => { + return Err(JsError::new(&format!( + "Unknown fundamental matrix method: {method}" + ))) + } + }; + + let mut mask_vec = Vec::new(); + let f_mat = find_fundamental_mat( + &points1.inner, + &points2.inner, + method_enum, + ransac_reproj_threshold, + confidence, + max_iters, + Some(&mut mask_vec), + ) + .map_err(|e| JsError::new(&format!("{e}")))?; + + let obj = js_sys::Object::new(); + let js_f = Mat { + inner: purecv::core::dynamic::DynamicMatrix { + data: purecv::core::dynamic::DynamicData::F64(f_mat), + }, + }; + js_sys::Reflect::set( + &obj, + &JsValue::from_str("fundamental"), + &JsValue::from(js_f), + ) + .map_err(|_| JsError::new("Failed to set fundamental"))?; + + let js_mask = js_sys::Uint8Array::from(mask_vec.as_slice()); + js_sys::Reflect::set(&obj, &JsValue::from_str("mask"), &js_mask) + .map_err(|_| JsError::new("Failed to set mask"))?; + + Ok(obj.into()) +} + // --------------------------------------------------------------------------- // features2d — KeyPoint, FAST, and ORB bindings // --------------------------------------------------------------------------- diff --git a/crates/wasm/www/example_calibration.html b/crates/wasm/www/example_calibration.html new file mode 100644 index 0000000..31038ca --- /dev/null +++ b/crates/wasm/www/example_calibration.html @@ -0,0 +1,340 @@ + + + + + + + PureCV - Calibration & Warping + + + +
+ +

Initializing PureCV WASM...

+
+ +
+
+

Calibration & Warping

+

Interactive Camera Lens Undistortion and Geometric Perspective Warp

+
+ +
+ + +
+ +
+ +
+

Controls

+ +
+

Click or drag custom image

+ +
+ + +
+
+ + +
+
+ + +
+
+ + +
+
+ + +
+
+ + +
+ +
+ + + +
+ + +
+
+
+

Source Image

+
+ +
+
+
+

Processed Result

+
+ +
+
+
+
+
+
+ + + + diff --git a/crates/wasm/www/example_calibration.js b/crates/wasm/www/example_calibration.js new file mode 100644 index 0000000..7cfe42a --- /dev/null +++ b/crates/wasm/www/example_calibration.js @@ -0,0 +1,403 @@ +/* + * example_calibration.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'; + +// --------------------------------------------------------------------------- +// State & Constants +// --------------------------------------------------------------------------- + +let cv = null; +let sourceImage = null; +let currentMode = 'undistort'; // 'undistort' or 'warp' + +const CANVAS_SIZE = 400; +const GRAB_RADIUS = 14; + +// Draggable corners for Perspective Warp mode +let corners = [ + { x: 80, y: 80 }, // TL + { x: 320, y: 70 }, // TR + { x: 330, y: 330 }, // BR + { x: 70, y: 320 } // BL +]; +let draggingIndex = -1; + +// Elements +const canvasSrc = document.getElementById('canvas-src'); +const ctxSrc = canvasSrc.getContext('2d', { willReadFrequently: true }); +const canvasDst = document.getElementById('canvas-dst'); +const ctxDst = canvasDst.getContext('2d', { willReadFrequently: true }); + +const tabUndistort = document.getElementById('tab-undistort'); +const tabWarp = document.getElementById('tab-warp'); +const panelUndistort = document.getElementById('panel-undistort'); +const panelWarp = document.getElementById('panel-warp'); +const labelCanvasSrc = document.getElementById('label-canvas-src'); + +const paramK1 = document.getElementById('param-k1'); +const paramK2 = document.getElementById('param-k2'); +const paramP1 = document.getElementById('param-p1'); +const paramP2 = document.getElementById('param-p2'); +const paramF = document.getElementById('param-f'); + +const valK1 = document.getElementById('val-k1'); +const valK2 = document.getElementById('val-k2'); +const valP1 = document.getElementById('val-p1'); +const valP2 = document.getElementById('val-p2'); +const valF = document.getElementById('val-f'); + +// --------------------------------------------------------------------------- +// Main logic +// --------------------------------------------------------------------------- + +async function start() { + try { + cv = await initWasm(); + document.getElementById('loader').classList.add('hidden'); + + // Load the default butterfly image + sourceImage = await loadImage('https://raw.githubusercontent.com/opencv/opencv/master/samples/data/butterfly.jpg'); + + setupEventListeners(); + updateSliderLabels(); + render(); + } catch (err) { + console.error("WASM Initialization failed:", err); + document.getElementById('loader').innerHTML = `

Error loading WASM: ${err.message}

`; + } +} + +function render() { + if (!sourceImage || !cv) return; + + if (currentMode === 'undistort') { + runUndistortion(); + } else { + runPerspectiveWarp(); + } +} + +// --------------------------------------------------------------------------- +// Mode 1: Camera Undistortion & Remap +// --------------------------------------------------------------------------- + +function runUndistortion() { + labelCanvasSrc.textContent = "Original Image"; + + // Setup source canvas with scaled image + const { width, height } = getScaledDimensions(sourceImage, CANVAS_SIZE); + canvasSrc.width = width; + canvasSrc.height = height; + ctxSrc.drawImage(sourceImage, 0, 0, width, height); + + // Convert source to PureCV Mat + const srcMat = canvasToMat(cv, canvasSrc, ctxSrc); + + // Read parameters + const k1 = parseFloat(paramK1.value); + const k2 = parseFloat(paramK2.value); + const p1 = parseFloat(paramP1.value); + const p2 = parseFloat(paramP2.value); + const f = parseFloat(paramF.value); + + try { + // Construct standard camera matrix + // [ fx, 0, cx ] + // [ 0, fy, cy ] + // [ 0, 0, 1 ] + const cx = width / 2.0; + const cy = height / 2.0; + const fx = width * f; + const fy = height * f; + const camMat = cv.Mat.fromF64Data(3, 3, 1, new Float64Array([ + fx, 0, cx, + 0, fy, cy, + 0, 0, 1.0 + ])); + + // Distortion coefficients: [k1, k2, p1, p2, k3, k4, k5, k6] + const distCoeffs = cv.Mat.fromF64Data(1, 5, 1, new Float64Array([ + k1, k2, p1, p2, 0.0 + ])); + + // Optional rectification (Identity) + const rMat = cv.Mat.fromF64Data(3, 3, 1, new Float64Array([ + 1.0, 0, 0, + 0, 1.0, 0, + 0, 0, 1.0 + ])); + + // New camera matrix (Same as original for simplicity) + const newCamMat = cv.Mat.fromF64Data(3, 3, 1, new Float64Array([ + fx, 0, cx, + 0, fy, cy, + 0, 0, 1.0 + ])); + + // Generate maps + const maps = cv.initUndistortRectifyMap(camMat, distCoeffs, rMat, newCamMat, width, height); + const map1 = maps.map1; + const map2 = maps.map2; + + // Apply remap (Linear interpolation, Constant border) + const dstMat = cv.remap(srcMat, map1, map2, 1, cv.BORDER_CONSTANT(), cv.Scalar.all(0)); + + // Render destination canvas + matToCanvas(dstMat, canvasDst); + + // Deallocate WASM matrices + camMat.free(); + distCoeffs.free(); + rMat.free(); + newCamMat.free(); + map1.free(); + map2.free(); + dstMat.free(); + } catch (e) { + console.error("Undistortion remapping failed:", e); + } + + srcMat.free(); +} + +// --------------------------------------------------------------------------- +// Mode 2: Perspective Warping +// --------------------------------------------------------------------------- + +function runPerspectiveWarp() { + labelCanvasSrc.textContent = "Source Image (Drag green corners)"; + + const { width, height } = getScaledDimensions(sourceImage, CANVAS_SIZE); + canvasSrc.width = width; + canvasSrc.height = height; + + // Draw the image first + ctxSrc.drawImage(sourceImage, 0, 0, width, height); + + // Draw lines connecting the corners + ctxSrc.strokeStyle = '#10b981'; + ctxSrc.lineWidth = 2; + ctxSrc.beginPath(); + ctxSrc.moveTo(corners[0].x, corners[0].y); + for (let i = 1; i < corners.length; i++) { + ctxSrc.lineTo(corners[i].x, corners[i].y); + } + ctxSrc.closePath(); + ctxSrc.stroke(); + + // Draw corners + corners.forEach((p, idx) => { + ctxSrc.fillStyle = (idx === draggingIndex) ? '#38bdf8' : '#10b981'; + ctxSrc.beginPath(); + ctxSrc.arc(p.x, p.y, 8, 0, 2 * Math.PI); + ctxSrc.fill(); + ctxSrc.lineWidth = 2; + ctxSrc.strokeStyle = '#ffffff'; + ctxSrc.stroke(); + }); + + // Run warp perspective + const srcMat = canvasToMat(cv, canvasSrc, ctxSrc); + + try { + const srcPts = new cv.Point2fVector(); + const dstPts = new cv.Point2fVector(); + + // Target quad is the entire output canvas + srcPts.push(0, 0); + srcPts.push(width, 0); + srcPts.push(width, height); + srcPts.push(0, height); + + // Source quad is the draggable polygon + corners.forEach(p => dstPts.push(p.x, p.y)); + + // Find homography (Source to Destination) + const hRes = cv.findHomography(srcPts, dstPts, 0, 3.0); + const hMat = hRes.homography; + + // Warp image onto the draggable quadrilateral + const dstMat = cv.warpPerspective(srcMat, hMat, width, height, 1, cv.BORDER_CONSTANT(), new cv.Scalar(0, 0, 0, 0)); + + matToCanvas(dstMat, canvasDst); + + srcPts.free(); + dstPts.free(); + hMat.free(); + hRes.mask.free(); + dstMat.free(); + } catch (e) { + console.error("Warp perspective failed:", e); + } + + srcMat.free(); +} + +// --------------------------------------------------------------------------- +// Interaction & Event Listeners +// --------------------------------------------------------------------------- + +function setupEventListeners() { + // Tab Switching + tabUndistort.onclick = () => { + tabUndistort.classList.add('active'); + tabWarp.classList.remove('active'); + panelUndistort.classList.remove('hidden'); + panelWarp.classList.add('hidden'); + currentMode = 'undistort'; + render(); + }; + + tabWarp.onclick = () => { + tabWarp.classList.add('active'); + tabUndistort.classList.remove('active'); + panelWarp.classList.remove('hidden'); + panelUndistort.classList.add('hidden'); + currentMode = 'warp'; + render(); + }; + + // Sliders + [paramK1, paramK2, paramP1, paramP2, paramF].forEach(slider => { + slider.oninput = () => { + updateSliderLabels(); + render(); + }; + }); + + // Resets + document.getElementById('reset-undistort').onclick = () => { + paramK1.value = -0.40; + paramK2.value = 0.10; + paramP1.value = 0.00; + paramP2.value = 0.00; + paramF.value = 0.95; + updateSliderLabels(); + render(); + }; + + document.getElementById('reset-warp').onclick = () => { + const { width, height } = getScaledDimensions(sourceImage, CANVAS_SIZE); + corners = [ + { x: width * 0.2, y: height * 0.2 }, + { x: width * 0.8, y: height * 0.18 }, + { x: width * 0.82, y: height * 0.82 }, + { x: width * 0.18, y: height * 0.8 } + ]; + render(); + }; + + // Custom Image Uploader + const fileInput = document.getElementById('file-input'); + const dropZone = document.getElementById('drop-zone'); + dropZone.onclick = () => fileInput.click(); + fileInput.onchange = async (e) => { + if (e.target.files.length > 0) { + const url = URL.createObjectURL(e.target.files[0]); + sourceImage = await loadImage(url); + + // Reset corners position for the new image size + const { width, height } = getScaledDimensions(sourceImage, CANVAS_SIZE); + corners = [ + { x: width * 0.2, y: height * 0.2 }, + { x: width * 0.8, y: height * 0.18 }, + { x: width * 0.82, y: height * 0.82 }, + { x: width * 0.18, y: height * 0.8 } + ]; + render(); + } + }; + + // Dragging corners interaction on canvasSrc + canvasSrc.addEventListener('pointerdown', (e) => { + if (currentMode !== 'warp') return; + const pos = getMousePos(e); + for (let i = 0; i < corners.length; i++) { + const dx = pos.x - corners[i].x; + const dy = pos.y - corners[i].y; + if (Math.sqrt(dx * dx + dy * dy) < GRAB_RADIUS) { + draggingIndex = i; + canvasSrc.setPointerCapture(e.pointerId); + render(); + break; + } + } + }); + + canvasSrc.addEventListener('pointermove', (e) => { + if (draggingIndex < 0) return; + const pos = getMousePos(e); + const { width, height } = getScaledDimensions(sourceImage, CANVAS_SIZE); + corners[draggingIndex].x = Math.max(10, Math.min(width - 10, pos.x)); + corners[draggingIndex].y = Math.max(10, Math.min(height - 10, pos.y)); + render(); + }); + + const stopDragging = () => { + if (draggingIndex >= 0) { + draggingIndex = -1; + render(); + } + }; + canvasSrc.addEventListener('pointerup', stopDragging); + canvasSrc.addEventListener('pointercancel', stopDragging); +} + +function getMousePos(e) { + const rect = canvasSrc.getBoundingClientRect(); + const { width, height } = getScaledDimensions(sourceImage, CANVAS_SIZE); + const scaleX = width / rect.width; + const scaleY = height / rect.height; + return { + x: (e.clientX - rect.left) * scaleX, + y: (e.clientY - rect.top) * scaleY, + }; +} + +function updateSliderLabels() { + valK1.textContent = parseFloat(paramK1.value).toFixed(2); + valK2.textContent = parseFloat(paramK2.value).toFixed(2); + valP1.textContent = parseFloat(paramP1.value).toFixed(2); + valP2.textContent = parseFloat(paramP2.value).toFixed(2); + valF.textContent = parseFloat(paramF.value).toFixed(2); +} + +// --------------------------------------------------------------------------- +// Bootstrap +// --------------------------------------------------------------------------- + +start(); diff --git a/crates/wasm/www/index.html b/crates/wasm/www/index.html index d0af3d7..b52b2c8 100644 --- a/crates/wasm/www/index.html +++ b/crates/wasm/www/index.html @@ -149,6 +149,11 @@

Optical Flow

Pose Estimation

Drag marker corners to compute camera pose (solvePnP, rodrigues) and homography in real-time.

+ + Calib3d & Imgproc +

Calibration & Warping

+

Interact with lens distortion parameters (initUndistortRectifyMap, remap) and warp perspective views in real-time.

+ diff --git a/examples/rectification.rs b/examples/rectification.rs new file mode 100644 index 0000000..81380cc --- /dev/null +++ b/examples/rectification.rs @@ -0,0 +1,152 @@ +/* + * rectification.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, Rgb}; +use purecv::core::Matrix; +use purecv::prelude::*; +use purecv::version; +use std::path::Path; +use std::time::Instant; + +fn main() -> Result<(), Box> { + println!("--- purecv Rectification & Warp Example ---"); + println!("purecv v{}", version::get_version()); + + // 1. Load the image + let img_path = "examples/data/butterfly.jpg"; + if !Path::new(img_path).exists() { + eprintln!("Error: {} not found. Run from the project root.", img_path); + return Ok(()); + } + + let img = image::open(img_path)?; + let (width, height) = img.dimensions(); + println!("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()); + + std::fs::create_dir_all("examples/data/out")?; + + // 2. Camera Undistortion + println!("Computing camera undistortion map..."); + let camera_matrix = Matrix::from_vec( + 3, + 3, + 1, + vec![ + width as f64 * 0.8, + 0.0, + width as f64 * 0.5, + 0.0, + height as f64 * 0.8, + height as f64 * 0.5, + 0.0, + 0.0, + 1.0, + ], + ); + // Strong barrel distortion (k1 = -0.5, k2 = 0.1) + let dist_coeffs = Matrix::from_vec(1, 4, 1, vec![-0.5, 0.1, 0.0, 0.0]); + let new_camera_matrix = camera_matrix.clone(); + + let t = Instant::now(); + let (map1, map2) = init_undistort_rectify_map( + &camera_matrix, + &dist_coeffs, + None, + &new_camera_matrix, + Size2i::new(width as i32, height as i32), + )?; + println!(" Map initialized in {:.2?}", t.elapsed()); + + println!("Remapping image (camera undistortion)..."); + let t_remap = Instant::now(); + let undistorted = remap( + &mat_rgb, + &map1, + &map2, + InterpolationFlags::Linear, + BorderTypes::Constant, + Scalar::new(0, 0, 0, 0), + )?; + println!(" Remapped in {:.2?}", t_remap.elapsed()); + save_matrix_rgb(&undistorted, "examples/data/out/butterfly_undistorted.png")?; + println!(" Saved: examples/data/out/butterfly_undistorted.png"); + + // 3. Perspective Warp + println!("Applying perspective warp..."); + // Dynamic perspective homography matrix (slight rotation, skew, and zoom) + let m = Matrix::from_vec( + 3, + 3, + 1, + vec![ + 0.8, + 0.1, + width as f64 * 0.1, + -0.1, + 0.8, + height as f64 * 0.1, + 0.0005, + 0.0005, + 1.0, + ], + ); + + let t_warp = Instant::now(); + let warped = warp_perspective( + &mat_rgb, + &m, + Size2i::new(width as i32, height as i32), + InterpolationFlags::Linear, + BorderTypes::Constant, + Scalar::new(0, 0, 0, 0), + )?; + println!(" Perspective warp completed in {:.2?}", t_warp.elapsed()); + save_matrix_rgb(&warped, "examples/data/out/butterfly_warped.png")?; + println!(" Saved: examples/data/out/butterfly_warped.png"); + + println!("\nRectification and Warp operations applied successfully!"); + Ok(()) +} + +fn save_matrix_rgb(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::ImageRgb8(img).save(filename) +} diff --git a/src/calib3d.rs b/src/calib3d.rs index 8f1775f..cef5ff7 100644 --- a/src/calib3d.rs +++ b/src/calib3d.rs @@ -54,10 +54,12 @@ //! //! [ocv]: https://docs.opencv.org/4.10.0/d9/d0c/group__calib3d.html +pub mod fundamental; pub mod geometry; pub mod homography; pub(crate) mod linalg; pub mod pose; +pub mod undistort; #[cfg(feature = "simd")] pub(crate) mod simd; @@ -69,6 +71,8 @@ mod tests; // Re-exports // --------------------------------------------------------------------------- +pub use fundamental::{find_fundamental_mat, FundamentalMatMethod}; pub use geometry::rodrigues; pub use homography::{find_homography, HomographyMethod}; pub use pose::{solve_pnp, solve_pnp_ransac, SolvePnPMethod}; +pub use undistort::init_undistort_rectify_map; diff --git a/src/calib3d/fundamental.rs b/src/calib3d/fundamental.rs new file mode 100644 index 0000000..67bc198 --- /dev/null +++ b/src/calib3d/fundamental.rs @@ -0,0 +1,378 @@ +/* + * fundamental.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 super::linalg::{mat3_mul, null_space_vector, svd_3x3, Lcg}; +use crate::core::error::{PureCvError, Result}; +use crate::core::types::Point2f; +use crate::core::Matrix; + +/// Method flags for computing the fundamental matrix. +#[allow(non_camel_case_types)] +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum FundamentalMatMethod { + FM_8POINT = 2, + FM_RANSAC = 8, +} + +/// Calculates a fundamental matrix from the corresponding points in two images. +/// +/// # Arguments +/// +/// * `points1` - Array of N points from the first image. +/// * `points2` - Array of N points from the second image. +/// * `method` - Method for computing the fundamental matrix (e.g., `FundamentalMatMethod::FM_RANSAC`). +/// * `ransac_reproj_threshold` - Parameter used only for RANSAC. It denotes the maximum +/// distance from a point to its epipolar line in pixels. +/// * `confidence` - Parameter used only for RANSAC. It denotes the confidence level +/// that the estimated matrix is correct. +/// * `max_iters` - Maximum number of iterations for RANSAC. +/// * `mask` - Optional output mask of inliers/outliers (1 for inliers, 0 for outliers). +/// +/// # Returns +/// +/// Returns a `Result>` containing the computed 3x3 fundamental matrix. +/// +/// # Errors +/// +/// Returns an error if: +/// * `points1` and `points2` have different lengths. +/// * There are fewer than 8 point correspondences. +/// * The RANSAC algorithm fails to find a valid fundamental matrix model. +/// * The null space vector computation fails (8-point algorithm). +/// +/// # Examples +/// +/// ``` +/// use purecv::core::types::Point2f; +/// use purecv::calib3d::fundamental::{find_fundamental_mat, FundamentalMatMethod}; +/// +/// let points1 = vec![ +/// Point2f::new(0.0, 0.0), Point2f::new(1.0, 0.0), +/// Point2f::new(0.0, 1.0), Point2f::new(1.0, 1.0), +/// Point2f::new(0.5, 0.5), Point2f::new(0.2, 0.8), +/// Point2f::new(0.8, 0.2), Point2f::new(0.3, 0.3), +/// ]; +/// let points2 = points1.clone(); // In practice, these would be matched points +/// +/// let mut mask = Vec::new(); +/// let f_mat = find_fundamental_mat( +/// &points1, +/// &points2, +/// FundamentalMatMethod::FM_8POINT, +/// 3.0, +/// 0.99, +/// 1000, +/// Some(&mut mask), +/// ).unwrap(); +/// ``` +pub fn find_fundamental_mat( + points1: &[Point2f], + points2: &[Point2f], + method: FundamentalMatMethod, + ransac_reproj_threshold: f64, + confidence: f64, + max_iters: usize, + mask: Option<&mut Vec>, +) -> Result> { + if points1.len() != points2.len() { + return Err(PureCvError::InvalidInput( + "points1 and points2 must have the same length".to_string(), + )); + } + if points1.len() < 8 { + return Err(PureCvError::InvalidInput( + "find_fundamental_mat requires at least 8 point correspondences".to_string(), + )); + } + + match method { + FundamentalMatMethod::FM_8POINT => { + let f = solve_8point(points1, points2)?; + if let Some(m) = mask { + m.clear(); + m.resize(points1.len(), 1); + } + Ok(f) + } + FundamentalMatMethod::FM_RANSAC => ransac_fundamental( + points1, + points2, + ransac_reproj_threshold, + confidence, + max_iters, + mask, + ), + } +} + +/// Normalized 8-Point Algorithm to solve for the fundamental matrix. +fn solve_8point(pts1: &[Point2f], pts2: &[Point2f]) -> Result> { + let n = pts1.len(); + let (cx1, cy1, s1) = compute_normalization(pts1); + let (cx2, cy2, s2) = compute_normalization(pts2); + + let mut a = vec![0.0f64; n * 9]; + for i in 0..n { + let u1 = s1 * (pts1[i].x as f64 - cx1); + let v1 = s1 * (pts1[i].y as f64 - cy1); + let u2 = s2 * (pts2[i].x as f64 - cx2); + let v2 = s2 * (pts2[i].y as f64 - cy2); + + a[i * 9] = u2 * u1; + a[i * 9 + 1] = u2 * v1; + a[i * 9 + 2] = u2; + a[i * 9 + 3] = v2 * u1; + a[i * 9 + 4] = v2 * v1; + a[i * 9 + 5] = v2; + a[i * 9 + 6] = u1; + a[i * 9 + 7] = v1; + a[i * 9 + 8] = 1.0; + } + + let f = null_space_vector(&a, n, 9); + if f.len() != 9 { + return Err(PureCvError::InternalError( + "Null space vector computation failed".to_string(), + )); + } + + let f_mat = [f[0], f[1], f[2], f[3], f[4], f[5], f[6], f[7], f[8]]; + + // Enforce rank-2 constraint: set the smallest singular value to 0 + let (u, sigma, vt) = svd_3x3(&f_mat); + let sigma_diag = [sigma[0], sigma[1], 0.0]; + + let mut u_sigma = [0.0f64; 9]; + for row in 0..3 { + u_sigma[row * 3] = u[row * 3] * sigma_diag[0]; + u_sigma[row * 3 + 1] = u[row * 3 + 1] * sigma_diag[1]; + u_sigma[row * 3 + 2] = 0.0; + } + + let f_constrained = mat3_mul(&u_sigma, &vt); + + // De-normalize: F = T2^T * F_constrained * T1 + let t1 = [s1, 0.0, -s1 * cx1, 0.0, s1, -s1 * cy1, 0.0, 0.0, 1.0]; + let t2_t = [s2, 0.0, 0.0, 0.0, s2, 0.0, -s2 * cx2, -s2 * cy2, 1.0]; + + let temp = mat3_mul(&t2_t, &f_constrained); + let f_final = mat3_mul(&temp, &t1); + + Ok(Matrix::from_vec(3, 3, 1, f_final.to_vec())) +} + +/// Compute center of mass and scale factor to normalize a set of 2D points. +fn compute_normalization(pts: &[Point2f]) -> (f64, f64, f64) { + let n = pts.len() as f64; + let mut cx = 0.0; + let mut cy = 0.0; + for p in pts { + cx += p.x as f64; + cy += p.y as f64; + } + cx /= n; + cy /= n; + + let mut mean_dist = 0.0; + for p in pts { + let dx = p.x as f64 - cx; + let dy = p.y as f64 - cy; + mean_dist += (dx * dx + dy * dy).sqrt(); + } + mean_dist /= n; + + let scale = if mean_dist > 1e-10 { + 2.0f64.sqrt() / mean_dist + } else { + 1.0 + }; + + (cx, cy, scale) +} + +/// Compute Sampson distance (first-order geometric error approximation) for fundamental matrix. +fn sampson_distance(pt1: Point2f, pt2: Point2f, f: &[f64; 9]) -> f64 { + let u = pt1.x as f64; + let v = pt1.y as f64; + let u_prime = pt2.x as f64; + let v_prime = pt2.y as f64; + + let f0 = f[0]; + let f1 = f[1]; + let f2 = f[2]; + let f3 = f[3]; + let f4 = f[4]; + let f5 = f[5]; + let f6 = f[6]; + let f7 = f[7]; + let f8 = f[8]; + + let e = u_prime * (f0 * u + f1 * v + f2) + + v_prime * (f3 * u + f4 * v + f5) + + (f6 * u + f7 * v + f8); + + let fx0 = f0 * u + f1 * v + f2; + let fx1 = f3 * u + f4 * v + f5; + + let ftx0 = f0 * u_prime + f3 * v_prime + f6; + let ftx1 = f1 * u_prime + f4 * v_prime + f7; + + let den = fx0 * fx0 + fx1 * fx1 + ftx0 * ftx0 + ftx1 * ftx1; + if den.abs() > 1e-10 { + (e * e) / den + } else { + 0.0 + } +} + +/// Sample `k` distinct indices from `[0, n)` without replacement. +fn sample_no_replace(rng: &mut Lcg, n: usize, k: usize) -> Vec { + let mut used = vec![false; n]; + let mut out = Vec::with_capacity(k); + while out.len() < k { + let i = rng.next_usize(n); + if !used[i] { + used[i] = true; + out.push(i); + } + } + out +} + +/// Robust RANSAC-based estimation of the fundamental matrix. +fn ransac_fundamental( + pts1: &[Point2f], + pts2: &[Point2f], + threshold: f64, + confidence: f64, + max_iters: usize, + mask: Option<&mut Vec>, +) -> Result> { + let n = pts1.len(); + let threshold_sq = threshold * threshold; + + let mut rng = Lcg::new(0x9e37_79b9_7f4a_7c15); + let mut best_count = 0; + let mut best_inliers = vec![0u8; n]; + let mut best_f: Option<[f64; 9]> = None; + + for _ in 0..max_iters { + let idx = sample_no_replace(&mut rng, n, 8); + let sample1: Vec = idx.iter().map(|&i| pts1[i]).collect(); + let sample2: Vec = idx.iter().map(|&i| pts2[i]).collect(); + + let f = match solve_8point(&sample1, &sample2) { + Ok(mat) => mat, + Err(_) => continue, + }; + + let f_arr: [f64; 9] = match f.data.clone().try_into() { + Ok(arr) => arr, + Err(_) => continue, + }; + + let mut count = 0; + let mut inliers = vec![0u8; n]; + for i in 0..n { + let err = sampson_distance(pts1[i], pts2[i], &f_arr); + if err <= threshold_sq { + inliers[i] = 1; + count += 1; + } + } + + if count > best_count { + best_count = count; + best_inliers = inliers; + best_f = Some(f_arr); + } + + // Early exit check + if best_count >= (n as f64 * confidence) as usize { + break; + } + } + + let best_f_arr = match best_f { + Some(arr) => arr, + None => { + return Err(PureCvError::InvalidInput( + "RANSAC: No valid fundamental matrix model found".to_string(), + )) + } + }; + + if best_count < 8 { + return Err(PureCvError::InvalidInput( + "RANSAC: Less than 8 inliers found".to_string(), + )); + } + + // Refit using all inliers + let inliers_pts1: Vec = pts1 + .iter() + .zip(best_inliers.iter()) + .filter_map(|(&p, &m)| if m == 1 { Some(p) } else { None }) + .collect(); + let inliers_pts2: Vec = pts2 + .iter() + .zip(best_inliers.iter()) + .filter_map(|(&p, &m)| if m == 1 { Some(p) } else { None }) + .collect(); + + let refined = if inliers_pts1.len() >= 8 { + match solve_8point(&inliers_pts1, &inliers_pts2) { + Ok(mat) => mat, + Err(_) => Matrix::from_vec(3, 3, 1, best_f_arr.to_vec()), + } + } else { + Matrix::from_vec(3, 3, 1, best_f_arr.to_vec()) + }; + + if let Some(m) = mask { + m.clear(); + m.resize(n, 0); + let refined_arr: [f64; 9] = refined.data.clone().try_into().unwrap_or(best_f_arr); + for i in 0..n { + let err = sampson_distance(pts1[i], pts2[i], &refined_arr); + if err <= threshold_sq { + m[i] = 1; + } + } + } + + Ok(refined) +} diff --git a/src/calib3d/tests.rs b/src/calib3d/tests.rs index ab9770b..b353dfe 100644 --- a/src/calib3d/tests.rs +++ b/src/calib3d/tests.rs @@ -408,4 +408,161 @@ mod calib3d_tests { tvec.data[2] ); } + + #[test] + fn test_init_undistort_rectify_map_identity() { + use crate::calib3d::init_undistort_rectify_map; + use crate::core::types::Size2i; + let cam = Matrix::from_vec(3, 3, 1, vec![1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]); + let dist = Matrix::from_vec(1, 4, 1, vec![0.0, 0.0, 0.0, 0.0]); + let (map1, map2) = + init_undistort_rectify_map(&cam, &dist, None, &cam, Size2i::new(2, 2)).unwrap(); + + assert_eq!(map1.data, vec![0.0f32, 1.0, 0.0, 1.0]); + assert_eq!(map2.data, vec![0.0f32, 0.0, 1.0, 1.0]); + } + + fn generate_synthetic_fundamental_data() -> (Vec, Vec) { + use crate::calib3d::rodrigues; + // 12 synthetic 3D points in general positions (not coplanar) + let pts3d = vec![ + [0.1, -0.2, 2.5], + [-0.3, 0.4, 3.0], + [0.5, 0.2, 2.0], + [-0.2, -0.5, 3.5], + [0.4, -0.4, 2.8], + [-0.1, 0.3, 2.2], + [0.3, 0.5, 3.2], + [-0.4, -0.2, 2.7], + [0.2, -0.1, 2.4], + [-0.2, 0.1, 3.1], + [0.1, 0.4, 2.9], + [-0.5, 0.5, 2.6], + ]; + + // Camera 1: Identity projection + let mut pts1 = Vec::new(); + for p in &pts3d { + pts1.push(Point2f { + x: (p[0] / p[2]) as f32, + y: (p[1] / p[2]) as f32, + }); + } + + // Camera 2: Rotated using a general rotation vector and translated + let rvec = Matrix::from_vec(3, 1, 1, vec![0.15, -0.2, 0.1]); + let mut rmat = Matrix::::new(3, 3, 1); + rodrigues(&rvec, &mut rmat).unwrap(); + + let tx = 0.3; + let ty = -0.2; + let tz = 0.5; + + let mut pts2 = Vec::new(); + for p in &pts3d { + let x = p[0]; + let y = p[1]; + let z = p[2]; + let x2 = rmat.data[0] * x + rmat.data[1] * y + rmat.data[2] * z + tx; + let y2 = rmat.data[3] * x + rmat.data[4] * y + rmat.data[5] * z + ty; + let z2 = rmat.data[6] * x + rmat.data[7] * y + rmat.data[8] * z + tz; + pts2.push(Point2f { + x: (x2 / z2) as f32, + y: (y2 / z2) as f32, + }); + } + + (pts1, pts2) + } + + #[test] + fn test_find_fundamental_mat_8point() { + use crate::calib3d::{find_fundamental_mat, FundamentalMatMethod}; + let (pts1, pts2) = generate_synthetic_fundamental_data(); + + let mut mask = Vec::new(); + let f = find_fundamental_mat( + &pts1, + &pts2, + FundamentalMatMethod::FM_8POINT, + 1.0, + 0.99, + 1000, + Some(&mut mask), + ) + .unwrap(); + + assert_eq!(f.rows, 3); + assert_eq!(f.cols, 3); + assert_eq!(mask.len(), 12); + assert!(mask.iter().all(|&m| m == 1)); + + // Verify epipolar constraint x'^T * F * x = 0 (tolerance 1e-4) + for i in 0..12 { + let u1 = pts1[i].x as f64; + let v1 = pts1[i].y as f64; + let u2 = pts2[i].x as f64; + let v2 = pts2[i].y as f64; + + let x2 = [u2, v2, 1.0]; + let x1 = [u1, v1, 1.0]; + + let fx1 = [ + f.data[0] * x1[0] + f.data[1] * x1[1] + f.data[2] * x1[2], + f.data[3] * x1[0] + f.data[4] * x1[1] + f.data[5] * x1[2], + f.data[6] * x1[0] + f.data[7] * x1[1] + f.data[8] * x1[2], + ]; + let err = x2[0] * fx1[0] + x2[1] * fx1[1] + x2[2] * fx1[2]; + assert!(err.abs() < 1e-4, "Error at point {i} is {err}"); + } + } + + #[test] + fn test_find_fundamental_mat_ransac() { + use crate::calib3d::{find_fundamental_mat, FundamentalMatMethod}; + let (pts1, mut pts2) = generate_synthetic_fundamental_data(); + + // Add 2 outliers (perturb the last two points) + pts2[10].y += 5.0f32; + pts2[11].y += 5.0f32; + + let mut mask = Vec::new(); + let f = find_fundamental_mat( + &pts1, + &pts2, + FundamentalMatMethod::FM_RANSAC, + 0.01, + 0.99, + 1000, + Some(&mut mask), + ) + .unwrap(); + + // Verify epipolar constraint for inliers (tolerance 1e-3) + for i in 0..10 { + let u1 = pts1[i].x as f64; + let v1 = pts1[i].y as f64; + let u2 = pts2[i].x as f64; + let v2 = pts2[i].y as f64; + + let x2 = [u2, v2, 1.0]; + let x1 = [u1, v1, 1.0]; + + let fx1 = [ + f.data[0] * x1[0] + f.data[1] * x1[1] + f.data[2] * x1[2], + f.data[3] * x1[0] + f.data[4] * x1[1] + f.data[5] * x1[2], + f.data[6] * x1[0] + f.data[7] * x1[1] + f.data[8] * x1[2], + ]; + let err = x2[0] * fx1[0] + x2[1] * fx1[1] + x2[2] * fx1[2]; + assert!(err.abs() < 1e-3, "Error at inlier {i} is {err}"); + } + + assert_eq!(mask.len(), 12); + // Clean points should be inliers (1), outliers should be 0 + for (i, &val) in mask.iter().enumerate().take(10) { + assert_eq!(val, 1, "Expected inlier at {i}"); + } + assert_eq!(mask[10], 0, "Expected outlier at 10"); + assert_eq!(mask[11], 0, "Expected outlier at 11"); + } } diff --git a/src/calib3d/undistort.rs b/src/calib3d/undistort.rs new file mode 100644 index 0000000..16ddd06 --- /dev/null +++ b/src/calib3d/undistort.rs @@ -0,0 +1,268 @@ +/* + * undistort.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 crate::core::error::{PureCvError, Result}; +use crate::core::types::Size2i; +use crate::core::Matrix; + +#[cfg(feature = "parallel")] +use rayon::prelude::*; + +/// Computes the undistortion and rectification transformation map. +/// +/// The function computes the joint undistortion and rectification transformation map +/// and stores the result in `map1` and `map2`. +/// +/// # Arguments +/// +/// * `camera_matrix` - Input camera matrix (3x3). +/// * `dist_coeffs` - Input vector of distortion coefficients (k1, k2, p1, p2, k3, k4, k5, k6). +/// * `r` - Optional rectification transformation in object space (3x3 matrix). If `None`, identity is used. +/// * `new_camera_matrix` - New camera matrix (3x3). +/// * `size` - Undistorted image size. +/// +/// # Returns +/// +/// Returns a `Result<(Matrix, Matrix)>` containing the computed maps `map1` and `map2`. +/// +/// # Errors +/// +/// Returns an error if: +/// * `camera_matrix`, `new_camera_matrix`, or `r` (if provided) are not 3x3 single-channel matrices. +/// +/// # Examples +/// +/// ``` +/// use purecv::core::{Matrix, types::Size2i}; +/// use purecv::calib3d::undistort::init_undistort_rectify_map; +/// +/// let mut camera_matrix = Matrix::::new(3, 3, 1); +/// camera_matrix.data[0] = 800.0; // fx +/// camera_matrix.data[4] = 800.0; // fy +/// camera_matrix.data[2] = 320.0; // cx +/// camera_matrix.data[5] = 240.0; // cy +/// camera_matrix.data[8] = 1.0; +/// +/// let dist_coeffs = Matrix::::new(1, 5, 1); // e.g., k1, k2, p1, p2, k3 +/// let size = Size2i::new(640, 480); +/// +/// let (map1, map2) = init_undistort_rectify_map( +/// &camera_matrix, +/// &dist_coeffs, +/// None, +/// &camera_matrix, +/// size +/// ).unwrap(); +/// ``` +pub fn init_undistort_rectify_map( + camera_matrix: &Matrix, + dist_coeffs: &Matrix, + r: Option<&Matrix>, + new_camera_matrix: &Matrix, + size: Size2i, +) -> Result<(Matrix, Matrix)> { + if camera_matrix.rows != 3 || camera_matrix.cols != 3 || camera_matrix.channels != 1 { + return Err(PureCvError::InvalidDimensions( + "camera_matrix must be 3x3 single-channel".to_string(), + )); + } + if new_camera_matrix.rows != 3 || new_camera_matrix.cols != 3 || new_camera_matrix.channels != 1 + { + return Err(PureCvError::InvalidDimensions( + "new_camera_matrix must be 3x3 single-channel".to_string(), + )); + } + if let Some(r_mat) = r { + if r_mat.rows != 3 || r_mat.cols != 3 || r_mat.channels != 1 { + return Err(PureCvError::InvalidDimensions( + "R must be 3x3 single-channel".to_string(), + )); + } + } + + let fx = camera_matrix.data[0]; + let cx = camera_matrix.data[2]; + let fy = camera_matrix.data[4]; + let cy = camera_matrix.data[5]; + + let fx_prime = new_camera_matrix.data[0]; + let cx_prime = new_camera_matrix.data[2]; + let fy_prime = new_camera_matrix.data[4]; + let cy_prime = new_camera_matrix.data[5]; + + let mut k1 = 0.0; + let mut k2 = 0.0; + let mut p1 = 0.0; + let mut p2 = 0.0; + let mut k3 = 0.0; + let mut k4 = 0.0; + let mut k5 = 0.0; + let mut k6 = 0.0; + + let len = dist_coeffs.data.len(); + if len >= 1 { + k1 = dist_coeffs.data[0]; + } + if len >= 2 { + k2 = dist_coeffs.data[1]; + } + if len >= 3 { + p1 = dist_coeffs.data[2]; + } + if len >= 4 { + p2 = dist_coeffs.data[3]; + } + if len >= 5 { + k3 = dist_coeffs.data[4]; + } + if len >= 6 { + k4 = dist_coeffs.data[5]; + } + if len >= 7 { + k5 = dist_coeffs.data[6]; + } + if len >= 8 { + k6 = dist_coeffs.data[7]; + } + + let (r00, r01, r02, r10, r11, r12, r20, r21, r22) = match r { + Some(r_mat) => ( + r_mat.data[0], + r_mat.data[3], + r_mat.data[6], // Row 0 of transpose (Col 0 of R) + r_mat.data[1], + r_mat.data[4], + r_mat.data[7], // Row 1 of transpose (Col 1 of R) + r_mat.data[2], + r_mat.data[5], + r_mat.data[8], // Row 2 of transpose (Col 2 of R) + ), + None => (1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0), + }; + + let mut map1 = Matrix::::new(size.height as usize, size.width as usize, 1); + let mut map2 = Matrix::::new(size.height as usize, size.width as usize, 1); + + #[cfg(feature = "parallel")] + { + map1.data + .par_chunks_exact_mut(size.width as usize) + .zip(map2.data.par_chunks_exact_mut(size.width as usize)) + .enumerate() + .for_each(|(v, (map1_row, map2_row))| { + let y_val = v as f64; + for u in 0..size.width as usize { + let x_val = u as f64; + let x = (x_val - cx_prime) / fx_prime; + let y = (y_val - cy_prime) / fy_prime; + + let rx = r00 * x + r01 * y + r02; + let ry = r10 * x + r11 * y + r12; + let rz = r20 * x + r21 * y + r22; + + let (xp, yp) = if rz.abs() > 1e-10 { + (rx / rz, ry / rz) + } else { + (x, y) + }; + + let r2 = xp * xp + yp * yp; + let radial_num = 1.0 + k1 * r2 + k2 * r2 * r2 + k3 * r2 * r2 * r2; + let radial_den = 1.0 + k4 * r2 + k5 * r2 * r2 + k6 * r2 * r2 * r2; + let radial = if radial_den.abs() > 1e-10 { + radial_num / radial_den + } else { + radial_num + }; + + let dx = 2.0 * p1 * xp * yp + p2 * (r2 + 2.0 * xp * xp); + let dy = p1 * (r2 + 2.0 * yp * yp) + 2.0 * p2 * xp * yp; + + let x_dist = xp * radial + dx; + let y_dist = yp * radial + dy; + + map1_row[u] = (fx * x_dist + cx) as f32; + map2_row[u] = (fy * y_dist + cy) as f32; + } + }); + } + + #[cfg(not(feature = "parallel"))] + { + let width = size.width as usize; + let height = size.height as usize; + for v in 0..height { + let map1_row = &mut map1.data[v * width..(v + 1) * width]; + let map2_row = &mut map2.data[v * width..(v + 1) * width]; + let y_val = v as f64; + for u in 0..width { + let x_val = u as f64; + let x = (x_val - cx_prime) / fx_prime; + let y = (y_val - cy_prime) / fy_prime; + + let rx = r00 * x + r01 * y + r02; + let ry = r10 * x + r11 * y + r12; + let rz = r20 * x + r21 * y + r22; + + let (xp, yp) = if rz.abs() > 1e-10 { + (rx / rz, ry / rz) + } else { + (x, y) + }; + + let r2 = xp * xp + yp * yp; + let radial_num = 1.0 + k1 * r2 + k2 * r2 * r2 + k3 * r2 * r2 * r2; + let radial_den = 1.0 + k4 * r2 + k5 * r2 * r2 + k6 * r2 * r2 * r2; + let radial = if radial_den.abs() > 1e-10 { + radial_num / radial_den + } else { + radial_num + }; + + let dx = 2.0 * p1 * xp * yp + p2 * (r2 + 2.0 * xp * xp); + let dy = p1 * (r2 + 2.0 * yp * yp) + 2.0 * p2 * xp * yp; + + let x_dist = xp * radial + dx; + let y_dist = yp * radial + dy; + + map1_row[u] = (fx * x_dist + cx) as f32; + map2_row[u] = (fy * y_dist + cy) as f32; + } + } + } + + Ok((map1, map2)) +} diff --git a/src/core/simd.rs b/src/core/simd.rs index 181c432..b20133d 100644 --- a/src/core/simd.rs +++ b/src/core/simd.rs @@ -52,6 +52,8 @@ //! WASM `simd128`). No target-specific `#[cfg(target_arch)]` is needed — //! the same source compiles for native and `wasm32` targets. +use crate::core::types::{BorderTypes, Scalar}; + // --------------------------------------------------------------------------- // SimdElement trait — compile-time type routing (always available) // --------------------------------------------------------------------------- @@ -211,6 +213,38 @@ pub trait SimdElement: Copy + Send + Sync + 'static { fn simd_min_max_col(_dst: &mut [Self], _rows: &[&[Self]], _is_erode: bool) -> bool { false } + + /// Remap a row with bilinear interpolation. + #[allow(clippy::too_many_arguments)] + fn simd_remap_bilinear_row( + _dst_row: &mut [Self], + _src: &[Self], + _src_cols: usize, + _src_rows: usize, + _channels: usize, + _map1_row: &[f32], + _map2_row: &[f32], + _border_mode: BorderTypes, + _border_value: Scalar, + ) -> bool { + false + } + + /// Remap a row with nearest neighbor interpolation. + #[allow(clippy::too_many_arguments)] + fn simd_remap_nearest_row( + _dst_row: &mut [Self], + _src: &[Self], + _src_cols: usize, + _src_rows: usize, + _channels: usize, + _map1_row: &[f32], + _map2_row: &[f32], + _border_mode: BorderTypes, + _border_value: Scalar, + ) -> bool { + false + } } // --------------------------------------------------------------------------- @@ -241,6 +275,8 @@ impl SimdElement for u64 {} #[cfg(feature = "simd")] mod simd_impls { use super::SimdElement; + use crate::core::types::{BorderTypes, Scalar}; + use crate::core::utils::border_interpolate; use pulp::Arch; // ----------------------------------------------------------------------- @@ -536,6 +572,123 @@ mod simd_impls { }); true } + + fn simd_remap_bilinear_row( + dst_row: &mut [Self], + src: &[Self], + src_cols: usize, + src_rows: usize, + channels: usize, + map1_row: &[f32], + map2_row: &[f32], + border_mode: BorderTypes, + border_value: Scalar, + ) -> bool { + let arch = pulp::Arch::new(); + arch.dispatch(|| { + let dst_cols = dst_row.len() / channels; + for col in 0..dst_cols { + let x = map1_row[col]; + let y = map2_row[col]; + + let x1 = x.floor() as i32; + let y1 = y.floor() as i32; + let x2 = x1 + 1; + let y2 = y1 + 1; + + let wx = x - x.floor(); + let wy = y - y.floor(); + + let w00 = (1.0 - wx) * (1.0 - wy); + let w10 = wx * (1.0 - wy); + let w01 = (1.0 - wx) * wy; + let w11 = wx * wy; + + let dst_idx = col * channels; + + if x1 >= 0 && x2 < src_cols as i32 && y1 >= 0 && y2 < src_rows as i32 { + let idx00 = (y1 as usize * src_cols + x1 as usize) * channels; + let idx10 = (y1 as usize * src_cols + x2 as usize) * channels; + let idx01 = (y2 as usize * src_cols + x1 as usize) * channels; + let idx11 = (y2 as usize * src_cols + x2 as usize) * channels; + + for c in 0..channels { + let v00 = src[idx00 + c]; + let v10 = src[idx10 + c]; + let v01 = src[idx01 + c]; + let v11 = src[idx11 + c]; + dst_row[dst_idx + c] = w00 * v00 + w10 * v10 + w01 * v01 + w11 * v11; + } + } else { + for c in 0..channels { + let get_pixel = |ix: i32, iy: i32| -> f32 { + let ix_interp = + border_interpolate(ix, src_cols as i32, border_mode); + let iy_interp = + border_interpolate(iy, src_rows as i32, border_mode); + if ix_interp >= 0 && iy_interp >= 0 { + src[(iy_interp as usize * src_cols + ix_interp as usize) + * channels + + c] + } else { + border_value.v[c] + } + }; + + let v00 = get_pixel(x1, y1); + let v10 = get_pixel(x2, y1); + let v01 = get_pixel(x1, y2); + let v11 = get_pixel(x2, y2); + dst_row[dst_idx + c] = w00 * v00 + w10 * v10 + w01 * v01 + w11 * v11; + } + } + } + }); + true + } + + fn simd_remap_nearest_row( + dst_row: &mut [Self], + src: &[Self], + src_cols: usize, + src_rows: usize, + channels: usize, + map1_row: &[f32], + map2_row: &[f32], + border_mode: BorderTypes, + border_value: Scalar, + ) -> bool { + let arch = pulp::Arch::new(); + arch.dispatch(|| { + let dst_cols = dst_row.len() / channels; + for col in 0..dst_cols { + let x = map1_row[col]; + let y = map2_row[col]; + let ix = x.round() as i32; + let iy = y.round() as i32; + let dst_idx = col * channels; + + if ix >= 0 && ix < src_cols as i32 && iy >= 0 && iy < src_rows as i32 { + let src_idx = (iy as usize * src_cols + ix as usize) * channels; + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&src[src_idx..src_idx + channels]); + } else { + let ix_interp = border_interpolate(ix, src_cols as i32, border_mode); + let iy_interp = border_interpolate(iy, src_rows as i32, border_mode); + if ix_interp >= 0 && iy_interp >= 0 { + let src_idx = + (iy_interp as usize * src_cols + ix_interp as usize) * channels; + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&src[src_idx..src_idx + channels]); + } else { + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&border_value.v[..channels]); + } + } + } + }); + true + } } // ----------------------------------------------------------------------- @@ -826,6 +979,123 @@ mod simd_impls { }); true } + + fn simd_remap_bilinear_row( + dst_row: &mut [Self], + src: &[Self], + src_cols: usize, + src_rows: usize, + channels: usize, + map1_row: &[f32], + map2_row: &[f32], + border_mode: BorderTypes, + border_value: Scalar, + ) -> bool { + let arch = pulp::Arch::new(); + arch.dispatch(|| { + let dst_cols = dst_row.len() / channels; + for col in 0..dst_cols { + let x = map1_row[col] as f64; + let y = map2_row[col] as f64; + + let x1 = x.floor() as i32; + let y1 = y.floor() as i32; + let x2 = x1 + 1; + let y2 = y1 + 1; + + let wx = x - x.floor(); + let wy = y - y.floor(); + + let w00 = (1.0 - wx) * (1.0 - wy); + let w10 = wx * (1.0 - wy); + let w01 = (1.0 - wx) * wy; + let w11 = wx * wy; + + let dst_idx = col * channels; + + if x1 >= 0 && x2 < src_cols as i32 && y1 >= 0 && y2 < src_rows as i32 { + let idx00 = (y1 as usize * src_cols + x1 as usize) * channels; + let idx10 = (y1 as usize * src_cols + x2 as usize) * channels; + let idx01 = (y2 as usize * src_cols + x1 as usize) * channels; + let idx11 = (y2 as usize * src_cols + x2 as usize) * channels; + + for c in 0..channels { + let v00 = src[idx00 + c]; + let v10 = src[idx10 + c]; + let v01 = src[idx01 + c]; + let v11 = src[idx11 + c]; + dst_row[dst_idx + c] = w00 * v00 + w10 * v10 + w01 * v01 + w11 * v11; + } + } else { + for c in 0..channels { + let get_pixel = |ix: i32, iy: i32| -> f64 { + let ix_interp = + border_interpolate(ix, src_cols as i32, border_mode); + let iy_interp = + border_interpolate(iy, src_rows as i32, border_mode); + if ix_interp >= 0 && iy_interp >= 0 { + src[(iy_interp as usize * src_cols + ix_interp as usize) + * channels + + c] + } else { + border_value.v[c] + } + }; + + let v00 = get_pixel(x1, y1); + let v10 = get_pixel(x2, y1); + let v01 = get_pixel(x1, y2); + let v11 = get_pixel(x2, y2); + dst_row[dst_idx + c] = w00 * v00 + w10 * v10 + w01 * v01 + w11 * v11; + } + } + } + }); + true + } + + fn simd_remap_nearest_row( + dst_row: &mut [Self], + src: &[Self], + src_cols: usize, + src_rows: usize, + channels: usize, + map1_row: &[f32], + map2_row: &[f32], + border_mode: BorderTypes, + border_value: Scalar, + ) -> bool { + let arch = pulp::Arch::new(); + arch.dispatch(|| { + let dst_cols = dst_row.len() / channels; + for col in 0..dst_cols { + let x = map1_row[col]; + let y = map2_row[col]; + let ix = x.round() as i32; + let iy = y.round() as i32; + let dst_idx = col * channels; + + if ix >= 0 && ix < src_cols as i32 && iy >= 0 && iy < src_rows as i32 { + let src_idx = (iy as usize * src_cols + ix as usize) * channels; + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&src[src_idx..src_idx + channels]); + } else { + let ix_interp = border_interpolate(ix, src_cols as i32, border_mode); + let iy_interp = border_interpolate(iy, src_rows as i32, border_mode); + if ix_interp >= 0 && iy_interp >= 0 { + let src_idx = + (iy_interp as usize * src_cols + ix_interp as usize) * channels; + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&src[src_idx..src_idx + channels]); + } else { + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&border_value.v[..channels]); + } + } + } + }); + true + } } // ----------------------------------------------------------------------- @@ -1064,6 +1334,125 @@ mod simd_impls { }); true } + + fn simd_remap_bilinear_row( + dst_row: &mut [Self], + src: &[Self], + src_cols: usize, + src_rows: usize, + channels: usize, + map1_row: &[f32], + map2_row: &[f32], + border_mode: BorderTypes, + border_value: Scalar, + ) -> bool { + let arch = pulp::Arch::new(); + arch.dispatch(|| { + let dst_cols = dst_row.len() / channels; + for col in 0..dst_cols { + let x = map1_row[col]; + let y = map2_row[col]; + + let x1 = x.floor() as i32; + let y1 = y.floor() as i32; + let x2 = x1 + 1; + let y2 = y1 + 1; + + let wx = x - x.floor(); + let wy = y - y.floor(); + + let w00 = (1.0 - wx) * (1.0 - wy); + let w10 = wx * (1.0 - wy); + let w01 = (1.0 - wx) * wy; + let w11 = wx * wy; + + let dst_idx = col * channels; + + if x1 >= 0 && x2 < src_cols as i32 && y1 >= 0 && y2 < src_rows as i32 { + let idx00 = (y1 as usize * src_cols + x1 as usize) * channels; + let idx10 = (y1 as usize * src_cols + x2 as usize) * channels; + let idx01 = (y2 as usize * src_cols + x1 as usize) * channels; + let idx11 = (y2 as usize * src_cols + x2 as usize) * channels; + + for c in 0..channels { + let v00 = src[idx00 + c] as f32; + let v10 = src[idx10 + c] as f32; + let v01 = src[idx01 + c] as f32; + let v11 = src[idx11 + c] as f32; + let val = w00 * v00 + w10 * v10 + w01 * v01 + w11 * v11; + dst_row[dst_idx + c] = val.clamp(0.0, 255.0).round() as u8; + } + } else { + for c in 0..channels { + let get_pixel = |ix: i32, iy: i32| -> f32 { + let ix_interp = + border_interpolate(ix, src_cols as i32, border_mode); + let iy_interp = + border_interpolate(iy, src_rows as i32, border_mode); + if ix_interp >= 0 && iy_interp >= 0 { + src[(iy_interp as usize * src_cols + ix_interp as usize) + * channels + + c] as f32 + } else { + border_value.v[c] as f32 + } + }; + + let v00 = get_pixel(x1, y1); + let v10 = get_pixel(x2, y1); + let v01 = get_pixel(x1, y2); + let v11 = get_pixel(x2, y2); + let val = w00 * v00 + w10 * v10 + w01 * v01 + w11 * v11; + dst_row[dst_idx + c] = val.clamp(0.0, 255.0).round() as u8; + } + } + } + }); + true + } + + fn simd_remap_nearest_row( + dst_row: &mut [Self], + src: &[Self], + src_cols: usize, + src_rows: usize, + channels: usize, + map1_row: &[f32], + map2_row: &[f32], + border_mode: BorderTypes, + border_value: Scalar, + ) -> bool { + let arch = pulp::Arch::new(); + arch.dispatch(|| { + let dst_cols = dst_row.len() / channels; + for col in 0..dst_cols { + let x = map1_row[col]; + let y = map2_row[col]; + let ix = x.round() as i32; + let iy = y.round() as i32; + let dst_idx = col * channels; + + if ix >= 0 && ix < src_cols as i32 && iy >= 0 && iy < src_rows as i32 { + let src_idx = (iy as usize * src_cols + ix as usize) * channels; + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&src[src_idx..src_idx + channels]); + } else { + let ix_interp = border_interpolate(ix, src_cols as i32, border_mode); + let iy_interp = border_interpolate(iy, src_rows as i32, border_mode); + if ix_interp >= 0 && iy_interp >= 0 { + let src_idx = + (iy_interp as usize * src_cols + ix_interp as usize) * channels; + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&src[src_idx..src_idx + channels]); + } else { + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&border_value.v[..channels]); + } + } + } + }); + true + } } } diff --git a/src/imgproc.rs b/src/imgproc.rs index 0f4cd58..d9482d8 100644 --- a/src/imgproc.rs +++ b/src/imgproc.rs @@ -39,6 +39,7 @@ pub mod derivatives; pub mod edge; pub mod feature; pub mod filter; +pub mod geometric; pub mod hough; pub mod morph; pub mod pyramid; @@ -60,6 +61,7 @@ pub use derivatives::*; pub use edge::*; pub use feature::*; pub use filter::*; +pub use geometric::{remap, warp_perspective, InterpolationFlags}; 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/geometric.rs b/src/imgproc/geometric.rs new file mode 100644 index 0000000..f4742a8 --- /dev/null +++ b/src/imgproc/geometric.rs @@ -0,0 +1,543 @@ +/* + * geometric.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 crate::core::arithm::{invert, DecompTypes}; +use crate::core::error::{PureCvError, Result}; +use crate::core::simd::SimdElement; +use crate::core::types::{BorderTypes, Scalar, Size2i}; +use crate::core::utils::border_interpolate; +use crate::core::Matrix; +use num_traits::{FromPrimitive, ToPrimitive}; + +#[cfg(feature = "parallel")] +use rayon::prelude::*; + +/// Interpolation methods. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum InterpolationFlags { + Nearest = 0, + Linear = 1, +} + +/// Applies a generic geometrical transformation to an image. +/// +/// The function `remap` transforms the source image using the specified map: +/// `dst(x,y) = src(map1(x,y), map2(x,y))` +/// +/// # Arguments +/// +/// * `src` - Source image as a `Matrix`. +/// * `map1` - The first map of either `(x,y)` points or just `x` values, with type `f32`. +/// * `map2` - The second map of `y` values with type `f32`. +/// * `interpolation` - Interpolation method to use (e.g., `InterpolationFlags::Linear`). +/// * `border_mode` - Pixel extrapolation method (e.g., `BorderTypes::Constant`). +/// * `border_value` - Value used in case of a constant border. +/// +/// # Returns +/// +/// Returns a `Result>` containing the transformed image. +/// +/// # Errors +/// +/// Returns an error if: +/// * `map1` and `map2` do not have the same dimensions. +/// * `map1` or `map2` are not single-channel matrices. +/// +/// # Examples +/// +/// ``` +/// use purecv::core::{Matrix, types::{BorderTypes, Scalar}}; +/// use purecv::imgproc::geometric::{remap, InterpolationFlags}; +/// +/// let src = Matrix::::new(10, 10, 1); +/// let map1 = Matrix::::new(10, 10, 1); +/// let map2 = Matrix::::new(10, 10, 1); +/// +/// let result = remap( +/// &src, +/// &map1, +/// &map2, +/// InterpolationFlags::Linear, +/// BorderTypes::Constant, +/// Scalar::all(0) +/// ).unwrap(); +/// ``` +pub fn remap( + src: &Matrix, + map1: &Matrix, + map2: &Matrix, + interpolation: InterpolationFlags, + border_mode: BorderTypes, + border_value: Scalar, +) -> Result> +where + T: Default + Clone + Copy + ToPrimitive + FromPrimitive + Send + Sync + SimdElement, +{ + if map1.rows != map2.rows || map1.cols != map2.cols { + return Err(PureCvError::InvalidInput( + "map1 and map2 must have the same dimensions".to_string(), + )); + } + if map1.channels != 1 || map2.channels != 1 { + return Err(PureCvError::InvalidInput( + "map1 and map2 must be single-channel matrices".to_string(), + )); + } + + let rows = map1.rows; + let cols = map1.cols; + let channels = src.channels; + + let mut dst = Matrix::::new(rows, cols, channels); + + #[cfg(feature = "parallel")] + { + dst.data + .par_chunks_exact_mut(cols * channels) + .enumerate() + .for_each(|(r, row_data)| { + let map1_row = &map1.data[r * cols..(r + 1) * cols]; + let map2_row = &map2.data[r * cols..(r + 1) * cols]; + + let handled = match interpolation { + InterpolationFlags::Nearest => T::simd_remap_nearest_row( + row_data, + &src.data, + src.cols, + src.rows, + channels, + map1_row, + map2_row, + border_mode, + border_value, + ), + InterpolationFlags::Linear => T::simd_remap_bilinear_row( + row_data, + &src.data, + src.cols, + src.rows, + channels, + map1_row, + map2_row, + border_mode, + border_value, + ), + }; + + if !handled { + remap_row_scalar( + row_data, + &src.data, + src.cols, + src.rows, + channels, + map1_row, + map2_row, + interpolation, + border_mode, + border_value, + ); + } + }); + } + + #[cfg(not(feature = "parallel"))] + { + for r in 0..rows { + let row_data = &mut dst.data[r * cols * channels..(r + 1) * cols * channels]; + let map1_row = &map1.data[r * cols..(r + 1) * cols]; + let map2_row = &map2.data[r * cols..(r + 1) * cols]; + + let handled = match interpolation { + InterpolationFlags::Nearest => T::simd_remap_nearest_row( + row_data, + &src.data, + src.cols, + src.rows, + channels, + map1_row, + map2_row, + border_mode, + border_value, + ), + InterpolationFlags::Linear => T::simd_remap_bilinear_row( + row_data, + &src.data, + src.cols, + src.rows, + channels, + map1_row, + map2_row, + border_mode, + border_value, + ), + }; + + if !handled { + remap_row_scalar( + row_data, + &src.data, + src.cols, + src.rows, + channels, + map1_row, + map2_row, + interpolation, + border_mode, + border_value, + ); + } + } + } + + Ok(dst) +} + +/// Applies a perspective transformation to an image. +/// +/// The function `warp_perspective` transforms the source image using the specified matrix: +/// `dst(x,y) = src((M_00*x + M_01*y + M_02)/(M_20*x + M_21*y + M_22), (M_10*x + M_11*y + M_12)/(M_20*x + M_21*y + M_22))` +/// +/// # Arguments +/// +/// * `src` - Source image. +/// * `m` - 3x3 perspective transformation matrix (`f64`). +/// * `dsize` - Size of the destination image. +/// * `flags` - Interpolation method to use (e.g., `InterpolationFlags::Linear`). +/// * `border_mode` - Pixel extrapolation method (e.g., `BorderTypes::Constant`). +/// * `border_value` - Value used in case of a constant border. +/// +/// # Returns +/// +/// Returns a `Result>` containing the perspective-transformed image. +/// +/// # Errors +/// +/// Returns an error if the homography matrix `m` is not a 3x3 single-channel matrix. +/// +/// # Examples +/// +/// ``` +/// use purecv::core::{Matrix, types::{BorderTypes, Scalar, Size2i}}; +/// use purecv::imgproc::geometric::{warp_perspective, InterpolationFlags}; +/// +/// let src = Matrix::::new(10, 10, 1); +/// let mut m = Matrix::::new(3, 3, 1); +/// // Set identity matrix for example +/// m.data[0] = 1.0; m.data[4] = 1.0; m.data[8] = 1.0; +/// +/// let result = warp_perspective( +/// &src, +/// &m, +/// Size2i::new(10, 10), +/// InterpolationFlags::Linear, +/// BorderTypes::Constant, +/// Scalar::all(0) +/// ).unwrap(); +/// ``` +pub fn warp_perspective( + src: &Matrix, + m: &Matrix, + dsize: Size2i, + flags: InterpolationFlags, + border_mode: BorderTypes, + border_value: Scalar, +) -> Result> +where + T: Default + Clone + Copy + ToPrimitive + FromPrimitive + Send + Sync + SimdElement, +{ + if m.rows != 3 || m.cols != 3 || m.channels != 1 { + return Err(PureCvError::InvalidDimensions( + "Homography matrix must be 3x3 single-channel".to_string(), + )); + } + + let mut m_inv = Matrix::::new(3, 3, 1); + invert(m, &mut m_inv, DecompTypes::DECOMP_LU)?; + + let m00 = m_inv.data[0]; + let m01 = m_inv.data[1]; + let m02 = m_inv.data[2]; + let m10 = m_inv.data[3]; + let m11 = m_inv.data[4]; + let m12 = m_inv.data[5]; + let m20 = m_inv.data[6]; + let m21 = m_inv.data[7]; + let m22 = m_inv.data[8]; + + let channels = src.channels; + let mut dst = Matrix::::new(dsize.height as usize, dsize.width as usize, channels); + + #[cfg(feature = "parallel")] + { + dst.data + .par_chunks_exact_mut(dsize.width as usize * channels) + .enumerate() + .for_each(|(r, row_data)| { + let dst_y = r as f64; + let mut map_x = vec![0.0f32; dsize.width as usize]; + let mut map_y = vec![0.0f32; dsize.width as usize]; + + for col in 0..dsize.width as usize { + let dst_x = col as f64; + let w = m20 * dst_x + m21 * dst_y + m22; + if w.abs() > 1e-10 { + map_x[col] = ((m00 * dst_x + m01 * dst_y + m02) / w) as f32; + map_y[col] = ((m10 * dst_x + m11 * dst_y + m12) / w) as f32; + } else { + map_x[col] = -1.0; + map_y[col] = -1.0; + } + } + + let handled = match flags { + InterpolationFlags::Nearest => T::simd_remap_nearest_row( + row_data, + &src.data, + src.cols, + src.rows, + channels, + &map_x, + &map_y, + border_mode, + border_value, + ), + InterpolationFlags::Linear => T::simd_remap_bilinear_row( + row_data, + &src.data, + src.cols, + src.rows, + channels, + &map_x, + &map_y, + border_mode, + border_value, + ), + }; + + if !handled { + remap_row_scalar( + row_data, + &src.data, + src.cols, + src.rows, + channels, + &map_x, + &map_y, + flags, + border_mode, + border_value, + ); + } + }); + } + + #[cfg(not(feature = "parallel"))] + { + let width = dsize.width as usize; + let height = dsize.height as usize; + for r in 0..height { + let row_data = &mut dst.data[r * width * channels..(r + 1) * width * channels]; + let dst_y = r as f64; + let mut map_x = vec![0.0f32; width]; + let mut map_y = vec![0.0f32; width]; + + for col in 0..width { + let dst_x = col as f64; + let w = m20 * dst_x + m21 * dst_y + m22; + if w.abs() > 1e-10 { + map_x[col] = ((m00 * dst_x + m01 * dst_y + m02) / w) as f32; + map_y[col] = ((m10 * dst_x + m11 * dst_y + m12) / w) as f32; + } else { + map_x[col] = -1.0; + map_y[col] = -1.0; + } + } + + let handled = match flags { + InterpolationFlags::Nearest => T::simd_remap_nearest_row( + row_data, + &src.data, + src.cols, + src.rows, + channels, + &map_x, + &map_y, + border_mode, + border_value, + ), + InterpolationFlags::Linear => T::simd_remap_bilinear_row( + row_data, + &src.data, + src.cols, + src.rows, + channels, + &map_x, + &map_y, + border_mode, + border_value, + ), + }; + + if !handled { + remap_row_scalar( + row_data, + &src.data, + src.cols, + src.rows, + channels, + &map_x, + &map_y, + flags, + border_mode, + border_value, + ); + } + } + } + + Ok(dst) +} + +/// Fallback scalar implementation for remapping a row. +#[allow(clippy::too_many_arguments)] +fn remap_row_scalar( + dst_row: &mut [T], + src: &[T], + src_cols: usize, + src_rows: usize, + channels: usize, + map1_row: &[f32], + map2_row: &[f32], + interpolation: InterpolationFlags, + border_mode: BorderTypes, + border_value: Scalar, +) where + T: Default + Clone + Copy + ToPrimitive + FromPrimitive, +{ + let dst_cols = dst_row.len() / channels; + + match interpolation { + InterpolationFlags::Nearest => { + for col in 0..dst_cols { + let x = map1_row[col]; + let y = map2_row[col]; + let ix = x.round() as i32; + let iy = y.round() as i32; + let dst_idx = col * channels; + + if ix >= 0 && ix < src_cols as i32 && iy >= 0 && iy < src_rows as i32 { + let src_idx = (iy as usize * src_cols + ix as usize) * channels; + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&src[src_idx..src_idx + channels]); + } else { + let ix_interp = border_interpolate(ix, src_cols as i32, border_mode); + let iy_interp = border_interpolate(iy, src_rows as i32, border_mode); + if ix_interp >= 0 && iy_interp >= 0 { + let src_idx = + (iy_interp as usize * src_cols + ix_interp as usize) * channels; + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&src[src_idx..src_idx + channels]); + } else { + dst_row[dst_idx..dst_idx + channels] + .copy_from_slice(&border_value.v[..channels]); + } + } + } + } + InterpolationFlags::Linear => { + for col in 0..dst_cols { + let x = map1_row[col]; + let y = map2_row[col]; + + let x1 = x.floor() as i32; + let y1 = y.floor() as i32; + let x2 = x1 + 1; + let y2 = y1 + 1; + + let wx = x - x.floor(); + let wy = y - y.floor(); + + let w00 = (1.0 - wx) * (1.0 - wy); + let w10 = wx * (1.0 - wy); + let w01 = (1.0 - wx) * wy; + let w11 = wx * wy; + + let dst_idx = col * channels; + + if x1 >= 0 && x2 < src_cols as i32 && y1 >= 0 && y2 < src_rows as i32 { + let idx00 = (y1 as usize * src_cols + x1 as usize) * channels; + let idx10 = (y1 as usize * src_cols + x2 as usize) * channels; + let idx01 = (y2 as usize * src_cols + x1 as usize) * channels; + let idx11 = (y2 as usize * src_cols + x2 as usize) * channels; + + for c in 0..channels { + let v00 = src[idx00 + c].to_f32().unwrap_or(0.0); + let v10 = src[idx10 + c].to_f32().unwrap_or(0.0); + let v01 = src[idx01 + c].to_f32().unwrap_or(0.0); + let v11 = src[idx11 + c].to_f32().unwrap_or(0.0); + + let val = w00 * v00 + w10 * v10 + w01 * v01 + w11 * v11; + dst_row[dst_idx + c] = T::from_f32(val).unwrap_or_default(); + } + } else { + for c in 0..channels { + let get_pixel = |ix: i32, iy: i32| -> f32 { + let ix_interp = border_interpolate(ix, src_cols as i32, border_mode); + let iy_interp = border_interpolate(iy, src_rows as i32, border_mode); + if ix_interp >= 0 && iy_interp >= 0 { + src[(iy_interp as usize * src_cols + ix_interp as usize) * channels + + c] + .to_f32() + .unwrap_or(0.0) + } else { + border_value.v[c].to_f32().unwrap_or(0.0) + } + }; + + let v00 = get_pixel(x1, y1); + let v10 = get_pixel(x2, y1); + let v01 = get_pixel(x1, y2); + let v11 = get_pixel(x2, y2); + + let val = w00 * v00 + w10 * v10 + w01 * v01 + w11 * v11; + dst_row[dst_idx + c] = T::from_f32(val).unwrap_or_default(); + } + } + } + } + } +} diff --git a/src/imgproc/tests.rs b/src/imgproc/tests.rs index 8e92ed3..b758584 100644 --- a/src/imgproc/tests.rs +++ b/src/imgproc/tests.rs @@ -1229,4 +1229,91 @@ mod imgproc_tests { assert!(resize(&src, Size::new(0, 4)).is_err()); assert!(resize(&src, Size::new(4, 0)).is_err()); } + + #[test] + fn test_remap_nearest() { + let src = Matrix::from_vec(3, 3, 1, vec![10u8, 20, 30, 40, 50, 60, 70, 80, 90]); + let mut map1 = Matrix::::new(3, 3, 1); + let mut map2 = Matrix::::new(3, 3, 1); + for y in 0..3 { + for x in 0..3 { + *map1.at_mut(y, x, 0).unwrap() = (x as f32) + 1.0; + *map2.at_mut(y, x, 0).unwrap() = y as f32; + } + } + + let res = remap( + &src, + &map1, + &map2, + InterpolationFlags::Nearest, + BorderTypes::Constant, + Scalar::all(0u8), + ) + .unwrap(); + + assert_eq!(res.data, vec![20, 30, 0, 50, 60, 0, 80, 90, 0,]); + } + + #[test] + fn test_remap_bilinear() { + let src = Matrix::from_vec(2, 2, 1, vec![10.0f32, 20.0, 30.0, 40.0]); + let mut map1 = Matrix::::new(2, 2, 1); + let mut map2 = Matrix::::new(2, 2, 1); + for y in 0..2 { + for x in 0..2 { + *map1.at_mut(y, x, 0).unwrap() = 0.5; + *map2.at_mut(y, x, 0).unwrap() = 0.5; + } + } + + let res = remap( + &src, + &map1, + &map2, + InterpolationFlags::Linear, + BorderTypes::Constant, + Scalar::all(0.0f32), + ) + .unwrap(); + + for val in res.data { + assert!((val - 25.0).abs() < 1e-4); + } + } + + #[test] + fn test_warp_perspective_identity() { + let src = Matrix::from_vec(2, 2, 1, vec![10u8, 20, 30, 40]); + let m = Matrix::from_vec(3, 3, 1, vec![1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]); + let res = warp_perspective( + &src, + &m, + Size2i::new(2, 2), + InterpolationFlags::Nearest, + BorderTypes::Constant, + Scalar::all(0u8), + ) + .unwrap(); + + assert_eq!(res.data, vec![10, 20, 30, 40]); + } + + #[test] + fn test_warp_perspective_translation() { + let src = Matrix::from_vec(2, 2, 1, vec![10f32, 20.0, 30.0, 40.0]); + // Translate by dx=1.0, dy=0.0 + let m = Matrix::from_vec(3, 3, 1, vec![1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]); + let res = warp_perspective( + &src, + &m, + Size2i::new(2, 2), + InterpolationFlags::Nearest, + BorderTypes::Constant, + Scalar::all(0.0f32), + ) + .unwrap(); + + assert_eq!(res.data, vec![0.0, 10.0, 0.0, 30.0]); + } } diff --git a/src/lib.rs b/src/lib.rs index 1e901e7..065302c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -46,7 +46,8 @@ pub mod video; /// Prelude to easily import common structures pub mod prelude { pub use crate::calib3d::{ - find_homography, rodrigues, solve_pnp, solve_pnp_ransac, HomographyMethod, SolvePnPMethod, + find_fundamental_mat, find_homography, init_undistort_rectify_map, rodrigues, solve_pnp, + solve_pnp_ransac, FundamentalMatMethod, HomographyMethod, SolvePnPMethod, }; pub use crate::core::types::{ BorderTypes, Point2f, Point2i, Point3f, Rect2f, Rect2i, Scalar, Size2f, Size2i, @@ -66,7 +67,9 @@ pub mod prelude { }; pub use crate::imgproc::filter::{bilateral_filter, box_filter, gaussian_blur}; pub use crate::imgproc::threshold::{threshold, ThresholdTypes}; - pub use crate::imgproc::{cvt_color, ColorConversionCode}; + pub use crate::imgproc::{ + cvt_color, remap, warp_perspective, ColorConversionCode, InterpolationFlags, + }; pub use crate::video::optical_flow::{ build_optical_flow_pyramid, calc_optical_flow_pyramid_lk, OpticalFlowPyramid, OPTFLOW_LK_GET_MIN_EIGENVALS, OPTFLOW_USE_INITIAL_FLOW,