oxigeo_algorithms/simd/

texture_simd.rs

//! Texture analysis using Gray-Level Co-occurrence Matrix (GLCM)
//!
//! This module provides implementations of GLCM computation and Haralick feature
//! extraction.
//!
//! # Implementation status
//!
//! GLCM construction is a scatter/accumulate operation; these routines are
//! currently scalar and do NOT contain hand-written SIMD intrinsics. They are
//! correct and unit-tested; vectorization is planned future work.
//!
//! # Supported Operations
//!
//! - **glcm_construct_simd**: SIMD-optimized GLCM matrix construction
//! - **glcm_normalize_simd**: Fast SIMD normalization
//! - **haralick_features_simd**: SIMD-accelerated feature computation
//! - **texture_contrast_simd**: Fast contrast feature extraction
//! - **texture_energy_simd**: Energy/ASM computation with SIMD
//!
//! # Example
//!
//! ```rust
//! # fn main() -> Result<(), Box<dyn std::error::Error>> {
//! use oxigeo_algorithms::simd::texture_simd::glcm_construct_simd;
//!
//! let quantized = vec![0_u8; 1000];
//! let mut glcm = vec![0.0_f32; 256 * 256];
//!
//! glcm_construct_simd(&quantized, &mut glcm, 100, 10, 256, 1, 0)?;
//! # Ok(())
//! # }
//! ```

use crate::error::{AlgorithmError, Result};

/// SIMD-accelerated GLCM construction
///
/// Constructs a Gray-Level Co-occurrence Matrix from quantized image data
/// using SIMD-optimized histogram updates.
///
/// # Arguments
///
/// * `quantized` - Quantized image data (values 0..gray_levels-1)
/// * `glcm` - Output GLCM matrix (gray_levels x gray_levels, row-major)
/// * `width` - Image width
/// * `height` - Image height
/// * `gray_levels` - Number of gray levels
/// * `dx` - X offset for co-occurrence
/// * `dy` - Y offset for co-occurrence
///
/// # Errors
///
/// Returns an error if parameters are invalid
#[allow(clippy::too_many_arguments)]
pub fn glcm_construct_simd(
    quantized: &[u8],
    glcm: &mut [f32],
    width: usize,
    height: usize,
    gray_levels: usize,
    dx: i64,
    dy: i64,
) -> Result<()> {
    if width == 0 || height == 0 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "dimensions",
            message: "Width and height must be greater than zero".to_string(),
        });
    }

    if gray_levels == 0 || gray_levels > 256 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "gray_levels",
            message: "Gray levels must be between 1 and 256".to_string(),
        });
    }

    if quantized.len() != width * height {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "quantized",
            message: "Quantized data size must match width * height".to_string(),
        });
    }

    if glcm.len() != gray_levels * gray_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM size must be gray_levels * gray_levels".to_string(),
        });
    }

    // Initialize GLCM to zero
    const LANES: usize = 8;
    let chunks = glcm.len() / LANES;

    for i in 0..chunks {
        let start = i * LANES;
        let end = start + LANES;
        for j in start..end {
            glcm[j] = 0.0;
        }
    }

    let remainder_start = chunks * LANES;
    for i in remainder_start..glcm.len() {
        glcm[i] = 0.0;
    }

    // Build co-occurrence matrix
    for y in 0..height {
        let ny = (y as i64 + dy) as usize;
        if ny >= height {
            continue;
        }

        for x in 0..width {
            let nx = (x as i64 + dx) as usize;
            if nx >= width {
                continue;
            }

            let i = quantized[y * width + x] as usize;
            let j = quantized[ny * width + nx] as usize;

            if i < gray_levels && j < gray_levels {
                glcm[i * gray_levels + j] += 1.0;
            }
        }
    }

    Ok(())
}

/// SIMD-accelerated GLCM normalization
///
/// Normalizes a GLCM matrix so that all entries sum to 1.0.
///
/// # Arguments
///
/// * `glcm` - GLCM matrix to normalize (modified in-place)
/// * `gray_levels` - Number of gray levels
///
/// # Errors
///
/// Returns an error if the GLCM size is invalid
pub fn glcm_normalize_simd(glcm: &mut [f32], gray_levels: usize) -> Result<()> {
    if glcm.len() != gray_levels * gray_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM size must be gray_levels * gray_levels".to_string(),
        });
    }

    // Compute sum with SIMD
    let mut sum = 0.0_f32;
    const LANES: usize = 8;
    let chunks = glcm.len() / LANES;

    for i in 0..chunks {
        let start = i * LANES;
        let end = start + LANES;

        for j in start..end {
            sum += glcm[j];
        }
    }

    let remainder_start = chunks * LANES;
    for i in remainder_start..glcm.len() {
        sum += glcm[i];
    }

    if sum == 0.0 {
        return Ok(()); // Empty GLCM, nothing to normalize
    }

    // Normalize with SIMD
    for i in 0..chunks {
        let start = i * LANES;
        let end = start + LANES;

        for j in start..end {
            glcm[j] /= sum;
        }
    }

    for i in remainder_start..glcm.len() {
        glcm[i] /= sum;
    }

    Ok(())
}

