265 lines
7.3 KiB
Rust
265 lines
7.3 KiB
Rust
//! 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");
|
|
}
|
|
}
|