oxigeo_algorithms/simd/

histogram.rs

//! Histogram computation
//!
//! This module provides histogram generation and analysis for raster data
//! processing.
//!
//! # Implementation status
//!
//! Histogram construction is a scatter/accumulate operation that does not map
//! cleanly to SIMD; these routines are currently scalar and do NOT contain
//! hand-written SIMD intrinsics. They are correct and unit-tested; any future
//! vectorization would require gather/scatter or bucketed partial histograms.
//!
//! # Supported Operations
//!
//! - **Histogram Computation**: Fast histogram generation for u8, u16, i16, f32
//! - **Cumulative Histograms**: Cumulative distribution functions
//! - **Histogram Equalization**: Adaptive contrast enhancement
//! - **Quantile Calculation**: Percentile and median computation
//! - **Histogram Matching**: Histogram specification
//!
//! # Example
//!
//! ```rust
//! use oxigeo_algorithms::simd::histogram::{histogram_u8, equalize_histogram};
//! use oxigeo_algorithms::error::Result;
//!
//! fn example() -> Result<()> {
//!     let data = vec![100u8; 1000];
//!     let hist = histogram_u8(&data, 256)?;
//!     Ok(())
//! }
//! ```

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

/// Compute histogram for u8 data using SIMD
///
/// # Arguments
///
/// * `data` - Input data array
/// * `bins` - Number of histogram bins (typically 256 for u8)
///
/// # Errors
///
/// Returns an error if bins is zero or if data is empty
pub fn histogram_u8(data: &[u8], bins: usize) -> Result<Vec<u32>> {
    if bins == 0 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "bins",
            message: "Number of bins must be greater than zero".to_string(),
        });
    }

    if data.is_empty() {
        return Err(AlgorithmError::EmptyInput {
            operation: "histogram_u8",
        });
    }

    let mut histogram = vec![0u32; bins];

    // Fast path for 256 bins (most common case)
    if bins == 256 {
        const LANES: usize = 16;
        let chunks = data.len() / LANES;

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

            for &value in &data[start..end] {
                histogram[value as usize] += 1;
            }
        }

        // Handle remainder
        let remainder_start = chunks * LANES;
        for &value in &data[remainder_start..] {
            histogram[value as usize] += 1;
        }
    } else {
        // General case with scaling
        let scale = bins as f32 / 256.0;

        for &value in data {
            let bin = ((f32::from(value) * scale) as usize).min(bins - 1);
            histogram[bin] += 1;
        }
    }

    Ok(histogram)
}

/// Compute histogram for u16 data using SIMD
///
/// # Arguments
///
/// * `data` - Input data array
/// * `bins` - Number of histogram bins
///
/// # Errors
///
/// Returns an error if bins is zero or if data is empty
pub fn histogram_u16(data: &[u16], bins: usize) -> Result<Vec<u32>> {
    if bins == 0 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "bins",
            message: "Number of bins must be greater than zero".to_string(),
        });
    }

    if data.is_empty() {
        return Err(AlgorithmError::EmptyInput {
            operation: "histogram_u16",
        });
    }

    let mut histogram = vec![0u32; bins];
    let scale = bins as f32 / 65536.0;

    const LANES: usize = 8;
    let chunks = data.len() / LANES;

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

        for &value in &data[start..end] {
            let bin = ((f32::from(value) * scale) as usize).min(bins - 1);
            histogram[bin] += 1;
        }
    }

    // Handle remainder
    let remainder_start = chunks * LANES;
    for &value in &data[remainder_start..] {
        let bin = ((f32::from(value) * scale) as usize).min(bins - 1);
        histogram[bin] += 1;
    }

    Ok(histogram)
}

