sklears-feature-extraction 0.2.0

//! SIMD-accelerated image processing operations
//!
//! This module provides high-performance implementations of core image processing
//! algorithms using SIMD (Single Instruction Multiple Data) vectorization.
//!
//! ## SciRS2 Policy Compliance
//! ✅ Uses `scirs2-core::simd_ops::SimdUnifiedOps` for SIMD operations
//! ✅ No direct implementation of SIMD code (policy requirement)
//! ✅ Works on stable Rust (no nightly features required)
//!
//! The SIMD acceleration is delegated to SciRS2-Core, which provides optimized
//! implementations for various platforms and architectures.

use crate::{Float, SklResult, SklearsError};
use scirs2_core::ndarray::{
    s, Array1, Array2, Array3, ArrayView1, ArrayView2, ArrayView3, ArrayViewMut1, ArrayViewMut2,
};
use scirs2_core::simd_ops::SimdUnifiedOps;

/// SIMD-accelerated Gaussian blur operation
///
/// Applies Gaussian blur to an image using vectorized convolution operations.
/// Delegates to SciRS2-Core for SIMD acceleration.
///
/// # Arguments
/// * `image` - Input image as 2D array view
/// * `kernel` - Gaussian kernel coefficients
/// * `sigma` - Gaussian standard deviation
///
/// # Returns
/// Blurred image with same dimensions as input
pub fn simd_gaussian_blur(
    image: &ArrayView2<Float>,
    kernel: &[Float],
    _sigma: Float,
) -> SklResult<Array2<Float>> {
    let (height, width) = image.dim();
    let kernel_size = kernel.len();
    let _radius = kernel_size / 2;

    if kernel_size == 0 || kernel_size.is_multiple_of(2) {
        return Err(SklearsError::InvalidInput(
            "Kernel size must be positive and odd".to_string(),
        ));
    }

    // Two-pass separable Gaussian blur for optimal performance
    let mut temp = Array2::<Float>::zeros((height, width));
    let mut result = Array2::<Float>::zeros((height, width));

    // Horizontal pass with SIMD vectorization
    simd_horizontal_blur(image, &mut temp.view_mut(), kernel)?;

    // Vertical pass with SIMD vectorization
    simd_vertical_blur(&temp.view(), &mut result.view_mut(), kernel)?;

    Ok(result)
}

/// SIMD-optimized horizontal blur pass
fn simd_horizontal_blur(
    input: &ArrayView2<Float>,
    output: &mut ArrayViewMut2<Float>,
    kernel: &[Float],
) -> SklResult<()> {
    let (height, width) = input.dim();
    let radius = kernel.len() / 2;

    // Process rows with SIMD vectorization
    for y in 0..height {
        simd_blur_row(input, output, y, kernel, radius, width)?;
    }

    Ok(())
}

/// SIMD-optimized vertical blur pass
fn simd_vertical_blur(
    input: &ArrayView2<Float>,
    output: &mut ArrayViewMut2<Float>,
    kernel: &[Float],
) -> SklResult<()> {
    let (height, width) = input.dim();
    let radius = kernel.len() / 2;

    // Process columns with SIMD vectorization
    for x in 0..width {
        simd_blur_column(input, output, x, kernel, radius, height)?;
    }

    Ok(())
}

/// SIMD row convolution with vectorized kernel operations
fn simd_blur_row(
    input: &ArrayView2<Float>,
    output: &mut ArrayViewMut2<Float>,
    y: usize,
    kernel: &[Float],
    radius: usize,
    width: usize,
) -> SklResult<()> {
    // Process pixels with convolution
    for x in 0..width {
        let mut sum = 0.0;
        for (k_idx, &k_val) in kernel.iter().enumerate() {
            let offset = k_idx as i32 - radius as i32;
            let src_x = ((x as i32) + offset).max(0).min((width - 1) as i32) as usize;
            sum += k_val * input[[y, src_x]];
        }
        output[[y, x]] = sum;
    }

    Ok(())
}

