neuroformats 0.5.0

//! Functions for reading FreeSurfer label files.
//!
//! A label groups a number of vertices (for surface label) or voxels (for volume labels) together. E.g., all
//! vertices which are part of a certain brain region can be stored in a label. Note though that nothing requires that the
//! vertices of a label form a spatially adjacent patch. Each vertex or voxel that is part of the label can be assigned a scalar value.


use std::fs::File;
use std::io::{BufRead, BufReader, Write, LineWriter};
use std::path::{Path};
use std::fmt;


use crate::error::{NeuroformatsError, Result};
use crate::util::vec32minmax;
use crate::config;

#[derive(Debug, Clone, PartialEq)]
pub struct FsLabel {
    pub vertexes: Vec<FsLabelVertex>,
}

impl FsLabel {

    /// Determine whether this is a binary label.
    ///
    /// A binary label assigns the same value (typically `0.0`) to all its vertices.
    /// Such a label is typically used to define a region of some sort, e.g., a single brain region extracted from a brain
    /// surface parcellation (see FsAnnot). Whether or not the label is intended as a binary inside/outside region definition
    /// cannot be known, so treat the return value as an educated guess.
    ///
    /// # Panics
    ///
    /// * If the label is empty, i.e., contains no values.
    pub fn is_binary(&self) -> bool {
        let mut values_iter = self.vertexes.iter().map(|x| x.value);
        let first_val = values_iter.next().expect("Empty label");
        values_iter.all(|val| val == first_val)
    }


    /// Determine for each vertex whether it is part of this label.
    ///
    /// This is a simple convenience function. Note that you need to supply the total number of vertices of
    /// the respective surface, as that number is not stored in the label.
    ///
    /// # Panics
    ///
    /// * If `num_surface_verts` is smaller than the max index stored in the label. If this happens, the label cannot belong to the respective surface.
    pub fn is_surface_vertex_in_label(&self, num_surface_verts: usize) -> Vec<bool> {
        if num_surface_verts < self.vertexes.len() {
            // In this case, num_surface_verts is definitely wrong, but we do not check the max index here, which means stuff can still go wrong below.
            panic!("Invalid vertex count 'num_surface_verts' for surface: label contains {} vertices, surface cannot contain only {}.", self.vertexes.len(), num_surface_verts);
        }
        let mut data_bin = vec![false; num_surface_verts];
        for label_vert in self.vertexes.iter() {
            data_bin[label_vert.index as usize] = true;
        }
        data_bin
    }


    /// Generate data for the whole surface from this label.
    ///
    /// This is a simple convenience function that creates a data vector with the specified length and fills it with the label
    /// value for vertices which are part of this label and sets the rest to the `not_in_label_value` (typically `f32::NAN`).
    ///
    /// # Panics
    ///
    /// * If `num_surface_verts` is smaller than the max index stored in the label. If this happens, the label cannot belong to the respective surface.
    pub fn as_surface_data(&self, num_surface_verts : usize, not_in_label_value : f32) -> Vec<f32> {
        let mut surface_data : Vec<f32> = vec![not_in_label_value; num_surface_verts];
        for surface_vert in self.vertexes.iter() {
            surface_data[surface_vert.index as usize] = surface_vert.value;
        }
        surface_data
    }
}

impl fmt::Display for FsLabel {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        let (min, max) = vec32minmax(self.vertexes.iter().map(|v| v.value), false);
        write!(f, "Label for {} vertices/voxels, with label values in range {} to {}.", self.vertexes.len(), min, max)
    }
}

#[derive(Debug, Clone, PartialEq)]
pub struct FsLabelVertex {
    pub index: i32,
    pub coord1: f32,
    pub coord2: f32,
    pub coord3: f32,
    pub value: f32,
}