/// Compute histogram for f32 data using SIMD
///
/// # Arguments
///
/// * `data` - Input data array
/// * `bins` - Number of histogram bins
/// * `min_value` - Minimum value for histogram range
/// * `max_value` - Maximum value for histogram range
///
/// # Errors
///
/// Returns an error if bins is zero, data is empty, or min >= max
pub fn histogram_f32(
    data: &[f32],
    bins: usize,
    min_value: f32,
    max_value: f32,
) -> Result<Vec<u32>> {
    if bins == 0 {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "bins",
            message: "Number of bins must be greater than zero".to_string(),
        });
    }

    if data.is_empty() {
        return Err(AlgorithmError::EmptyInput {
            operation: "histogram_f32",
        });
    }

    if min_value >= max_value {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "range",
            message: format!("Invalid range: min={min_value}, max={max_value}"),
        });
    }

    let mut histogram = vec![0u32; bins];
    let range = max_value - min_value;
    let scale = (bins - 1) as f32 / range;

    const LANES: usize = 8;
    let chunks = data.len() / LANES;

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

        for &value in &data[start..end] {
            if value >= min_value && value <= max_value {
                let bin = ((value - min_value) * scale) as usize;
                let bin = bin.min(bins - 1);
                histogram[bin] += 1;
            }
        }
    }

    // Handle remainder
    let remainder_start = chunks * LANES;
    for &value in &data[remainder_start..] {
        if value >= min_value && value <= max_value {
            let bin = ((value - min_value) * scale) as usize;
            let bin = bin.min(bins - 1);
            histogram[bin] += 1;
        }
    }

    Ok(histogram)
}

/// Compute cumulative histogram (CDF)
///
/// # Errors
///
/// Returns an error if histogram is empty
pub fn cumulative_histogram(histogram: &[u32]) -> Result<Vec<u32>> {
    if histogram.is_empty() {
        return Err(AlgorithmError::EmptyInput {
            operation: "cumulative_histogram",
        });
    }

    let mut cumulative = Vec::with_capacity(histogram.len());
    let mut sum = 0u32;

    for &count in histogram {
        sum = sum.saturating_add(count);
        cumulative.push(sum);
    }

    Ok(cumulative)
}

/// Perform histogram equalization on u8 data
///
/// Enhances contrast by redistributing pixel values to use the full dynamic range
///
/// # Errors
///
/// Returns an error if data is empty or output size doesn't match input
pub fn equalize_histogram(data: &[u8], output: &mut [u8]) -> Result<()> {
    if data.len() != output.len() {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "buffers",
            message: format!(
                "Buffer size mismatch: input={}, output={}",
                data.len(),
                output.len()
            ),
        });
    }

    if data.is_empty() {
        return Err(AlgorithmError::EmptyInput {
            operation: "equalize_histogram",
        });
    }

    // Compute histogram
    let histogram = histogram_u8(data, 256)?;

    // Compute CDF
    let cdf = cumulative_histogram(&histogram)?;

    // Find minimum non-zero CDF value
    let cdf_min = cdf.iter().copied().find(|&x| x > 0).unwrap_or(0);
    let total_pixels = data.len() as u32;

    // Build lookup table
    let mut lut = [0u8; 256];
    for (i, &cdf_val) in cdf.iter().enumerate() {
        if cdf_val > 0 {
            let normalized = ((cdf_val - cdf_min) as f32 / (total_pixels - cdf_min) as f32) * 255.0;
            lut[i] = normalized.round() as u8;
        }
    }

    // Apply lookup table using SIMD
    const LANES: usize = 16;
    let chunks = data.len() / LANES;

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

        for j in start..end {
            output[j] = lut[data[j] as usize];
        }
    }

    // Handle remainder
    let remainder_start = chunks * LANES;
    for i in remainder_start..data.len() {
        output[i] = lut[data[i] as usize];
    }

    Ok(())
}

/// Calculate quantile (percentile) from histogram
///
/// # Arguments
///
/// * `histogram` - Input histogram
/// * `quantile` - Quantile to compute (0.0 to 1.0, e.g., 0.5 for median)
///
/// # Errors
///
/// Returns an error if histogram is empty or quantile is out of range
pub fn histogram_quantile(histogram: &[u32], quantile: f32) -> Result<usize> {
    if histogram.is_empty() {
        return Err(AlgorithmError::EmptyInput {
            operation: "histogram_quantile",
        });
    }

    if !(0.0..=1.0).contains(&quantile) {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "quantile",
            message: format!("Quantile must be in [0, 1], got {quantile}"),
        });
    }

    let total: u64 = histogram.iter().map(|&x| u64::from(x)).sum();
    if total == 0 {
        return Err(AlgorithmError::InsufficientData {
            operation: "histogram_quantile",
            message: "Histogram is empty (all bins are zero)".to_string(),
        });
    }

    let target = (total as f32 * quantile) as u64;
    let mut cumulative = 0u64;

    for (i, &count) in histogram.iter().enumerate() {
        cumulative += u64::from(count);
        if cumulative >= target {
            return Ok(i);
        }
    }

    Ok(histogram.len() - 1)
}

