Skip to content

Commit

Permalink
Merge pull request #4595 from cookpa/nifti_xyzt_header_units
Browse files Browse the repository at this point in the history
BUG: Set NIFTI header xyzt_units field correctly in output
  • Loading branch information
thewtex authored May 2, 2024
2 parents 66ce55f + d9029dc commit 0ffa6d5
Show file tree
Hide file tree
Showing 6 changed files with 218 additions and 14 deletions.
2 changes: 1 addition & 1 deletion Modules/IO/NIFTI/include/itkNiftiImageIO.h
Original file line number Diff line number Diff line change
Expand Up @@ -256,7 +256,7 @@ class ITKIONIFTI_EXPORT NiftiImageIO : public ImageIOBase
SetNIfTIOrientationFromImageIO(unsigned short origdims, unsigned short dims);

void
SetImageIOOrientationFromNIfTI(unsigned short dims);
SetImageIOOrientationFromNIfTI(unsigned short dims, double spacingscale, double timingscale);

void
SetImageIOMetadataFromNIfTI();
Expand Down
34 changes: 21 additions & 13 deletions Modules/IO/NIFTI/src/itkNiftiImageIO.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -1309,7 +1309,7 @@ NiftiImageIO::ReadImageInformation()
//
// set up the dimension stuff
double spacingscale = 1.0; // default to mm
switch (XYZT_TO_SPACE(this->m_NiftiImage->xyz_units))
switch (this->m_NiftiImage->xyz_units)
{
case NIFTI_UNITS_METER:
spacingscale = 1e3;
Expand All @@ -1322,7 +1322,7 @@ NiftiImageIO::ReadImageInformation()
break;
}
double timingscale = 1.0; // Default to seconds
switch (XYZT_TO_TIME(this->m_NiftiImage->xyz_units))
switch (this->m_NiftiImage->time_units)
{
case NIFTI_UNITS_SEC:
timingscale = 1.0;
Expand Down Expand Up @@ -1393,7 +1393,7 @@ NiftiImageIO::ReadImageInformation()
EncapsulateMetaData<std::string>(thisDic, ITK_InputFilterName, classname);

// set the image orientation
this->SetImageIOOrientationFromNIfTI(dims);
this->SetImageIOOrientationFromNIfTI(dims, spacingscale, timingscale);

// Set the metadata.
this->SetImageIOMetadataFromNIfTI();
Expand Down Expand Up @@ -1503,13 +1503,14 @@ NiftiImageIO::WriteImageInformation()
// dim[2] is required for 2-D volumes,
// dim[3] for 3-D volumes, etc.
this->m_NiftiImage->nvox = 1;
// Spacial dims in ITK are given in mm.
// Spatial dims in ITK are given in mm.
// If 4D assume 4thD is in SECONDS, for all of ITK.
// NOTE: Due to an ambiguity in the nifti specification, some developers
// external tools believe that the time units must be set, even if there
// is only one dataset. Having the time specified for a purly spatial
// is only one dataset. Having the time specified for a purely spatial
// image has no consequence, so go ahead and set it to seconds.
this->m_NiftiImage->xyz_units = static_cast<int>(NIFTI_UNITS_MM | NIFTI_UNITS_SEC);
this->m_NiftiImage->xyz_units = static_cast<int>(NIFTI_UNITS_MM);
this->m_NiftiImage->time_units = static_cast<int>(NIFTI_UNITS_SEC);
this->m_NiftiImage->dim[7] = this->m_NiftiImage->nw = 1;
this->m_NiftiImage->dim[6] = this->m_NiftiImage->nv = 1;
this->m_NiftiImage->dim[5] = this->m_NiftiImage->nu = 1;
Expand Down Expand Up @@ -1537,6 +1538,9 @@ NiftiImageIO::WriteImageInformation()
case 4:
this->m_NiftiImage->dim[4] = this->m_NiftiImage->nt = this->GetDimensions(3);
this->m_NiftiImage->pixdim[4] = this->m_NiftiImage->dt = static_cast<float>(this->GetSpacing(3));
// Add time origin here because it is not set with the spatial origin in
// SetNIfTIOrientationFromImageIO
this->m_NiftiImage->toffset = static_cast<float>(this->GetOrigin(3));
this->m_NiftiImage->nvox *= this->m_NiftiImage->dim[4];
[[fallthrough]];
case 3:
Expand Down Expand Up @@ -1861,7 +1865,7 @@ IsAffine(const mat44 & nifti_mat)
} // namespace