/// SIMD-accelerated contrast feature computation
///
/// Computes the contrast Haralick feature from a normalized GLCM.
///
/// # Arguments
///
/// * `glcm` - Normalized GLCM matrix
/// * `gray_levels` - Number of gray levels
///
/// # Errors
///
/// Returns an error if the GLCM size is invalid
pub fn texture_contrast_simd(glcm: &[f32], gray_levels: usize) -> Result<f32> {
    if glcm.len() != gray_levels * gray_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM size must be gray_levels * gray_levels".to_string(),
        });
    }

    let mut contrast = 0.0_f32;

    // Compute contrast: sum of (i-j)^2 * P(i,j)
    for i in 0..gray_levels {
        let row_offset = i * gray_levels;

        // SIMD-friendly inner loop
        const LANES: usize = 8;
        let chunks = gray_levels / LANES;

        for chunk in 0..chunks {
            let j_start = chunk * LANES;
            let j_end = j_start + LANES;

            for j in j_start..j_end {
                let diff = (i as i64 - j as i64) as f32;
                contrast += diff * diff * glcm[row_offset + j];
            }
        }

        // Scalar remainder
        let remainder_start = chunks * LANES;
        for j in remainder_start..gray_levels {
            let diff = (i as i64 - j as i64) as f32;
            contrast += diff * diff * glcm[row_offset + j];
        }
    }

    Ok(contrast)
}

/// SIMD-accelerated energy (Angular Second Moment) feature computation
///
/// Computes the energy/ASM Haralick feature from a normalized GLCM.
///
/// # Arguments
///
/// * `glcm` - Normalized GLCM matrix
/// * `gray_levels` - Number of gray levels
///
/// # Errors
///
/// Returns an error if the GLCM size is invalid
pub fn texture_energy_simd(glcm: &[f32], gray_levels: usize) -> Result<f32> {
    if glcm.len() != gray_levels * gray_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM size must be gray_levels * gray_levels".to_string(),
        });
    }

    let mut energy = 0.0_f32;

    // Compute energy: sum of P(i,j)^2
    const LANES: usize = 8;
    let chunks = glcm.len() / LANES;

    for i in 0..chunks {
        let start = i * LANES;
        let end = start + LANES;

        for j in start..end {
            energy += glcm[j] * glcm[j];
        }
    }

    let remainder_start = chunks * LANES;
    for i in remainder_start..glcm.len() {
        energy += glcm[i] * glcm[i];
    }

    Ok(energy)
}

/// SIMD-accelerated entropy feature computation
///
/// Computes the entropy Haralick feature from a normalized GLCM.
///
/// # Arguments
///
/// * `glcm` - Normalized GLCM matrix
/// * `gray_levels` - Number of gray levels
///
/// # Errors
///
/// Returns an error if the GLCM size is invalid
pub fn texture_entropy_simd(glcm: &[f32], gray_levels: usize) -> Result<f32> {
    if glcm.len() != gray_levels * gray_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM size must be gray_levels * gray_levels".to_string(),
        });
    }

    let mut entropy = 0.0_f32;

    // Compute entropy: -sum of P(i,j) * log(P(i,j))
    const LANES: usize = 8;
    let chunks = glcm.len() / LANES;

    for i in 0..chunks {
        let start = i * LANES;
        let end = start + LANES;

        for j in start..end {
            let p = glcm[j];
            if p > 0.0 {
                entropy -= p * p.ln();
            }
        }
    }

    let remainder_start = chunks * LANES;
    for i in remainder_start..glcm.len() {
        let p = glcm[i];
        if p > 0.0 {
            entropy -= p * p.ln();
        }
    }

    Ok(entropy)
}

/// SIMD-accelerated homogeneity (Inverse Difference Moment) feature computation
///
/// Computes the homogeneity Haralick feature from a normalized GLCM.
///
/// # Arguments
///
/// * `glcm` - Normalized GLCM matrix
/// * `gray_levels` - Number of gray levels
///
/// # Errors
///
/// Returns an error if the GLCM size is invalid
pub fn texture_homogeneity_simd(glcm: &[f32], gray_levels: usize) -> Result<f32> {
    if glcm.len() != gray_levels * gray_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM size must be gray_levels * gray_levels".to_string(),
        });
    }

    let mut homogeneity = 0.0_f32;

    // Compute homogeneity: sum of P(i,j) / (1 + (i-j)^2)
    for i in 0..gray_levels {
        let row_offset = i * gray_levels;

        const LANES: usize = 8;
        let chunks = gray_levels / LANES;

        for chunk in 0..chunks {
            let j_start = chunk * LANES;
            let j_end = j_start + LANES;

            for j in j_start..j_end {
                let diff = (i as i64 - j as i64) as f32;
                homogeneity += glcm[row_offset + j] / (1.0 + diff * diff);
            }
        }

        let remainder_start = chunks * LANES;
        for j in remainder_start..gray_levels {
            let diff = (i as i64 - j as i64) as f32;
            homogeneity += glcm[row_offset + j] / (1.0 + diff * diff);
        }
    }

    Ok(homogeneity)
}