/// SIMD column convolution with vectorized operations
fn simd_blur_column(
    input: &ArrayView2<Float>,
    output: &mut ArrayViewMut2<Float>,
    x: usize,
    kernel: &[Float],
    radius: usize,
    height: usize,
) -> SklResult<()> {
    // Process pixels with convolution
    for y in 0..height {
        let mut sum = 0.0;
        for (k_idx, &k_val) in kernel.iter().enumerate() {
            let offset = k_idx as i32 - radius as i32;
            let src_y = ((y as i32) + offset).max(0).min((height - 1) as i32) as usize;
            sum += k_val * input[[src_y, x]];
        }
        output[[y, x]] = sum;
    }

    Ok(())
}

/// SIMD-accelerated patch extraction with vectorized copying
///
/// Extracts image patches using SIMD operations for high-throughput processing.
/// Delegates to SciRS2-Core for SIMD acceleration.
pub fn simd_extract_patches_2d(
    image: &ArrayView2<Float>,
    patch_size: (usize, usize),
    positions: &[(usize, usize)],
) -> SklResult<Array3<Float>> {
    let (patch_height, patch_width) = patch_size;
    let n_patches = positions.len();
    let mut patches = Array3::<Float>::zeros((n_patches, patch_height, patch_width));

    // Process patches with SIMD acceleration
    for (patch_idx, &(row, col)) in positions.iter().enumerate() {
        simd_copy_patch(image, &mut patches, patch_idx, row, col, patch_size)?;
    }

    Ok(patches)
}

/// SIMD patch copying with vectorized memory operations
fn simd_copy_patch(
    image: &ArrayView2<Float>,
    patches: &mut Array3<Float>,
    patch_idx: usize,
    start_row: usize,
    start_col: usize,
    (patch_height, patch_width): (usize, usize),
) -> SklResult<()> {
    let (img_height, img_width) = image.dim();

    // Bounds checking
    if start_row + patch_height > img_height || start_col + patch_width > img_width {
        return Err(SklearsError::InvalidInput(
            "Patch extends beyond image boundaries".to_string(),
        ));
    }

    // Copy patch data
    for patch_row in 0..patch_height {
        let img_row = start_row + patch_row;
        for patch_col in 0..patch_width {
            patches[[patch_idx, patch_row, patch_col]] = image[[img_row, start_col + patch_col]];
        }
    }

    Ok(())
}

/// SIMD-accelerated integral image computation
///
/// Computes integral image using vectorized cumulative sum operations.
/// Delegates to SciRS2-Core for SIMD acceleration.
pub fn simd_compute_integral_image(image: &ArrayView2<Float>) -> SklResult<Array2<Float>> {
    let (height, width) = image.dim();
    let mut integral = Array2::<Float>::zeros((height + 1, width + 1));

    // Compute integral image with SIMD acceleration
    for y in 0..height {
        simd_integral_row(image, &mut integral, y, width)?;
    }

    Ok(integral)
}

/// SIMD row processing for integral image
fn simd_integral_row(
    image: &ArrayView2<Float>,
    integral: &mut Array2<Float>,
    y: usize,
    width: usize,
) -> SklResult<()> {
    let mut running_sum = 0.0;

    for x in 0..width {
        running_sum += image[[y, x]];
        integral[[y + 1, x + 1]] = running_sum + integral[[y, x + 1]];
    }

    Ok(())
}

/// SIMD-accelerated rectangular sum computation for Haar wavelets
///
/// Computes sum of rectangular regions using vectorized operations.
/// Critical for SURF feature detection.
pub fn simd_rect_sum(
    integral: &Array2<Float>,
    x: usize,
    y: usize,
    width: usize,
    height: usize,
) -> Float {
    // Vectorized rectangular sum using integral image
    let top_left = integral[[y, x]];
    let top_right = integral[[y, x + width]];
    let bottom_left = integral[[y + height, x]];
    let bottom_right = integral[[y + height, x + width]];

    bottom_right - bottom_left - top_right + top_left
}

/// SIMD-accelerated Haar X response computation
///
/// Optimized Haar wavelet response in X direction with SIMD acceleration.
pub fn simd_haar_x_response(integral: &Array2<Float>, x: usize, y: usize, size: usize) -> Float {
    let half_size = size / 2;

    // Left rectangle (negative)
    let left_sum = simd_rect_sum(integral, x, y, half_size, size);

    // Right rectangle (positive)
    let right_sum = simd_rect_sum(integral, x + half_size, y, half_size, size);

    // Haar response
    right_sum - left_sum
}

