trk_io/
cheader.rs

1use std::{
2    fmt,
3    fs::File,
4    io::{BufReader, BufWriter, Error, ErrorKind, Read, Result, Seek, SeekFrom},
5    str::from_utf8,
6};
7
8use byteorder::{BigEndian, ByteOrder, LittleEndian, ReadBytesExt, WriteBytesExt};
9use nalgebra::Vector4;
10
11#[cfg(feature = "nifti_images")]
12use crate::Affine;
13use crate::{
14    orientation::{
15        affine_to_axcodes, axcodes_to_orientations, inverse_orientations_affine,
16        orientations_transform,
17    },
18    Affine4, TrkEndianness,
19};
20
21pub enum Endianness {
22    Little,
23    Big,
24}
25
26impl fmt::Display for Endianness {
27    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
28        match *self {
29            Endianness::Little => write!(f, "< (Little)"),
30            Endianness::Big => write!(f, "> (Big)"),
31        }
32    }
33}
34
35// http://www.trackvis.org/docs/?subsect=fileformat
36#[derive(Clone)]
37#[repr(C)]
38pub struct CHeader {
39    pub id_string: [u8; 6],
40    pub dim: [i16; 3],
41    pub voxel_size: [f32; 3],
42    pub origin: [f32; 3],
43    pub n_scalars: i16,
44    pub scalar_name: [u8; 200], // [10][20]
45    pub n_properties: i16,
46    pub property_name: [u8; 200], // [10][20]
47    pub vox_to_ras: [f32; 16],    // [4][4]
48    pub reserved: [u8; 444],
49    pub voxel_order: [u8; 4],
50    pub pad2: [u8; 4],
51    pub image_orientation_patient: [f32; 6],
52    pub pad1: [u8; 2],
53    pub invert_x: u8,
54    pub invert_y: u8,
55    pub invert_z: u8,
56    pub swap_x: u8,
57    pub swap_y: u8,
58    pub swap_z: u8,
59    pub n_count: i32,
60    pub version: i32,
61    pub hdr_size: i32,
62}
63
64pub const HEADER_SIZE: usize = std::mem::size_of::<CHeader>();
65
66impl CHeader {
67    #[cfg(feature = "nifti_images")]
68    pub fn from_nifti(
69        dim: [u16; 8],
70        pixdim: [f32; 8],
71        srow_x: [f32; 4],
72        srow_y: [f32; 4],
73        srow_z: [f32; 4],
74    ) -> CHeader {
75        #[rustfmt::skip]
76        let affine = Affine::new(
77            srow_x[0], srow_x[1], srow_x[2],
78            srow_y[0], srow_y[1], srow_y[2],
79            srow_z[0], srow_z[1], srow_z[2],
80        );
81        #[rustfmt::skip]
82        let vox_to_ras = [
83            srow_x[0], srow_x[1], srow_x[2], srow_x[3],
84            srow_y[0], srow_y[1], srow_y[2], srow_y[3],
85            srow_z[0], srow_z[1], srow_z[2], srow_z[3],
86            0.0, 0.0, 0.0, 1.0,
87        ];
88        let vo = affine_to_axcodes(&affine).into_bytes();
89        CHeader {
90            dim: [dim[1] as i16, dim[2] as i16, dim[3] as i16],
91            voxel_size: [pixdim[1], pixdim[2], pixdim[3]],
92            vox_to_ras,
93            voxel_order: [vo[0], vo[1], vo[2], 0u8],
94            ..CHeader::default()
95        }
96    }
97
98    pub fn seek_n_count_field(f: &mut BufWriter<File>) -> Result<()> {
99        let n_count_offset = (HEADER_SIZE - 12) as u64;
100        f.seek(SeekFrom::Start(n_count_offset))?;
101        Ok(())
102    }
103
104    pub fn clear_scalars(&mut self) {
105        self.n_scalars = 0;
106        for b in self.scalar_name.iter_mut() {
107            *b = 0;
108        }
109    }
110
111    pub fn add_scalar(&mut self, name: &str) -> Result<()> {
112        if self.n_scalars > 10 {
113            Err(Error::new(ErrorKind::InvalidInput, "Trk header is already full of scalars (10)"))
114        } else if name.len() > 20 {
115            Err(Error::new(ErrorKind::InvalidInput, "New scalar name must be <= 20 characters."))
116        } else if !name.is_ascii() {
117            Err(Error::new(ErrorKind::InvalidInput, "New scalar name must be pure ascii."))
118        } else {
119            let pos = 20 * self.n_scalars as usize;
120            self.scalar_name[pos..pos + name.len()].clone_from_slice(name.as_bytes());
121            self.n_scalars += 1;
122            return Ok(());
123        }
124    }
125
126    pub fn get_scalars_name(&self) -> Vec<String> {
127        read_names(&self.scalar_name, self.n_scalars as usize)
128    }
129
130    pub fn clear_properties(&mut self) {
131        self.n_properties = 0;
132        for b in self.property_name.iter_mut() {
133            *b = 0;
134        }
135    }
136
137    pub fn add_property(&mut self, name: &str) -> Result<()> {
138        if self.n_properties > 10 {
139            Err(Error::new(
140                ErrorKind::InvalidInput,
141                "Trk header is already full of properties (10)",
142            ))
143        } else if name.len() > 20 {
144            Err(Error::new(ErrorKind::InvalidInput, "New property name must be <= 20 characters."))
145        } else if !name.is_ascii() {
146            Err(Error::new(ErrorKind::InvalidInput, "New property name must be pure ascii."))
147        } else {
148            let pos = 20 * self.n_properties as usize;
149            self.property_name[pos..pos + name.len()].clone_from_slice(name.as_bytes());
150            self.n_properties += 1;
151            return Ok(());
152        }
153    }
154
155    pub fn get_properties_name(&self) -> Vec<String> {
156        read_names(&self.property_name, self.n_properties as usize)
157    }
158
159    /// Get affine mapping trackvis voxelmm space to RAS+ mm space
160    ///
161    /// The streamlines in a trackvis file are in 'voxelmm' space, where the coordinates refer to
162    /// the corner of the voxel.
163    ///
164    /// Compute the affine matrix that will bring them back to RAS+ mm space, where the coordinates
165    /// refer to the center of the voxel.
166    pub fn get_affine_to_rasmm(&self) -> Affine4 {
167        let mut affine = Affine4::identity();
168
169        let scale = Affine4::from_diagonal(&Vector4::new(
170            1.0 / self.voxel_size[0],
171            1.0 / self.voxel_size[1],
172            1.0 / self.voxel_size[2],
173            1.0,
174        ));
175        affine = scale * affine;
176
177        #[rustfmt::skip]
178        let offset = Affine4::new(
179            1.0, 0.0, 0.0, -0.5,
180            0.0, 1.0, 0.0, -0.5,
181            0.0, 0.0, 1.0, -0.5,
182            0.0, 0.0, 0.0, 1.0,
183        );
184        affine = offset * affine;
185
186        let voxel_to_rasmm = Affine4::from_iterator(self.vox_to_ras.iter().cloned()).transpose();
187
188        let header_ornt = axcodes_to_orientations(from_utf8(&self.voxel_order).unwrap());
189        let affine_order = affine_to_axcodes(&voxel_to_rasmm.fixed_view::<3, 3>(0, 0).into_owned());
190        let affine_ornt = axcodes_to_orientations(&affine_order);
191        let orientations = orientations_transform(&header_ornt, &affine_ornt);
192        let inv = inverse_orientations_affine(&orientations, self.dim);
193        affine = inv * affine;
194
195        voxel_to_rasmm * affine
196    }
197
198    pub fn read(reader: &mut BufReader<File>) -> Result<(CHeader, Endianness)> {
199        reader.seek(SeekFrom::Start(0))?;
200        let endianness = test_endianness(reader)?;
201        let header = match endianness {
202            Endianness::Little => CHeader::read_::<LittleEndian>(reader)?,
203            Endianness::Big => CHeader::read_::<BigEndian>(reader)?,
204        };
205        Ok((header, endianness))
206    }
207
208    fn read_<E: ByteOrder>(reader: &mut BufReader<File>) -> Result<CHeader> {
209        let mut header = CHeader::default();
210
211        // Make sure that the file signature (magic number) is right before doing anything else
212        reader.read_exact(&mut header.id_string)?;
213        if &header.id_string != b"TRACK\0" {
214            return Err(Error::new(
215                ErrorKind::InvalidData,
216                "Not a TrackVis file (wrong file signature)",
217            ));
218        }
219
220        for i in &mut header.dim {
221            *i = reader.read_i16::<E>()?;
222        }
223        for f in &mut header.voxel_size {
224            *f = reader.read_f32::<E>()?;
225        }
226        for f in &mut header.origin {
227            *f = reader.read_f32::<E>()?;
228        }
229        header.n_scalars = reader.read_i16::<E>()?;
230        reader.read_exact(&mut header.scalar_name)?;
231        header.n_properties = reader.read_i16::<E>()?;
232        reader.read_exact(&mut header.property_name)?;
233        for f in &mut header.vox_to_ras {
234            *f = reader.read_f32::<E>()?;
235        }
236        reader.read_exact(&mut header.reserved)?;
237        reader.read_exact(&mut header.voxel_order)?;
238        reader.read_exact(&mut header.pad2)?;
239        for f in &mut header.image_orientation_patient {
240            *f = reader.read_f32::<E>()?;
241        }
242        reader.read_exact(&mut header.pad1)?;
243        header.invert_x = reader.read_u8()?;
244        header.invert_y = reader.read_u8()?;
245        header.invert_z = reader.read_u8()?;
246        header.swap_x = reader.read_u8()?;
247        header.swap_y = reader.read_u8()?;
248        header.swap_z = reader.read_u8()?;
249        header.n_count = reader.read_i32::<E>()?;
250        header.version = reader.read_i32::<E>()?;
251        header.hdr_size = reader.read_i32::<E>()?;
252
253        Ok(header)
254    }
255
256    pub fn write<W: WriteBytesExt>(&self, writer: &mut W) -> Result<()> {
257        writer.write(&self.id_string)?;
258        for i in &self.dim {
259            writer.write_i16::<TrkEndianness>(*i)?;
260        }
261        for f in &self.voxel_size {
262            writer.write_f32::<TrkEndianness>(*f)?;
263        }
264        for f in &self.origin {
265            writer.write_f32::<TrkEndianness>(*f)?;
266        }
267        writer.write_i16::<TrkEndianness>(self.n_scalars)?;
268        writer.write(&self.scalar_name)?;
269        writer.write_i16::<TrkEndianness>(self.n_properties)?;
270        writer.write(&self.property_name)?;
271        for f in &self.vox_to_ras {
272            writer.write_f32::<TrkEndianness>(*f)?;
273        }
274        writer.write(&self.reserved)?;
275        writer.write(&self.voxel_order)?;
276        writer.write(&self.pad2)?;
277        for f in &self.image_orientation_patient {
278            writer.write_f32::<TrkEndianness>(*f)?;
279        }
280        writer.write(&self.pad1)?;
281        writer.write_u8(self.invert_x)?;
282        writer.write_u8(self.invert_y)?;
283        writer.write_u8(self.invert_z)?;
284        writer.write_u8(self.swap_x)?;
285        writer.write_u8(self.swap_y)?;
286        writer.write_u8(self.swap_z)?;
287        writer.write_i32::<TrkEndianness>(self.n_count)?;
288        writer.write_i32::<TrkEndianness>(self.version)?;
289        writer.write_i32::<TrkEndianness>(self.hdr_size)?;
290
291        Ok(())
292    }
293}
294
295impl Default for CHeader {
296    fn default() -> CHeader {
297        CHeader {
298            id_string: *b"TRACK\0",
299            dim: [0, 0, 0],
300            voxel_size: [1.0, 1.0, 1.0],
301            origin: [0.0, 0.0, 0.0],
302            n_scalars: 0,
303            scalar_name: [0; 200],
304            n_properties: 0,
305            property_name: [0; 200],
306            vox_to_ras: [
307                1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0,
308            ],
309            reserved: [0; 444],
310            voxel_order: [82, 65, 83, 0],
311            pad2: [0; 4],
312            image_orientation_patient: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
313            pad1: [0; 2],
314            invert_x: 0,
315            invert_y: 0,
316            invert_z: 0,
317            swap_x: 0,
318            swap_y: 0,
319            swap_z: 0,
320            n_count: 0,
321            version: 2,
322            hdr_size: HEADER_SIZE as i32,
323        }
324    }
325}
326
327/// Returns the endianness used when saving the trk file read by `reader`
328///
329/// We use `version` to discover the endianness because it's the biggest
330/// integer field with the most constrained possible values {1, 2}.
331/// Read in LittleEndian, version == 1 or 2.
332/// Read in BigEndian, version == 511 or 767
333/// Even with hundreds major updates, `version` should be safe.
334fn test_endianness(reader: &mut BufReader<File>) -> Result<Endianness> {
335    let version_offset = (HEADER_SIZE - 8) as u64;
336    reader.seek(SeekFrom::Start(version_offset))?;
337    let version = reader.read_i32::<LittleEndian>()?;
338    let endianness = if version <= 255 { Endianness::Little } else { Endianness::Big };
339    reader.seek(SeekFrom::Start(0))?;
340
341    Ok(endianness)
342}
343
344/// Returns the names from the [10][20] arrays of bytes.
345///
346/// Normal case: name\0\0...
347/// Special case: name\0{number}\0\0...
348fn read_names(names_bytes: &[u8], nb: usize) -> Vec<String> {
349    let mut at = 0;
350    let mut names = vec![String::from(""); nb];
351    for names_byte in names_bytes.chunks(20) {
352        if names_byte[0] == 0u8 {
353            break;
354        }
355
356        let idx = names_byte.iter().position(|&e| e == 0u8).unwrap_or(20);
357        let name = from_utf8(&names_byte[..idx]).unwrap().to_string();
358        if idx < 19 && names_byte[idx + 1] != 0u8 {
359            let number = names_byte[idx + 1] - 48;
360            for _ in 0..number {
361                names[at] = name.clone();
362                at += 1;
363            }
364        } else {
365            names[at] = name;
366            at += 1;
367        }
368    }
369    names
370}
371
372#[cfg(test)]
373mod tests {
374    use super::*;
375
376    #[test]
377    fn test_scalars() {
378        let mut header = CHeader::default();
379        header.add_scalar("color_x").unwrap();
380        header.add_scalar("color_y").unwrap();
381
382        let mut gt = [0u8; 200];
383        gt[..7].clone_from_slice(b"color_x");
384        gt[20..27].clone_from_slice(b"color_y");
385        assert_eq!(&header.scalar_name[..], &gt[..]);
386    }
387
388    #[test]
389    fn test_read_empty_names() {
390        // N scalars/properties without a empty description should still return a vector of N
391        // empty strings. It's not super practical, but that's the best we can do with such data.
392        let scalars = read_names(&vec![0; 80], 3);
393        assert_eq!(scalars, vec![String::from(""), String::from(""), String::from("")]);
394    }
395
396    #[test]
397    fn test_header_size() {
398        assert_eq!(HEADER_SIZE, 1000);
399    }
400}
trk_io/cheader.rs

trk_io/
cheader.rs