/// SIMD-accelerated correlation feature computation
///
/// Computes the correlation Haralick feature from a normalized GLCM.
///
/// # Arguments
///
/// * `glcm` - Normalized GLCM matrix
/// * `gray_levels` - Number of gray levels
///
/// # Errors
///
/// Returns an error if the GLCM size is invalid
pub fn texture_correlation_simd(glcm: &[f32], gray_levels: usize) -> Result<f32> {
    if glcm.len() != gray_levels * gray_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM size must be gray_levels * gray_levels".to_string(),
        });
    }

    // Compute marginal probabilities with SIMD
    let mut px = vec![0.0_f32; gray_levels];
    let mut py = vec![0.0_f32; gray_levels];

    for i in 0..gray_levels {
        let row_offset = i * gray_levels;

        const LANES: usize = 8;
        let chunks = gray_levels / LANES;

        for chunk in 0..chunks {
            let j_start = chunk * LANES;
            let j_end = j_start + LANES;

            for j in j_start..j_end {
                let val = glcm[row_offset + j];
                px[i] += val;
                py[j] += val;
            }
        }

        let remainder_start = chunks * LANES;
        for j in remainder_start..gray_levels {
            let val = glcm[row_offset + j];
            px[i] += val;
            py[j] += val;
        }
    }

    // Compute means
    let mut mu_x = 0.0_f32;
    let mut mu_y = 0.0_f32;

    for i in 0..gray_levels {
        mu_x += i as f32 * px[i];
        mu_y += i as f32 * py[i];
    }

    // Compute standard deviations
    let mut sigma_x = 0.0_f32;
    let mut sigma_y = 0.0_f32;

    for i in 0..gray_levels {
        let dx = i as f32 - mu_x;
        let dy = i as f32 - mu_y;
        sigma_x += dx * dx * px[i];
        sigma_y += dy * dy * py[i];
    }

    sigma_x = sigma_x.sqrt();
    sigma_y = sigma_y.sqrt();

    if sigma_x == 0.0 || sigma_y == 0.0 {
        return Ok(0.0);
    }

    // Compute correlation
    let mut correlation = 0.0_f32;

    for i in 0..gray_levels {
        let row_offset = i * gray_levels;

        const LANES: usize = 8;
        let chunks = gray_levels / LANES;

        for chunk in 0..chunks {
            let j_start = chunk * LANES;
            let j_end = j_start + LANES;

            for j in j_start..j_end {
                let term = ((i as f32 - mu_x) * (j as f32 - mu_y) * glcm[row_offset + j])
                    / (sigma_x * sigma_y);
                correlation += term;
            }
        }

        let remainder_start = chunks * LANES;
        for j in remainder_start..gray_levels {
            let term = ((i as f32 - mu_x) * (j as f32 - mu_y) * glcm[row_offset + j])
                / (sigma_x * sigma_y);
            correlation += term;
        }
    }

    Ok(correlation)
}

/// SIMD-accelerated dissimilarity feature computation
///
/// Computes the dissimilarity Haralick feature from a normalized GLCM.
///
/// # Arguments
///
/// * `glcm` - Normalized GLCM matrix
/// * `gray_levels` - Number of gray levels
///
/// # Errors
///
/// Returns an error if the GLCM size is invalid
pub fn texture_dissimilarity_simd(glcm: &[f32], gray_levels: usize) -> Result<f32> {
    if glcm.len() != gray_levels * gray_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM size must be gray_levels * gray_levels".to_string(),
        });
    }

    let mut dissimilarity = 0.0_f32;

    // Compute dissimilarity: sum of |i-j| * P(i,j)
    for i in 0..gray_levels {
        let row_offset = i * gray_levels;

        const LANES: usize = 8;
        let chunks = gray_levels / LANES;

        for chunk in 0..chunks {
            let j_start = chunk * LANES;
            let j_end = j_start + LANES;

            for j in j_start..j_end {
                let diff = (i as i64 - j as i64).abs() as f32;
                dissimilarity += diff * glcm[row_offset + j];
            }
        }

        let remainder_start = chunks * LANES;
        for j in remainder_start..gray_levels {
            let diff = (i as i64 - j as i64).abs() as f32;
            dissimilarity += diff * glcm[row_offset + j];
        }
    }

    Ok(dissimilarity)
}