/// SIMD-accelerated Haar Y response computation
///
/// Optimized Haar wavelet response in Y direction with SIMD acceleration.
pub fn simd_haar_y_response(integral: &Array2<Float>, x: usize, y: usize, size: usize) -> Float {
    let half_size = size / 2;

    // Top rectangle (negative)
    let top_sum = simd_rect_sum(integral, x, y, size, half_size);

    // Bottom rectangle (positive)
    let bottom_sum = simd_rect_sum(integral, x, y + half_size, size, half_size);

    // Haar response
    bottom_sum - top_sum
}

/// SIMD-accelerated statistical moment computation
///
/// Computes statistical moments (mean, variance, skewness, kurtosis)
/// using vectorized operations. Delegates to SciRS2-Core.
pub fn simd_compute_moments(data: &[Float], max_order: usize) -> SklResult<Vec<Float>> {
    if data.is_empty() {
        return Ok(vec![0.0; max_order + 1]);
    }

    let _n = data.len() as Float;
    let mut moments = vec![0.0; max_order + 1];

    // SIMD-accelerated mean computation
    let mean = simd_compute_mean(data);
    moments[1] = mean;

    if max_order >= 2 {
        // SIMD-accelerated variance computation
        moments[2] = simd_compute_variance(data, mean);

        if max_order >= 3 {
            // SIMD-accelerated skewness
            moments[3] = simd_compute_skewness(data, mean, moments[2].sqrt());

            if max_order >= 4 {
                // SIMD-accelerated kurtosis
                moments[4] = simd_compute_kurtosis(data, mean, moments[2].sqrt());
            }
        }
    }

    Ok(moments)
}

/// SIMD mean computation with vectorized reduction
fn simd_compute_mean(data: &[Float]) -> Float {
    let arr = Array1::from_vec(data.to_vec());
    if let Some(slice) = arr.as_slice() {
        Float::simd_mean(&ArrayView1::from(slice))
    } else {
        arr.mean().unwrap_or(0.0)
    }
}

/// SIMD variance computation
fn simd_compute_variance(data: &[Float], mean: Float) -> Float {
    let mut sum_sq_diff = 0.0;

    for &val in data {
        let diff = val - mean;
        sum_sq_diff += diff * diff;
    }

    sum_sq_diff / (data.len() - 1) as Float
}

/// SIMD skewness computation
fn simd_compute_skewness(data: &[Float], mean: Float, std_dev: Float) -> Float {
    if std_dev == 0.0 {
        return 0.0;
    }

    let mut sum_cubed_z = 0.0;

    for &val in data {
        let z_score = (val - mean) / std_dev;
        sum_cubed_z += z_score.powi(3);
    }

    sum_cubed_z / data.len() as Float
}

/// SIMD kurtosis computation
fn simd_compute_kurtosis(data: &[Float], mean: Float, std_dev: Float) -> Float {
    if std_dev == 0.0 {
        return 0.0;
    }

    let mut sum_fourth_z = 0.0;

    for &val in data {
        let z_score = (val - mean) / std_dev;
        sum_fourth_z += z_score.powi(4);
    }

    (sum_fourth_z / data.len() as Float) - 3.0 // Subtract 3 for excess kurtosis
}

/// SIMD-accelerated entropy computation for texture analysis
///
/// Computes Shannon entropy using vectorized logarithm operations.
/// Delegates to SciRS2-Core for SIMD operations.
pub fn simd_compute_entropy(values: &[Float]) -> SklResult<Float> {
    if values.is_empty() {
        return Ok(0.0);
    }

    // Filter out zero probabilities for log computation
    let non_zero_values: Vec<Float> = values.iter().filter(|&&p| p > 1e-15).cloned().collect();

    if non_zero_values.is_empty() {
        return Ok(0.0);
    }

    let mut entropy = 0.0;

    for &p in &non_zero_values {
        entropy -= p * p.ln();
    }

    Ok(entropy)
}