// Read a FsLabelVertex from a line.
impl std::str::FromStr for FsLabelVertex {
    type Err = NeuroformatsError;
    fn from_str(s: &str) -> std::result::Result<Self, Self::Err> {
        let mut iter = s.split_whitespace();

        let index_str = iter
            .next()
            .ok_or_else(|| NeuroformatsError::InvalidVertexValue("Missing vertex index".to_string()))?;
        let index = index_str
            .parse::<i32>()
            .map_err(|_| NeuroformatsError::InvalidVertexValue(format!("Invalid vertex index: '{}'", index_str)))?;

        let parse_f32 = |field: &str, label: &str| -> std::result::Result<f32, NeuroformatsError> {
            let val = field
                .parse::<f32>()
                .map_err(|_| NeuroformatsError::InvalidVertexValue(format!("Invalid {} value: '{}'", label, field)))?;
            if !val.is_finite() {
                return Err(NeuroformatsError::InvalidVertexValue(format!(
                    "{} is not finite (value: {})",
                    label, val
                )));
            }
            Ok(val)
        };

        let coord1_str = iter
            .next()
            .ok_or_else(|| NeuroformatsError::InvalidVertexValue("Missing coord1".to_string()))?;
        let coord1 = parse_f32(coord1_str, "coord1")?;

        let coord2_str = iter
            .next()
            .ok_or_else(|| NeuroformatsError::InvalidVertexValue("Missing coord2".to_string()))?;
        let coord2 = parse_f32(coord2_str, "coord2")?;

        let coord3_str = iter
            .next()
            .ok_or_else(|| NeuroformatsError::InvalidVertexValue("Missing coord3".to_string()))?;
        let coord3 = parse_f32(coord3_str, "coord3")?;

        let value_str = iter
            .next()
            .ok_or_else(|| NeuroformatsError::InvalidVertexValue("Missing vertex value".to_string()))?;
        let value = parse_f32(value_str, "vertex value")?;

        Ok(FsLabelVertex {
            index,
            coord1,
            coord2,
            coord3,
            value,
        })
    }
}

/// Read a surface label or volume label from a file in FreeSurfer label format.
///
/// A label groups a number of vertices (for surface label) or voxels (for volume labels) together. It can
/// also assign a scalar value to each element.
///
/// # Examples
///
/// ```no_run
/// let label = neuroformats::read_label("/path/to/subjects_dir/subject1/label/lh.entorhinal_exvivo.label").unwrap();
/// let first = &label.vertexes[0];
/// println!("Vertex #{} has coordinates {} {} {} and is assigned value {}.", first.index, first.coord1, first.coord2, first.coord3, first.value);
/// ```
pub fn read_label<P: AsRef<Path>>(path: P) -> Result<FsLabel> {

    let reader = BufReader::new(File::open(path)?);

    // Read the file line by line using the lines() iterator from std::io::BufRead.
    let mut lines = reader.lines();
    // We ignore the first line at index 0: it is a comment line.
    let _comment_line = lines.next().transpose()?;
    // The line 1 (after comment) is the header
    let hdr_num_entries_str = lines
        .next()
        .transpose()?
        .ok_or_else(|| NeuroformatsError::InvalidFsLabelFormat)?;
    let hdr_num_entries: i32 = hdr_num_entries_str
        .parse::<i32>()
        .map_err(|_| NeuroformatsError::InvalidFsLabelFormat)?;

    if hdr_num_entries < 0 {
        return Err(NeuroformatsError::InvalidHeaderValue(format!(
            "Negative label entry count: {}",
            hdr_num_entries
        )));
    }
    let num_entries = hdr_num_entries as usize;
    if num_entries > config::max_label_entries() {
        return Err(NeuroformatsError::AllocationTooLarge);
    }

    let mut vertexes = Vec::with_capacity(num_entries);
    for line in lines {
        let line = line?;
        let vertex = line.parse()?;
        vertexes.push(vertex);
    }

    if num_entries != vertexes.len() {
        Err(NeuroformatsError::InvalidFsLabelFormat)
    } else {
        Ok(FsLabel { vertexes })
    }
}