/// Complete Haralick features computed with SIMD
#[derive(Debug, Clone, Default)]
pub struct HaralickFeaturesSIMD {
    /// Contrast (variance of differences)
    pub contrast: f32,
    /// Correlation (linear dependency)
    pub correlation: f32,
    /// Energy/Angular Second Moment (uniformity)
    pub energy: f32,
    /// Homogeneity/Inverse Difference Moment (smoothness)
    pub homogeneity: f32,
    /// Entropy (randomness)
    pub entropy: f32,
    /// Dissimilarity (linear contrast)
    pub dissimilarity: f32,
}

/// Compute all major Haralick features with SIMD acceleration
///
/// # Arguments
///
/// * `glcm` - Normalized GLCM matrix
/// * `gray_levels` - Number of gray levels
///
/// # Errors
///
/// Returns an error if parameters are invalid
pub fn compute_haralick_features_simd(
    glcm: &[f32],
    gray_levels: usize,
) -> Result<HaralickFeaturesSIMD> {
    Ok(HaralickFeaturesSIMD {
        contrast: texture_contrast_simd(glcm, gray_levels)?,
        correlation: texture_correlation_simd(glcm, gray_levels)?,
        energy: texture_energy_simd(glcm, gray_levels)?,
        homogeneity: texture_homogeneity_simd(glcm, gray_levels)?,
        entropy: texture_entropy_simd(glcm, gray_levels)?,
        dissimilarity: texture_dissimilarity_simd(glcm, gray_levels)?,
    })
}

/// SIMD-accelerated GLCM computation from u8 data
///
/// Convenience function that computes a normalized GLCM from raw image data.
///
/// # Arguments
///
/// * `data` - Input image data (grayscale 0-255)
/// * `glcm` - Output GLCM matrix (num_levels x num_levels)
/// * `width` - Image width
/// * `height` - Image height
/// * `distance` - Pixel distance for co-occurrence
/// * `angle` - Angle in radians (0 = horizontal)
/// * `num_levels` - Number of gray levels in GLCM
///
/// # Errors
///
/// Returns an error if parameters are invalid
#[allow(clippy::too_many_arguments)]
pub fn compute_glcm_simd(
    data: &[u8],
    glcm: &mut [f32],
    width: usize,
    height: usize,
    distance: i64,
    angle: i64,
    num_levels: usize,
) -> Result<()> {
    if data.len() != width * height {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "data",
            message: "Data length must match width * height".to_string(),
        });
    }

    if glcm.len() != num_levels * num_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM must be num_levels x num_levels".to_string(),
        });
    }

    if num_levels == 0 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "num_levels",
            message: "Number of levels must be positive".to_string(),
        });
    }

    // Convert angle to dx, dy
    let (dx, dy) = match angle {
        0 => (distance, 0),
        45 => (distance, distance),
        90 => (0, distance),
        135 => (-distance, distance),
        _ => (distance, 0), // Default to horizontal
    };

    let scale = 256.0 / num_levels as f32;

    // Initialize GLCM
    for val in glcm.iter_mut() {
        *val = 0.0;
    }

    // Build co-occurrence matrix
    for y in 0..height as i64 {
        for x in 0..width as i64 {
            let nx = x + dx;
            let ny = y + dy;

            if nx >= 0 && nx < width as i64 && ny >= 0 && ny < height as i64 {
                let i = (data[(y as usize) * width + (x as usize)] as f32 / scale) as usize;
                let j = (data[(ny as usize) * width + (nx as usize)] as f32 / scale) as usize;

                let i = i.min(num_levels - 1);
                let j = j.min(num_levels - 1);

                glcm[i * num_levels + j] += 1.0;
            }
        }
    }

    // Normalize GLCM
    let sum: f32 = glcm.iter().sum();
    if sum > 0.0 {
        for val in glcm.iter_mut() {
            *val /= sum;
        }
    }

    Ok(())
}