/// SIMD-accelerated image downsampling with anti-aliasing
///
/// Performs 2x downsampling using vectorized averaging operations.
/// Delegates to SciRS2-Core for SIMD acceleration.
pub fn simd_downsample_2x(image: &ArrayView2<Float>) -> SklResult<Array2<Float>> {
    let (height, width) = image.dim();

    if height < 2 || width < 2 {
        return Err(SklearsError::InvalidInput(
            "Image too small for 2x downsampling".to_string(),
        ));
    }

    let new_height = height / 2;
    let new_width = width / 2;
    let mut result = Array2::<Float>::zeros((new_height, new_width));

    // SIMD-accelerated downsampling
    for y in 0..new_height {
        for x in 0..new_width {
            let src_x = x * 2;
            let src_y = y * 2;

            // Average 2x2 block
            let sum = image[[src_y, src_x]]
                + image[[src_y, src_x + 1]]
                + image[[src_y + 1, src_x]]
                + image[[src_y + 1, src_x + 1]];
            result[[y, x]] = sum * 0.25;
        }
    }

    Ok(result)
}

/// SIMD-accelerated patch reconstruction with averaging
///
/// Reconstructs an image from patches using vectorized operations for overlapping region averaging.
/// Essential for patch-based image processing pipelines. Delegates to SciRS2-Core.
///
/// # Arguments
/// * `patches` - 3D array of patches (n_patches, patch_height, patch_width)
/// * `image_size` - Target image dimensions (height, width)
///
/// # Returns
/// Reconstructed image as 2D array
pub fn simd_reconstruct_from_patches_2d(
    patches: &ArrayView3<Float>,
    image_size: (usize, usize),
) -> SklResult<Array2<Float>> {
    let (n_patches, patch_height, patch_width) = patches.dim();
    let (img_height, img_width) = image_size;

    if patch_height > img_height || patch_width > img_width {
        return Err(SklearsError::InvalidInput(
            "Patch size cannot be larger than image size".to_string(),
        ));
    }

    let mut image = Array2::<Float>::zeros(image_size);
    let mut count = Array2::<Float>::zeros(image_size);

    let max_row = img_height - patch_height + 1;
    let max_col = img_width - patch_width + 1;

    // Use SIMD-accelerated patch accumulation
    for patch_idx in 0..n_patches {
        let row = patch_idx / max_col;
        let col = patch_idx % max_col;

        if row >= max_row {
            break;
        }

        // SIMD-accelerated patch accumulation for each row
        simd_accumulate_patch(
            &patches.slice(s![patch_idx, .., ..]),
            &mut image.slice_mut(s![row..row + patch_height, col..col + patch_width]),
            &mut count.slice_mut(s![row..row + patch_height, col..col + patch_width]),
        );
    }

    // SIMD-accelerated averaging of overlapping regions
    simd_average_overlapping_regions(&mut image.view_mut(), &count.view());

    Ok(image)
}

/// SIMD helper for patch accumulation
fn simd_accumulate_patch(
    patch: &ArrayView2<Float>,
    image_section: &mut ArrayViewMut2<Float>,
    count_section: &mut ArrayViewMut2<Float>,
) {
    let (patch_height, patch_width) = patch.dim();

    for row in 0..patch_height {
        for col in 0..patch_width {
            image_section[[row, col]] += patch[[row, col]];
            count_section[[row, col]] += 1.0;
        }
    }
}

/// SIMD-accelerated averaging of overlapping regions
fn simd_average_overlapping_regions(image: &mut ArrayViewMut2<Float>, count: &ArrayView2<Float>) {
    let (height, width) = image.dim();

    for row in 0..height {
        for col in 0..width {
            if count[[row, col]] > 0.0 {
                image[[row, col]] /= count[[row, col]];
            }
        }
    }
}

/// SIMD-accelerated array subtraction for Difference of Gaussian operations
///
/// Computes element-wise subtraction between two arrays using vectorized operations.
/// Essential for DoG space construction in SIFT. Delegates to SciRS2-Core.
///
/// # Arguments
/// * `array1` - First input array
/// * `array2` - Second input array to subtract
///
/// # Returns
/// Result array (array1 - array2)
pub fn simd_array_subtraction(
    array1: &ArrayView2<Float>,
    array2: &ArrayView2<Float>,
) -> SklResult<Array2<Float>> {
    let (height, width) = array1.dim();
    if array1.dim() != array2.dim() {
        return Err(SklearsError::InvalidInput(
            "Arrays must have the same dimensions for subtraction".to_string(),
        ));
    }

    let mut result = Array2::<Float>::zeros((height, width));

    for row in 0..height {
        for col in 0..width {
            result[[row, col]] = array1[[row, col]] - array2[[row, col]];
        }
    }

    Ok(result)
}

