diff --git a/Modules/IO/MINC/include/itkMINCImageIO.h b/Modules/IO/MINC/include/itkMINCImageIO.h index 67a82966ed4..8850f913def 100644 --- a/Modules/IO/MINC/include/itkMINCImageIO.h +++ b/Modules/IO/MINC/include/itkMINCImageIO.h @@ -85,11 +85,12 @@ class ITKIOMINC_EXPORT MINCImageIO : public ImageIOBase /** \see LightObject::GetNameOfClass() */ itkOverrideGetNameOfClassMacro(MINCImageIO); - /** Right now MINC supports up to 3D with multiple components */ + /** MINC supports up to 3 spatial dimensions plus a time dimension (4D), + * each with multiple components (MINC vector_dimension). */ bool SupportsDimension(unsigned long dim) override { - return dim < 4; + return dim < 5; } /*-------- This part of the interface deals with reading data. ------ */ diff --git a/Modules/IO/MINC/src/itkMINCImageIO.cxx b/Modules/IO/MINC/src/itkMINCImageIO.cxx index 07c5c5bf0bb..233954a82d4 100644 --- a/Modules/IO/MINC/src/itkMINCImageIO.cxx +++ b/Modules/IO/MINC/src/itkMINCImageIO.cxx @@ -440,12 +440,14 @@ MINCImageIO::ReadImageInformation() itkExceptionStringMacro(" minc files without spatial dimensions are not supported!"); } - if (m_MINCPImpl->m_DimensionIndices[0] != -1 && m_MINCPImpl->m_DimensionIndices[4] != -1) - { - itkExceptionStringMacro(" 4D minc files vector dimension are not supported currently"); - } + const bool haveTimeDimension = (m_MINCPImpl->m_DimensionIndices[4] != -1); - this->SetNumberOfDimensions(spatial_dimension_count); + // The MINC time dimension becomes an additional (highest) ITK dimension, while + // the MINC vector_dimension is mapped to ITK components. The two are + // independent and may coexist (e.g. a 4D vector image). + const unsigned int itkDimensionCount = spatial_dimension_count + (haveTimeDimension ? 1 : 0); + + this->SetNumberOfDimensions(itkDimensionCount); int numberOfComponents = 1; unsigned int usableDimensions = 0; @@ -457,11 +459,37 @@ MINCImageIO::ReadImageInformation() auto RAStofromLPS = Matrix::GetIdentity(); RAStofromLPS(0, 0) = -1.0; RAStofromLPS(1, 1) = -1.0; - std::vector dir_cos_temp(3); + std::vector dir_cos_temp(itkDimensionCount, 0.0); Vector origin{}; Vector oOrigin{}; + // Build the MINC apparent dimension order slowest-varying first. The time + // dimension is the slowest, matching the ITK buffer layout used by + // Read()/Write() where it is the highest ITK dimension. + if (haveTimeDimension) + { + m_MINCPImpl->m_MincApparentDims[usableDimensions] = m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[4]]; + // always use positive + miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_POSITIVE); + misize_t _sz = 0; + miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_sz); + double _sep = NAN; + miget_dimension_separation(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_ORDER_APPARENT, &_sep); + double _start = NAN; + miget_dimension_start(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_ORDER_APPARENT, &_start); + + this->SetDimensions(spatial_dimension_count, static_cast(_sz)); + this->SetSpacing(spatial_dimension_count, _sep); + this->SetOrigin(spatial_dimension_count, _start); + + std::vector time_dir(itkDimensionCount, 0.0); + time_dir[spatial_dimension_count] = 1.0; + this->SetDirection(spatial_dimension_count, time_dir); + + ++usableDimensions; + } + // minc api uses inverse order of dimensions , fastest varying are last Vector sep; for (int i = 3; i > 0; i--) @@ -511,7 +539,7 @@ MINCImageIO::ReadImageInformation() for (int i = 0; i < spatial_dimension_count; ++i) { this->SetOrigin(i, oOrigin[i]); - for (unsigned int j = 0; j < 3; j++) + for (unsigned int j = 0; j < static_cast(spatial_dimension_count); ++j) { dir_cos_temp[j] = dir_cos[j][i]; } @@ -530,18 +558,6 @@ MINCImageIO::ReadImageInformation() ++usableDimensions; } - if (m_MINCPImpl->m_DimensionIndices[4] != -1) // have time dimension - { - // micopy_dimension(hdim[m_MINCPImpl->m_DimensionIndices[4]],&apparent_dimension_order[usable_dimensions]); - m_MINCPImpl->m_MincApparentDims[usableDimensions] = m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[4]]; - // always use positive - miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_POSITIVE); - misize_t _sz = 0; - miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_sz); - numberOfComponents = _sz; - ++usableDimensions; - } - // Set apparent dimension order to the MINC2 api if (miset_apparent_dimension_order(m_MINCPImpl->m_Volume, usableDimensions, m_MINCPImpl->m_MincApparentDims) < 0) { @@ -893,33 +909,9 @@ MINCImageIO::WriteImageInformation() MetaDataDictionary & thisDic = GetMetaDataDictionary(); unsigned int minc_dimensions = 0; - if (nComp > 3) // last dimension will be either vector or time - { - micreate_dimension(MItime, - MI_DIMCLASS_TIME, - MI_DIMATTR_REGULARLY_SAMPLED, - nComp, - &m_MINCPImpl->m_MincApparentDims[m_MINCPImpl->m_NDims - minc_dimensions - 1]); - - double tstart = 0.0; - if (!ExposeMetaData(thisDic, "tstart", tstart)) - { - tstart = 0.0; - } - - miset_dimension_start(m_MINCPImpl->m_MincApparentDims[m_MINCPImpl->m_NDims - minc_dimensions - 1], tstart); - - double tstep = 1.0; - if (!ExposeMetaData(thisDic, "tstep", tstep)) - { - tstep = 1.0; - } - - miset_dimension_separation(m_MINCPImpl->m_MincApparentDims[m_MINCPImpl->m_NDims - minc_dimensions - 1], tstep); - - ++minc_dimensions; - } - else if (nComp > 1) + // ITK components are always written as a MINC vector_dimension (fastest + // varying). A time series is an ITK dimension, not components (see below). + if (nComp > 1) { micreate_dimension(MIvector_dimension, MI_DIMCLASS_RECORD, @@ -956,10 +948,20 @@ MINCImageIO::WriteImageInformation() ++minc_dimensions; } - if (nDims > 3) + if (nDims > 4) { MINCIOFreeTmpDimHandle(minc_dimensions, m_MINCPImpl->m_MincApparentDims); - itkExceptionStringMacro("Unfortunately, only up to 3D volume are supported now."); + itkExceptionStringMacro("Unfortunately, only up to 4D volumes (3 spatial + time) are supported."); + } + + if (nDims > 3) // the 4th ITK dimension is written as a MINC time dimension (slowest varying) + { + micreate_dimension(MItime, + MI_DIMCLASS_TIME, + MI_DIMATTR_REGULARLY_SAMPLED, + this->GetDimensions(3), + &m_MINCPImpl->m_MincApparentDims[m_MINCPImpl->m_NDims - minc_dimensions - 1]); + ++minc_dimensions; } // allocating dimensions @@ -994,21 +996,18 @@ MINCImageIO::WriteImageInformation() for (unsigned int i = 0; i < nDims; ++i) { const unsigned int j = i + (nComp > 1 ? 1 : 0); - double dir_cos[3]; - for (unsigned int k = 0; k < 3; ++k) + const unsigned int apparentIndex = minc_dimensions - j - 1; + miset_dimension_separation(m_MINCPImpl->m_MincApparentDims[apparentIndex], this->GetSpacing(i)); + miset_dimension_start(m_MINCPImpl->m_MincApparentDims[apparentIndex], origin[i]); + if (i < 3) // spatial dimension: also set the direction cosines (time has none) { - if (k < nDims) + double dir_cos[3]; + for (unsigned int k = 0; k < 3; ++k) { - dir_cos[k] = directionCosineMatrix[i][k]; - } - else - { - dir_cos[k] = 0.0; + dir_cos[k] = (k < nDims) ? directionCosineMatrix[i][k] : 0.0; } + miset_dimension_cosines(m_MINCPImpl->m_MincApparentDims[apparentIndex], dir_cos); } - miset_dimension_separation(m_MINCPImpl->m_MincApparentDims[minc_dimensions - j - 1], this->GetSpacing(i)); - miset_dimension_start(m_MINCPImpl->m_MincApparentDims[minc_dimensions - j - 1], origin[i]); - miset_dimension_cosines(m_MINCPImpl->m_MincApparentDims[minc_dimensions - j - 1], dir_cos); } // TODO: fix this to appropriate @@ -1107,14 +1106,14 @@ MINCImageIO::WriteImageInformation() break; case 't': case 'T': - if (nComp <= 1) + if (nDims <= 3) // no time dimension present { itkDebugMacro("Dimension order is incorrect " << dimension_order); dimorder_good = false; } else { - j = m_MINCPImpl->m_NDims - 1; + j = 0; // time is the slowest-varying apparent dimension } break; case 'x': diff --git a/Modules/IO/MINC/test/CMakeLists.txt b/Modules/IO/MINC/test/CMakeLists.txt index e0c78181660..3467cf06a21 100644 --- a/Modules/IO/MINC/test/CMakeLists.txt +++ b/Modules/IO/MINC/test/CMakeLists.txt @@ -129,6 +129,7 @@ itk_add_test( ${ITK_TEST_OUTPUT_DIR}/dti_sample.mnc ITK_REMOVE_TEMPORARY_TEST_FILES ${ITK_TEST_OUTPUT_DIR}/dti_sample.mnc + ${ITK_TEST_OUTPUT_DIR}/dti_sample.mnc_synthetic.mnc ) itk_add_test( diff --git a/Modules/IO/MINC/test/itkMINCImageIOTest2.cxx b/Modules/IO/MINC/test/itkMINCImageIOTest2.cxx index ab0061d0695..1d898be08d7 100644 --- a/Modules/IO/MINC/test/itkMINCImageIOTest2.cxx +++ b/Modules/IO/MINC/test/itkMINCImageIOTest2.cxx @@ -45,7 +45,7 @@ itkMINCImageIOTest2(int argc, char * argv[]) ITK_EXERCISE_BASIC_OBJECT_METHODS(mincIO1, MINCImageIO, ImageIOBase); - constexpr unsigned int supportedDimCount{ 4 }; // includes the degenerate 0-dimensional case + constexpr unsigned int supportedDimCount{ 5 }; // 3 spatial + time, plus the degenerate 0-dimensional case std::vector supportedDims(supportedDimCount); std::iota(std::begin(supportedDims), std::end(supportedDims), 0); for (const auto & value : supportedDims) diff --git a/Modules/IO/MINC/test/itkMINCImageIOTest_4D.cxx b/Modules/IO/MINC/test/itkMINCImageIOTest_4D.cxx index 9b708523829..f71a2d9b18b 100644 --- a/Modules/IO/MINC/test/itkMINCImageIOTest_4D.cxx +++ b/Modules/IO/MINC/test/itkMINCImageIOTest_4D.cxx @@ -16,12 +16,14 @@ * *=========================================================================*/ -#include "itkImageAlgorithm.h" #include #include "itkMINCImageIOFactory.h" #include "itkImageFileReader.h" #include "itkImageFileWriter.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionConstIterator.h" +#include "itkMath.h" #include "itkTestingMacros.h" @@ -38,7 +40,10 @@ itkMINCImageIOTest_4D(int argc, char * argv[]) itk::MINCImageIOFactory::RegisterOneFactory(); - using ImageType = itk::VectorImage; + // A MINC volume with x/y/z spatial dimensions plus a time dimension is read as + // a genuine 4D itk::Image, with the time step/start mapped to the 4th + // dimension spacing/origin. + using ImageType = itk::Image; using ReaderType = itk::ImageFileReader; using WriterType = itk::ImageFileWriter; @@ -52,11 +57,72 @@ itkMINCImageIOTest_4D(int argc, char * argv[]) ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); - const ImageType::ConstPointer image = reader->GetOutput(); - image->Print(std::cout); + // Round-trip a synthetic 4D volume with a non-trivial time step/start and + // verify the time geometry survives, together with the voxel values. + const std::string synthetic = std::string(argv[2]) + "_synthetic.mnc"; + + auto synthImage = ImageType::New(); + constexpr float timeStep = 2.5F; + constexpr float timeStart = 10.0F; + + ImageType::RegionType region; + region.SetSize({ { 3, 4, 5, 2 } }); + synthImage->SetRegions(region); + + ImageType::SpacingType spacing; + spacing[0] = 1.0; + spacing[1] = 1.0; + spacing[2] = 1.0; + spacing[3] = timeStep; + synthImage->SetSpacing(spacing); + + ImageType::PointType origin; + origin[0] = 0.0; + origin[1] = 0.0; + origin[2] = 0.0; + origin[3] = timeStart; + synthImage->SetOrigin(origin); + + synthImage->Allocate(); + + { + float v = 0.0F; + itk::ImageRegionIterator it(synthImage, synthImage->GetLargestPossibleRegion()); + for (; !it.IsAtEnd(); ++it, ++v) + { + it.Set(v); + } + } + + auto synthWriter = WriterType::New(); + synthWriter->SetFileName(synthetic); + synthWriter->SetInput(synthImage); + ITK_TRY_EXPECT_NO_EXCEPTION(synthWriter->Update()); + + auto synthReader = ReaderType::New(); + synthReader->SetFileName(synthetic); + ITK_TRY_EXPECT_NO_EXCEPTION(synthReader->Update()); + + const ImageType::ConstPointer synthBack = synthReader->GetOutput(); + + ITK_TEST_EXPECT_EQUAL(synthBack->GetLargestPossibleRegion().GetSize()[3], 2); + ITK_TEST_EXPECT_TRUE(itk::Math::FloatAlmostEqual(static_cast(synthBack->GetSpacing()[3]), timeStep)); + ITK_TEST_EXPECT_TRUE(itk::Math::FloatAlmostEqual(static_cast(synthBack->GetOrigin()[3]), timeStart)); + + itk::ImageRegionConstIterator inIt(synthImage, synthImage->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator outIt(synthBack, synthBack->GetLargestPossibleRegion()); + for (; !inIt.IsAtEnd(); ++inIt, ++outIt) + { + if (itk::Math::NotAlmostEquals(inIt.Get(), outIt.Get())) + { + std::cerr << "Voxel mismatch after 4D round-trip at " << inIt.ComputeIndex() << ": " << inIt.Get() + << " != " << outIt.Get() << std::endl; + return EXIT_FAILURE; + } + } std::cout << "Test finished." << std::endl; return EXIT_SUCCESS;