/// SIMD-accelerated multi-directional GLCM computation
///
/// Computes GLCM averaging over multiple directions for rotation invariance.
///
/// # Arguments
///
/// * `data` - Input image data (grayscale)
/// * `glcm` - Output GLCM matrix (num_levels x num_levels)
/// * `width` - Image width
/// * `height` - Image height
/// * `distance` - Pixel distance for co-occurrence
/// * `num_levels` - Number of gray levels in GLCM
///
/// # Errors
///
/// Returns an error if parameters are invalid
pub fn compute_glcm_multidirectional_simd(
    data: &[u8],
    glcm: &mut [f32],
    width: usize,
    height: usize,
    distance: i64,
    num_levels: usize,
) -> Result<()> {
    if data.len() != width * height {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "data",
            message: "Data length must match width * height".to_string(),
        });
    }

    if glcm.len() != num_levels * num_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM must be num_levels x num_levels".to_string(),
        });
    }

    if num_levels == 0 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "num_levels",
            message: "Number of levels must be positive".to_string(),
        });
    }

    // Initialize GLCM
    for val in glcm.iter_mut() {
        *val = 0.0;
    }

    // Compute GLCM for 4 directions: 0, 45, 90, 135 degrees
    let directions: [(i64, i64); 4] = [
        (distance, 0),         // 0 degrees
        (distance, distance),  // 45 degrees
        (0, distance),         // 90 degrees
        (-distance, distance), // 135 degrees
    ];

    let scale = 256.0 / num_levels as f32;

    for (dx, dy) in &directions {
        for y in 0..height as i64 {
            for x in 0..width as i64 {
                let nx = x + dx;
                let ny = y + dy;

                if nx >= 0 && nx < width as i64 && ny >= 0 && ny < height as i64 {
                    let i = (data[(y as usize) * width + (x as usize)] as f32 / scale) as usize;
                    let j = (data[(ny as usize) * width + (nx as usize)] as f32 / scale) as usize;

                    let i = i.min(num_levels - 1);
                    let j = j.min(num_levels - 1);

                    // Symmetric GLCM
                    glcm[i * num_levels + j] += 1.0;
                    glcm[j * num_levels + i] += 1.0;
                }
            }
        }
    }

    // Normalize GLCM
    let sum: f32 = glcm.iter().sum();
    if sum > 0.0 {
        for val in glcm.iter_mut() {
            *val /= sum;
        }
    }

    Ok(())
}

/// All Haralick texture features
#[derive(Debug, Clone, Copy, Default)]
pub struct TextureFeatures {
    /// Contrast - measure of local variations
    pub contrast: f32,
    /// Dissimilarity - similar to contrast but linear
    pub dissimilarity: f32,
    /// Homogeneity - local homogeneity (inverse difference moment)
    pub homogeneity: f32,
    /// ASM - Angular Second Moment (uniformity)
    pub asm: f32,
    /// Energy - square root of ASM
    pub energy: f32,
    /// Entropy - randomness measure
    pub entropy: f32,
    /// Correlation - linear dependency of gray levels
    pub correlation: f32,
    /// Mean - average gray level
    pub mean: f32,
    /// Variance - gray level variance
    pub variance: f32,
}

/// SIMD-accelerated texture feature extraction from GLCM
///
/// Computes all Haralick texture features from a pre-computed GLCM.
///
/// # Arguments
///
/// * `glcm` - Pre-computed normalized GLCM
/// * `num_levels` - Number of gray levels in GLCM
///
/// # Returns
///
/// A struct containing all texture features
///
/// # Errors
///
/// Returns an error if parameters are invalid
pub fn compute_all_texture_features_simd(
    glcm: &[f32],
    num_levels: usize,
) -> Result<TextureFeatures> {
    if glcm.len() != num_levels * num_levels {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "glcm",
            message: "GLCM must be num_levels x num_levels".to_string(),
        });
    }

    let mut contrast = 0.0_f32;
    let mut dissimilarity = 0.0_f32;
    let mut homogeneity = 0.0_f32;
    let mut asm = 0.0_f32;
    let mut entropy = 0.0_f32;

    // Compute marginal means and standard deviations
    let mut mean_i = 0.0_f32;
    let mut mean_j = 0.0_f32;

    for i in 0..num_levels {
        for j in 0..num_levels {
            let p = glcm[i * num_levels + j];
            mean_i += p * (i as f32);
            mean_j += p * (j as f32);
        }
    }

    let mut std_i = 0.0_f32;
    let mut std_j = 0.0_f32;

    for i in 0..num_levels {
        for j in 0..num_levels {
            let p = glcm[i * num_levels + j];
            std_i += p * (i as f32 - mean_i).powi(2);
            std_j += p * (j as f32 - mean_j).powi(2);
        }
    }

    std_i = std_i.sqrt();
    std_j = std_j.sqrt();

    // Compute correlation
    let mut correlation = 0.0_f32;

    // Compute features
    for i in 0..num_levels {
        for j in 0..num_levels {
            let p = glcm[i * num_levels + j];
            let diff = (i as i32 - j as i32).abs() as f32;

            contrast += p * diff * diff;
            dissimilarity += p * diff;
            homogeneity += p / (1.0 + diff);
            asm += p * p;

            if p > 0.0 {
                entropy -= p * p.ln();
            }

            if std_i > 0.0 && std_j > 0.0 {
                correlation += p * (i as f32 - mean_i) * (j as f32 - mean_j) / (std_i * std_j);
            }
        }
    }

    Ok(TextureFeatures {
        contrast,
        dissimilarity,
        homogeneity,
        asm,
        energy: asm.sqrt(),
        entropy,
        correlation,
        mean: (mean_i + mean_j) / 2.0,
        variance: (std_i * std_i + std_j * std_j) / 2.0,
    })
}

/// Texture feature type for computing feature images
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum TextureFeatureType {
    /// Contrast
    Contrast,
    /// Dissimilarity
    Dissimilarity,
    /// Homogeneity
    Homogeneity,
    /// ASM (Angular Second Moment)
    ASM,
    /// Energy (sqrt of ASM)
    Energy,
    /// Entropy
    Entropy,
    /// Correlation
    Correlation,
}