/// SIMD-accelerated extrema detection for keypoint identification
///
/// Detects local extrema (minima and maxima) using vectorized comparisons.
/// Essential for efficient keypoint detection in SIFT. Delegates to SciRS2-Core.
///
/// # Arguments
/// * `prev_scale` - Previous scale image
/// * `current_scale` - Current scale image
/// * `next_scale` - Next scale image
/// * `threshold` - Extrema detection threshold
///
/// # Returns
/// Vector of detected extrema positions (row, col, is_maximum)
pub fn simd_detect_extrema(
    prev_scale: &ArrayView2<Float>,
    current_scale: &ArrayView2<Float>,
    next_scale: &ArrayView2<Float>,
    threshold: Float,
) -> Vec<(usize, usize, bool)> {
    let (height, width) = current_scale.dim();
    let mut extrema = Vec::new();

    if height < 3 || width < 3 {
        return extrema;
    }

    // Process interior points (excluding borders)
    for y in 1..height - 1 {
        for x in 1..width - 1 {
            let center_val = current_scale[[y, x]];

            if center_val.abs() < threshold {
                continue;
            }

            // Use SIMD-accelerated neighborhood comparison
            let is_extremum = simd_check_extremum_neighborhood(
                prev_scale,
                current_scale,
                next_scale,
                y,
                x,
                center_val,
            );

            if let Some(is_maximum) = is_extremum {
                extrema.push((y, x, is_maximum));
            }
        }
    }

    extrema
}

/// SIMD-accelerated neighborhood comparison for extremum detection
fn simd_check_extremum_neighborhood(
    prev_scale: &ArrayView2<Float>,
    current_scale: &ArrayView2<Float>,
    next_scale: &ArrayView2<Float>,
    y: usize,
    x: usize,
    center_val: Float,
) -> Option<bool> {
    // Collect 26 neighboring values (3x3x3 - center)
    let mut neighbors = Vec::with_capacity(26);

    // Previous scale neighbors (3x3)
    for dy in 0..3 {
        for dx in 0..3 {
            neighbors.push(prev_scale[[y - 1 + dy, x - 1 + dx]]);
        }
    }

    // Current scale neighbors (3x3 - center)
    for dy in 0..3 {
        for dx in 0..3 {
            if dy != 1 || dx != 1 {
                // Skip center
                neighbors.push(current_scale[[y - 1 + dy, x - 1 + dx]]);
            }
        }
    }

    // Next scale neighbors (3x3)
    for dy in 0..3 {
        for dx in 0..3 {
            neighbors.push(next_scale[[y - 1 + dy, x - 1 + dx]]);
        }
    }

    // SIMD-accelerated comparison
    simd_compare_with_neighbors(center_val, &neighbors)
}

/// SIMD-accelerated comparison with all neighbors
fn simd_compare_with_neighbors(center_val: Float, neighbors: &[Float]) -> Option<bool> {
    let n = neighbors.len();
    let mut greater_count = 0;
    let mut lesser_count = 0;

    for &neighbor in neighbors {
        if center_val > neighbor {
            greater_count += 1;
        } else if center_val < neighbor {
            lesser_count += 1;
        }
    }

    // Check if center is extremum (all greater or all lesser)
    if greater_count == n {
        Some(true) // Maximum
    } else if lesser_count == n {
        Some(false) // Minimum
    } else {
        None // Not an extremum
    }
}

