From a8811a05e99bc3ea0bdc838361e84e087120482b Mon Sep 17 00:00:00 2001 From: "Gabriel A. Devenyi" Date: Fri, 27 Sep 2024 16:22:21 -0400 Subject: [PATCH] BUG: Fix TransformMINC convert MINC PositiveCoordinateOrientation RAS coordinate system to ITK PositiveCoordinateOrientation LPS --- CMakeLists.txt | 4 + Modules/IO/MINC/include/itkMINCImageIO.h | 6 + Modules/IO/MINC/src/CMakeLists.txt | 3 + Modules/IO/MINC/src/itkMINCImageIO.cxx | 34 ++-- .../src/itkMINCImageIOConfigurePrivate.h.in | 27 +++ Modules/IO/MINC/test/CMakeLists.txt | 79 +++++--- Modules/IO/MINC/test/itkMINCImageIOTest4.cxx | 46 +++-- .../include/itkMINCTransformAdapter.h | 1 + .../include/itkMINCTransformIO.h | 10 +- Modules/IO/TransformMINC/src/CMakeLists.txt | 3 + .../src/itkMINCImageIOConfigurePrivate.h.in | 27 +++ .../TransformMINC/src/itkMINCTransformIO.cxx | 145 ++++++++++++-- .../test/itkIOTransformMINCTest.cxx | 178 ++++++++++++------ .../test/itkMINCTransformAdapterTest.cxx | 11 +- 14 files changed, 452 insertions(+), 122 deletions(-) create mode 100644 Modules/IO/MINC/src/itkMINCImageIOConfigurePrivate.h.in create mode 100644 Modules/IO/TransformMINC/src/itkMINCImageIOConfigurePrivate.h.in diff --git a/CMakeLists.txt b/CMakeLists.txt index cc05e3deb5d..520ca0984ca 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -372,6 +372,10 @@ mark_as_advanced(ITK_NIFTI_IO_ANALYZE_FLAVOR) option(ITK_NIFTI_IO_SFORM_PERMISSIVE_DEFAULT "Allow non-orthogonal rotation matrix in NIFTI sform by default" OFF) mark_as_advanced(ITK_NIFTI_IO_SFORM_PERMISSIVE_DEFAULT) +# Automatically convert MINC coordinate system from PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS +option(ITK_MINC_IO_RAS_TO_LPS "Convert MINC coordinate system from PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS by default" ON) + +mark_as_advanced(ITK_MINC_IO_RAS_TO_LPS) #----------------------------------------------------------------------------- # ITK build classes that are in the review process # ITK_USE_REVIEW is deprecated, only kept for backward compatibility diff --git a/Modules/IO/MINC/include/itkMINCImageIO.h b/Modules/IO/MINC/include/itkMINCImageIO.h index f8c8b46a536..822b4535567 100644 --- a/Modules/IO/MINC/include/itkMINCImageIO.h +++ b/Modules/IO/MINC/include/itkMINCImageIO.h @@ -123,6 +123,11 @@ class ITKIOMINC_EXPORT MINCImageIO : public ImageIOBase void Write(const void * buffer) override; + /** Set to automatically convert from PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS*/ + itkSetMacro(RAStoLPS, bool); + itkGetConstMacro(RAStoLPS, bool); + itkBooleanMacro(RAStoLPS); + protected: MINCImageIO(); ~MINCImageIO() override; @@ -153,6 +158,7 @@ class ITKIOMINC_EXPORT MINCImageIO : public ImageIOBase // complex type images, composed of complex numbers // int m_Complex; + bool m_RAStoLPS; }; } // end namespace itk diff --git a/Modules/IO/MINC/src/CMakeLists.txt b/Modules/IO/MINC/src/CMakeLists.txt index db5afa07594..e3d20fed5e1 100644 --- a/Modules/IO/MINC/src/CMakeLists.txt +++ b/Modules/IO/MINC/src/CMakeLists.txt @@ -1,3 +1,6 @@ +configure_file(itkMINCImageIOConfigurePrivate.h.in itkMINCImageIOConfigurePrivate.h @ONLY) +include_directories("${CMAKE_CURRENT_BINARY_DIR}") + set(ITKIOMINC_SRC itkMINCImageIO.cxx itkMINCImageIOFactory.cxx) add_library(ITKIOMINC ${ITK_LIBRARY_BUILD_TYPE} ${ITKIOMINC_SRC}) diff --git a/Modules/IO/MINC/src/itkMINCImageIO.cxx b/Modules/IO/MINC/src/itkMINCImageIO.cxx index 90f2df8a8dd..5a2fccb1223 100644 --- a/Modules/IO/MINC/src/itkMINCImageIO.cxx +++ b/Modules/IO/MINC/src/itkMINCImageIO.cxx @@ -29,6 +29,9 @@ #include // For unique_ptr. +#include "itkMINCImageIOConfigurePrivate.h" + + extern "C" { void @@ -235,6 +238,7 @@ MINCImageIO::CloseVolume() MINCImageIO::MINCImageIO() : m_MINCPImpl(std::make_unique()) + , m_RAStoLPS(ITK_MINC_IO_RAS_TO_LPS) { for (auto & dimensionIndex : m_MINCPImpl->m_DimensionIndices) { @@ -283,6 +287,7 @@ MINCImageIO::PrintSelf(std::ostream & os, Indent indent) const os << indent << "MINCPImpl: " << m_MINCPImpl.get() << std::endl; os << indent << "DirectionCosines: " << m_DirectionCosines << std::endl; + os << indent << "RAStoLPS: " << m_RAStoLPS << std::endl; } void @@ -511,9 +516,11 @@ MINCImageIO::ReadImageInformation() } } + // Transform MINC PositiveCoordinateOrientation RAS coordinates to // internal ITK PositiveCoordinateOrientation LPS Coordinates - dir_cos = RAS_tofrom_LPS * dir_cos; + if (this->m_RAStoLPS) + dir_cos = RAS_tofrom_LPS * dir_cos; // Transform origin coordinates o_origin = dir_cos * origin; @@ -688,6 +695,9 @@ MINCImageIO::ReadImageInformation() std::string classname(GetNameOfClass()); // EncapsulateMetaData(thisDic,ITK_InputFilterName, // classname); + // preserve information if the volume was PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS + // converted + EncapsulateMetaData(thisDic, "RAStoLPS", m_RAStoLPS); // store history size_t minc_history_length = 0; @@ -997,8 +1007,10 @@ MINCImageIO::WriteImageInformation() const vnl_matrix inverseDirectionCosines{ vnl_matrix_inverse(dircosmatrix).as_matrix() }; origin *= inverseDirectionCosines; // transform to minc convention + // Convert ITK direction cosines from PositiveCoordinateOrientation LPS to PositiveCoordinateOrientation RAS - dircosmatrix *= RAS_tofrom_LPS; + if (this->m_RAStoLPS) + dircosmatrix *= RAS_tofrom_LPS; for (unsigned int i = 0; i < nDims; ++i) { @@ -1028,27 +1040,21 @@ MINCImageIO::WriteImageInformation() { case IOComponentEnum::UCHAR: m_MINCPImpl->m_Volume_type = MI_TYPE_UBYTE; - // m_MINCPImpl->m_Volume_class=MI_CLASS_INT; break; case IOComponentEnum::CHAR: m_MINCPImpl->m_Volume_type = MI_TYPE_BYTE; - // m_MINCPImpl->m_Volume_class=MI_CLASS_INT; break; case IOComponentEnum::USHORT: m_MINCPImpl->m_Volume_type = MI_TYPE_USHORT; - // m_MINCPImpl->m_Volume_class=MI_CLASS_INT; break; case IOComponentEnum::SHORT: m_MINCPImpl->m_Volume_type = MI_TYPE_SHORT; - // m_MINCPImpl->m_Volume_class=MI_CLASS_INT; break; case IOComponentEnum::UINT: m_MINCPImpl->m_Volume_type = MI_TYPE_UINT; - // m_MINCPImpl->m_Volume_class=MI_CLASS_INT; break; case IOComponentEnum::INT: m_MINCPImpl->m_Volume_type = MI_TYPE_INT; - // m_MINCPImpl->m_Volume_class=MI_CLASS_INT; break; // case IOComponentEnum::ULONG://TODO: make sure we are cross-platform here! // volume_data_type=MI_TYPE_ULONG; @@ -1056,10 +1062,10 @@ MINCImageIO::WriteImageInformation() // case IOComponentEnum::LONG://TODO: make sure we are cross-platform here! // volume_data_type=MI_TYPE_LONG; // break; - case IOComponentEnum::FLOAT: // TODO: make sure we are cross-platform here! + case IOComponentEnum::FLOAT: m_MINCPImpl->m_Volume_type = MI_TYPE_FLOAT; break; - case IOComponentEnum::DOUBLE: // TODO: make sure we are cross-platform here! + case IOComponentEnum::DOUBLE: m_MINCPImpl->m_Volume_type = MI_TYPE_DOUBLE; break; default: @@ -1099,7 +1105,6 @@ MINCImageIO::WriteImageInformation() if (ExposeMetaData(thisDic, "dimension_order", dimension_order)) { // the format should be ((+|-)(X|Y|Z|V|T))* - // std::cout<<"Restoring original dimension order:"<m_RAStoLPS; + miset_attr_values(m_MINCPImpl->m_Volume, MI_TYPE_INT, "itk", "RAStoLPS", 1, &tmp); + } + mifree_volume_props(hprops); } diff --git a/Modules/IO/MINC/src/itkMINCImageIOConfigurePrivate.h.in b/Modules/IO/MINC/src/itkMINCImageIOConfigurePrivate.h.in new file mode 100644 index 00000000000..cd43c9181ea --- /dev/null +++ b/Modules/IO/MINC/src/itkMINCImageIOConfigurePrivate.h.in @@ -0,0 +1,27 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ +#ifndef itkMINCImageIOConfigurePrivate_h +#define itkMINCImageIOConfigurePrivate_h + +// This file is intended to hold preprocessor macros only used internally by +// IO/MINC module and should not be installed. + +// Convert from PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS by default +#cmakedefine01 ITK_MINC_IO_RAS_TO_LPS + +#endif // itkMINCImageIOConfigurePrivate_h diff --git a/Modules/IO/MINC/test/CMakeLists.txt b/Modules/IO/MINC/test/CMakeLists.txt index 6f2b2f8551e..050e6d78488 100644 --- a/Modules/IO/MINC/test/CMakeLists.txt +++ b/Modules/IO/MINC/test/CMakeLists.txt @@ -1,5 +1,9 @@ itk_module_test() +# to know system default MINC PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS mode +include_directories("${CMAKE_CURRENT_BINARY_DIR}/../src") + + if(NOT ITK_USE_SYSTEM_HDF5) include_directories("${ITKHDF5_SOURCE_DIR}/src/itkhdf5" "${ITKHDF5_SOURCE_DIR}/src/itkhdf5/src" "${ITKHDF5_BINARY_DIR}/src/itkhdf5") @@ -42,7 +46,9 @@ itk_add_test( ${ITK_TEST_OUTPUT_DIR}/t1_z+_byte_cor_2.mnc itkMINCImageIOTest4 DATA{Input/t1_z+_byte_cor.mnc} - ${ITK_TEST_OUTPUT_DIR}/t1_z+_byte_cor_2.mnc) + ${ITK_TEST_OUTPUT_DIR}/t1_z+_byte_cor_2.mnc + 0 + ) itk_add_test( NAME @@ -54,7 +60,9 @@ itk_add_test( ${ITK_TEST_OUTPUT_DIR}/HeadMRVolume.mnc itkMINCImageIOTest4 DATA{${ITK_DATA_ROOT}/Input/HeadMRVolume.mhd,HeadMRVolume.raw} - ${ITK_TEST_OUTPUT_DIR}/HeadMRVolume.mnc) + ${ITK_TEST_OUTPUT_DIR}/HeadMRVolume.mnc + -1 + ) # Test to make sure that inter-slice normalization was properly dealt with itk_add_test( @@ -67,7 +75,8 @@ itk_add_test( ${ITK_TEST_OUTPUT_DIR}/t1_z+_ubyte_yxz_nonorm_noParams.mnc itkMINCImageIOTest4 DATA{Input/t1_z+_ubyte_yxz_nonorm.mnc} - ${ITK_TEST_OUTPUT_DIR}/t1_z+_ubyte_yxz_nonorm_noParams.mnc) + ${ITK_TEST_OUTPUT_DIR}/t1_z+_ubyte_yxz_nonorm_noParams.mnc + 0) itk_add_test( NAME @@ -79,7 +88,8 @@ itk_add_test( ${ITK_TEST_OUTPUT_DIR}/t1_z+_float_yxz_nonorm_noParams.mnc itkMINCImageIOTest4 DATA{Input/t1_z+_float_yxz_nonorm.mnc} - ${ITK_TEST_OUTPUT_DIR}/t1_z+_float_yxz_nonorm_noParams.mnc) + ${ITK_TEST_OUTPUT_DIR}/t1_z+_float_yxz_nonorm_noParams.mnc + 0) itk_add_test( NAME @@ -130,6 +140,7 @@ itk_add_test( itkMINCImageIOTest4 DATA{Input/t1_z+_float_yxz_nonorm.mnc} ${ITK_TEST_OUTPUT_DIR}/t1_z+_float_yxz_nonorm.mnc + 0 # RAStoLPS 427621.7839 -8.195741583 72.45998819 @@ -146,6 +157,7 @@ itk_add_test( itkMINCImageIOTest4 DATA{Input/t1_z+_float_yxz_norm.mnc} ${ITK_TEST_OUTPUT_DIR}/t1_z+_float_yxz_norm.mnc + 0 # RAStoLPS 427621.7839 -8.195741583 72.45998819 @@ -162,47 +174,51 @@ itk_add_test( itkMINCImageIOTest4 DATA{Input/t1_z+_ubyte_yxz_nonorm.mnc} ${ITK_TEST_OUTPUT_DIR}/t1_z+_ubyte_yxz_nonorm.mnc + 0 # RAStoLPS 427621.7838 -8.195741583 72.45998819 -3.148534512) # multiple loops because of different numerical parameters +foreach(rastoeps 0;1) + foreach(type byte;short;ubyte) + foreach(axis cor;sag;trans) + foreach(plusMinus -;+) + set(imageName t1_z${plusMinus}_${type}_${axis}) + set(outImage ${ITK_TEST_OUTPUT_DIR}/${imageName}-${rastoeps}.mnc) -foreach(type byte;short;ubyte) - foreach(axis cor;sag;trans) - foreach(plusMinus -;+) - set(imageName t1_z${plusMinus}_${type}_${axis}) - set(outImage ${ITK_TEST_OUTPUT_DIR}/${imageName}.mnc) - - itk_add_test( - NAME - itkMINCImageIOTest-COM-${imageName} - COMMAND - ITKIOMINCTestDriver - --compare - DATA{Input/${imageName}.mnc} - ${outImage} - itkMINCImageIOTest4 - DATA{Input/${imageName}.mnc} - ${outImage} - 427620.3115 - -8.195582241 - 72.46002233 - -3.148594157) # this line is different + itk_add_test( + NAME + itkMINCImageIOTest-COM-${imageName}-${rastoeps} + COMMAND + ITKIOMINCTestDriver + --compare + DATA{Input/${imageName}.mnc} + ${outImage} + itkMINCImageIOTest4 + DATA{Input/${imageName}.mnc} + ${outImage} + ${rastoeps} # RAStoLPS + 427620.3115 + -8.195582241 + 72.46002233 + -3.148594157) # this line is different + endforeach() endforeach() endforeach() endforeach() +foreach(rastoeps 0;1) foreach(type double;float;long;ulong) foreach(axis cor;sag;trans) foreach(plusMinus -;+) set(imageName t1_z${plusMinus}_${type}_${axis}) - set(outImage ${ITK_TEST_OUTPUT_DIR}/${imageName}.mnc) + set(outImage ${ITK_TEST_OUTPUT_DIR}/${imageName}_${rastoeps}.mnc) itk_add_test( NAME - itkMINCImageIOTest-COM-${imageName} + itkMINCImageIOTest-COM-${imageName}-${rastoeps} COMMAND ITKIOMINCTestDriver --compare @@ -211,6 +227,7 @@ foreach(type double;float;long;ulong) itkMINCImageIOTest4 DATA{Input/${imageName}.mnc} ${outImage} + ${rastoeps} # RAStoLPS 427590.7631 -8.195995507 72.45943584 @@ -218,16 +235,18 @@ foreach(type double;float;long;ulong) endforeach() endforeach() endforeach() +endforeach() +foreach(rastoeps 0;1) foreach(type ushort) foreach(axis cor;sag;trans) foreach(plusMinus -;+) set(imageName t1_z${plusMinus}_${type}_${axis}) - set(outImage ${ITK_TEST_OUTPUT_DIR}/${imageName}.mnc) + set(outImage ${ITK_TEST_OUTPUT_DIR}/${imageName}_${rastoeps}.mnc) itk_add_test( NAME - itkMINCImageIOTest-COM-${imageName} + itkMINCImageIOTest-COM-${imageName}-${rastoeps} COMMAND ITKIOMINCTestDriver --compare @@ -236,6 +255,7 @@ foreach(type ushort) itkMINCImageIOTest4 DATA{Input/${imageName}.mnc} ${outImage} + ${rastoeps} # RAStoLPS 427590.7957 -8.195997123 72.45943721 @@ -243,3 +263,4 @@ foreach(type ushort) endforeach() endforeach() endforeach() +endforeach() diff --git a/Modules/IO/MINC/test/itkMINCImageIOTest4.cxx b/Modules/IO/MINC/test/itkMINCImageIOTest4.cxx index da1197e7b3d..abca2e3ff2e 100644 --- a/Modules/IO/MINC/test/itkMINCImageIOTest4.cxx +++ b/Modules/IO/MINC/test/itkMINCImageIOTest4.cxx @@ -23,11 +23,14 @@ #include "itkMINCImageIOFactory.h" #include "itkImageFileReader.h" #include "itkImageFileWriter.h" +#include "itksys/SystemTools.hxx" #include "itkImageMomentsCalculator.h" #include "itkStdStreamStateSave.h" #include "itkTestingMacros.h" +#include "itkMINCImageIOConfigurePrivate.h" + template int test_image_moments(const char * input_image, @@ -36,9 +39,17 @@ test_image_moments(const char * input_image, double mx, double my, double mz, - double epsilon) + double epsilon, + bool ras_to_lps) { - // itk::MINCImageIO::Pointer mincIO1 = itk::MINCImageIO::New(); + if (ras_to_lps) + { // need to flip expected moments to match flipped coordinate system + mx = -mx; + my = -my; + } + + itk::MINCImageIO::Pointer mincIO1 = itk::MINCImageIO::New(); + mincIO1->SetRAStoLPS(ras_to_lps); using ReaderType = itk::ImageFileReader; @@ -50,7 +61,8 @@ test_image_moments(const char * input_image, auto calculator = MomentsCalculatorType::New(); - // reader->SetImageIO( mincIO1 ); + if (itksys::SystemTools::StringEndsWith(input_image, ".mnc")) + reader->SetImageIO(mincIO1); reader->SetFileName(input_image); @@ -90,6 +102,11 @@ test_image_moments(const char * input_image, { auto writer = WriterType::New(); writer->SetFileName(output_image); + // writer should use default RAS to LPS flag, to satisfy comparison after + mincIO1->SetRAStoLPS(ras_to_lps); + if (itksys::SystemTools::StringEndsWith(output_image, ".mnc")) // HACK to enable use .mhd files + writer->SetImageIO(mincIO1); + writer->SetInput(reader->GetOutput()); writer->Update(); } @@ -105,16 +122,18 @@ itkMINCImageIOTest4(int argc, char * argv[]) // scope. itk::StdStreamStateSave coutState(std::cout); - if (argc < 3) + if (argc < 4) { std::cerr << "Missing parameters." << std::endl; std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv); - std::cerr << " inputfile outputfile [sum mx my mz ]" << std::endl; + std::cerr << " inputfile outputfile LPS> [sum mx my mz ]" << std::endl; return EXIT_FAILURE; } const char * input = argv[1]; const char * output = argv[2]; + int ras_to_lps_test = atoi(argv[3]); + bool ras_to_lps = ras_to_lps_test < 0 ? ITK_MINC_IO_RAS_TO_LPS : ras_to_lps_test == 1; double total = 0.0; double mx = 0.0; @@ -123,14 +142,14 @@ itkMINCImageIOTest4(int argc, char * argv[]) itk::MINCImageIOFactory::RegisterOneFactory(); - if (argc > 3) + if (argc > 4) { - if (argc == 7) + if (argc == 8) { - total = std::stod(argv[3]); - mx = std::stod(argv[4]); - my = std::stod(argv[5]); - mz = std::stod(argv[6]); + total = std::stod(argv[4]); + mx = std::stod(argv[5]); + my = std::stod(argv[6]); + mz = std::stod(argv[7]); } else { @@ -148,12 +167,13 @@ itkMINCImageIOTest4(int argc, char * argv[]) int ret = EXIT_SUCCESS; std::cout.precision(10); - if (test_image_moments>(input, nullptr, total, mx, my, mz, epsilon) != EXIT_SUCCESS) + if (test_image_moments>(input, nullptr, total, mx, my, mz, epsilon, ras_to_lps) != + EXIT_SUCCESS) { ret = EXIT_FAILURE; } // write out only float image - if (test_image_moments>(input, output, total, mx, my, mz, epsilon) != EXIT_SUCCESS) + if (test_image_moments>(input, output, total, mx, my, mz, epsilon, ras_to_lps) != EXIT_SUCCESS) { ret = EXIT_FAILURE; } diff --git a/Modules/IO/TransformMINC/include/itkMINCTransformAdapter.h b/Modules/IO/TransformMINC/include/itkMINCTransformAdapter.h index 428ad144a79..e55578718d6 100644 --- a/Modules/IO/TransformMINC/include/itkMINCTransformAdapter.h +++ b/Modules/IO/TransformMINC/include/itkMINCTransformAdapter.h @@ -39,6 +39,7 @@ namespace itk /** \class MINCTransformAdapter * \ingroup ITKIOTransformMINC * \brief ITK wrapper around MINC general transform functions, supports all the transformations that MINC XFM supports + * note, this wrapper does not take into account RAS to LPS conversion flag, and always assumes MINC convention * * \author Vladimir S. FONOV * Brain Imaging Center, Montreal Neurological Institute, McGill University, Montreal Canada 2012 diff --git a/Modules/IO/TransformMINC/include/itkMINCTransformIO.h b/Modules/IO/TransformMINC/include/itkMINCTransformIO.h index a0e945e0972..1579d436296 100644 --- a/Modules/IO/TransformMINC/include/itkMINCTransformIO.h +++ b/Modules/IO/TransformMINC/include/itkMINCTransformIO.h @@ -33,6 +33,8 @@ namespace itk /** \class MINCTransformIOTemplate * * \brief Read and write transforms in MINC format (.xfm). + * Takes into account PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS conversion + * flag to convert from internal MINC to ITK conventions and back, if enabled. * * \author Vladimir S. FONOV * Brain Imaging Center, Montreal Neurological Institute, McGill University, Montreal Canada 2012 @@ -80,6 +82,11 @@ class ITK_TEMPLATE_EXPORT MINCTransformIOTemplate : public TransformIOBaseTempla void Write() override; + /** Set to automatically convert from RAS PositiveCoordinateOrientation to LPS PositiveCoordinateOrientation*/ + itkSetMacro(RAStoLPS, bool); + itkGetConstMacro(RAStoLPS, bool); + itkBooleanMacro(RAStoLPS); + protected: MINCTransformIOTemplate(); ~MINCTransformIOTemplate() override; @@ -98,7 +105,8 @@ class ITK_TEMPLATE_EXPORT MINCTransformIOTemplate : public TransformIOBaseTempla int & serial); void - ReadOneTransform(VIO_General_transform * xfm); + ReadOneTransform(VIO_General_transform * xfm); + bool m_RAStoLPS; }; /** This helps to meet backward compatibility */ diff --git a/Modules/IO/TransformMINC/src/CMakeLists.txt b/Modules/IO/TransformMINC/src/CMakeLists.txt index bcced1a399c..9ff584a2a0e 100644 --- a/Modules/IO/TransformMINC/src/CMakeLists.txt +++ b/Modules/IO/TransformMINC/src/CMakeLists.txt @@ -1,3 +1,6 @@ +configure_file(itkMINCImageIOConfigurePrivate.h.in itkMINCImageIOConfigurePrivate.h @ONLY) +include_directories("${CMAKE_CURRENT_BINARY_DIR}") + set(ITKIOTransformMINC_SRC itkMINCTransformIO.cxx itkMINCTransformIOFactory.cxx) add_library(ITKIOTransformMINC ${ITK_LIBRARY_BUILD_TYPE} ${ITKIOTransformMINC_SRC}) diff --git a/Modules/IO/TransformMINC/src/itkMINCImageIOConfigurePrivate.h.in b/Modules/IO/TransformMINC/src/itkMINCImageIOConfigurePrivate.h.in new file mode 100644 index 00000000000..6034b0df09c --- /dev/null +++ b/Modules/IO/TransformMINC/src/itkMINCImageIOConfigurePrivate.h.in @@ -0,0 +1,27 @@ +/*========================================================================= + * + * 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. + * + *=========================================================================*/ +#ifndef itkMINCImageIOConfigurePrivate_h +#define itkMINCImageIOConfigurePrivate_h + +// This file is intended to hold preprocessor macros only used internally by +// IO/MINC module and should not be installed. + +// Convert from RAS PositiveCoordinateOrientation to LPS PositiveCoordinateOrientation by default +#cmakedefine01 ITK_MINC_IO_RAS_TO_LPS + +#endif // itkMINCImageIOConfigurePrivate_h diff --git a/Modules/IO/TransformMINC/src/itkMINCTransformIO.cxx b/Modules/IO/TransformMINC/src/itkMINCTransformIO.cxx index d87a1fe117c..ad39250f6e8 100644 --- a/Modules/IO/TransformMINC/src/itkMINCTransformIO.cxx +++ b/Modules/IO/TransformMINC/src/itkMINCTransformIO.cxx @@ -31,15 +31,20 @@ #include "itkVector.h" #include "itkDisplacementFieldTransform.h" #include "itkMetaDataObject.h" +#include "itkImageRegionIterator.h" +#include "vnl/vnl_vector_fixed.h" + +#include "itkMINCImageIOConfigurePrivate.h" + namespace itk { template MINCTransformIOTemplate::MINCTransformIOTemplate() -{ - m_XFM_initialized = false; -} + : m_XFM_initialized(false) + , m_RAStoLPS(ITK_MINC_IO_RAS_TO_LPS) +{} template MINCTransformIOTemplate::~MINCTransformIOTemplate() @@ -93,13 +98,38 @@ MINCTransformIOTemplate::ReadOneTransform(VIO_General_tran ParametersType parameterArray; parameterArray.SetSize(12); + + Matrix _affine_transform; + + _affine_transform.SetIdentity(); + + // MINC stores transforms in PositiveCoordinateOrientation RAS + // need to convert to PositiveCoordinateOrientation LPS for ITK + Matrix RAS_tofrom_LPS; + RAS_tofrom_LPS.SetIdentity(); + RAS_tofrom_LPS(0, 0) = -1.0; + RAS_tofrom_LPS(1, 1) = -1.0; + for (int j = 0; j < 3; ++j) { for (int i = 0; i < 3; ++i) { - parameterArray.SetElement(i + j * 3, Transform_elem(*lin, j, i)); + _affine_transform(i, j) = Transform_elem(*lin, j, i); } - parameterArray.SetElement(j + 9, Transform_elem(*lin, j, 3)); + // shifts + _affine_transform(3, j) = Transform_elem(*lin, j, 3); + } + + if (this->m_RAStoLPS) // flip RAS PositiveCoordinateOrientation to LPS PositiveCoordinateOrientation + _affine_transform = RAS_tofrom_LPS * _affine_transform * RAS_tofrom_LPS; + + for (int j = 0; j < 3; ++j) + { + for (int i = 0; i < 3; ++i) + { + parameterArray.SetElement(i + j * 3, _affine_transform(i, j)); + } + parameterArray.SetElement(j + 9, _affine_transform(3, j)); } if (xfm->inverse_flag) @@ -138,14 +168,40 @@ MINCTransformIOTemplate::ReadOneTransform(VIO_General_tran using DisplacementFieldTransformType = DisplacementFieldTransform; using GridImageType = typename DisplacementFieldTransformType::DisplacementFieldType; using MincReaderType = ImageFileReader; + using OutputPixelType = vnl_vector_fixed; + const OutputPixelType RAS_tofrom_LPS_vector = { -1, -1, 1 }; auto mincIO = MINCImageIO::New(); + mincIO->SetRAStoLPS(this->m_RAStoLPS); + auto reader = MincReaderType::New(); reader->SetImageIO(mincIO); reader->SetFileName(xfm->displacement_volume_file); reader->Update(); - typename GridImageType::Pointer grid = reader->GetOutput(); + typename GridImageType::Pointer LPSgrid = GridImageType::New(); + + if (this->m_RAStoLPS) + { + LPSgrid->CopyInformation(grid); + LPSgrid->SetRegions(grid->GetBufferedRegion()); + LPSgrid->Allocate(true); + + itk::MultiThreaderBase::Pointer mt = itk::MultiThreaderBase::New(); + mt->ParallelizeImageRegion<3>( + grid->GetBufferedRegion(), + [grid, LPSgrid, RAS_tofrom_LPS_vector](const typename GridImageType::RegionType & region) { + itk::Vector p; + itk::ImageRegionConstIterator iIt(grid, region); + itk::ImageRegionIterator oIt(LPSgrid, region); + for (; !iIt.IsAtEnd(); ++iIt, ++oIt) + { + p.SetVnlVector(element_product(iIt.Get().GetVnlVector(), RAS_tofrom_LPS_vector)); + oIt.Set(p); + } + }, + nullptr); + } TransformPointer transform; std::string transformTypeName = "DisplacementFieldTransform_"; @@ -155,11 +211,17 @@ MINCTransformIOTemplate::ReadOneTransform(VIO_General_tran auto * gridTransform = static_cast(transform.GetPointer()); if (xfm->inverse_flag) // TODO: invert grid transform? { - gridTransform->SetInverseDisplacementField(grid); + if (this->m_RAStoLPS) + gridTransform->SetInverseDisplacementField(LPSgrid); + else + gridTransform->SetInverseDisplacementField(grid); } else { - gridTransform->SetDisplacementField(grid); + if (this->m_RAStoLPS) + gridTransform->SetDisplacementField(LPSgrid); + else + gridTransform->SetDisplacementField(grid); } this->GetReadTransformList().push_back(transform); @@ -224,13 +286,36 @@ MINCTransformIOTemplate::WriteOneTransform(const int MatrixType matrix = matrixOffsetTransform->GetMatrix(); OffsetType offset = matrixOffsetTransform->GetOffset(); + // MINC stores everything in PositiveCoordinateOrientation RAS + // need to convert from PositiveCoordinateOrientation LPS + Matrix RAS_tofrom_LPS; + Matrix _affine_transform; + + RAS_tofrom_LPS.SetIdentity(); + RAS_tofrom_LPS(0, 0) = -1.0; + RAS_tofrom_LPS(1, 1) = -1.0; + _affine_transform.SetIdentity(); + + + for (int j = 0; j < 3; ++j) + { + for (int i = 0; i < 3; ++i) + { + _affine_transform(j, i) = matrix(j, i); + } + _affine_transform(3, j) = offset[j]; + } + + if (this->m_RAStoLPS) + _affine_transform = RAS_tofrom_LPS * _affine_transform * RAS_tofrom_LPS; + for (int j = 0; j < 3; ++j) { for (int i = 0; i < 3; ++i) { - Transform_elem(lin, j, i) = matrix(j, i); + Transform_elem(lin, j, i) = _affine_transform(j, i); } - Transform_elem(lin, j, 3) = offset[j]; + Transform_elem(lin, j, 3) = _affine_transform(3, j); } // add 4th normalization row (not stored) Transform_elem(lin, 3, 3) = 1.0; @@ -245,29 +330,65 @@ MINCTransformIOTemplate::WriteOneTransform(const int using DisplacementFieldTransformType = DisplacementFieldTransform; using GridImageType = typename DisplacementFieldTransformType::DisplacementFieldType; using MincWriterType = ImageFileWriter; + typename GridImageType::Pointer grid = GridImageType::New(); + using OutputPixelType = vnl_vector_fixed; + const OutputPixelType RAS_tofrom_LPS_vector = { -1, -1, 1 }; auto * _grid_transform = static_cast(const_cast(curTransform)); char tmp[1024]; snprintf(tmp, sizeof(tmp), "%s_grid_%d.mnc", xfm_file_base, serial); ++serial; auto mincIO = MINCImageIO::New(); + // writer should follow the same notation, in case we manually switched the coordinate system + mincIO->SetRAStoLPS(this->m_RAStoLPS); + auto writer = MincWriterType::New(); + writer->SetImageIO(mincIO); writer->SetFileName(tmp); if (_grid_transform->GetDisplacementField()) { - writer->SetInput(_grid_transform->GetModifiableDisplacementField()); + grid = _grid_transform->GetModifiableDisplacementField(); } else if (_grid_transform->GetInverseDisplacementField()) { - writer->SetInput(_grid_transform->GetModifiableInverseDisplacementField()); + grid = _grid_transform->GetModifiableInverseDisplacementField(); _inverse_grid = true; } else { itkExceptionMacro("Trying to write-out displacement transform without displacement field"); } + + if (this->m_RAStoLPS) + { + typename GridImageType::Pointer ras_grid = GridImageType::New(); // flipped grid for RAS->LPS conversion + ras_grid->CopyInformation(grid); + ras_grid->SetRegions(grid->GetBufferedRegion()); + ras_grid->Allocate(true); + + itk::MultiThreaderBase::Pointer mt = itk::MultiThreaderBase::New(); + mt->ParallelizeImageRegion<3>( + grid->GetBufferedRegion(), + [grid, ras_grid, RAS_tofrom_LPS_vector](const typename GridImageType::RegionType & region) { + itk::Vector p; + itk::ImageRegionConstIterator iIt(grid, region); + itk::ImageRegionIterator oIt(ras_grid, region); + for (; !iIt.IsAtEnd(); ++iIt, ++oIt) + { + p.SetVnlVector(element_product(iIt.Get().GetVnlVector(), RAS_tofrom_LPS_vector)); + oIt.Set(p); + } + }, + nullptr); + + writer->SetInput(ras_grid); + } + else + { + writer->SetInput(grid); + } writer->Update(); xfm.emplace_back(); diff --git a/Modules/IO/TransformMINC/test/itkIOTransformMINCTest.cxx b/Modules/IO/TransformMINC/test/itkIOTransformMINCTest.cxx index dbc92baf3bf..64c7a946472 100644 --- a/Modules/IO/TransformMINC/test/itkIOTransformMINCTest.cxx +++ b/Modules/IO/TransformMINC/test/itkIOTransformMINCTest.cxx @@ -20,6 +20,8 @@ #include #include "itkMINCTransformIOFactory.h" #include "itkMINCImageIOFactory.h" +#include "itkMINCTransformIO.h" +#include "itkMINCImageIO.h" #include "itkTransformFileWriter.h" #include "itkTransformFileReader.h" #include "itkAffineTransform.h" @@ -34,6 +36,7 @@ static constexpr int point_counter = 1000; using TransformFileReader = itk::TransformFileReaderTemplate; using TransformFileWriter = itk::TransformFileWriterTemplate; +using MINCTransformIO = itk::MINCTransformIO; template void @@ -57,8 +60,10 @@ RandomPoint(vnl_random & randgen, itk::Point & pix, double _max = itk::Num static int -check_linear(const char * linear_transform) +check_linear(const char * linear_transform, bool ras_to_lps) { + std::cout << "check_linear, ras_to_lps=" << ras_to_lps << std::endl << std::endl; + int testStatus = EXIT_SUCCESS; using AffineTransformType = itk::AffineTransform; @@ -88,32 +93,34 @@ check_linear(const char * linear_transform) TransformFileWriter::Pointer writer; TransformFileReader::Pointer reader; + MINCTransformIO::Pointer mincIO = MINCTransformIO::New(); + mincIO->SetRAStoLPS(ras_to_lps); reader = TransformFileReader::New(); writer = TransformFileWriter::New(); + writer->SetTransformIO(mincIO); + reader->SetTransformIO(mincIO); + writer->AddTransform(affine); writer->SetFileName(linear_transform); reader->SetFileName(linear_transform); - // Testing writing std::cout << "Testing write : "; - affine->Print(std::cout); + // affine->Print(std::cout); ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); - std::cout << "Testing read : " << std::endl; ITK_TRY_EXPECT_NO_EXCEPTION(reader->Update()); - TransformFileReader::TransformListType list = *reader->GetTransformList(); ITK_TEST_EXPECT_EQUAL("AffineTransform_double_3_3", list.front()->GetTransformTypeAsString()); AffineTransformType::Pointer affine2 = static_cast(list.front().GetPointer()); std::cout << "Read transform : " << std::endl; - affine2->Print(std::cout); + // affine2->Print(std::cout); vnl_random randgen(12345678); @@ -133,7 +140,8 @@ check_linear(const char * linear_transform) double expectedDiff = 0.0; if (!itk::Math::FloatAlmostEqual(expectedDiff, (v1 - v2).GetSquaredNorm(), 10, tolerance)) { - std::cerr << "Test failed!" << std::endl; + std::cerr << "Test failed!" + << " In " __FILE__ ", line " << __LINE__ << std::endl; std::cerr << "Error in pixel value at point ( " << pnt << ")" << std::endl; std::cerr << "Expected value " << v1; std::cerr << " differs from " << v2; @@ -147,8 +155,9 @@ check_linear(const char * linear_transform) } static int -check_nonlinear_double(const char * nonlinear_transform) +check_nonlinear_double(const char * nonlinear_transform, bool ras_to_lps) { + std::cout << "check_nonlinear_double, ras_to_lps=" << ras_to_lps << std::endl << std::endl; int testStatus = EXIT_SUCCESS; const double tolerance = 1e-5; @@ -198,13 +207,18 @@ check_nonlinear_double(const char * nonlinear_transform) disp->SetDisplacementField(field); - disp->Print(std::cout); + // disp->Print(std::cout); TransformFileWriter::Pointer nlwriter; TransformFileReader::Pointer nlreader; + MINCTransformIO::Pointer mincIO = MINCTransformIO::New(); + mincIO->SetRAStoLPS(ras_to_lps); nlreader = TransformFileReader::New(); nlwriter = TransformFileWriter::New(); + nlreader->SetTransformIO(mincIO); + nlwriter->SetTransformIO(mincIO); + nlwriter->AddTransform(disp); nlwriter->SetFileName(nonlinear_transform); nlreader->SetFileName(nonlinear_transform); @@ -220,7 +234,6 @@ check_nonlinear_double(const char * nonlinear_transform) ITK_TRY_EXPECT_NO_EXCEPTION(nlreader->Update()); - std::cout << "[PASSED]" << std::endl; std::cout << "Comparing of non linear transform (double) : " << std::endl; @@ -240,7 +253,8 @@ check_nonlinear_double(const char * nonlinear_transform) { if (it.Value() != it2.Value()) { - std::cerr << "Test failed!" << std::endl; + std::cerr << "Test failed!" + << " In " __FILE__ ", line " << __LINE__ << std::endl; std::cerr << "Error in pixel value" << std::endl; std::cerr << "Expected value " << it.Value(); std::cerr << " differs from " << it2.Value() << std::endl; @@ -255,7 +269,8 @@ check_nonlinear_double(const char * nonlinear_transform) double expectedDiff = 0.0; if (!itk::Math::FloatAlmostEqual(expectedDiff, (it.Value() - it2.Value()).GetSquaredNorm(), 10, tolerance)) { - std::cerr << "Test failed!" << std::endl; + std::cerr << "Test failed!" + << " In " __FILE__ ", line " << __LINE__ << std::endl; std::cerr << "Error in pixel value" << std::endl; std::cerr << "Expected value " << it.Value(); std::cerr << " differs from " << it2.Value(); @@ -270,14 +285,18 @@ check_nonlinear_double(const char * nonlinear_transform) static int -check_nonlinear_float(const char * nonlinear_transform) +check_nonlinear_float(const char * nonlinear_transform, bool ras_to_lps) { + + std::cout << "check_nonlinear_float, ras_to_lps=" << ras_to_lps << std::endl << std::endl; + int testStatus = EXIT_SUCCESS; double tolerance = 1e-5; using TransformFileWriterFloat = itk::TransformFileWriterTemplate; using TransformFileReaderFloat = itk::TransformFileReaderTemplate; + using MINCTransformIOFloat = itk::MINCTransformIOTemplate; using DisplacementFieldTransform = itk::DisplacementFieldTransform; using DisplacementFieldType = DisplacementFieldTransform::DisplacementFieldType; @@ -324,13 +343,19 @@ check_nonlinear_float(const char * nonlinear_transform) disp->SetDisplacementField(field); - disp->Print(std::cout); + // disp->Print(std::cout); TransformFileWriterFloat::Pointer nlwriter; TransformFileReaderFloat::Pointer nlreader; nlreader = TransformFileReaderFloat::New(); nlwriter = TransformFileWriterFloat::New(); + + MINCTransformIOFloat::Pointer mincIO = MINCTransformIOFloat::New(); + mincIO->SetRAStoLPS(ras_to_lps); + nlreader->SetTransformIO(mincIO); + nlwriter->SetTransformIO(mincIO); + nlwriter->AddTransform(disp); nlwriter->SetFileName(nonlinear_transform); nlreader->SetFileName(nonlinear_transform); @@ -366,7 +391,8 @@ check_nonlinear_float(const char * nonlinear_transform) { if (it.Value() != it2.Value()) { - std::cerr << "Test failed!" << std::endl; + std::cerr << "Test failed!" + << " In " __FILE__ ", line " << __LINE__ << std::endl; std::cerr << "Error in pixel value" << std::endl; std::cerr << "Expected value " << it.Value(); std::cerr << " differs from " << it2.Value() << std::endl; @@ -381,7 +407,8 @@ check_nonlinear_float(const char * nonlinear_transform) double expectedDiff = 0.0; if (!itk::Math::FloatAlmostEqual(expectedDiff, (it.Value() - it2.Value()).GetSquaredNorm(), 10, tolerance)) { - std::cerr << "Test failed!" << std::endl; + std::cerr << "Test failed!" + << " In " __FILE__ ", line " << __LINE__ << std::endl; std::cerr << "Error in pixel value" << std::endl; std::cerr << "Expected value " << it.Value(); std::cerr << " differs from " << it2.Value(); @@ -395,8 +422,11 @@ check_nonlinear_float(const char * nonlinear_transform) } static int -secondTest() +secondTest(bool ras_to_lps) { + std::cout << "secondTest, ras_to_lps=" << ras_to_lps << std::endl << std::endl; + + // rotation of 30 degrees around X axis std::filebuf fb; fb.open("Rotation.xfm", std::ios::out); std::ostream os(&fb); @@ -410,24 +440,25 @@ secondTest() TransformFileReader::Pointer reader; reader = TransformFileReader::New(); + MINCTransformIO::Pointer mincIO = MINCTransformIO::New(); + mincIO->SetRAStoLPS(ras_to_lps); + reader->SetTransformIO(mincIO); reader->SetFileName("Rotation.xfm"); - reader->Update(); const TransformFileReader::TransformListType * list = reader->GetTransformList(); - auto lit = list->begin(); - while (lit != list->end()) - { - (*lit)->Print(std::cout); - lit++; - } + + ITK_TEST_EXPECT_EQUAL(1, list->size()); + return EXIT_SUCCESS; } static int -check_composite(const char * transform_file) +check_composite(const char * transform_file, bool ras_to_lps) { + std::cout << "check_composite, ras_to_lps=" << ras_to_lps << std::endl << std::endl; + int testStatus = EXIT_SUCCESS; using AffineTransformType = itk::AffineTransform; @@ -464,23 +495,27 @@ check_composite(const char * transform_file) TransformFileWriter::Pointer writer; TransformFileReader::Pointer reader; + MINCTransformIO::Pointer mincIO = MINCTransformIO::New(); + mincIO->SetRAStoLPS(ras_to_lps); + reader = TransformFileReader::New(); writer = TransformFileWriter::New(); + reader->SetTransformIO(mincIO); + writer->SetTransformIO(mincIO); + writer->AddTransform(compositeTransform); writer->SetFileName(transform_file); reader->SetFileName(transform_file); - compositeTransform->Print(std::cout); + // compositeTransform->Print(std::cout); ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); - std::cout << "Testing read : " << std::endl; ITK_TRY_EXPECT_NO_EXCEPTION(reader->Update()); - TransformFileReader::TransformListType list = *reader->GetTransformList(); // MINC XFM internally collapeses two concatenated linear transforms into one @@ -489,7 +524,7 @@ check_composite(const char * transform_file) AffineTransformType::Pointer affine_xfm = static_cast(list.front().GetPointer()); std::cout << "Read transform : " << std::endl; - affine_xfm->Print(std::cout); + // affine_xfm->(std::cout); vnl_random randgen(12345678); @@ -509,7 +544,8 @@ check_composite(const char * transform_file) double expectedDiff = 0.0; if (!itk::Math::FloatAlmostEqual(expectedDiff, (v1 - v2).GetSquaredNorm(), 10, tolerance)) { - std::cerr << "Test failed!" << std::endl; + std::cerr << "Test failed!" + << " In " __FILE__ ", line " << __LINE__ << std::endl; std::cerr << "Error in pixel value at point ( " << pnt << ")" << std::endl; std::cerr << "Expected value " << v1; std::cerr << " differs from " << v2; @@ -523,8 +559,10 @@ check_composite(const char * transform_file) } static int -check_composite2(const char * transform_file, const char * transform_grid_file) +check_composite2(const char * transform_file, const char * transform_grid_file, bool ras_to_lps) { + std::cout << "check_composite2, ras_to_lps=" << ras_to_lps << std::endl << std::endl; + int testStatus = EXIT_SUCCESS; const double tolerance = 1e-5; @@ -532,7 +570,11 @@ check_composite2(const char * transform_file, const char * transform_grid_file) std::filebuf fb; fb.open(transform_file, std::ios::out); std::ostream os(&fb); - std::cout << "Testing reading of composite transform" << std::endl << std::endl << std::endl << std::endl; + std::cout << "Testing reading of composite transform, ras_to_lps=" << ras_to_lps << std::endl + << std::endl + << std::endl + << std::endl; + // create concatenation of two transforms: // nonlinear grid with shift by 1 // rotation by 45 deg @@ -555,7 +597,7 @@ check_composite2(const char * transform_file, const char * transform_grid_file) auto field = DisplacementFieldType::New(); - // create zero displacement field + // create displacement field in MINC space (always) DisplacementFieldType::SizeType imageSize3D = { { 10, 10, 10 } }; DisplacementFieldType::IndexType startIndex3D = { { 0, 0, 0 } }; @@ -579,27 +621,35 @@ check_composite2(const char * transform_file, const char * transform_grid_file) using MincWriterType = itk::ImageFileWriter; auto writer = MincWriterType::New(); + auto mincImageIO = itk::MINCImageIO::New(); + mincImageIO->SetRAStoLPS(false); + writer->SetImageIO(mincImageIO); + + // expecting .mnc here writer->SetFileName(transform_grid_file); - writer->SetInput(field); - ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); + // NOW, read back with RAS to LPS conversion + TransformFileReader::Pointer reader; using CompositeTransformType = itk::CompositeTransform; reader = TransformFileReader::New(); + + MINCTransformIO::Pointer mincIO = MINCTransformIO::New(); + mincIO->SetRAStoLPS(ras_to_lps); + reader->SetTransformIO(mincIO); reader->SetFileName(transform_file); ITK_TRY_EXPECT_NO_EXCEPTION(reader->Update()); - const TransformFileReader::TransformListType * list = reader->GetTransformList(); ITK_TEST_EXPECT_EQUAL(2, list->size()); - std::cout << "Reading back transforms : " << list->size() << std::endl << std::endl; + std::cout << "Reading back transforms : " << list->size() << std::endl; using CompositeTransformType = itk::CompositeTransform; using TransformType = itk::Transform; @@ -607,35 +657,45 @@ check_composite2(const char * transform_file, const char * transform_grid_file) auto _xfm = CompositeTransformType::New(); for (const auto & it : *list) { - it->Print(std::cout); - std::cout << std::endl; _xfm->AddTransform(static_cast(it.GetPointer())); } std::cout << std::endl; std::cout << std::endl; std::cout << std::endl; - _xfm->Print(std::cout); - std::cout << "Testing that transformations is expected ..." << std::endl; CompositeTransformType::InputPointType pnt; CompositeTransformType::OutputPointType v, v2; - pnt[0] = 1.0; + if (ras_to_lps) + { + pnt[0] = -1.0; + } + else + { + pnt[0] = 1.0; + } pnt[1] = pnt[2] = 0.0; // expected transform: shift by 1 , rotate by 45 deg - v2[0] = v2[1] = sqrt(2); + if (ras_to_lps) + { + v2[0] = v2[1] = -sqrt(2); + } + else + { + v2[0] = v2[1] = sqrt(2); + } v2[2] = 0.0; v = _xfm->TransformPoint(pnt); - std::cout << "In:" << pnt << " Out:" << v << std::endl; double expectedDiff = 0.0; if (!itk::Math::FloatAlmostEqual(expectedDiff, (v - v2).GetSquaredNorm(), 10, tolerance)) { - std::cerr << "Test failed!" << std::endl; + std::cerr << "Test failed!" + << " In " __FILE__ ", line " << __LINE__ << std::endl; std::cerr << "Error in pixel value at point ( " << pnt << ")" << std::endl; std::cerr << "Expected value " << v2; std::cerr << " differs from " << v; @@ -654,7 +714,6 @@ itkIOTransformMINCTest(int argc, char * argv[]) { itk::ObjectFactoryBase::RegisterFactory(itk::MINCImageIOFactory::New(), itk::ObjectFactoryBase::InsertionPositionEnum::INSERT_AT_FRONT); - if (argc > 1) { itksys::SystemTools::ChangeDirectory(argv[1]); @@ -662,14 +721,23 @@ itkIOTransformMINCTest(int argc, char * argv[]) itk::TransformFactory>::RegisterTransform(); itk::TransformFactory>::RegisterTransform(); - int result1 = check_linear("itkIOTransformMINCTestTransformLinear.xfm"); - int result2 = check_nonlinear_double("itkIOTransformMINCTestTransformNonLinear.xfm"); - int result3 = check_nonlinear_float("itkIOTransformMINCTestTransformNonLinear_float.xfm"); - int result4 = secondTest(); - int result5 = check_composite("itkIOTransformMINCTestTransformComposite.xfm"); - int result6 = check_composite2("itkIOTransformMINCTestTransformComposite2.xfm", - "itkIOTransformMINCTestTransformComposite2_grid_0.mnc"); + bool result = true; + for (bool ras_to_lps : { false, true }) + { + std::cout << "RAS to LPS=" << ras_to_lps << std::endl << std::endl << std::endl; + + bool result1 = check_linear("itkIOTransformMINCTestTransformLinear.xfm", ras_to_lps) == EXIT_SUCCESS; + bool result2 = check_nonlinear_double("itkIOTransformMINCTestTransformNonLinear.xfm", ras_to_lps) == EXIT_SUCCESS; + bool result3 = + check_nonlinear_float("itkIOTransformMINCTestTransformNonLinear_float.xfm", ras_to_lps) == EXIT_SUCCESS; + bool result4 = secondTest(ras_to_lps) == EXIT_SUCCESS; + bool result5 = check_composite("itkIOTransformMINCTestTransformComposite.xfm", ras_to_lps) == EXIT_SUCCESS; + bool result6 = check_composite2("itkIOTransformMINCTestTransformComposite2.xfm", + "itkIOTransformMINCTestTransformComposite2_grid_0.mnc", + ras_to_lps) == EXIT_SUCCESS; + + result = result && result1 && result2 && result3 && result4 && result5 && result6; + } - return !(result1 == EXIT_SUCCESS && result2 == EXIT_SUCCESS && result3 == EXIT_SUCCESS && result4 == EXIT_SUCCESS && - result5 == EXIT_SUCCESS && result6 == EXIT_SUCCESS); + return result ? EXIT_SUCCESS : EXIT_FAILURE; } diff --git a/Modules/IO/TransformMINC/test/itkMINCTransformAdapterTest.cxx b/Modules/IO/TransformMINC/test/itkMINCTransformAdapterTest.cxx index 799fcc0d74a..493dd8ce04f 100644 --- a/Modules/IO/TransformMINC/test/itkMINCTransformAdapterTest.cxx +++ b/Modules/IO/TransformMINC/test/itkMINCTransformAdapterTest.cxx @@ -27,6 +27,7 @@ #include "itkDisplacementFieldTransform.h" #include "itkIOTestHelper.h" #include "itkMINCTransformAdapter.h" +#include "itkMINCTransformIO.h" #include "itkMath.h" #include "itkTestingMacros.h" @@ -82,8 +83,12 @@ compare_linear(const char * linear_transform) affine->Scale(1.2); itk::TransformFileWriter::Pointer writer; + itk::MINCTransformIO::Pointer mincIO = itk::MINCTransformIO::New(); + // MINC standard is always LPS + mincIO->SetRAStoLPS(false); writer = itk::TransformFileWriter::New(); + writer->SetTransformIO(mincIO); writer->AddTransform(affine); writer->SetFileName(linear_transform); @@ -112,7 +117,7 @@ compare_linear(const char * linear_transform) if ((v1 - v2).GetSquaredNorm() > tolerance) { - std::cout << "Original Pixel (" << v1 << ") doesn't match read-in Pixel (" << v2 << " ) " + std::cout << "Original Coordinate (" << v1 << ") doesn't match read-in Coordinate (" << v2 << " ) " << " in " << linear_transform << " at " << pnt << std::endl; return EXIT_FAILURE; } @@ -179,8 +184,12 @@ compare_nonlinear_double(const char * nonlinear_transform) disp->SetDisplacementField(field); itk::TransformFileWriter::Pointer nlwriter; + itk::MINCTransformIO::Pointer mincIO = itk::MINCTransformIO::New(); + // MINC standard is always LPS + mincIO->SetRAStoLPS(false); nlwriter = itk::TransformFileWriter::New(); + nlwriter->SetTransformIO(mincIO); nlwriter->AddTransform(disp); nlwriter->SetFileName(nonlinear_transform);