1#[cfg(feature = "nalgebra_affine")]
5use crate::affine::*;
6use crate::error::{NiftiError, Result};
7use crate::typedef::*;
8use crate::util::{is_gz_file, validate_dim, validate_dimensionality};
9use byteordered::{ByteOrdered, Endian, Endianness};
10use flate2::bufread::GzDecoder;
11#[cfg(feature = "nalgebra_affine")]
12use nalgebra::{Matrix3, Matrix4, Quaternion, RealField, Vector3};
13use num_traits::FromPrimitive;
14#[cfg(feature = "nalgebra_affine")]
15use num_traits::ToPrimitive;
16#[cfg(feature = "nalgebra_affine")]
17use simba::scalar::SubsetOf;
18use std::fs::File;
19use std::io::{BufReader, Read};
20use std::ops::Deref;
21use std::path::Path;
22
23pub const MAGIC_CODE_NI1: &[u8; 4] = b"ni1\0";
25pub const MAGIC_CODE_NIP1: &[u8; 4] = b"n+1\0";
27
28#[derive(Debug, Clone, PartialEq)]
60pub struct NiftiHeader {
61 pub sizeof_hdr: i32,
63 pub data_type: [u8; 10],
65 pub db_name: [u8; 18],
67 pub extents: i32,
69 pub session_error: i16,
71 pub regular: u8,
73 pub dim_info: u8,
75 pub dim: [u16; 8],
77 pub intent_p1: f32,
79 pub intent_p2: f32,
81 pub intent_p3: f32,
83 pub intent_code: i16,
85 pub datatype: i16,
87 pub bitpix: i16,
89 pub slice_start: i16,
91 pub pixdim: [f32; 8],
93 pub vox_offset: f32,
98 pub scl_slope: f32,
100 pub scl_inter: f32,
102 pub slice_end: i16,
104 pub slice_code: u8,
106 pub xyzt_units: u8,
108 pub cal_max: f32,
110 pub cal_min: f32,
112 pub slice_duration: f32,
114 pub toffset: f32,
116 pub glmax: i32,
118 pub glmin: i32,
120
121 pub descrip: Vec<u8>,
123 pub aux_file: [u8; 24],
125 pub qform_code: i16,
127 pub sform_code: i16,
129 pub quatern_b: f32,
131 pub quatern_c: f32,
133 pub quatern_d: f32,
135 pub quatern_x: f32,
137 pub quatern_y: f32,
139 pub quatern_z: f32,
141
142 pub srow_x: [f32; 4],
144 pub srow_y: [f32; 4],
146 pub srow_z: [f32; 4],
148
149 pub intent_name: [u8; 16],
151
152 pub magic: [u8; 4],
154
155 pub endianness: Endianness,
157}
158
159impl Default for NiftiHeader {
160 fn default() -> NiftiHeader {
161 NiftiHeader {
162 sizeof_hdr: 348,
163 data_type: [0; 10],
164 db_name: [0; 18],
165 extents: 0,
166 session_error: 0,
167 regular: 0,
168 dim_info: 0,
169 dim: [1, 0, 0, 0, 0, 0, 0, 0],
170 intent_p1: 0.,
171 intent_p2: 0.,
172 intent_p3: 0.,
173 intent_code: 0,
174 datatype: 0,
175 bitpix: 0,
176 slice_start: 0,
177 pixdim: [1.; 8],
178 vox_offset: 352.,
179 scl_slope: 0.,
180 scl_inter: 0.,
181 slice_end: 0,
182 slice_code: 0,
183 xyzt_units: 0,
184 cal_max: 0.,
185 cal_min: 0.,
186 slice_duration: 0.,
187 toffset: 0.,
188 glmax: 0,
189 glmin: 0,
190
191 descrip: vec![0; 80],
192 aux_file: [0; 24],
193 qform_code: 1,
194 sform_code: 1,
195 quatern_b: 0.,
196 quatern_c: 0.,
197 quatern_d: 0.,
198 quatern_x: 0.,
199 quatern_y: 0.,
200 quatern_z: 0.,
201
202 srow_x: [1., 0., 0., 0.],
203 srow_y: [0., 1., 0., 0.],
204 srow_z: [0., 0., 1., 0.],
205
206 intent_name: [0; 16],
207
208 magic: *MAGIC_CODE_NI1,
209
210 endianness: Endianness::native(),
211 }
212 }
213}
214
215impl NiftiHeader {
216 pub fn from_file<P: AsRef<Path>>(path: P) -> Result<NiftiHeader> {
219 let gz = is_gz_file(&path);
220 let file = BufReader::new(File::open(path)?);
221 if gz {
222 NiftiHeader::from_reader(GzDecoder::new(file))
223 } else {
224 NiftiHeader::from_reader(file)
225 }
226 }
227
228 pub fn from_reader<S>(input: S) -> Result<NiftiHeader>
232 where
233 S: Read,
234 {
235 parse_header_1(input)
236 }
237
238 pub fn fix(&mut self) {
243 if !self.is_pixdim_0_valid() {
244 self.pixdim[0] = 1.0;
245 }
246 }
247
248 pub fn dim(&self) -> Result<&[u16]> {
257 validate_dim(&self.dim)
258 }
259
260 pub fn dimensionality(&self) -> Result<usize> {
268 validate_dimensionality(&self.dim)
269 }
270
271 pub fn data_type(&self) -> Result<NiftiType> {
273 FromPrimitive::from_i16(self.datatype)
274 .ok_or(NiftiError::InvalidCode("datatype", self.datatype))
275 }
276
277 pub fn xyzt_to_space(&self) -> Result<Unit> {
279 let space_code = self.xyzt_units & 0o0007;
280 FromPrimitive::from_u8(space_code).ok_or(NiftiError::InvalidCode(
281 "xyzt units (space)",
282 space_code as i16,
283 ))
284 }
285
286 pub fn xyzt_to_time(&self) -> Result<Unit> {
288 let time_code = self.xyzt_units & 0o0070;
289 FromPrimitive::from_u8(time_code).ok_or(NiftiError::InvalidCode(
290 "xyzt units (time)",
291 time_code as i16,
292 ))
293 }
294
295 pub fn xyzt_units(&self) -> Result<(Unit, Unit)> {
297 Ok((self.xyzt_to_space()?, self.xyzt_to_time()?))
298 }
299
300 pub fn slice_order(&self) -> Result<SliceOrder> {
302 FromPrimitive::from_u8(self.slice_code).ok_or(NiftiError::InvalidCode(
303 "slice order",
304 self.slice_code as i16,
305 ))
306 }
307
308 pub fn intent(&self) -> Result<Intent> {
310 FromPrimitive::from_i16(self.intent_code)
311 .ok_or(NiftiError::InvalidCode("intent", self.intent_code))
312 }
313
314 pub fn qform(&self) -> Result<XForm> {
316 FromPrimitive::from_i16(self.qform_code)
317 .ok_or(NiftiError::InvalidCode("qform", self.qform_code))
318 }
319
320 pub fn sform(&self) -> Result<XForm> {
322 FromPrimitive::from_i16(self.sform_code)
323 .ok_or(NiftiError::InvalidCode("sform", self.sform_code))
324 }
325
326 pub fn validate_description(&mut self) -> Result<()> {
330 let len = self.descrip.len();
331 if len > 80 {
332 Err(NiftiError::IncorrectDescriptionLength(len))
333 } else {
334 if len < 80 {
335 self.descrip.extend((len..80).map(|_| 0));
336 }
337 Ok(())
338 }
339 }
340
341 pub fn set_description<D>(&mut self, description: D) -> Result<()>
343 where
344 D: Into<Vec<u8>>,
345 D: Deref<Target = [u8]>,
346 {
347 let len = description.len();
348 match len.cmp(&80) {
349 std::cmp::Ordering::Less => {
350 let mut descrip = vec![0; 80];
351 descrip[..len].copy_from_slice(&description);
352 self.descrip = descrip;
353 Ok(())
354 }
355 std::cmp::Ordering::Equal => {
356 self.descrip = description.into();
357 Ok(())
358 }
359 _ => Err(NiftiError::IncorrectDescriptionLength(len)),
360 }
361 }
362
363 pub fn set_description_str<T>(&mut self, description: T) -> Result<()>
365 where
366 T: Into<String>,
367 {
368 self.set_description(description.into().as_bytes())
369 }
370
371 #[inline]
373 fn is_pixdim_0_valid(&self) -> bool {
374 (self.pixdim[0].abs() - 1.).abs() < 1e-11
375 }
376}
377
378#[cfg(feature = "nalgebra_affine")]
379impl NiftiHeader {
380 pub fn affine<T>(&self) -> Matrix4<T>
389 where
390 T: RealField,
391 f32: SubsetOf<T>,
392 {
393 if self.sform_code != 0 {
394 self.sform_affine::<T>()
395 } else if self.qform_code != 0 {
396 self.qform_affine::<T>()
397 } else {
398 self.base_affine::<T>()
399 }
400 }
401
402 pub fn sform_affine<T>(&self) -> Matrix4<T>
404 where
405 T: RealField,
406 f32: SubsetOf<T>,
407 {
408 #[rustfmt::skip]
409 let affine = Matrix4::new(
410 self.srow_x[0], self.srow_x[1], self.srow_x[2], self.srow_x[3],
411 self.srow_y[0], self.srow_y[1], self.srow_y[2], self.srow_y[3],
412 self.srow_z[0], self.srow_z[1], self.srow_z[2], self.srow_z[3],
413 0.0, 0.0, 0.0, 1.0,
414 );
415 nalgebra::convert(affine)
416 }
417
418 pub fn qform_affine<T>(&self) -> Matrix4<T>
420 where
421 T: RealField,
422 {
423 if self.pixdim[1] < 0.0 || self.pixdim[2] < 0.0 || self.pixdim[3] < 0.0 {
424 panic!("All spacings (pixdim) should be positive");
425 }
426 if !self.is_pixdim_0_valid() {
427 panic!("qfac (pixdim[0]) should be 1 or -1");
428 }
429
430 let quaternion = self.qform_quaternion();
431 let r = quaternion_to_affine(quaternion);
432 let s = Matrix3::from_diagonal(&Vector3::new(
433 self.pixdim[1] as f64,
434 self.pixdim[2] as f64,
435 self.pixdim[3] as f64 * self.pixdim[0] as f64,
436 ));
437 let m = r * s;
438 #[rustfmt::skip]
439 let affine = Matrix4::new(
440 m[0], m[3], m[6], self.quatern_x as f64,
441 m[1], m[4], m[7], self.quatern_y as f64,
442 m[2], m[5], m[8], self.quatern_z as f64,
443 0.0, 0.0, 0.0, 1.0,
444 );
445 nalgebra::convert(affine)
446 }
447
448 fn base_affine<T>(&self) -> Matrix4<T>
452 where
453 T: RealField,
454 {
455 let d = self.dim[0] as usize;
456 let affine = shape_zoom_affine(&self.dim[1..d + 1], &self.pixdim[1..d + 1]);
457 nalgebra::convert(affine)
458 }
459
460 fn qform_quaternion(&self) -> Quaternion<f64> {
464 let xyz = Vector3::new(
465 self.quatern_b as f64,
466 self.quatern_c as f64,
467 self.quatern_d as f64,
468 );
469 fill_positive(xyz)
470 }
471
472 pub fn set_affine<T>(&mut self, affine: &Matrix4<T>)
480 where
481 T: RealField,
482 T: SubsetOf<f64>,
483 T: ToPrimitive,
484 {
485 self.set_sform(affine, XForm::AlignedAnat);
487
488 self.set_qform(affine, XForm::Unknown);
490 }
491
492 pub fn set_sform<T>(&mut self, affine: &Matrix4<T>, code: XForm)
494 where
495 T: RealField,
496 T: ToPrimitive,
497 {
498 self.sform_code = code as i16;
499 self.srow_x[0] = affine[0].to_f32().unwrap();
500 self.srow_x[1] = affine[4].to_f32().unwrap();
501 self.srow_x[2] = affine[8].to_f32().unwrap();
502 self.srow_x[3] = affine[12].to_f32().unwrap();
503 self.srow_y[0] = affine[1].to_f32().unwrap();
504 self.srow_y[1] = affine[5].to_f32().unwrap();
505 self.srow_y[2] = affine[9].to_f32().unwrap();
506 self.srow_y[3] = affine[13].to_f32().unwrap();
507 self.srow_z[0] = affine[2].to_f32().unwrap();
508 self.srow_z[1] = affine[6].to_f32().unwrap();
509 self.srow_z[2] = affine[10].to_f32().unwrap();
510 self.srow_z[3] = affine[14].to_f32().unwrap();
511 }
512
513 pub fn set_qform<T>(&mut self, affine4: &Matrix4<T>, code: XForm)
520 where
521 T: RealField,
522 T: SubsetOf<f64>,
523 T: ToPrimitive,
524 {
525 let affine4: Matrix4<f64> = nalgebra::convert(affine4.clone());
526 let (affine, translation) = affine_and_translation(&affine4);
527 let aff2 = affine.component_mul(&affine);
528 let spacing = (
529 (aff2[0] + aff2[1] + aff2[2]).sqrt(),
530 (aff2[3] + aff2[4] + aff2[5]).sqrt(),
531 (aff2[6] + aff2[7] + aff2[8]).sqrt(),
532 );
533 #[rustfmt::skip]
534 let mut r = Matrix3::new(
535 affine[0] / spacing.0, affine[3] / spacing.1, affine[6] / spacing.2,
536 affine[1] / spacing.0, affine[4] / spacing.1, affine[7] / spacing.2,
537 affine[2] / spacing.0, affine[5] / spacing.1, affine[8] / spacing.2,
538 );
539
540 let qfac = if r.determinant() > 0.0 {
542 1.0
543 } else {
544 r[6] *= -1.0;
545 r[7] *= -1.0;
546 r[8] *= -1.0;
547 -1.0
548 };
549
550 let svd = r.svd(true, true);
555 let pr = svd.u.unwrap() * svd.v_t.unwrap();
556 let quaternion = affine_to_quaternion(&pr);
557
558 self.qform_code = code as i16;
559 self.pixdim[0] = qfac;
560 self.pixdim[1] = spacing.0 as f32;
561 self.pixdim[2] = spacing.1 as f32;
562 self.pixdim[3] = spacing.2 as f32;
563 self.quatern_b = quaternion[1] as f32;
564 self.quatern_c = quaternion[2] as f32;
565 self.quatern_d = quaternion[3] as f32;
566 self.quatern_x = translation[0] as f32;
567 self.quatern_y = translation[1] as f32;
568 self.quatern_z = translation[2] as f32;
569 }
570}
571
572fn parse_header_1<S>(input: S) -> Result<NiftiHeader>
573where
574 S: Read,
575{
576 let mut h = NiftiHeader::default();
577
578 let mut input = ByteOrdered::native(input);
580
581 h.sizeof_hdr = input.read_i32()?;
582 input.read_exact(&mut h.data_type)?;
583 input.read_exact(&mut h.db_name)?;
584 h.extents = input.read_i32()?;
585 h.session_error = input.read_i16()?;
586 h.regular = input.read_u8()?;
587 h.dim_info = input.read_u8()?;
588 h.dim[0] = input.read_u16()?;
589
590 if h.dim[0] > 7 {
591 h.endianness = Endianness::native().to_opposite();
592
593 h.sizeof_hdr = h.sizeof_hdr.swap_bytes();
595 h.extents = h.extents.swap_bytes();
596 h.session_error = h.session_error.swap_bytes();
597 h.dim[0] = h.dim[0].swap_bytes();
598 parse_header_2(h, input.into_opposite())
599 } else {
600 h.endianness = Endianness::native();
602 parse_header_2(h, input)
603 }
604}
605
606fn parse_header_2<S, E>(mut h: NiftiHeader, mut input: ByteOrdered<S, E>) -> Result<NiftiHeader>
608where
609 S: Read,
610 E: Endian,
611{
612 for v in &mut h.dim[1..] {
613 *v = input.read_u16()?;
614 }
615 h.intent_p1 = input.read_f32()?;
616 h.intent_p2 = input.read_f32()?;
617 h.intent_p3 = input.read_f32()?;
618 h.intent_code = input.read_i16()?;
619 h.datatype = input.read_i16()?;
620 h.bitpix = input.read_i16()?;
621 h.slice_start = input.read_i16()?;
622 for v in &mut h.pixdim {
623 *v = input.read_f32()?;
624 }
625 h.vox_offset = input.read_f32()?;
626 h.scl_slope = input.read_f32()?;
627 h.scl_inter = input.read_f32()?;
628 h.slice_end = input.read_i16()?;
629 h.slice_code = input.read_u8()?;
630 h.xyzt_units = input.read_u8()?;
631 h.cal_max = input.read_f32()?;
632 h.cal_min = input.read_f32()?;
633 h.slice_duration = input.read_f32()?;
634 h.toffset = input.read_f32()?;
635 h.glmax = input.read_i32()?;
636 h.glmin = input.read_i32()?;
637
638 input.read_exact(h.descrip.as_mut_slice())?;
640 input.read_exact(&mut h.aux_file)?;
641 h.qform_code = input.read_i16()?;
642 h.sform_code = input.read_i16()?;
643 h.quatern_b = input.read_f32()?;
644 h.quatern_c = input.read_f32()?;
645 h.quatern_d = input.read_f32()?;
646 h.quatern_x = input.read_f32()?;
647 h.quatern_y = input.read_f32()?;
648 h.quatern_z = input.read_f32()?;
649 for v in &mut h.srow_x {
650 *v = input.read_f32()?;
651 }
652 for v in &mut h.srow_y {
653 *v = input.read_f32()?;
654 }
655 for v in &mut h.srow_z {
656 *v = input.read_f32()?;
657 }
658 input.read_exact(&mut h.intent_name)?;
659 input.read_exact(&mut h.magic)?;
660
661 debug_assert_eq!(h.descrip.len(), 80);
662
663 if &h.magic != MAGIC_CODE_NI1 && &h.magic != MAGIC_CODE_NIP1 {
664 Err(NiftiError::InvalidFormat)
665 } else {
666 Ok(h)
667 }
668}