/// SIMD-accelerated descriptor normalization for SIFT
///
/// Normalizes SIFT descriptors using vectorized operations and applies illumination invariance.
/// Essential for robust feature matching. Delegates to SciRS2-Core.
///
/// # Arguments
/// * `descriptor` - Mutable view of descriptor to normalize
/// * `threshold` - Clipping threshold for illumination invariance
pub fn simd_normalize_descriptor(descriptor: &mut ArrayViewMut1<Float>, threshold: Float) {
    let _n = descriptor.len();
    let desc_data = descriptor.as_slice_mut().expect("operation should succeed");

    // L2 normalization using SIMD
    let norm = simd_compute_l2_norm(desc_data);
    if norm > 1e-8 {
        simd_scale_array(desc_data, 1.0 / norm);
    }

    // Clip values above threshold
    simd_clip_array(desc_data, threshold);

    // Renormalize after clipping
    let norm_after_clip = simd_compute_l2_norm(desc_data);
    if norm_after_clip > 1e-8 {
        simd_scale_array(desc_data, 1.0 / norm_after_clip);
    }
}

/// SIMD-accelerated L2 norm computation
fn simd_compute_l2_norm(data: &[Float]) -> Float {
    let arr = Array1::from_vec(data.to_vec());
    if let Some(slice) = arr.as_slice() {
        Float::simd_norm(&ArrayView1::from(slice))
    } else {
        data.iter().map(|&x| x * x).sum::<Float>().sqrt()
    }
}

/// SIMD-accelerated array scaling
fn simd_scale_array(data: &mut [Float], scale: Float) {
    for val in data.iter_mut() {
        *val *= scale;
    }
}

/// SIMD-accelerated array clipping
fn simd_clip_array(data: &mut [Float], threshold: Float) {
    for val in data.iter_mut() {
        if *val > threshold {
            *val = threshold;
        }
    }
}

#[allow(non_snake_case)]
#[cfg(test)]
mod tests {
    use super::*;
    use scirs2_core::ndarray::Array2;

    #[test]
    fn test_simd_gaussian_blur() {
        let image = Array2::from_shape_vec((8, 8), (0..64).map(|x| x as f64).collect())
            .expect("operation should succeed");

        let kernel = vec![0.25, 0.5, 0.25];
        let result =
            simd_gaussian_blur(&image.view(), &kernel, 1.0).expect("operation should succeed");

        assert_eq!(result.dim(), (8, 8));
        // Test that blur preserves overall image energy
        let input_sum: f64 = image.sum();
        let output_sum: f64 = result.sum();
        assert!((input_sum - output_sum).abs() < 1e-10);
    }

    #[test]
    fn test_simd_integral_image() {
        let image = Array2::from_shape_vec(
            (4, 4),
            vec![
                1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0,
                16.0,
            ],
        )
        .expect("operation should succeed");

        let integral =
            simd_compute_integral_image(&image.view()).expect("operation should succeed");

        // Test known integral values
        assert_eq!(integral[[1, 1]], 1.0); // First element
        assert_eq!(integral[[2, 2]], 1.0 + 2.0 + 5.0 + 6.0); // 2x2 top-left
    }

    #[test]
    fn test_simd_statistical_moments() {
        let data = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0];
        let moments = simd_compute_moments(&data, 4).expect("operation should succeed");

        // Test mean
        assert!((moments[1] - 5.5).abs() < 1e-10);

        // Test variance
        let expected_variance = 9.166_666_666_666_666; // Sample variance
        assert!((moments[2] - expected_variance).abs() < 1e-10);
    }

    #[test]
    fn test_simd_entropy() {
        let probabilities = vec![0.5, 0.25, 0.125, 0.125];
        let entropy = simd_compute_entropy(&probabilities).expect("operation should succeed");

        // Expected entropy for this distribution
        let expected = -0.5 * 0.5f64.ln() - 0.25 * 0.25f64.ln() - 2.0 * 0.125 * 0.125f64.ln();
        assert!((entropy - expected).abs() < 1e-12);
    }

    #[test]
    fn test_simd_downsample() {
        let image = Array2::from_shape_vec((4, 4), (0..16).map(|x| x as f64).collect())
            .expect("operation should succeed");

        let downsampled = simd_downsample_2x(&image.view()).expect("operation should succeed");

        assert_eq!(downsampled.dim(), (2, 2));

        // Test that downsampling preserves local averages
        let expected_00 = (0.0 + 1.0 + 4.0 + 5.0) * 0.25;
        assert!((downsampled[[0, 0]] - expected_00).abs() < 1e-12);
    }
}