From c5d5b7bf74a99a3a491932cd8ccbed48ee7d76d1 Mon Sep 17 00:00:00 2001 From: J-Donald Tournier Date: Wed, 17 Apr 2019 00:13:15 +0100 Subject: [PATCH] NIfTI: simplify NIfTI 1 & 2 handling into single template implementation --- core/file/mgh.cpp | 1 - core/file/nifti1_utils.cpp | 474 --------------------------------- core/file/nifti1_utils.h | 42 --- core/file/nifti2_utils.cpp | 422 ----------------------------- core/file/nifti2_utils.h | 43 --- core/file/nifti_utils.cpp | 524 ++++++++++++++++++++++++++++++++++++- core/file/nifti_utils.h | 14 + core/formats/analyse.cpp | 5 +- core/formats/nifti1.cpp | 11 +- core/formats/nifti1_gz.cpp | 23 +- core/formats/nifti2.cpp | 7 +- core/formats/nifti2_gz.cpp | 23 +- 12 files changed, 572 insertions(+), 1017 deletions(-) delete mode 100644 core/file/nifti1_utils.cpp delete mode 100644 core/file/nifti1_utils.h delete mode 100644 core/file/nifti2_utils.cpp delete mode 100644 core/file/nifti2_utils.h diff --git a/core/file/mgh.cpp b/core/file/mgh.cpp index ad0c661ca4..a9db7a76dc 100644 --- a/core/file/mgh.cpp +++ b/core/file/mgh.cpp @@ -19,7 +19,6 @@ #include "header.h" #include "file/ofstream.h" #include "file/mgh.h" -#include "file/nifti1_utils.h" namespace MR { diff --git a/core/file/nifti1_utils.cpp b/core/file/nifti1_utils.cpp deleted file mode 100644 index f411df3d7a..0000000000 --- a/core/file/nifti1_utils.cpp +++ /dev/null @@ -1,474 +0,0 @@ -/* Copyright (c) 2008-2019 the MRtrix3 contributors. - * - * This Source Code Form is subject to the terms of the Mozilla Public - * License, v. 2.0. If a copy of the MPL was not distributed with this - * file, You can obtain one at http://mozilla.org/MPL/2.0/. - * - * Covered Software is provided under this License on an "as is" - * basis, without warranty of any kind, either expressed, implied, or - * statutory, including, without limitation, warranties that the - * Covered Software is free of defects, merchantable, fit for a - * particular purpose or non-infringing. - * See the Mozilla Public License v. 2.0 for more details. - * - * For more details, see http://www.mrtrix.org/. - */ - -#include "header.h" -#include "raw.h" -#include "file/config.h" -#include "file/json_utils.h" -#include "file/nifti_utils.h" -#include "file/nifti1_utils.h" - -namespace MR -{ - namespace File - { - namespace NIfTI1 - { - - - - size_t read (Header& H, const nifti_1_header& NH) - { - bool is_BE = false; - if (Raw::fetch_ (&NH.sizeof_hdr, is_BE) != header_size) { - is_BE = true; - if (Raw::fetch_ (&NH.sizeof_hdr, is_BE) != header_size) - throw Exception ("image \"" + H.name() + "\" is not in NIfTI-1.1 format (sizeof_hdr != " + str(header_size) + ")"); - } - - bool is_nifti = true; - if (memcmp (NH.magic, "n+1\0", 4) && memcmp (NH.magic, "ni1\0", 4)) { - is_nifti = false; - DEBUG ("assuming image \"" + H.name() + "\" is in AnalyseAVW format."); - } - - char db_name[19]; - strncpy (db_name, NH.db_name, 18); - if (db_name[0]) { - db_name[18] = '\0'; - add_line (H.keyval()["comments"], db_name); - } - - // data set dimensions: - int ndim = Raw::fetch_ (&NH.dim, is_BE); - if (ndim < 1) - throw Exception ("too few dimensions specified in NIfTI-1.1 image \"" + H.name() + "\""); - if (ndim > 7) - throw Exception ("too many dimensions specified in NIfTI-1.1 image \"" + H.name() + "\""); - H.ndim() = ndim; - - - for (int i = 0; i < ndim; i++) { - H.size(i) = Raw::fetch_ (&NH.dim[i+1], is_BE); - if (H.size (i) < 0) { - INFO ("dimension along axis " + str (i) + " specified as negative in NIfTI-1.1 image \"" + H.name() + "\" - taking absolute value"); - H.size(i) = abs (H.size (i)); - } - if (!H.size (i)) - H.size(i) = 1; - H.stride(i) = i+1; - } - - // data type: - DataType dtype; - switch (Raw::fetch_ (&NH.datatype, is_BE)) { - case DT_BINARY: - dtype = DataType::Bit; - break; - case DT_INT8: - dtype = DataType::Int8; - break; - case DT_UINT8: - dtype = DataType::UInt8; - break; - case DT_INT16: - dtype = DataType::Int16; - break; - case DT_UINT16: - dtype = DataType::UInt16; - break; - case DT_INT32: - dtype = DataType::Int32; - break; - case DT_UINT32: - dtype = DataType::UInt32; - break; - case DT_INT64: - dtype = DataType::Int64; - break; - case DT_UINT64: - dtype = DataType::UInt64; - break; - case DT_FLOAT32: - dtype = DataType::Float32; - break; - case DT_FLOAT64: - dtype = DataType::Float64; - break; - case DT_COMPLEX64: - dtype = DataType::CFloat32; - break; - case DT_COMPLEX128: - dtype = DataType::CFloat64; - break; - default: - throw Exception ("unknown data type for NIfTI-1.1 image \"" + H.name() + "\""); - } - - if (! (dtype.is (DataType::Bit) || dtype.is (DataType::UInt8) || dtype.is (DataType::Int8))) { - if (is_BE) - dtype.set_flag (DataType::BigEndian); - else - dtype.set_flag (DataType::LittleEndian); - } - - if (Raw::fetch_ (&NH.bitpix, is_BE) != (int16_t) dtype.bits()) - WARN ("bitpix field does not match data type in NIfTI-1.1 image \"" + H.name() + "\" - ignored"); - - H.datatype() = dtype; - - // voxel sizes: - for (int i = 0; i < ndim; i++) { - H.spacing(i) = Raw::fetch_ (&NH.pixdim[i+1], is_BE); - if (H.spacing (i) < 0.0) { - INFO ("voxel size along axis " + str (i) + " specified as negative in NIfTI-1.1 image \"" + H.name() + "\" - taking absolute value"); - H.spacing(i) = abs (H.spacing (i)); - } - } - - - // offset and scale: - H.intensity_scale() = Raw::fetch_ (&NH.scl_slope, is_BE); - if (std::isfinite (H.intensity_scale()) && H.intensity_scale() != 0.0) { - H.intensity_offset() = Raw::fetch_ (&NH.scl_inter, is_BE); - if (!std::isfinite (H.intensity_offset())) - H.intensity_offset() = 0.0; - } - else - H.reset_intensity_scaling(); - - size_t data_offset = (size_t) Raw::fetch_ (&NH.vox_offset, is_BE); - - char descrip[81]; - strncpy (descrip, NH.descrip, 80); - if (descrip[0]) { - descrip[80] = '\0'; - if (strncmp (descrip, "MRtrix version: ", 16) == 0) - H.keyval()["mrtrix_version"] = descrip+16; - else - add_line (H.keyval()["comments"], descrip); - } - - if (is_nifti) { - bool sform_code = Raw::fetch_ (&NH.sform_code, is_BE); - if (sform_code) { - auto& M (H.transform().matrix()); - - M(0,0) = Raw::fetch_ (&NH.srow_x[0], is_BE); - M(0,1) = Raw::fetch_ (&NH.srow_x[1], is_BE); - M(0,2) = Raw::fetch_ (&NH.srow_x[2], is_BE); - M(0,3) = Raw::fetch_ (&NH.srow_x[3], is_BE); - - M(1,0) = Raw::fetch_ (&NH.srow_y[0], is_BE); - M(1,1) = Raw::fetch_ (&NH.srow_y[1], is_BE); - M(1,2) = Raw::fetch_ (&NH.srow_y[2], is_BE); - M(1,3) = Raw::fetch_ (&NH.srow_y[3], is_BE); - - M(2,0) = Raw::fetch_ (&NH.srow_z[0], is_BE); - M(2,1) = Raw::fetch_ (&NH.srow_z[1], is_BE); - M(2,2) = Raw::fetch_ (&NH.srow_z[2], is_BE); - M(2,3) = Raw::fetch_ (&NH.srow_z[3], is_BE); - - // check voxel sizes: - for (size_t axis = 0; axis != 3; ++axis) { - if (size_t(ndim) > axis) - if (abs(H.spacing(axis) - M.col(axis).head<3>().norm()) > 1e-4) { - WARN ("voxel spacings inconsistent between NIFTI s-form and header field pixdim"); - break; - } - } - - // normalize each transform axis: - for (size_t axis = 0; axis != 3; ++axis) { - if (size_t(ndim) > axis) - M.col(axis).normalize(); - } - - } - - if (Raw::fetch_ (&NH.qform_code, is_BE)) { - transform_type M_qform; - - Eigen::Quaterniond Q (0.0, Raw::fetch_ (&NH.quatern_b, is_BE), Raw::fetch_ (&NH.quatern_c, is_BE), Raw::fetch_ (&NH.quatern_d, is_BE)); - const double w = 1.0 - Q.squaredNorm(); - if (w < 1.0e-7) - Q.normalize(); - else - Q.w() = std::sqrt (w); - M_qform.matrix().topLeftCorner<3,3>() = Q.matrix(); - - M_qform.translation()[0] = Raw::fetch_ (&NH.qoffset_x, is_BE); - M_qform.translation()[1] = Raw::fetch_ (&NH.qoffset_y, is_BE); - M_qform.translation()[2] = Raw::fetch_ (&NH.qoffset_z, is_BE); - - // qfac: - float qfac = Raw::fetch_ (&NH.pixdim[0], is_BE) >= 0.0 ? 1.0 : -1.0; - if (qfac < 0.0) - M_qform.matrix().col(2) *= qfac; - - if (sform_code) { - Header header2 (H); - header2.transform() = M_qform; - if (!voxel_grids_match_in_scanner_space (H, header2, 0.1)) { - //CONF option: NIfTIUseSform - //CONF default: 1 (true) - //CONF A boolean value to control whether, in cases where both - //CONF the sform and qform transformations are defined in an - //CONF input NIfTI image, but those transformations differ, the - //CONF sform transformation should be used in preference to the - //CONF qform matrix. - const bool use_sform = File::Config::get_bool ("NIfTIUseSform", true); - WARN ("qform and sform are inconsistent in NIfTI image \"" + H.name() + "\" - using " + (use_sform ? "sform" : "qform")); - if (!use_sform) - H.transform() = M_qform; - } - } - else - H.transform() = M_qform; - } - - //CONF option: NIfTIAutoLoadJSON - //CONF default: 0 (false) - //CONF A boolean value to indicate whether, when opening NIfTI images, - //CONF any corresponding JSON file should be automatically loaded. - if (File::Config::get_bool ("NIfTIAutoLoadJSON", false)) { - std::string json_path = H.name(); - if (Path::has_suffix (json_path, ".nii.gz")) - json_path = json_path.substr (0, json_path.size()-7); - else if (Path::has_suffix (json_path, ".nii")) - json_path = json_path.substr (0, json_path.size()-4); - else - assert (0); - json_path += ".json"; - if (Path::exists (json_path)) - File::JSON::load (H, json_path); - } - - } - else { - H.transform()(0,0) = std::numeric_limits::quiet_NaN(); - //CONF option: AnalyseLeftToRight - //CONF default: 0 (false) - //CONF A boolean value to indicate whether images in Analyse format - //CONF should be assumed to be in LAS orientation (default) or RAS - //CONF (when this is option is turned on). - if (!File::Config::get_bool ("AnalyseLeftToRight", false)) - H.stride(0) = -H.stride (0); - if (!File::NIfTI::right_left_warning_issued) { - INFO ("assuming Analyse images are encoded " + std::string (H.stride (0) > 0 ? "left to right" : "right to left")); - File::NIfTI::right_left_warning_issued = true; - } - } - - return data_offset; - } - - - - - - - - void write (nifti_1_header& NH, const Header& H, const bool single_file) - { - if (H.ndim() > 7) - throw Exception ("NIfTI-1.1 format cannot support more than 7 dimensions for image \"" + H.name() + "\""); - - bool is_BE = H.datatype().is_big_endian(); - - vector axes; - auto M = File::NIfTI::adjust_transform (H, axes); - - memset (&NH, 0, sizeof (NH)); - - // magic number: - Raw::store (header_size, &NH.sizeof_hdr, is_BE); - - const auto hit = H.keyval().find("comments"); - auto comments = split_lines (hit == H.keyval().end() ? std::string() : hit->second); -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wpragmas" -#pragma GCC diagnostic ignored "-Wunknown-warning-option" -#pragma GCC diagnostic ignored "-Wstringop-truncation" - strncpy ( (char*) &NH.db_name, comments.size() ? comments[0].c_str() : "untitled\0\0\0\0\0\0\0\0\0\0\0", 18); -#pragma GCC diagnostic pop - Raw::store (16384, &NH.extents, is_BE); - NH.regular = 'r'; - NH.dim_info = 0; - - // data set dimensions: - Raw::store (H.ndim(), &NH.dim[0], is_BE); - { - size_t i = 0; - for (; i < 3; i++) - Raw::store (H.size (axes[i]), &NH.dim[i+1], is_BE); - for (; i < H.ndim(); i++) - Raw::store (H.size (i), &NH.dim[i+1], is_BE); - - // pad out the other dimensions with 1, fix for fslview - ++i; - for (; i < 8; i++) - Raw::store (1, &NH.dim[i], is_BE); - } - - - // data type: - int16_t dt = 0; - switch (H.datatype()()) { - case DataType::Bit: - dt = DT_BINARY; - break; - case DataType::Int8: - dt = DT_INT8; - break; - case DataType::UInt8: - dt = DT_UINT8; - break; - case DataType::Int16LE: - case DataType::Int16BE: - dt = DT_INT16; - break; - case DataType::UInt16LE: - case DataType::UInt16BE: - dt = DT_UINT16; - break; - case DataType::Int32LE: - case DataType::Int32BE: - dt = DT_INT32; - break; - case DataType::UInt32LE: - case DataType::UInt32BE: - dt = DT_UINT32; - break; - case DataType::Int64LE: - case DataType::Int64BE: - dt = DT_INT64; - break; - case DataType::UInt64LE: - case DataType::UInt64BE: - dt = DT_UINT64; - break; - case DataType::Float32LE: - case DataType::Float32BE: - dt = DT_FLOAT32; - break; - case DataType::Float64LE: - case DataType::Float64BE: - dt = DT_FLOAT64; - break; - case DataType::CFloat32LE: - case DataType::CFloat32BE: - dt = DT_COMPLEX64; - break; - case DataType::CFloat64LE: - case DataType::CFloat64BE: - dt = DT_COMPLEX128; - break; - default: - throw Exception ("unknown data type for NIfTI-1.1 image \"" + H.name() + "\""); - } - - Raw::store (dt, &NH.datatype, is_BE); - Raw::store (H.datatype().bits(), &NH.bitpix, is_BE); - - Raw::store (1.0, &NH.pixdim[0], is_BE); - // voxel sizes: - for (size_t i = 0; i < 3; ++i) - Raw::store (H.spacing (axes[i]), &NH.pixdim[i+1], is_BE); - for (size_t i = 3; i < H.ndim(); ++i) - Raw::store (H.spacing (i), &NH.pixdim[i+1], is_BE); - - Raw::store (float32(header_with_ext_size), &NH.vox_offset, is_BE); - - // offset and scale: - Raw::store (H.intensity_scale(), &NH.scl_slope, is_BE); - Raw::store (H.intensity_offset(), &NH.scl_inter, is_BE); - - NH.xyzt_units = SPACE_TIME_TO_XYZT (NIFTI_UNITS_MM, NIFTI_UNITS_SEC); - - memset ((char*) &NH.descrip, 0, 80); - std::string version_string = std::string("MRtrix version: ") + App::mrtrix_version; - if (App::project_version) - version_string += std::string(", project version: ") + App::project_version; - strncpy ( (char*) &NH.descrip, version_string.c_str(), 79); - - Raw::store (NIFTI_XFORM_SCANNER_ANAT, &NH.qform_code, is_BE); - Raw::store (NIFTI_XFORM_SCANNER_ANAT, &NH.sform_code, is_BE); - - // qform: - Eigen::Matrix3d R = M.matrix().topLeftCorner<3,3>(); - if (R.determinant() < 0.0) { - R.col(2) = -R.col(2); - Raw::store (-1.0, &NH.pixdim[0], is_BE); - } - Eigen::Quaterniond Q (R); - - if (Q.w() < 0.0) - Q.vec() = -Q.vec(); - - Raw::store (Q.x(), &NH.quatern_b, is_BE); - Raw::store (Q.y(), &NH.quatern_c, is_BE); - Raw::store (Q.z(), &NH.quatern_d, is_BE); - - - // sform: - - Raw::store (M(0,3), &NH.qoffset_x, is_BE); - Raw::store (M(1,3), &NH.qoffset_y, is_BE); - Raw::store (M(2,3), &NH.qoffset_z, is_BE); - - Raw::store (H.spacing (axes[0]) * M(0,0), &NH.srow_x[0], is_BE); - Raw::store (H.spacing (axes[1]) * M(0,1), &NH.srow_x[1], is_BE); - Raw::store (H.spacing (axes[2]) * M(0,2), &NH.srow_x[2], is_BE); - Raw::store (M (0,3), &NH.srow_x[3], is_BE); - - Raw::store (H.spacing (axes[0]) * M(1,0), &NH.srow_y[0], is_BE); - Raw::store (H.spacing (axes[1]) * M(1,1), &NH.srow_y[1], is_BE); - Raw::store (H.spacing (axes[2]) * M(1,2), &NH.srow_y[2], is_BE); - Raw::store (M (1,3), &NH.srow_y[3], is_BE); - - Raw::store (H.spacing (axes[0]) * M(2,0), &NH.srow_z[0], is_BE); - Raw::store (H.spacing (axes[1]) * M(2,1), &NH.srow_z[1], is_BE); - Raw::store (H.spacing (axes[2]) * M(2,2), &NH.srow_z[2], is_BE); - Raw::store (M (2,3), &NH.srow_z[3], is_BE); - - strncpy ( (char*) &NH.magic, single_file ? "n+1\0" : "ni1\0", 4); - - //CONF option: NIfTIAutoSaveJSON - //CONF default: 0 (false) - //CONF A boolean value to indicate whether, when writing NIfTI images, - //CONF a corresponding JSON file should be automatically created in order - //CONF to save any header entries that cannot be stored in the NIfTI - //CONF header. - if (single_file && File::Config::get_bool ("NIfTIAutoSaveJSON", false)) { - std::string json_path = H.name(); - if (Path::has_suffix (json_path, ".nii.gz")) - json_path = json_path.substr (0, json_path.size()-7); - else if (Path::has_suffix (json_path, ".nii")) - json_path = json_path.substr (0, json_path.size()-4); - else - assert (0); - json_path += ".json"; - File::JSON::save (H, json_path); - } - } - - - - } - } -} - diff --git a/core/file/nifti1_utils.h b/core/file/nifti1_utils.h deleted file mode 100644 index 20d4bb8a15..0000000000 --- a/core/file/nifti1_utils.h +++ /dev/null @@ -1,42 +0,0 @@ -/* Copyright (c) 2008-2019 the MRtrix3 contributors. - * - * This Source Code Form is subject to the terms of the Mozilla Public - * License, v. 2.0. If a copy of the MPL was not distributed with this - * file, You can obtain one at http://mozilla.org/MPL/2.0/. - * - * Covered Software is provided under this License on an "as is" - * basis, without warranty of any kind, either expressed, implied, or - * statutory, including, without limitation, warranties that the - * Covered Software is free of defects, merchantable, fit for a - * particular purpose or non-infringing. - * See the Mozilla Public License v. 2.0 for more details. - * - * For more details, see http://www.mrtrix.org/. - */ - -#ifndef __file_nifti1_utils_h__ -#define __file_nifti1_utils_h__ - -#include "file/nifti1.h" - -namespace MR -{ - class Header; - - namespace File - { - namespace NIfTI1 - { - - constexpr size_t header_size = 348; - constexpr size_t header_with_ext_size = 352; - - size_t read (Header& H, const nifti_1_header& NH); - void write (nifti_1_header& NH, const Header& H, const bool single_file); - - } - } -} - -#endif - diff --git a/core/file/nifti2_utils.cpp b/core/file/nifti2_utils.cpp deleted file mode 100644 index 6a9dea201d..0000000000 --- a/core/file/nifti2_utils.cpp +++ /dev/null @@ -1,422 +0,0 @@ -/* Copyright (c) 2008-2019 the MRtrix3 contributors. - * - * This Source Code Form is subject to the terms of the Mozilla Public - * License, v. 2.0. If a copy of the MPL was not distributed with this - * file, You can obtain one at http://mozilla.org/MPL/2.0/. - * - * Covered Software is provided under this License on an "as is" - * basis, without warranty of any kind, either expressed, implied, or - * statutory, including, without limitation, warranties that the - * Covered Software is free of defects, merchantable, fit for a - * particular purpose or non-infringing. - * See the Mozilla Public License v. 2.0 for more details. - * - * For more details, see http://www.mrtrix.org/. - */ - -#include "header.h" -#include "raw.h" -#include "file/config.h" -#include "file/json_utils.h" -#include "file/nifti1.h" -#include "file/nifti_utils.h" -#include "file/nifti2_utils.h" - -namespace MR -{ - namespace File - { - namespace NIfTI2 - { - - - - size_t read (Header& H, const nifti_2_header& NH) - { - bool is_BE = false; - if (Raw::fetch_ (&NH.sizeof_hdr, is_BE) != header_size) { - is_BE = true; - if (Raw::fetch_ (&NH.sizeof_hdr, is_BE) != header_size) - throw Exception ("image \"" + H.name() + "\" is not in NIfTI-2 format (sizeof_hdr != " + str(header_size) + ")"); - } - - if (memcmp (NH.magic, "n+2\0", 4) && memcmp (NH.magic, "ni2\0", 4)) - throw Exception ("image \"" + H.name() + "\" is not in NIfTI-2 format (invalid magic signature)"); - - if (memcmp (NH.magic+4, signature_extra, 4)) - WARN ("possible file transfer corruption of file \"" + H.name() + "\" (invalid magic signature)"); - - // data type: - DataType dtype; - switch (Raw::fetch_ (&NH.datatype, is_BE)) { - case DT_BINARY: - dtype = DataType::Bit; - break; - case DT_INT8: - dtype = DataType::Int8; - break; - case DT_UINT8: - dtype = DataType::UInt8; - break; - case DT_INT16: - dtype = DataType::Int16; - break; - case DT_UINT16: - dtype = DataType::UInt16; - break; - case DT_INT32: - dtype = DataType::Int32; - break; - case DT_UINT32: - dtype = DataType::UInt32; - break; - case DT_INT64: - dtype = DataType::Int64; - break; - case DT_UINT64: - dtype = DataType::UInt64; - break; - case DT_FLOAT32: - dtype = DataType::Float32; - break; - case DT_FLOAT64: - dtype = DataType::Float64; - break; - case DT_COMPLEX64: - dtype = DataType::CFloat32; - break; - case DT_COMPLEX128: - dtype = DataType::CFloat64; - break; - default: - throw Exception ("unknown data type for NIfTI-2 image \"" + H.name() + "\""); - } - - if (! (dtype.is (DataType::Bit) || dtype.is (DataType::UInt8) || dtype.is (DataType::Int8))) { - if (is_BE) - dtype.set_flag (DataType::BigEndian); - else - dtype.set_flag (DataType::LittleEndian); - } - - H.datatype() = dtype; - - if (Raw::fetch_ (&NH.bitpix, is_BE) != (int16_t) dtype.bits()) - WARN ("bitpix field does not match data type in NIfTI-2 image \"" + H.name() + "\" - ignored"); - - // data set dimensions: - const int64_t ndim = Raw::fetch_ (&NH.dim, is_BE); - if (ndim < 1) - throw Exception ("too few dimensions specified in NIfTI-2 image \"" + H.name() + "\""); - if (ndim > 7) - throw Exception ("too many dimensions specified in NIfTI-2 image \"" + H.name() + "\""); - H.ndim() = ndim; - for (int64_t i = 0; i != ndim; i++) { - H.size(i) = Raw::fetch_ (&NH.dim[i+1], is_BE); - if (H.size (i) < 0) { - INFO ("dimension along axis " + str (i) + " specified as negative in NIfTI-2 image \"" + H.name() + "\" - taking absolute value"); - H.size(i) = abs (H.size (i)); - } - if (!H.size (i)) - H.size(i) = 1; - H.stride(i) = i+1; - } - - // voxel sizes: - for (int i = 0; i < ndim; i++) { - H.spacing(i) = Raw::fetch_ (&NH.pixdim[i+1], is_BE); - if (H.spacing (i) < 0.0) { - INFO ("voxel size along axis " + str (i) + " specified as negative in NIfTI-2 image \"" + H.name() + "\" - taking absolute value"); - H.spacing(i) = abs (H.spacing (i)); - } - } - - const int64_t data_offset = Raw::fetch_ (&NH.vox_offset, is_BE); - - // offset and scale: - H.intensity_scale() = Raw::fetch_ (&NH.scl_slope, is_BE); - if (std::isfinite (H.intensity_scale()) && H.intensity_scale() != 0.0) { - H.intensity_offset() = Raw::fetch_ (&NH.scl_inter, is_BE); - if (!std::isfinite (H.intensity_offset())) - H.intensity_offset() = 0.0; - } else { - H.reset_intensity_scaling(); - } - - char descrip[81]; - strncpy (descrip, NH.descrip, 80); - if (descrip[0]) { - descrip[80] = '\0'; - if (strncmp (descrip, "MRtrix version: ", 16) == 0) - H.keyval()["mrtrix_version"] = descrip+16; - else - add_line (H.keyval()["comments"], descrip); - } - - // Note: Unlike reading from a nifti_1_header class, here we - // don't have to worry about whether or not the file is in - // Analyse format; we can treat it as a NIfTI regardless of - // whether the hedaer & data are in the same file or not. - bool sform_code = Raw::fetch_ (&NH.sform_code, is_BE); - if (sform_code) { - auto& M (H.transform().matrix()); - - M(0,0) = Raw::fetch_ (&NH.srow_x[0], is_BE); - M(0,1) = Raw::fetch_ (&NH.srow_x[1], is_BE); - M(0,2) = Raw::fetch_ (&NH.srow_x[2], is_BE); - M(0,3) = Raw::fetch_ (&NH.srow_x[3], is_BE); - - M(1,0) = Raw::fetch_ (&NH.srow_y[0], is_BE); - M(1,1) = Raw::fetch_ (&NH.srow_y[1], is_BE); - M(1,2) = Raw::fetch_ (&NH.srow_y[2], is_BE); - M(1,3) = Raw::fetch_ (&NH.srow_y[3], is_BE); - - M(2,0) = Raw::fetch_ (&NH.srow_z[0], is_BE); - M(2,1) = Raw::fetch_ (&NH.srow_z[1], is_BE); - M(2,2) = Raw::fetch_ (&NH.srow_z[2], is_BE); - M(2,3) = Raw::fetch_ (&NH.srow_z[3], is_BE); - - // check voxel sizes: - for (size_t axis = 0; axis != 3; ++axis) { - if (size_t(ndim) > axis) - if (abs(H.spacing(axis) - M.col(axis).head<3>().norm()) > 1e-4) { - WARN ("voxel spacings inconsistent between NIFTI s-form and header field pixdim"); - break; - } - } - - // normalize each transform axis: - for (size_t axis = 0; axis != 3; ++axis) { - if (size_t(ndim) > axis) - M.col(axis).normalize(); - } - - } - - if (Raw::fetch_ (&NH.qform_code, is_BE)) { - transform_type M_qform; - - Eigen::Quaterniond Q (0.0, Raw::fetch_ (&NH.quatern_b, is_BE), Raw::fetch_ (&NH.quatern_c, is_BE), Raw::fetch_ (&NH.quatern_d, is_BE)); - const double w = 1.0 - Q.squaredNorm(); - if (w < 1.0e-15) - Q.normalize(); - else - Q.w() = std::sqrt (w); - M_qform.matrix().topLeftCorner<3,3>() = Q.matrix(); - - M_qform.translation()[0] = Raw::fetch_ (&NH.qoffset_x, is_BE); - M_qform.translation()[1] = Raw::fetch_ (&NH.qoffset_y, is_BE); - M_qform.translation()[2] = Raw::fetch_ (&NH.qoffset_z, is_BE); - - // qfac: - const float64 qfac = Raw::fetch_ (&NH.pixdim[0], is_BE) >= 0.0 ? 1.0 : -1.0; - if (qfac < 0.0) - M_qform.matrix().col(2) *= qfac; - - if (sform_code) { - Header header2 (H); - header2.transform() = M_qform; - if (!voxel_grids_match_in_scanner_space (H, header2, 0.1)) { - const bool use_sform = File::Config::get_bool ("NIfTIUseSform", true); - WARN ("qform and sform are inconsistent in NIfTI image \"" + H.name() + "\" - using " + (use_sform ? "sform" : "qform")); - if (!use_sform) - H.transform() = M_qform; - } - } - else - H.transform() = M_qform; - } - - if (File::Config::get_bool ("NIfTIAutoLoadJSON", false)) { - std::string json_path = H.name(); - if (Path::has_suffix (json_path, ".nii.gz")) - json_path = json_path.substr (0, json_path.size()-7); - else if (Path::has_suffix (json_path, ".nii")) - json_path = json_path.substr (0, json_path.size()-4); - else - assert (0); - json_path += ".json"; - if (Path::exists (json_path)) - File::JSON::load (H, json_path); - } - - return data_offset; - } - - - - - - - void write (nifti_2_header& NH, const Header& H, const bool single_file) - { - if (H.ndim() > 7) - throw Exception ("NIfTI-2 format cannot support more than 7 dimensions for image \"" + H.name() + "\""); - - bool is_BE = H.datatype().is_big_endian(); - - vector axes; - auto M = File::NIfTI::adjust_transform (H, axes); - - - memset (&NH, 0, sizeof (NH)); - - // magic number: - Raw::store (header_size, &NH.sizeof_hdr, is_BE); - - strncpy ( (char*) &NH.magic, single_file ? "n+2\0" : "ni2\0", 4); - strncpy ( (char*) &NH.magic+4, signature_extra, 4); - - // data type: - int16_t dt = 0; - switch (H.datatype()()) { - case DataType::Bit: - dt = DT_BINARY; - break; - case DataType::Int8: - dt = DT_INT8; - break; - case DataType::UInt8: - dt = DT_UINT8; - break; - case DataType::Int16LE: - case DataType::Int16BE: - dt = DT_INT16; - break; - case DataType::UInt16LE: - case DataType::UInt16BE: - dt = DT_UINT16; - break; - case DataType::Int32LE: - case DataType::Int32BE: - dt = DT_INT32; - break; - case DataType::UInt32LE: - case DataType::UInt32BE: - dt = DT_UINT32; - break; - case DataType::Int64LE: - case DataType::Int64BE: - dt = DT_INT64; - break; - case DataType::UInt64LE: - case DataType::UInt64BE: - dt = DT_UINT64; - break; - case DataType::Float32LE: - case DataType::Float32BE: - dt = DT_FLOAT32; - break; - case DataType::Float64LE: - case DataType::Float64BE: - dt = DT_FLOAT64; - break; - case DataType::CFloat32LE: - case DataType::CFloat32BE: - dt = DT_COMPLEX64; - break; - case DataType::CFloat64LE: - case DataType::CFloat64BE: - dt = DT_COMPLEX128; - break; - default: - throw Exception ("unknown data type for NIfTI-2 image \"" + H.name() + "\""); - } - Raw::store (dt, &NH.datatype, is_BE); - - Raw::store (H.datatype().bits(), &NH.bitpix, is_BE); - - // data set dimensions: - Raw::store (H.ndim(), &NH.dim[0], is_BE); - { - size_t i = 0; - for (; i < 3; i++) - Raw::store (H.size (axes[i]), &NH.dim[i+1], is_BE); - for (; i < H.ndim(); i++) - Raw::store (H.size (i), &NH.dim[i+1], is_BE); - - // pad out the other dimensions with 1, fix for fslview - ++i; - for (; i < 8; i++) - Raw::store (1, &NH.dim[i], is_BE); - } - - Raw::store (1.0, &NH.pixdim[0], is_BE); - // voxel sizes: - for (size_t i = 0; i < 3; ++i) - Raw::store (H.spacing (axes[i]), &NH.pixdim[i+1], is_BE); - for (size_t i = 3; i < H.ndim(); ++i) - Raw::store (H.spacing (i), &NH.pixdim[i+1], is_BE); - - Raw::store (int64_t(header_with_ext_size), &NH.vox_offset, is_BE); - - // offset and scale: - Raw::store (H.intensity_scale(), &NH.scl_slope, is_BE); - Raw::store (H.intensity_offset(), &NH.scl_inter, is_BE); - - std::string version_string = std::string("MRtrix version: ") + App::mrtrix_version; - if (App::project_version) - version_string += std::string(", project version: ") + App::project_version; - strncpy ( (char*) &NH.descrip, version_string.c_str(), 79); - - Raw::store (NIFTI_XFORM_SCANNER_ANAT, &NH.qform_code, is_BE); - Raw::store (NIFTI_XFORM_SCANNER_ANAT, &NH.sform_code, is_BE); - - // qform: - Eigen::Matrix3d R = M.matrix().topLeftCorner<3,3>(); - if (R.determinant() < 0.0) { - R.col(2) = -R.col(2); - Raw::store (-1.0, &NH.pixdim[0], is_BE); - } - Eigen::Quaterniond Q (R); - - if (Q.w() < 0.0) - Q.vec() = -Q.vec(); - - Raw::store (Q.x(), &NH.quatern_b, is_BE); - Raw::store (Q.y(), &NH.quatern_c, is_BE); - Raw::store (Q.z(), &NH.quatern_d, is_BE); - - Raw::store (M(0,3), &NH.qoffset_x, is_BE); - Raw::store (M(1,3), &NH.qoffset_y, is_BE); - Raw::store (M(2,3), &NH.qoffset_z, is_BE); - - // sform: - Raw::store (H.spacing (axes[0]) * M(0,0), &NH.srow_x[0], is_BE); - Raw::store (H.spacing (axes[1]) * M(0,1), &NH.srow_x[1], is_BE); - Raw::store (H.spacing (axes[2]) * M(0,2), &NH.srow_x[2], is_BE); - Raw::store (M (0,3), &NH.srow_x[3], is_BE); - - Raw::store (H.spacing (axes[0]) * M(1,0), &NH.srow_y[0], is_BE); - Raw::store (H.spacing (axes[1]) * M(1,1), &NH.srow_y[1], is_BE); - Raw::store (H.spacing (axes[2]) * M(1,2), &NH.srow_y[2], is_BE); - Raw::store (M (1,3), &NH.srow_y[3], is_BE); - - Raw::store (H.spacing (axes[0]) * M(2,0), &NH.srow_z[0], is_BE); - Raw::store (H.spacing (axes[1]) * M(2,1), &NH.srow_z[1], is_BE); - Raw::store (H.spacing (axes[2]) * M(2,2), &NH.srow_z[2], is_BE); - Raw::store (M (2,3), &NH.srow_z[3], is_BE); - - const char xyzt_units[4] { NIFTI_UNITS_MM, NIFTI_UNITS_MM, NIFTI_UNITS_MM, NIFTI_UNITS_SEC }; - const int32_t* const xyzt_units_as_int_ptr = reinterpret_cast(xyzt_units); - Raw::store (*xyzt_units_as_int_ptr, &NH.xyzt_units, is_BE); - - if (File::Config::get_bool ("NIfTI.AutoSaveJSON", false)) { - std::string json_path = H.name(); - if (Path::has_suffix (json_path, ".nii.gz")) - json_path = json_path.substr (0, json_path.size()-7); - else if (Path::has_suffix (json_path, ".nii")) - json_path = json_path.substr (0, json_path.size()-4); - else - assert (0); - json_path += ".json"; - File::JSON::save (H, json_path); - } - } - - - - } - } -} - diff --git a/core/file/nifti2_utils.h b/core/file/nifti2_utils.h deleted file mode 100644 index b2a6a8599e..0000000000 --- a/core/file/nifti2_utils.h +++ /dev/null @@ -1,43 +0,0 @@ -/* Copyright (c) 2008-2019 the MRtrix3 contributors. - * - * This Source Code Form is subject to the terms of the Mozilla Public - * License, v. 2.0. If a copy of the MPL was not distributed with this - * file, You can obtain one at http://mozilla.org/MPL/2.0/. - * - * Covered Software is provided under this License on an "as is" - * basis, without warranty of any kind, either expressed, implied, or - * statutory, including, without limitation, warranties that the - * Covered Software is free of defects, merchantable, fit for a - * particular purpose or non-infringing. - * See the Mozilla Public License v. 2.0 for more details. - * - * For more details, see http://www.mrtrix.org/. - */ - -#ifndef __file_nifti2_utils_h__ -#define __file_nifti2_utils_h__ - -#include "file/nifti2.h" - -namespace MR -{ - class Header; - - namespace File - { - namespace NIfTI2 - { - - constexpr size_t header_size = 540; - constexpr size_t header_with_ext_size = 544; - constexpr char signature_extra[4] { '\r', '\n', '\032', '\n' }; - - size_t read (Header& H, const nifti_2_header& NH); - void write (nifti_2_header& NH, const Header& H, const bool single_file); - - } - } -} - -#endif - diff --git a/core/file/nifti_utils.cpp b/core/file/nifti_utils.cpp index 40e27df43e..ec9af418c7 100644 --- a/core/file/nifti_utils.cpp +++ b/core/file/nifti_utils.cpp @@ -18,6 +18,7 @@ #include "raw.h" #include "file/config.h" #include "file/nifti_utils.h" +#include "file/json_utils.h" namespace MR { @@ -26,11 +27,530 @@ namespace MR namespace NIfTI { + namespace { + template struct Type { + using code_type = int16_t; + using float_type = float32; + using dim_type = int16_t; + using vox_offset_type = float32; + static constexpr bool is_version2 = false; + static const char* signature_extra() { return "\0\0\0\0"; } + static const char* magic1() { return "n+1\0"; } + static const char* magic2() { return "ni1\0"; } + static const std::string version() { return "NIFTI-1.1"; } + static const char* db_name (const NiftiHeaderType& NH) { return NH.db_name; } + static char* db_name (NiftiHeaderType& NH) { return NH.db_name; } + static int* extents (NiftiHeaderType& NH) { return &NH.extents; } + static char* regular (NiftiHeaderType& NH) { return &NH.regular; } + }; + + template <> struct Type { + using code_type = int32_t; + using float_type = float64; + using dim_type = int64_t; + using vox_offset_type = int64_t; + static constexpr bool is_version2 = true; + static const char* signature_extra() { return "\r\n\032\n"; } + static const char* magic1() { return "n+2\0"; } + static const char* magic2() { return "ni2\0"; } + static const std::string version() { return "NIFTI-2"; } + static const char* db_name (const nifti_2_header& NH) { return nullptr; } + static char* db_name (nifti_2_header& NH) { return nullptr; } + static int* extents (nifti_2_header& NH) { return nullptr; } + static char* regular (nifti_2_header& NH) { return nullptr; } + }; + + } + bool right_left_warning_issued = false; + + template + size_t read (Header& MH, const H& NH) + { + const std::string& version = Type::version(); + using dim_type = typename Type::dim_type; + using code_type = typename Type::code_type; + using float_type = typename Type::float_type; + + bool is_BE = false; + if (Raw::fetch_ (&NH.sizeof_hdr, is_BE) != header_size(NH)) { + is_BE = true; + if (Raw::fetch_ (&NH.sizeof_hdr, is_BE) != header_size(NH)) + throw Exception ("image \"" + MH.name() + "\" is not in " + version + + " format (sizeof_hdr != " + str(header_size(NH)) + ")"); + } + + bool is_nifti = true; + if (memcmp (NH.magic, Type::magic1(), 4) && memcmp (NH.magic, Type::magic2(), 4)) { + if (Type::is_version2) { + throw Exception ("image \"" + MH.name() + "\" is not in " + version + " format (invalid magic signature)"); + } + else { + is_nifti = false; + DEBUG ("assuming image \"" + MH.name() + "\" is in AnalyseAVW format."); + } + } + + if (Type::is_version2) { + if (memcmp (NH.magic+4, Type::signature_extra(), 4)) + WARN ("possible file transfer corruption of file \"" + MH.name() + "\" (invalid magic signature)"); + } + else{ + char db_name[19]; + strncpy (db_name, Type::db_name (NH), 18); + if (db_name[0]) { + db_name[18] = '\0'; + add_line (MH.keyval()["comments"], db_name); + } + } + + + + DataType dtype; + switch (Raw::fetch_ (&NH.datatype, is_BE)) { + case DT_BINARY: + dtype = DataType::Bit; + break; + case DT_INT8: + dtype = DataType::Int8; + break; + case DT_UINT8: + dtype = DataType::UInt8; + break; + case DT_INT16: + dtype = DataType::Int16; + break; + case DT_UINT16: + dtype = DataType::UInt16; + break; + case DT_INT32: + dtype = DataType::Int32; + break; + case DT_UINT32: + dtype = DataType::UInt32; + break; + case DT_INT64: + dtype = DataType::Int64; + break; + case DT_UINT64: + dtype = DataType::UInt64; + break; + case DT_FLOAT32: + dtype = DataType::Float32; + break; + case DT_FLOAT64: + dtype = DataType::Float64; + break; + case DT_COMPLEX64: + dtype = DataType::CFloat32; + break; + case DT_COMPLEX128: + dtype = DataType::CFloat64; + break; + default: + throw Exception ("unknown data type for " + version + " image \"" + MH.name() + "\""); + } + + if (! (dtype.is (DataType::Bit) || dtype.is (DataType::UInt8) || dtype.is (DataType::Int8))) { + if (is_BE) + dtype.set_flag (DataType::BigEndian); + else + dtype.set_flag (DataType::LittleEndian); + } + + if (Raw::fetch_ (&NH.bitpix, is_BE) != (int16_t) dtype.bits()) + WARN ("bitpix field does not match data type in " + version + " image \"" + MH.name() + "\" - ignored"); + + MH.datatype() = dtype; + + + + + const int ndim = Raw::fetch_ (&NH.dim, is_BE); + if (ndim < 1) + throw Exception ("too few dimensions specified in NIfTI image \"" + MH.name() + "\""); + if (ndim > 7) + throw Exception ("too many dimensions specified in NIfTI image \"" + MH.name() + "\""); + MH.ndim() = ndim; + for (int i = 0; i != ndim; i++) { + MH.size(i) = Raw::fetch_ (&NH.dim[i+1], is_BE); + if (MH.size (i) < 0) { + INFO ("dimension along axis " + str (i) + " specified as negative in NIfTI image \"" + MH.name() + "\" - taking absolute value"); + MH.size(i) = abs (MH.size (i)); + } + if (!MH.size (i)) + MH.size(i) = 1; + MH.stride(i) = i+1; + } + + // voxel sizes: + for (int i = 0; i < ndim; i++) { + MH.spacing(i) = Raw::fetch_ (&NH.pixdim[i+1], is_BE); + if (MH.spacing (i) < 0.0) { + INFO ("voxel size along axis " + str (i) + " specified as negative in NIfTI image \"" + MH.name() + "\" - taking absolute value"); + MH.spacing(i) = abs (MH.spacing (i)); + } + } + + + // offset and scale: + MH.intensity_scale() = Raw::fetch_ (&NH.scl_slope, is_BE); + if (std::isfinite (MH.intensity_scale()) && MH.intensity_scale() != 0.0) { + MH.intensity_offset() = Raw::fetch_ (&NH.scl_inter, is_BE); + if (!std::isfinite (MH.intensity_offset())) + MH.intensity_offset() = 0.0; + } else { + MH.reset_intensity_scaling(); + } + + + + const int64_t data_offset = Raw::fetch_ (&NH.vox_offset, is_BE); + + + + char descrip[81]; + strncpy (descrip, NH.descrip, 80); + if (descrip[0]) { + descrip[80] = '\0'; + if (strncmp (descrip, "MRtrix version: ", 16) == 0) + MH.keyval()["mrtrix_version"] = descrip+16; + else + add_line (MH.keyval()["comments"], descrip); + } + + if (is_nifti) { + bool sform_code = Raw::fetch_ (&NH.sform_code, is_BE); + if (sform_code) { + auto& M (MH.transform().matrix()); + + M(0,0) = Raw::fetch_ (&NH.srow_x[0], is_BE); + M(0,1) = Raw::fetch_ (&NH.srow_x[1], is_BE); + M(0,2) = Raw::fetch_ (&NH.srow_x[2], is_BE); + M(0,3) = Raw::fetch_ (&NH.srow_x[3], is_BE); + + M(1,0) = Raw::fetch_ (&NH.srow_y[0], is_BE); + M(1,1) = Raw::fetch_ (&NH.srow_y[1], is_BE); + M(1,2) = Raw::fetch_ (&NH.srow_y[2], is_BE); + M(1,3) = Raw::fetch_ (&NH.srow_y[3], is_BE); + + M(2,0) = Raw::fetch_ (&NH.srow_z[0], is_BE); + M(2,1) = Raw::fetch_ (&NH.srow_z[1], is_BE); + M(2,2) = Raw::fetch_ (&NH.srow_z[2], is_BE); + M(2,3) = Raw::fetch_ (&NH.srow_z[3], is_BE); + + // check voxel sizes: + for (size_t axis = 0; axis != 3; ++axis) { + if (size_t(MH.ndim()) > axis) + if (abs(MH.spacing(axis) - M.col(axis).head<3>().norm()) > 1e-4) { + WARN ("voxel spacings inconsistent between NIFTI s-form and header field pixdim"); + break; + } + } + + // normalize each transform axis: + for (size_t axis = 0; axis != 3; ++axis) { + if (size_t(MH.ndim()) > axis) + M.col(axis).normalize(); + } + + } + + if (Raw::fetch_ (&NH.qform_code, is_BE)) { + transform_type M_qform; + + Eigen::Quaterniond Q (0.0, + Raw::fetch_ (&NH.quatern_b, is_BE), + Raw::fetch_ (&NH.quatern_c, is_BE), + Raw::fetch_ (&NH.quatern_d, is_BE)); + const double w = 1.0 - Q.squaredNorm(); + if (w < 1.0e-7) + Q.normalize(); + else + Q.w() = std::sqrt (w); + M_qform.matrix().topLeftCorner<3,3>() = Q.matrix(); + + M_qform.translation()[0] = Raw::fetch_ (&NH.qoffset_x, is_BE); + M_qform.translation()[1] = Raw::fetch_ (&NH.qoffset_y, is_BE); + M_qform.translation()[2] = Raw::fetch_ (&NH.qoffset_z, is_BE); + + // qfac: + float qfac = Raw::fetch_ (&NH.pixdim[0], is_BE) >= 0.0 ? 1.0 : -1.0; + if (qfac < 0.0) + M_qform.matrix().col(2) *= qfac; + + if (sform_code) { + Header header2 (MH); + header2.transform() = M_qform; + if (!voxel_grids_match_in_scanner_space (MH, header2, 0.1)) { + //CONF option: NIfTIUseSform + //CONF default: 1 (true) + //CONF A boolean value to control whether, in cases where both + //CONF the sform and qform transformations are defined in an + //CONF input NIfTI image, but those transformations differ, the + //CONF sform transformation should be used in preference to the + //CONF qform matrix. + const bool use_sform = File::Config::get_bool ("NIfTIUseSform", true); + WARN ("qform and sform are inconsistent in NIfTI image \"" + MH.name() + "\" - using " + (use_sform ? "sform" : "qform")); + if (!use_sform) + MH.transform() = M_qform; + } + } + else + MH.transform() = M_qform; + } + + + + //CONF option: NIfTIAutoLoadJSON + //CONF default: 0 (false) + //CONF A boolean value to indicate whether, when opening NIfTI images, + //CONF any corresponding JSON file should be automatically loaded. + if (File::Config::get_bool ("NIfTIAutoLoadJSON", false)) { + std::string json_path = MH.name(); + if (Path::has_suffix (json_path, ".nii.gz")) + json_path = json_path.substr (0, json_path.size()-7); + else if (Path::has_suffix (json_path, ".nii")) + json_path = json_path.substr (0, json_path.size()-4); + else + assert (0); + json_path += ".json"; + if (Path::exists (json_path)) + File::JSON::load (MH, json_path); + } + } + else { + MH.transform()(0,0) = std::numeric_limits::quiet_NaN(); + //CONF option: AnalyseLeftToRight + //CONF default: 0 (false) + //CONF A boolean value to indicate whether images in Analyse format + //CONF should be assumed to be in LAS orientation (default) or RAS + //CONF (when this is option is turned on). + if (!File::Config::get_bool ("AnalyseLeftToRight", false)) + MH.stride(0) = -MH.stride (0); + if (!File::NIfTI::right_left_warning_issued) { + INFO ("assuming Analyse images are encoded " + std::string (MH.stride (0) > 0 ? "left to right" : "right to left")); + File::NIfTI::right_left_warning_issued = true; + } + } + + return data_offset; + } + + + + + + + template + void write (H& NH, const Header& MH, const bool single_file) + { + const std::string& version = Type::version(); + using dim_type = typename Type::dim_type; + using vox_offset_type = typename Type::vox_offset_type; + using code_type = typename Type::code_type; + using float_type = typename Type::float_type; + + if (MH.ndim() > 7) + throw Exception (version + " format cannot support more than 7 dimensions for image \"" + MH.name() + "\""); + + bool is_BE = MH.datatype().is_big_endian(); + + vector axes; + auto M = File::NIfTI::adjust_transform (MH, axes); + + + memset (&NH, 0, sizeof (NH)); + + // magic number: + Raw::store (header_size(NH), &NH.sizeof_hdr, is_BE); + + strncpy ( (char*) &NH.magic, single_file ? Type::magic1() : Type::magic2(), 4); + if (Type::is_version2) + strncpy ( (char*) &NH.magic+4, Type::signature_extra(), 4); + + if (!Type::is_version2) { + const auto hit = MH.keyval().find("comments"); + auto comments = split_lines (hit == MH.keyval().end() ? std::string() : hit->second); + strncpy ( Type::db_name (NH), comments.size() ? comments[0].c_str() : "untitled\0\0\0\0\0\0\0\0\0\0\0", 18); + Raw::store (16384, Type::extents (NH), is_BE); + *Type::regular(NH) = 'r'; + } + NH.dim_info = 0; + + // data type: + int16_t dt = 0; + switch (MH.datatype()()) { + case DataType::Bit: + dt = DT_BINARY; + break; + case DataType::Int8: + dt = DT_INT8; + break; + case DataType::UInt8: + dt = DT_UINT8; + break; + case DataType::Int16LE: + case DataType::Int16BE: + dt = DT_INT16; + break; + case DataType::UInt16LE: + case DataType::UInt16BE: + dt = DT_UINT16; + break; + case DataType::Int32LE: + case DataType::Int32BE: + dt = DT_INT32; + break; + case DataType::UInt32LE: + case DataType::UInt32BE: + dt = DT_UINT32; + break; + case DataType::Int64LE: + case DataType::Int64BE: + dt = DT_INT64; + break; + case DataType::UInt64LE: + case DataType::UInt64BE: + dt = DT_UINT64; + break; + case DataType::Float32LE: + case DataType::Float32BE: + dt = DT_FLOAT32; + break; + case DataType::Float64LE: + case DataType::Float64BE: + dt = DT_FLOAT64; + break; + case DataType::CFloat32LE: + case DataType::CFloat32BE: + dt = DT_COMPLEX64; + break; + case DataType::CFloat64LE: + case DataType::CFloat64BE: + dt = DT_COMPLEX128; + break; + default: + throw Exception ("unknown data type for " + version + " image \"" + MH.name() + "\""); + } + Raw::store (dt, &NH.datatype, is_BE); + + Raw::store (MH.datatype().bits(), &NH.bitpix, is_BE); + + // data set dimensions: + Raw::store (MH.ndim(), &NH.dim[0], is_BE); + { + size_t i = 0; + for (; i < 3; i++) + Raw::store (MH.size (axes[i]), &NH.dim[i+1], is_BE); + for (; i < MH.ndim(); i++) + Raw::store (MH.size (i), &NH.dim[i+1], is_BE); + + // pad out the other dimensions with 1, fix for fslview + ++i; + for (; i < 8; i++) + Raw::store (1, &NH.dim[i], is_BE); + } + + Raw::store (1.0, &NH.pixdim[0], is_BE); + // voxel sizes: + for (size_t i = 0; i < 3; ++i) + Raw::store (MH.spacing (axes[i]), &NH.pixdim[i+1], is_BE); + for (size_t i = 3; i < MH.ndim(); ++i) + Raw::store (MH.spacing (i), &NH.pixdim[i+1], is_BE); + + Raw::store (header_size(NH)+4, &NH.vox_offset, is_BE); + + // offset and scale: + Raw::store (MH.intensity_scale(), &NH.scl_slope, is_BE); + Raw::store (MH.intensity_offset(), &NH.scl_inter, is_BE); + + std::string version_string = std::string("MRtrix version: ") + App::mrtrix_version; + if (App::project_version) + version_string += std::string(", project version: ") + App::project_version; + strncpy ( (char*) &NH.descrip, version_string.c_str(), 79); + + Raw::store (NIFTI_XFORM_SCANNER_ANAT, &NH.qform_code, is_BE); + Raw::store (NIFTI_XFORM_SCANNER_ANAT, &NH.sform_code, is_BE); + + // qform: + Eigen::Matrix3d R = M.matrix().topLeftCorner<3,3>(); + if (R.determinant() < 0.0) { + R.col(2) = -R.col(2); + Raw::store (-1.0, &NH.pixdim[0], is_BE); + } + Eigen::Quaterniond Q (R); + + if (Q.w() < 0.0) + Q.vec() = -Q.vec(); + + Raw::store (Q.x(), &NH.quatern_b, is_BE); + Raw::store (Q.y(), &NH.quatern_c, is_BE); + Raw::store (Q.z(), &NH.quatern_d, is_BE); + + Raw::store (M(0,3), &NH.qoffset_x, is_BE); + Raw::store (M(1,3), &NH.qoffset_y, is_BE); + Raw::store (M(2,3), &NH.qoffset_z, is_BE); + + // sform: + Raw::store (MH.spacing (axes[0]) * M(0,0), &NH.srow_x[0], is_BE); + Raw::store (MH.spacing (axes[1]) * M(0,1), &NH.srow_x[1], is_BE); + Raw::store (MH.spacing (axes[2]) * M(0,2), &NH.srow_x[2], is_BE); + Raw::store (M (0,3), &NH.srow_x[3], is_BE); + + Raw::store (MH.spacing (axes[0]) * M(1,0), &NH.srow_y[0], is_BE); + Raw::store (MH.spacing (axes[1]) * M(1,1), &NH.srow_y[1], is_BE); + Raw::store (MH.spacing (axes[2]) * M(1,2), &NH.srow_y[2], is_BE); + Raw::store (M (1,3), &NH.srow_y[3], is_BE); + + Raw::store (MH.spacing (axes[0]) * M(2,0), &NH.srow_z[0], is_BE); + Raw::store (MH.spacing (axes[1]) * M(2,1), &NH.srow_z[1], is_BE); + Raw::store (MH.spacing (axes[2]) * M(2,2), &NH.srow_z[2], is_BE); + Raw::store (M (2,3), &NH.srow_z[3], is_BE); + + if (Type::is_version2) { + const char xyzt_units[4] { NIFTI_UNITS_MM, NIFTI_UNITS_MM, NIFTI_UNITS_MM, NIFTI_UNITS_SEC }; + const int32_t* const xyzt_units_as_int_ptr = reinterpret_cast(xyzt_units); + Raw::store (*xyzt_units_as_int_ptr, &NH.xyzt_units, is_BE); + } + else + NH.xyzt_units = SPACE_TIME_TO_XYZT (NIFTI_UNITS_MM, NIFTI_UNITS_SEC); + + if (single_file && File::Config::get_bool ("NIfTI.AutoSaveJSON", false)) { + std::string json_path = MH.name(); + if (Path::has_suffix (json_path, ".nii.gz")) + json_path = json_path.substr (0, json_path.size()-7); + else if (Path::has_suffix (json_path, ".nii")) + json_path = json_path.substr (0, json_path.size()-4); + else + assert (0); + json_path += ".json"; + File::JSON::save (MH, json_path); + } + } + + + + + + // force explicit instantiation: + template size_t read (Header& MH, const nifti_1_header& NH); + template size_t read (Header& MH, const nifti_2_header& NH); + template void write (nifti_1_header& NH, const Header& H, const bool single_file); + template void write (nifti_2_header& NH, const Header& H, const bool single_file); + + + + + + + + + transform_type adjust_transform (const Header& H, vector& axes) { Stride::List strides = Stride::get (H); @@ -83,7 +603,7 @@ namespace MR Stride::symbolise (H); // if .img, reset all strides to defaults, since it can't be assumed - // that downstream software will be able to parse the NIfTI transform + // that downstream software will be able to parse the NIfTI transform if (is_analyse) { for (size_t i = 0; i < H.ndim(); ++i) H.stride(i) = i+1; @@ -106,7 +626,7 @@ namespace MR //CONF data is permitted (most 3rd party software packages don't //CONF support bitwise data). If false (the default), data will be //CONF stored using more widely supported unsigned 8-bit integers. - if (H.datatype() == DataType::Bit) + if (H.datatype() == DataType::Bit) if (!File::Config::get_bool ("NIfTIAllowBitwise", false)) H.datatype() = DataType::UInt8; } diff --git a/core/file/nifti_utils.h b/core/file/nifti_utils.h index a12b7e5a26..9bdb1bc51f 100644 --- a/core/file/nifti_utils.h +++ b/core/file/nifti_utils.h @@ -18,6 +18,11 @@ #define __file_nifti_utils_h__ #include "types.h" +#include "raw.h" +#include "header.h" +#include "file/config.h" +#include "file/nifti1.h" +#include "file/nifti2.h" namespace MR { @@ -28,6 +33,15 @@ namespace MR namespace NIfTI { + template + size_t read (Header& H, const NiftiHeaderType& NH); + + template + void write (NiftiHeaderType& NH, const Header& H, const bool single_file); + + template inline int header_size (const NiftiHeaderType&) { return 348; } + template <> inline int header_size (const nifti_2_header&) { return 540; } + extern bool right_left_warning_issued; transform_type adjust_transform (const Header& H, vector& order); diff --git a/core/formats/analyse.cpp b/core/formats/analyse.cpp index 96afded5ac..b2c0689834 100644 --- a/core/formats/analyse.cpp +++ b/core/formats/analyse.cpp @@ -18,7 +18,6 @@ #include "file/utils.h" #include "file/entry.h" #include "file/nifti_utils.h" -#include "file/nifti1_utils.h" #include "header.h" #include "formats/list.h" #include "image_io/default.h" @@ -36,7 +35,7 @@ namespace MR File::MMap fmap (header_path); try { - File::NIfTI1::read (H, * ( (const nifti_1_header*) fmap.address())); + File::NIfTI::read (H, * ( (const nifti_1_header*) fmap.address())); std::unique_ptr io_handler (new ImageIO::Default (H)); io_handler->files.push_back (File::Entry (H.name())); return io_handler; @@ -75,7 +74,7 @@ namespace MR const std::string hdr_name (H.name().substr (0, H.name().size()-4) + ".hdr"); File::OFStream out (hdr_name); nifti_1_header NH; - File::NIfTI1::write (NH, H, false); + File::NIfTI::write (NH, H, false); out.write ( (char*) &NH, sizeof (nifti_1_header)); out.close(); diff --git a/core/formats/nifti1.cpp b/core/formats/nifti1.cpp index a380c11d82..a5ae8dde17 100644 --- a/core/formats/nifti1.cpp +++ b/core/formats/nifti1.cpp @@ -18,7 +18,6 @@ #include "file/path.h" #include "file/utils.h" #include "file/nifti_utils.h" -#include "file/nifti1_utils.h" #include "header.h" #include "image_io/default.h" #include "formats/list.h" @@ -30,12 +29,12 @@ namespace MR std::unique_ptr NIfTI1::read (Header& H) const { - if (!Path::has_suffix (H.name(), ".nii")) + if (!Path::has_suffix (H.name(), ".nii")) return std::unique_ptr(); try { File::MMap fmap (H.name()); - const size_t data_offset = File::NIfTI1::read (H, * ( (const nifti_1_header*) fmap.address())); + const size_t data_offset = File::NIfTI::read (H, * ( (const nifti_1_header*) fmap.address())); std::unique_ptr handler (new ImageIO::Default (H)); handler->files.push_back (File::Entry (H.name(), data_offset)); return std::move (handler); @@ -72,7 +71,7 @@ namespace MR throw Exception ("NIfTI-1.1 format cannot support more than 7 dimensions for image \"" + H.name() + "\""); nifti_1_header NH; - File::NIfTI1::write (NH, H, true); + File::NIfTI::write (NH, H, true); File::OFStream out (H.name(), std::ios::out | std::ios::binary); out.write ( (char*) &NH, sizeof (nifti_1_header)); nifti1_extender extender; @@ -80,10 +79,10 @@ namespace MR out.write (extender.extension, sizeof (nifti1_extender)); out.close(); - File::resize (H.name(), File::NIfTI1::header_with_ext_size + footprint(H)); + File::resize (H.name(), File::NIfTI::header_size(NH)+4 + footprint(H)); std::unique_ptr handler (new ImageIO::Default (H)); - handler->files.push_back (File::Entry (H.name(), File::NIfTI1::header_with_ext_size)); + handler->files.push_back (File::Entry (H.name(), File::NIfTI::header_size(NH)+4)); return std::move (handler); } diff --git a/core/formats/nifti1_gz.cpp b/core/formats/nifti1_gz.cpp index 7fc9dfd47f..f54912d255 100644 --- a/core/formats/nifti1_gz.cpp +++ b/core/formats/nifti1_gz.cpp @@ -18,7 +18,6 @@ #include "file/path.h" #include "file/gz.h" #include "file/nifti_utils.h" -#include "file/nifti1_utils.h" #include "header.h" #include "image_io/gz.h" #include "formats/list.h" @@ -31,19 +30,20 @@ namespace MR std::unique_ptr NIfTI1_GZ::read (Header& H) const { - if (!Path::has_suffix (H.name(), ".nii.gz")) + if (!Path::has_suffix (H.name(), ".nii.gz")) return std::unique_ptr(); nifti_1_header NH; + const size_t header_size = File::NIfTI::header_size (NH); File::GZ zf (H.name(), "rb"); - zf.read (reinterpret_cast (&NH), File::NIfTI1::header_size); + zf.read (reinterpret_cast (&NH), header_size); zf.close(); try { - const size_t data_offset = File::NIfTI1::read (H, NH); + const size_t data_offset = File::NIfTI::read (H, NH); std::unique_ptr io_handler (new ImageIO::GZ (H, data_offset)); - memcpy (io_handler.get()->header(), &NH, File::NIfTI1::header_size); - memset (io_handler.get()->header() + File::NIfTI1::header_size, 0, sizeof(nifti1_extender)); + memcpy (io_handler.get()->header(), &NH, header_size); + memset (io_handler.get()->header() + header_size, 0, sizeof(nifti1_extender)); io_handler->files.push_back (File::Entry (H.name(), data_offset)); return std::move (io_handler); } catch (...) { @@ -78,13 +78,16 @@ namespace MR if (H.ndim() > 7) throw Exception ("NIfTI-1.1 format cannot support more than 7 dimensions for image \"" + H.name() + "\""); - std::unique_ptr io_handler (new ImageIO::GZ (H, File::NIfTI1::header_with_ext_size)); + const size_t header_size = File::NIfTI::header_size (nifti_1_header()); - File::NIfTI1::write (*reinterpret_cast (io_handler->header()), H, true); - memset (io_handler->header()+File::NIfTI1::header_size, 0, sizeof(nifti1_extender)); + std::unique_ptr io_handler (new ImageIO::GZ (H, header_size+4)); + nifti_1_header& NH = *reinterpret_cast (io_handler->header()); + + File::NIfTI::write (*reinterpret_cast (io_handler->header()), H, true); + memset (io_handler->header()+header_size, 0, sizeof(nifti1_extender)); File::create (H.name()); - io_handler->files.push_back (File::Entry (H.name(), File::NIfTI1::header_with_ext_size)); + io_handler->files.push_back (File::Entry (H.name(), header_size+4)); return std::move (io_handler); } diff --git a/core/formats/nifti2.cpp b/core/formats/nifti2.cpp index e8d1226422..99623b0992 100644 --- a/core/formats/nifti2.cpp +++ b/core/formats/nifti2.cpp @@ -19,7 +19,6 @@ #include "file/utils.h" #include "file/nifti1.h" #include "file/nifti_utils.h" -#include "file/nifti2_utils.h" #include "header.h" #include "image_io/default.h" #include "formats/list.h" @@ -41,7 +40,7 @@ namespace MR File::MMap fmap (header_path); try { - const size_t data_offset = File::NIfTI2::read (H, * ( (const nifti_2_header*) fmap.address())); + const size_t data_offset = File::NIfTI::read (H, * ( (const nifti_2_header*) fmap.address())); std::unique_ptr handler (new ImageIO::Default (H)); handler->files.push_back (File::Entry (H.name(), data_offset)); return std::move (handler); @@ -83,7 +82,7 @@ namespace MR const std::string header_path = single_file ? H.name() : H.name().substr (0, H.name().size()-4) + ".hdr"; nifti_2_header NH; - File::NIfTI2::write (NH, H, true); + File::NIfTI::write (NH, H, true); File::OFStream out (header_path, std::ios::out | std::ios::binary); out.write ( (char*) &NH, sizeof (nifti_2_header)); nifti1_extender extender; @@ -91,7 +90,7 @@ namespace MR out.write (extender.extension, sizeof (nifti1_extender)); out.close(); - const size_t data_offset = single_file ? File::NIfTI2::header_with_ext_size : 0; + const size_t data_offset = single_file ? File::NIfTI::header_size(NH)+4 : 0; if (single_file) File::resize (H.name(), data_offset + footprint(H)); diff --git a/core/formats/nifti2_gz.cpp b/core/formats/nifti2_gz.cpp index 0d5bbc7676..8d98fa82c6 100644 --- a/core/formats/nifti2_gz.cpp +++ b/core/formats/nifti2_gz.cpp @@ -19,7 +19,6 @@ #include "file/gz.h" #include "file/nifti1.h" #include "file/nifti_utils.h" -#include "file/nifti2_utils.h" #include "header.h" #include "image_io/gz.h" #include "formats/list.h" @@ -32,18 +31,19 @@ namespace MR std::unique_ptr NIfTI2_GZ::read (Header& H) const { - if (!Path::has_suffix (H.name(), ".nii.gz")) + if (!Path::has_suffix (H.name(), ".nii.gz")) return std::unique_ptr(); nifti_2_header NH; + const size_t header_size = File::NIfTI::header_size (NH); File::GZ zf (H.name(), "rb"); - zf.read (reinterpret_cast (&NH), File::NIfTI2::header_size); + zf.read (reinterpret_cast (&NH), header_size); zf.close(); - const size_t data_offset = File::NIfTI2::read (H, NH); + const size_t data_offset = File::NIfTI::read (H, NH); std::unique_ptr io_handler (new ImageIO::GZ (H, data_offset)); - memcpy (io_handler.get()->header(), &NH, File::NIfTI2::header_size); - memset (io_handler.get()->header() + File::NIfTI2::header_size, 0, sizeof(nifti1_extender)); + memcpy (io_handler.get()->header(), &NH, header_size); + memset (io_handler.get()->header() + header_size, 0, sizeof(nifti1_extender)); io_handler->files.push_back (File::Entry (H.name(), data_offset)); return std::move (io_handler); @@ -76,13 +76,16 @@ namespace MR if (H.ndim() > 7) throw Exception ("NIfTI-2 format cannot support more than 7 dimensions for image \"" + H.name() + "\""); - std::unique_ptr io_handler (new ImageIO::GZ (H, File::NIfTI2::header_with_ext_size)); + const size_t header_size = File::NIfTI::header_size (nifti_2_header()); - File::NIfTI2::write (*reinterpret_cast (io_handler->header()), H, true); - memset (io_handler->header()+File::NIfTI2::header_size, 0, sizeof(nifti1_extender)); + std::unique_ptr io_handler (new ImageIO::GZ (H, header_size+4)); + nifti_2_header& NH = *reinterpret_cast (io_handler->header()); + + File::NIfTI::write (NH, H, true); + memset (io_handler->header()+header_size, 0, sizeof(nifti1_extender)); File::create (H.name()); - io_handler->files.push_back (File::Entry (H.name(), File::NIfTI2::header_with_ext_size)); + io_handler->files.push_back (File::Entry (H.name(), header_size+4)); return std::move (io_handler); }