/// SIMD-accelerated texture feature image computation
///
/// Computes a specific texture feature for each pixel using a sliding window GLCM.
///
/// # Arguments
///
/// * `data` - Input image data (grayscale)
/// * `output` - Output feature image
/// * `width` - Image width
/// * `height` - Image height
/// * `window_size` - Size of the sliding window (must be odd)
/// * `distance` - Pixel distance for co-occurrence
/// * `num_levels` - Number of gray levels in GLCM
/// * `feature` - Which texture feature to compute
///
/// # Errors
///
/// Returns an error if parameters are invalid
#[allow(clippy::too_many_arguments)]
pub fn compute_texture_feature_image_simd(
    data: &[u8],
    output: &mut [f32],
    width: usize,
    height: usize,
    window_size: usize,
    distance: i64,
    num_levels: usize,
    feature: TextureFeatureType,
) -> Result<()> {
    if data.len() != width * height {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "data",
            message: "Data length must match width * height".to_string(),
        });
    }

    if output.len() != width * height {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "output",
            message: "Output length must match width * height".to_string(),
        });
    }

    if window_size.is_multiple_of(2) || window_size < 3 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "window_size",
            message: "Window size must be odd and at least 3".to_string(),
        });
    }

    if num_levels == 0 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "num_levels",
            message: "Number of levels must be positive".to_string(),
        });
    }

    let half_window = window_size / 2;
    let scale = 256.0 / num_levels as f32;

    // Reusable GLCM buffer
    let mut glcm = vec![0.0_f32; num_levels * num_levels];

    // Process each pixel
    for y in 0..height {
        for x in 0..width {
            // Clear GLCM
            for val in glcm.iter_mut() {
                *val = 0.0;
            }

            // Compute GLCM for window
            let y_start = y.saturating_sub(half_window);
            let y_end = (y + half_window + 1).min(height);
            let x_start = x.saturating_sub(half_window);
            let x_end = (x + half_window + 1).min(width);

            let mut count = 0_usize;

            for wy in y_start..y_end {
                for wx in x_start..x_end {
                    let nx = wx as i64 + distance;
                    let ny = wy as i64;

                    if nx >= x_start as i64 && nx < x_end as i64 {
                        let i = (data[wy * width + wx] as f32 / scale) as usize;
                        let j = (data[ny as usize * width + nx as usize] as f32 / scale) as usize;

                        let i = i.min(num_levels - 1);
                        let j = j.min(num_levels - 1);

                        glcm[i * num_levels + j] += 1.0;
                        glcm[j * num_levels + i] += 1.0;
                        count += 2;
                    }
                }
            }

            // Normalize GLCM
            if count > 0 {
                for val in glcm.iter_mut() {
                    *val /= count as f32;
                }
            }

            // Compute feature
            output[y * width + x] = compute_single_texture_feature(&glcm, num_levels, feature);
        }
    }

    Ok(())
}

/// Compute a single texture feature from GLCM
fn compute_single_texture_feature(
    glcm: &[f32],
    num_levels: usize,
    feature: TextureFeatureType,
) -> f32 {
    match feature {
        TextureFeatureType::Contrast => {
            let mut val = 0.0_f32;
            for i in 0..num_levels {
                for j in 0..num_levels {
                    let diff = (i as i32 - j as i32).abs() as f32;
                    val += glcm[i * num_levels + j] * diff * diff;
                }
            }
            val
        }
        TextureFeatureType::Dissimilarity => {
            let mut val = 0.0_f32;
            for i in 0..num_levels {
                for j in 0..num_levels {
                    let diff = (i as i32 - j as i32).abs() as f32;
                    val += glcm[i * num_levels + j] * diff;
                }
            }
            val
        }
        TextureFeatureType::Homogeneity => {
            let mut val = 0.0_f32;
            for i in 0..num_levels {
                for j in 0..num_levels {
                    let diff = (i as i32 - j as i32).abs() as f32;
                    val += glcm[i * num_levels + j] / (1.0 + diff);
                }
            }
            val
        }
        TextureFeatureType::ASM => {
            let mut val = 0.0_f32;
            for p in glcm {
                val += p * p;
            }
            val
        }
        TextureFeatureType::Energy => {
            let asm: f32 = glcm.iter().map(|p| p * p).sum();
            asm.sqrt()
        }
        TextureFeatureType::Entropy => {
            let mut val = 0.0_f32;
            for &p in glcm {
                if p > 0.0 {
                    val -= p * p.ln();
                }
            }
            val
        }
        TextureFeatureType::Correlation => {
            let mut mean_i = 0.0_f32;
            let mut mean_j = 0.0_f32;

            for i in 0..num_levels {
                for j in 0..num_levels {
                    let p = glcm[i * num_levels + j];
                    mean_i += p * (i as f32);
                    mean_j += p * (j as f32);
                }
            }

            let mut std_i = 0.0_f32;
            let mut std_j = 0.0_f32;

            for i in 0..num_levels {
                for j in 0..num_levels {
                    let p = glcm[i * num_levels + j];
                    std_i += p * (i as f32 - mean_i).powi(2);
                    std_j += p * (j as f32 - mean_j).powi(2);
                }
            }

            std_i = std_i.sqrt();
            std_j = std_j.sqrt();

            if std_i > 0.0 && std_j > 0.0 {
                let mut correlation = 0.0_f32;
                for i in 0..num_levels {
                    for j in 0..num_levels {
                        let p = glcm[i * num_levels + j];
                        correlation +=
                            p * (i as f32 - mean_i) * (j as f32 - mean_j) / (std_i * std_j);
                    }
                }
                correlation
            } else {
                0.0
            }
        }
    }
}

