//! FreeSurfer curvature file reading
//!
//! Reads .curv files containing per-vertex curvature values.
//! Common curvature files:
//! - lh.curv / rh.curv - mean curvature
//! - lh.sulc / rh.sulc - sulcal depth
//! - lh.thickness / rh.thickness - cortical thickness

use std::fs::File;
use std::io::BufReader;
use std::path::Path;

use byteorder::{BigEndian, ReadBytesExt};
use serde::{Deserialize, Serialize};

use crate::error::{AnatomyError, Result};
use crate::surface::Hemisphere;

/// FreeSurfer new curvature file magic number
const NEW_CURV_MAGIC: i32 = -1;

/// Per-vertex curvature data
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Curvature {
    /// Per-vertex values
    pub values: Vec<f32>,
    /// Curvature type (e.g., "curv", "sulc", "thickness")
    pub curv_type: String,
    /// Hemisphere
    pub hemisphere: Hemisphere,
    /// Number of faces (from file header)
    pub n_faces: usize,
}

impl Curvature {
    /// Number of vertices
    pub fn n_vertices(&self) -> usize {
        self.values.len()
    }

    /// Get min and max values
    pub fn range(&self) -> (f32, f32) {
        let min = self.values.iter().copied().fold(f32::MAX, f32::min);
        let max = self.values.iter().copied().fold(f32::MIN, f32::max);
        (min, max)
    }

    /// Get mean value
    pub fn mean(&self) -> f32 {
        if self.values.is_empty() {
            return 0.0;
        }
        self.values.iter().sum::<f32>() / self.values.len() as f32
    }

    /// Normalize values to [0, 1] range
    pub fn normalize(&self) -> Vec<f32> {
        let (min, max) = self.range();
        let range = max - min;
        if range.abs() < 1e-10 {
            return vec![0.5; self.values.len()];
        }
        self.values.iter().map(|v| (v - min) / range).collect()
    }

    /// Get values suitable for colormapping (clamped and normalized)
    ///
    /// # Arguments
    ///
    /// * `vmin` - Minimum value for colormap
    /// * `vmax` - Maximum value for colormap
    pub fn colormap_values(&self, vmin: f32, vmax: f32) -> Vec<f32> {
        let range = vmax - vmin;
        if range.abs() < 1e-10 {
            return vec![0.5; self.values.len()];
        }
        self.values
            .iter()
            .map(|v| (v.clamp(vmin, vmax) - vmin) / range)
            .collect()
    }
}

/// Read a FreeSurfer curvature file
///
/// # Arguments
///
/// * `path` - Path to the curvature file
/// * `curv_type` - Type of curvature (e.g., "curv", "sulc")
/// * `hemisphere` - Hemisphere
///
/// # Returns
///
/// Curvature data with per-vertex values
///
/// # Example
///
/// ```ignore
/// use rtx_neuro_anatomy::curvature::{read_curvature, Hemisphere};
///
/// let curv = read_curvature(
///     "/path/to/subject/surf/lh.curv",
///     "curv",
///     Hemisphere::Left,
/// )?;
/// println!("Mean curvature: {}", curv.mean());
/// ```
pub fn read_curvature(
    path: impl AsRef<Path>,
    curv_type: &str,
    hemisphere: Hemisphere,
) -> Result<Curvature> {
    let path = path.as_ref();
    let file = File::open(path).map_err(|_| AnatomyError::CurvatureNotFound {
        path: path.to_path_buf(),
    })?;
    let mut reader = BufReader::new(file);

    // Read magic number (3 bytes as i24 big-endian)
    let b1 = reader.read_u8()? as i32;
    let b2 = reader.read_u8()? as i32;
    let b3 = reader.read_u8()? as i32;
    let magic = (b1 << 16) | (b2 << 8) | b3;

    // Check for new format
    if magic == (NEW_CURV_MAGIC & 0xFFFFFF) {
        // New format
        read_curvature_new(&mut reader, curv_type, hemisphere)
    } else {
        // Old format: magic is actually n_vertices
        let n_vertices = magic as usize;
        read_curvature_old(&mut reader, n_vertices, curv_type, hemisphere)
    }
}