/// Write an FsLabel struct to a new file.
pub fn write_label<P: AsRef<Path> + Copy>(path: P, label : &FsLabel) -> std::io::Result<()> {
    let file = File::create(path)?;
    let mut file = LineWriter::new(file);

    let header_lines = format!("# FreeSurfer label.\n{}\n", label.vertexes.len());
    let header_lines = header_lines.as_bytes();
    file.write_all(header_lines)?;

    for vertex in label.vertexes.iter() {
        let vline = format!("{} {} {} {} {}\n", vertex.index, vertex.coord1, vertex.coord2, vertex.coord3, vertex.value);
        let vline = vline.as_bytes();
        file.write_all(vline)?;
    }

    file.flush()?;

    Ok(())
}


#[cfg(test)]
mod test {
    use super::*;
    use tempfile::{tempdir};

    #[test]
    fn the_demo_surface_label_file_can_be_read() {
        const LABEL_FILE: &str = "resources/subjects_dir/subject1/label/lh.entorhinal_exvivo.label";
        let label = read_label(LABEL_FILE).unwrap();

        let expected_vertex_count: usize = 1085;
        assert_eq!(expected_vertex_count, label.vertexes.len());
    }

    #[test]
    fn the_label_utility_functions_work() {
        const LABEL_FILE: &str = "resources/subjects_dir/subject1/label/lh.entorhinal_exvivo.label";
        let label = read_label(LABEL_FILE).unwrap();

        let num_surface_verts: usize = 160_000;
        let label_mask = label.is_surface_vertex_in_label(num_surface_verts);
        assert_eq!(num_surface_verts, label_mask.len());

        let surface_data = label.as_surface_data(num_surface_verts, f32::NAN);
        assert_eq!(num_surface_verts, surface_data.len());

        assert_eq!(false, label.is_binary());
    }

    #[test]
    fn a_label_file_can_be_written_and_reread() {
        const LABEL_FILE: &str = "resources/subjects_dir/subject1/label/lh.entorhinal_exvivo.label";
        let label = read_label(LABEL_FILE).unwrap();

        let dir = tempdir().unwrap();

        let tfile_path = dir.path().join("temp-file.label");
        let tfile_path = tfile_path.to_str().unwrap();
        write_label(tfile_path, &label).unwrap();

        let label_re = read_label(tfile_path).unwrap();
        let expected_vertex_count: usize = 1085;
        assert_eq!(expected_vertex_count, label_re.vertexes.len());
    }

    #[test]
    fn label_vertex_from_str_rejects_nan_coord() {
        let line = "1 NaN 2.0 3.0 0.5";
        let result: std::result::Result<FsLabelVertex, NeuroformatsError> = line.parse();
        assert!(result.is_err());
    }

    #[test]
    fn label_vertex_from_str_rejects_inf_coord() {
        let line = "1 inf 2.0 3.0 0.5";
        let result: std::result::Result<FsLabelVertex, NeuroformatsError> = line.parse();
        assert!(result.is_err());
    }

    #[test]
    fn label_vertex_from_str_rejects_nan_value() {
        let line = "1 1.0 2.0 3.0 NaN";
        let result: std::result::Result<FsLabelVertex, NeuroformatsError> = line.parse();
        assert!(result.is_err());
    }

    #[test]
    fn label_vertex_from_str_rejects_invalid_index() {
        let line = "not_a_number 1.0 2.0 3.0 0.5";
        let result: std::result::Result<FsLabelVertex, NeuroformatsError> = line.parse();
        assert!(result.is_err());
    }

    #[test]
    fn label_vertex_from_str_rejects_missing_fields() {
        let line = "1 1.0 2.0";
        let result: std::result::Result<FsLabelVertex, NeuroformatsError> = line.parse();
        assert!(result.is_err());
    }

    #[test]
    fn label_vertex_from_str_parses_valid_line() {
        let line = "42 1.0 2.5 3.0 0.75";
        let vertex: FsLabelVertex = line.parse().unwrap();
        assert_eq!(vertex.index, 42);
        assert!((vertex.coord1 - 1.0).abs() < 1e-10);
        assert!((vertex.coord2 - 2.5).abs() < 1e-10);
        assert!((vertex.coord3 - 3.0).abs() < 1e-10);
        assert!((vertex.value - 0.75).abs() < 1e-10);
    }
}