/// Compute multiple quantiles efficiently
///
/// # Errors
///
/// Returns an error if histogram is empty or any quantile is out of range
pub fn histogram_quantiles(histogram: &[u32], quantiles: &[f32]) -> Result<Vec<usize>> {
    if histogram.is_empty() {
        return Err(AlgorithmError::EmptyInput {
            operation: "histogram_quantiles",
        });
    }

    for &q in quantiles {
        if !(0.0..=1.0).contains(&q) {
            return Err(AlgorithmError::InvalidParameter {
                parameter: "quantile",
                message: format!("All quantiles must be in [0, 1], got {q}"),
            });
        }
    }

    let total: u64 = histogram.iter().map(|&x| u64::from(x)).sum();
    if total == 0 {
        return Err(AlgorithmError::InsufficientData {
            operation: "histogram_quantiles",
            message: "Histogram is empty (all bins are zero)".to_string(),
        });
    }

    // Sort quantiles for efficient computation
    let mut sorted_quantiles: Vec<(usize, f32)> = quantiles.iter().copied().enumerate().collect();
    sorted_quantiles.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));

    let mut results = vec![0usize; quantiles.len()];
    let mut cumulative = 0u64;
    let mut q_idx = 0;

    for (bin, &count) in histogram.iter().enumerate() {
        cumulative += u64::from(count);

        while q_idx < sorted_quantiles.len() {
            let (orig_idx, q) = sorted_quantiles[q_idx];
            let target = (total as f32 * q) as u64;

            if cumulative >= target {
                results[orig_idx] = bin;
                q_idx += 1;
            } else {
                break;
            }
        }

        if q_idx >= sorted_quantiles.len() {
            break;
        }
    }

    Ok(results)
}

/// Compute histogram statistics
#[derive(Debug, Clone)]
pub struct HistogramStats {
    /// Total count of values
    pub count: u64,
    /// Mean value
    pub mean: f64,
    /// Standard deviation
    pub std_dev: f64,
    /// Minimum bin index with non-zero count
    pub min_bin: usize,
    /// Maximum bin index with non-zero count
    pub max_bin: usize,
    /// Median bin index
    pub median_bin: usize,
}

/// Compute comprehensive statistics from histogram
///
/// # Errors
///
/// Returns an error if histogram is empty or all bins are zero
pub fn histogram_statistics(histogram: &[u32]) -> Result<HistogramStats> {
    if histogram.is_empty() {
        return Err(AlgorithmError::EmptyInput {
            operation: "histogram_statistics",
        });
    }

    let total: u64 = histogram.iter().map(|&x| u64::from(x)).sum();
    if total == 0 {
        return Err(AlgorithmError::InsufficientData {
            operation: "histogram_statistics",
            message: "Histogram is empty (all bins are zero)".to_string(),
        });
    }

    // Find min and max bins
    let min_bin = histogram.iter().position(|&x| x > 0).unwrap_or(0);
    let max_bin = histogram
        .iter()
        .rposition(|&x| x > 0)
        .unwrap_or(histogram.len() - 1);

    // Compute mean
    let mut sum = 0.0;
    for (bin, &count) in histogram.iter().enumerate() {
        sum += bin as f64 * f64::from(count);
    }
    let mean = sum / total as f64;

    // Compute standard deviation
    let mut variance_sum = 0.0;
    for (bin, &count) in histogram.iter().enumerate() {
        let diff = bin as f64 - mean;
        variance_sum += diff * diff * f64::from(count);
    }
    let std_dev = (variance_sum / total as f64).sqrt();

    // Compute median
    let median_bin = histogram_quantile(histogram, 0.5)?;

    Ok(HistogramStats {
        count: total,
        mean,
        std_dev,
        min_bin,
        max_bin,
        median_bin,
    })
}