/// SIMD-accelerated local binary pattern (LBP) computation
///
/// Computes the Local Binary Pattern texture descriptor.
///
/// # Arguments
///
/// * `data` - Input image data (grayscale)
/// * `output` - Output LBP image
/// * `width` - Image width
/// * `height` - Image height
///
/// # Errors
///
/// Returns an error if parameters are invalid
pub fn compute_lbp_simd(data: &[u8], output: &mut [u8], width: usize, height: usize) -> Result<()> {
    if data.len() != width * height {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "data",
            message: "Data length must match width * height".to_string(),
        });
    }

    if output.len() != width * height {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "output",
            message: "Output length must match width * height".to_string(),
        });
    }

    if width < 3 || height < 3 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "dimensions",
            message: "Width and height must be at least 3".to_string(),
        });
    }

    // Initialize edges to 0
    for x in 0..width {
        output[x] = 0;
        output[(height - 1) * width + x] = 0;
    }

    for y in 0..height {
        output[y * width] = 0;
        output[y * width + width - 1] = 0;
    }

    // Process interior pixels
    for y in 1..(height - 1) {
        let prev_row = (y - 1) * width;
        let curr_row = y * width;
        let next_row = (y + 1) * width;

        for x in 1..(width - 1) {
            let center = data[curr_row + x];
            let mut lbp: u8 = 0;

            // 8-neighborhood pattern (clockwise from top-left)
            if data[prev_row + x - 1] >= center {
                lbp |= 1 << 0;
            }
            if data[prev_row + x] >= center {
                lbp |= 1 << 1;
            }
            if data[prev_row + x + 1] >= center {
                lbp |= 1 << 2;
            }
            if data[curr_row + x + 1] >= center {
                lbp |= 1 << 3;
            }
            if data[next_row + x + 1] >= center {
                lbp |= 1 << 4;
            }
            if data[next_row + x] >= center {
                lbp |= 1 << 5;
            }
            if data[next_row + x - 1] >= center {
                lbp |= 1 << 6;
            }
            if data[curr_row + x - 1] >= center {
                lbp |= 1 << 7;
            }

            output[curr_row + x] = lbp;
        }
    }

    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_abs_diff_eq;

    #[test]
    fn test_glcm_construct() {
        let quantized = vec![0_u8; 100]; // Uniform image
        let mut glcm = vec![0.0_f32; 4 * 4]; // 4 gray levels

        glcm_construct_simd(&quantized, &mut glcm, 10, 10, 4, 1, 0)
            .expect("Failed to construct GLCM for uniform image");

        // For uniform image, all co-occurrences should be at (0,0)
        assert!(glcm[0] > 0.0);
    }

    #[test]
    fn test_glcm_normalize() {
        let mut glcm = vec![2.0_f32; 16]; // 4x4 GLCM with all entries = 2.0

        glcm_normalize_simd(&mut glcm, 4).expect("Failed to normalize GLCM");

        // Sum should be 1.0 after normalization
        let sum: f32 = glcm.iter().sum();
        assert_abs_diff_eq!(sum, 1.0, epsilon = 0.001);
    }

    #[test]
    fn test_texture_energy() {
        let mut glcm = vec![0.0_f32; 16]; // 4x4 GLCM
        glcm[0] = 1.0; // Single entry with probability 1

        let energy = texture_energy_simd(&glcm, 4).expect("Failed to compute texture energy");

        // Energy should be 1.0 for single entry
        assert_abs_diff_eq!(energy, 1.0, epsilon = 0.001);
    }

    #[test]
    fn test_texture_contrast_uniform() {
        let mut glcm = vec![0.0_f32; 16]; // 4x4 GLCM

        // Diagonal matrix (perfect correlation, no contrast)
        glcm[0] = 0.25;
        glcm[5] = 0.25;
        glcm[10] = 0.25;
        glcm[15] = 0.25;

        let contrast = texture_contrast_simd(&glcm, 4).expect("Failed to compute texture contrast");

        // Contrast should be 0 for diagonal matrix
        assert_abs_diff_eq!(contrast, 0.0, epsilon = 0.001);
    }

    #[test]
    fn test_texture_entropy() {
        let glcm = vec![1.0 / 16.0_f32; 16]; // Uniform distribution

        let entropy = texture_entropy_simd(&glcm, 4).expect("Failed to compute texture entropy");

        // Uniform distribution should have maximum entropy
        // H = -sum(p * ln(p)) = -16 * (1/16 * ln(1/16)) = ln(16)
        assert_abs_diff_eq!(entropy, 16.0_f32.ln(), epsilon = 0.01);
    }

    #[test]
    fn test_haralick_features() {
        let mut glcm = vec![0.0_f32; 16]; // 4x4 GLCM
        glcm[0] = 0.5;
        glcm[5] = 0.3;
        glcm[10] = 0.2;

        let features =
            compute_haralick_features_simd(&glcm, 4).expect("Failed to compute Haralick features");

        // All features should be finite
        assert!(features.contrast.is_finite());
        assert!(features.correlation.is_finite());
        assert!(features.energy.is_finite());
        assert!(features.homogeneity.is_finite());
        assert!(features.entropy.is_finite());
        assert!(features.dissimilarity.is_finite());
    }

    #[test]
    fn test_invalid_gray_levels() {
        let quantized = vec![0_u8; 100];
        let mut glcm = vec![0.0_f32; 16];

        // Gray levels = 0 should fail
        let result = glcm_construct_simd(&quantized, &mut glcm, 10, 10, 0, 1, 0);
        assert!(result.is_err());

        // Gray levels > 256 should fail
        let result = glcm_construct_simd(&quantized, &mut glcm, 10, 10, 257, 1, 0);
        assert!(result.is_err());
    }

    #[test]
    fn test_texture_homogeneity() {
        let mut glcm = vec![0.0_f32; 16]; // 4x4 GLCM

        // Diagonal matrix should have maximum homogeneity
        glcm[0] = 0.25;
        glcm[5] = 0.25;
        glcm[10] = 0.25;
        glcm[15] = 0.25;

        let homogeneity =
            texture_homogeneity_simd(&glcm, 4).expect("Failed to compute texture homogeneity");

        // Homogeneity for diagonal should be 1.0
        assert_abs_diff_eq!(homogeneity, 1.0, epsilon = 0.001);
    }

    #[test]
    fn test_glcm_uniform() {
        let data = vec![5u8; 100]; // Uniform image
        let mut glcm = vec![0.0_f32; 256 * 256];

        compute_glcm_simd(&data, &mut glcm, 10, 10, 1, 0, 256)
            .expect("Failed to compute GLCM for uniform image");

        // Uniform image should have all weight at one cell
        let max_val = glcm.iter().copied().fold(f32::NEG_INFINITY, f32::max);
        assert!(max_val > 0.0);
    }

    #[test]
    fn test_glcm_multidirectional() {
        let data = vec![128u8; 100]; // Uniform image
        let mut glcm = vec![0.0_f32; 16 * 16];

        compute_glcm_multidirectional_simd(&data, &mut glcm, 10, 10, 1, 16)
            .expect("Failed to compute multidirectional GLCM");

        // Sum should be 1.0 after normalization
        let sum: f32 = glcm.iter().sum();
        assert_abs_diff_eq!(sum, 1.0, epsilon = 0.001);
    }

    #[test]
    fn test_all_texture_features() {
        let mut glcm = vec![0.0_f32; 16];
        glcm[0] = 0.5;
        glcm[5] = 0.3;
        glcm[10] = 0.2;

        let features = compute_all_texture_features_simd(&glcm, 4)
            .expect("Failed to compute all texture features");

        assert!(features.contrast.is_finite());
        assert!(features.energy.is_finite());
        assert!(features.homogeneity.is_finite());
    }

    #[test]
    fn test_texture_feature_image() {
        let data = vec![128u8; 100];
        let mut output = vec![0.0_f32; 100];

        compute_texture_feature_image_simd(
            &data,
            &mut output,
            10,
            10,
            3,
            1,
            8,
            TextureFeatureType::Contrast,
        )
        .expect("Failed to compute texture feature image");

        // All contrast values should be 0 for uniform image
        for &val in &output {
            assert_abs_diff_eq!(val, 0.0, epsilon = 0.01);
        }
    }

    #[test]
    fn test_lbp_uniform() {
        let data = vec![128u8; 100]; // Uniform image
        let mut output = vec![0u8; 100];

        compute_lbp_simd(&data, &mut output, 10, 10)
            .expect("Failed to compute LBP for uniform image");

        // All interior pixels should be 255 (all neighbors equal to center)
        for y in 1..9 {
            for x in 1..9 {
                assert_eq!(output[y * 10 + x], 255);
            }
        }
    }
}