void
NiftiImageIO::SetImageIOOrientationFromNIfTI(unsigned short dims)
NiftiImageIO::SetImageIOOrientationFromNIfTI(unsigned short dims, double spacingscale, double timingscale)
{
typedef SpatialOrientationAdapter OrientAdapterType;
// in the case of an Analyze75 file, use old analyze orient method.
Expand Down Expand Up @@ -1929,7 +1933,7 @@ NiftiImageIO::SetImageIOOrientationFromNIfTI(unsigned short dims)
// so check it first.
// Commonly the 4x4 double precision dicom information is stored in
// the 4x4 single precision sform fields, and that original representation
// is converted (with lossy conversoin) into the qform representation.
// is converted (with lossy conversion) into the qform representation.
const bool qform_sform_are_similar = [=]() -> bool {
vnl_matrix_fixed<float, 4, 4> sto_xyz{ &(this->m_NiftiImage->sto_xyz.m[0][0]) };
vnl_matrix_fixed<float, 4, 4> qto_xyz{ &(this->m_NiftiImage->qto_xyz.m[0][0]) };
Expand Down Expand Up @@ -2122,14 +2126,18 @@ NiftiImageIO::SetImageIOOrientationFromNIfTI(unsigned short dims)

//
// set origin
m_Origin[0] = -theMat.m[0][3];
m_Origin[0] = -theMat.m[0][3] * spacingscale;
if (dims > 1)
{
m_Origin[1] = -theMat.m[1][3];
m_Origin[1] = -theMat.m[1][3] * spacingscale;
}
if (dims > 2)
{
m_Origin[2] = theMat.m[2][3];
m_Origin[2] = theMat.m[2][3] * spacingscale;
}
if (dims > 3)
{
m_Origin[3] = this->m_NiftiImage->toffset * timingscale;
}

const int max_defined_orientation_dims = (dims > 3) ? 3 : dims;
Expand Down Expand Up @@ -2364,7 +2372,7 @@ NiftiImageIO::Write(const void * buffer)
// so will destructor of the image that really owns it.
if (nifti_write_status)
{
itkExceptionMacro("ERROR: nifti library failed to write image" << this->GetFileName());
itkExceptionMacro("ERROR: nifti library failed to write image: " << this->GetFileName());
}
}
else /// Image intent is vector image
Expand Down Expand Up @@ -2469,7 +2477,7 @@ NiftiImageIO::Write(const void * buffer)
this->m_NiftiImage->data = nullptr; // if left pointing to data buffer
if (nifti_write_status)
{
itkExceptionMacro("ERROR: nifti library failed to write image" << this->GetFileName());
itkExceptionMacro("ERROR: nifti library failed to write image: " << this->GetFileName());
}
}
}
Expand Down
23 changes: 23 additions & 0 deletions Modules/IO/NIFTI/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ set(ITKIONIFTITests
itkNiftiImageIOTest11.cxx
itkNiftiImageIOTest12.cxx
itkNiftiImageIOTest13.cxx
itkNiftiImageIOTest14.cxx
itkNiftiLargeImageRegionReadTest.cxx
itkNiftiReadAnalyzeTest.cxx
itkNiftiReadWriteDirectionTest.cxx
Expand Down Expand Up @@ -245,6 +246,28 @@ itk_add_test(
itkNiftiImageIOTest13
DATA{Input/SmallVoxels_AffinePrecision.nii.gz})

itk_add_test(
NAME
itkNiftiSpatialTemporalUnitsTest
COMMAND
ITKIONIFTITestDriver
itkNiftiImageIOTest14
${ITK_TEST_OUTPUT_DIR}
DATA{Input/xyzt_units_test_scl_mm_s.nii.gz}
DATA{Input/xyzt_units_test_scl.nii.gz}
)

itk_add_test(
NAME
itkNiftiSpatialTemporalUnitsTest2
COMMAND
ITKIONIFTITestDriver
itkNiftiImageIOTest14
${ITK_TEST_OUTPUT_DIR}
DATA{Input/xyzt_units_test_scl_mm_s.nii.gz}
DATA{Input/xyzt_units_test_scl_mm_s.nii.gz}
)

itk_add_test(
NAME
itkNiftiQFormSFormDifferentSpacingTest
Expand Down
1 change: 1 addition & 0 deletions Modules/IO/NIFTI/test/Input/xyzt_units_test_scl.nii.gz.cid
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
bafkreifgv6rke262bhoeaa47lu3fmrtcv4gdfm4wjynctbuhkr3w43ccta
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
bafkreiawqehmd75yk334c7pxcjumrdw4ide2dgg7x3ctopotkzyfz7cjva
171 changes: 171 additions & 0 deletions Modules/IO/NIFTI/test/itkNiftiImageIOTest14.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,171 @@
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/

#include "itkNiftiImageIOTest.h"
#include "itkMath.h"
#include "itkTestingMacros.h"

// Test scaling of different spatial or temporal units
// ITK uses mm for space and s for time. NIFTI allows other storage units
// These should be converted to mm and s on output
//
// NIFTI xyzt_units field should be 10 (mm, sec) for all output
int
itkNiftiImageIOTest14(int argc, char * argv[])
{
// usage: itkNiftiImageIOTest14 test_dir ref_image test_image
// images should have the same size, spacing, and origin, but may have different units
if (argc != 4)
{
std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv) << " test_dir ref_image test_image" << std::endl;
return EXIT_FAILURE;
}

// first arg is output dir
char * test_dir = argv[1];
itksys::SystemTools::ChangeDirectory(test_dir);

// second arg is reference image in mm and sec
char * ref_image_fn = argv[2];

// third arg is the test image in different units
char * test_image_fn = argv[3];

constexpr unsigned int Dimension = 4;

using PixelType = float;
using ImageType = itk::Image<PixelType, Dimension>;

auto ref_image = itk::IOTestHelper::ReadImage<ImageType>(std::string(ref_image_fn));
auto test_image = itk::IOTestHelper::ReadImage<ImageType>(std::string(test_image_fn));

// check spacing of ref and test are the same
const auto ref_spacing = ref_image->GetSpacing();
const auto test_spacing = test_image->GetSpacing();

// check size of bounding box of ref and test are the same
const auto ref_region_size = ref_image->GetLargestPossibleRegion().GetSize();
const auto test_region_size = test_image->GetLargestPossibleRegion().GetSize();

// check origin
const auto ref_origin = ref_image->GetOrigin();
const auto test_origin = test_image->GetOrigin();

for (unsigned int i = 0; i < Dimension; i++)
{
// check spacing is the same within tolerance
if (!itk::Math::FloatAlmostEqual(ref_spacing[i], test_spacing[i], 4, 1e-6))
{
std::cerr << "Spacing of reference and test images do not match" << std::endl;
return EXIT_FAILURE;
}

// check bounding box is the same
if (ref_region_size[i] != test_region_size[i])
{
std::cerr << "Size of reference and test images do not match" << std::endl;
return EXIT_FAILURE;
}

// check origin is the same within tolerance
if (!itk::Math::FloatAlmostEqual(ref_origin[i], test_origin[i], 4, 1e-6))
{
std::cerr << "Origin of reference and test images do not match" << std::endl;
return EXIT_FAILURE;
}
}

const char * output_test_fn = "xyzt_units_output_check.nii";

// set the origin of time to a non-zero value
auto newOrigin = ImageType::PointType();
newOrigin[0] = 1.0;
newOrigin[1] = 1.0;
newOrigin[2] = 1.0;
newOrigin[3] = 2.0;
test_image->SetOrigin(newOrigin);

itk::IOTestHelper::WriteImage<ImageType, itk::NiftiImageIO>(test_image, std::string(output_test_fn));

bool fileHasCorrectXYZTTUnits = false;
bool fileHasCorrectTimeOrigin = false;

try
{
// read the image back in
ImageType::Pointer image_from_disk = itk::IOTestHelper::ReadImage<ImageType>(std::string(output_test_fn));

// check the metadata for xyzt_units and toffset
itk::MetaDataDictionary & dictionary = image_from_disk->GetMetaDataDictionary();

std::string toffset;
if (!itk::ExposeMetaData<std::string>(dictionary, "toffset", toffset))
{
std::cerr << "toffset not found in metadata" << std::endl;
}
if (!itk::Math::FloatAlmostEqual(std::stod(toffset), 2.0, 4, 1e-6))
{
fileHasCorrectTimeOrigin = true;
}
else
{
std::cerr << "toffset not set correctly in metadata" << std::endl;
fileHasCorrectTimeOrigin = false;
}

auto read_origin = test_image->GetOrigin();

if (!itk::Math::FloatAlmostEqual(read_origin[3], 2.0, 4, 1e-6))
{
std::cerr << "Time origin not read back correctly from toffset field" << std::endl;
fileHasCorrectTimeOrigin = false;
}
else
{
fileHasCorrectTimeOrigin = true;
}


std::string xyzt_units;
if (!itk::ExposeMetaData<std::string>(dictionary, "xyzt_units", xyzt_units))
{
std::cerr << "xyzt_units not found in metadata" << std::endl;
}
if (xyzt_units == "10")
{
fileHasCorrectXYZTTUnits = true;
}
else
{
std::cerr << "xyzt_units not set correctly in metadata" << std::endl;
fileHasCorrectXYZTTUnits = false;
}
}
catch (...)
{
std::cerr << "Exception caught while reading image back in" << std::endl;
}

itk::IOTestHelper::Remove(output_test_fn);

if (fileHasCorrectXYZTTUnits && fileHasCorrectTimeOrigin)
{
return EXIT_SUCCESS;
}
return EXIT_FAILURE;
}

0 comments on commit 0ffa6d5

Please sign in to comment.