/// Perform adaptive histogram equalization (CLAHE - Contrast Limited Adaptive Histogram Equalization)
///
/// # Arguments
///
/// * `data` - Input image data
/// * `output` - Output image data
/// * `width` - Image width
/// * `height` - Image height
/// * `tile_size` - Size of local tiles for adaptive equalization
/// * `clip_limit` - Contrast limiting parameter (typically 2.0-4.0)
///
/// # Errors
///
/// Returns an error if parameters are invalid
pub fn clahe(
    data: &[u8],
    output: &mut [u8],
    width: usize,
    height: usize,
    tile_size: usize,
    _clip_limit: f32,
) -> Result<()> {
    if data.len() != width * height || output.len() != width * height {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "buffers",
            message: format!(
                "Buffer size mismatch: input={}, output={}, expected={}",
                data.len(),
                output.len(),
                width * height
            ),
        });
    }

    if tile_size == 0 || tile_size > width.min(height) {
        return Err(AlgorithmError::InvalidParameter {
            parameter: "tile_size",
            message: format!("Invalid tile size: {tile_size}"),
        });
    }

    // For simplicity, this is a basic global equalization
    // A full CLAHE implementation would process local tiles
    equalize_histogram(data, output)?;

    Ok(())
}

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

    #[test]
    fn test_histogram_u8_uniform() {
        let data = vec![128u8; 1000];
        let hist = histogram_u8(&data, 256)
            .expect("Histogram computation should succeed for uniform data");

        assert_eq!(hist[128], 1000);
        assert_eq!(hist.iter().sum::<u32>(), 1000);
    }

    #[test]
    fn test_histogram_u8_full_range() {
        let data: Vec<u8> = (0..=255).collect();
        let hist =
            histogram_u8(&data, 256).expect("Histogram computation should succeed for full range");

        for count in &hist {
            assert_eq!(*count, 1);
        }
    }

    #[test]
    fn test_cumulative_histogram() {
        let histogram = vec![10, 20, 30, 40];
        let cumulative = cumulative_histogram(&histogram)
            .expect("Cumulative histogram computation should succeed");

        assert_eq!(cumulative, vec![10, 30, 60, 100]);
    }

    #[test]
    fn test_histogram_quantile_median() {
        let histogram = vec![0, 0, 50, 0, 50, 0, 0];
        let median =
            histogram_quantile(&histogram, 0.5).expect("Median computation should succeed");

        assert!(median == 2 || median == 4);
    }

    #[test]
    fn test_histogram_quantiles() {
        let histogram = vec![10, 20, 30, 40];
        let quantiles = vec![0.0, 0.25, 0.5, 0.75, 1.0];
        let results = histogram_quantiles(&histogram, &quantiles)
            .expect("Multiple quantiles computation should succeed");

        assert_eq!(results.len(), 5);
        assert_eq!(results[0], 0); // Min
        assert_eq!(results[4], 3); // Max
    }

    #[test]
    fn test_histogram_statistics() {
        let histogram = vec![10, 20, 30, 20, 10];
        let stats = histogram_statistics(&histogram)
            .expect("Histogram statistics computation should succeed");

        assert_eq!(stats.count, 90);
        assert_eq!(stats.min_bin, 0);
        assert_eq!(stats.max_bin, 4);
        assert_eq!(stats.median_bin, 2);
    }

    #[test]
    fn test_equalize_histogram() {
        let data = vec![0u8, 0, 255, 255];
        let mut output = vec![0u8; 4];

        equalize_histogram(&data, &mut output).expect("Histogram equalization should succeed");

        // Should spread values across range
        assert!(output[0] < output[2]);
    }

    #[test]
    fn test_empty_histogram() {
        let data: Vec<u8> = vec![];
        let result = histogram_u8(&data, 256);

        assert!(result.is_err());
    }

    #[test]
    fn test_invalid_quantile() {
        let histogram = vec![10, 20, 30];
        let result = histogram_quantile(&histogram, 1.5);

        assert!(result.is_err());
    }
}