/// Read new format curvature file
fn read_curvature_new(
    reader: &mut BufReader<File>,
    curv_type: &str,
    hemisphere: Hemisphere,
) -> Result<Curvature> {
    // Read header
    let n_vertices = reader.read_i32::<BigEndian>()? as usize;
    let n_faces = reader.read_i32::<BigEndian>()? as usize;
    let _vals_per_vertex = reader.read_i32::<BigEndian>()?;

    // Read values
    let mut values = Vec::with_capacity(n_vertices);
    for _ in 0..n_vertices {
        values.push(reader.read_f32::<BigEndian>()?);
    }

    Ok(Curvature {
        values,
        curv_type: curv_type.to_string(),
        hemisphere,
        n_faces,
    })
}

/// Read old format curvature file
fn read_curvature_old(
    reader: &mut BufReader<File>,
    n_vertices: usize,
    curv_type: &str,
    hemisphere: Hemisphere,
) -> Result<Curvature> {
    // Read face count from remaining header
    let n_faces = reader.read_i24::<BigEndian>()? as usize;

    // Read values (stored as i16, scale by 100)
    let mut values = Vec::with_capacity(n_vertices);
    for _ in 0..n_vertices {
        let raw = reader.read_i16::<BigEndian>()? as f32;
        values.push(raw / 100.0);
    }

    Ok(Curvature {
        values,
        curv_type: curv_type.to_string(),
        hemisphere,
        n_faces,
    })
}

/// Get curvature filename
pub fn curvature_filename(hemisphere: &Hemisphere, curv_type: &str) -> String {
    format!("{}.{}", hemisphere.prefix(), curv_type)
}

#[cfg(test)]
mod tests {
    use super::*;
    use byteorder::WriteBytesExt;
    use std::io::Write;
    use tempfile::NamedTempFile;

    fn create_test_curvature_new() -> NamedTempFile {
        let mut file = NamedTempFile::new().unwrap();

        // Write magic (0xFFFFFF as 3 bytes)
        file.write_all(&[0xFF, 0xFF, 0xFF]).unwrap();

        // Write header
        file.write_i32::<BigEndian>(4).unwrap(); // n_vertices
        file.write_i32::<BigEndian>(2).unwrap(); // n_faces
        file.write_i32::<BigEndian>(1).unwrap(); // vals_per_vertex

        // Write values
        file.write_f32::<BigEndian>(0.1).unwrap();
        file.write_f32::<BigEndian>(0.2).unwrap();
        file.write_f32::<BigEndian>(0.3).unwrap();
        file.write_f32::<BigEndian>(0.4).unwrap();

        file.flush().unwrap();
        file
    }

    #[test]
    fn test_read_curvature_new() {
        let file = create_test_curvature_new();
        let curv = read_curvature(file.path(), "curv", Hemisphere::Left).unwrap();

        assert_eq!(curv.n_vertices(), 4);
        assert_eq!(curv.n_faces, 2);
        assert!((curv.values[0] - 0.1).abs() < 1e-6);
        assert!((curv.values[3] - 0.4).abs() < 1e-6);
    }

    #[test]
    fn test_curvature_range() {
        let curv = Curvature {
            values: vec![-0.5, 0.0, 0.5, 1.0],
            curv_type: "test".to_string(),
            hemisphere: Hemisphere::Left,
            n_faces: 0,
        };

        let (min, max) = curv.range();
        assert!((min - (-0.5)).abs() < 1e-6);
        assert!((max - 1.0).abs() < 1e-6);
    }

    #[test]
    fn test_curvature_normalize() {
        let curv = Curvature {
            values: vec![0.0, 0.5, 1.0],
            curv_type: "test".to_string(),
            hemisphere: Hemisphere::Left,
            n_faces: 0,
        };

        let normalized = curv.normalize();
        assert!((normalized[0] - 0.0).abs() < 1e-6);
        assert!((normalized[1] - 0.5).abs() < 1e-6);
        assert!((normalized[2] - 1.0).abs() < 1e-6);
    }

    #[test]
    fn test_curvature_filename() {
        assert_eq!(curvature_filename(&Hemisphere::Left, "curv"), "lh.curv");
        assert_eq!(curvature_filename(&Hemisphere::Right, "sulc"), "rh.sulc");
    }
}