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..b3b81726736 100644 --- a/Modules/IO/MINC/src/CMakeLists.txt +++ b/Modules/IO/MINC/src/CMakeLists.txt @@ -1,3 +1,4 @@ + 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 19c5fe38f1c..a13ccb29d8d 100644 --- a/Modules/IO/MINC/src/itkMINCImageIO.cxx +++ b/Modules/IO/MINC/src/itkMINCImageIO.cxx @@ -29,6 +29,7 @@ #include // For unique_ptr. + extern "C" { void @@ -235,6 +236,7 @@ MINCImageIO::CloseVolume() MINCImageIO::MINCImageIO() : m_MINCPImpl(std::make_unique()) + , m_RAStoLPS(false) { for (auto & dimensionIndex : m_MINCPImpl->m_DimensionIndices) { @@ -280,6 +282,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 @@ -446,14 +449,20 @@ MINCImageIO::ReadImageInformation() this->SetNumberOfDimensions(spatial_dimension_count); - int numberOfComponents = 1; - int usable_dimensions = 0; + int numberOfComponents = 1; + unsigned int usableDimensions = 0; + + auto dir_cos = Matrix::GetIdentity(); - Matrix dir_cos{}; - dir_cos.SetIdentity(); + // Conversion matrix for MINC PositiveCoordinateOrientation RAS (LeftToRight, PosteriorToAnterior, InferiorToSuperior) + // to ITK PositiveCoordinateOrientation LPS (RightToLeft, AnteriorToPosterior, InferiorToSuperior) + auto RAStofromLPS = Matrix::GetIdentity(); + RAStofromLPS(0, 0) = -1.0; + RAStofromLPS(1, 1) = -1.0; + std::vector dir_cos_temp(3); Vector origin{}; - Vector o_origin{}; + Vector oOrigin{}; // minc api uses inverse order of dimensions , fastest varying are last Vector sep; @@ -463,19 +472,19 @@ MINCImageIO::ReadImageInformation() { // MINC2: bad design! // micopy_dimension(hdim[m_MINCPImpl->m_DimensionIndices[i]],&apparent_dimension_order[usable_dimensions]); - m_MINCPImpl->m_MincApparentDims[usable_dimensions] = + m_MINCPImpl->m_MincApparentDims[usableDimensions] = m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[i]]; // always use positive - miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usable_dimensions], MI_POSITIVE); + miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_POSITIVE); misize_t _sz; - miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usable_dimensions], &_sz); + miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_sz); double _sep; - miget_dimension_separation(m_MINCPImpl->m_MincApparentDims[usable_dimensions], MI_ORDER_APPARENT, &_sep); + miget_dimension_separation(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_ORDER_APPARENT, &_sep); std::vector _dir(3); - miget_dimension_cosines(m_MINCPImpl->m_MincApparentDims[usable_dimensions], &_dir[0]); + miget_dimension_cosines(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_dir[0]); double _start; - miget_dimension_start(m_MINCPImpl->m_MincApparentDims[usable_dimensions], MI_ORDER_APPARENT, &_start); + miget_dimension_start(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_ORDER_APPARENT, &_start); for (int j = 0; j < 3; ++j) { @@ -486,52 +495,61 @@ MINCImageIO::ReadImageInformation() sep[i - 1] = _sep; this->SetDimensions(i - 1, static_cast(_sz)); - this->SetDirection(i - 1, _dir); this->SetSpacing(i - 1, _sep); - ++usable_dimensions; + ++usableDimensions; + } + } + + + // Transform MINC PositiveCoordinateOrientation RAS coordinates to + // internal ITK PositiveCoordinateOrientation LPS Coordinates + if (this->m_RAStoLPS) + dir_cos = RAStofromLPS * dir_cos; + + // Transform origin coordinates + oOrigin = dir_cos * origin; + + for (int i = 0; i < spatial_dimension_count; ++i) + { + this->SetOrigin(i, oOrigin[i]); + for (unsigned int j = 0; j < 3; j++) + { + dir_cos_temp[j] = dir_cos[j][i]; } + this->SetDirection(i, dir_cos_temp); } if (m_MINCPImpl->m_DimensionIndices[0] != -1) // have vector dimension { // micopy_dimension(m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[0]],&apparent_dimension_order[usable_dimensions]); - m_MINCPImpl->m_MincApparentDims[usable_dimensions] = - m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[0]]; + m_MINCPImpl->m_MincApparentDims[usableDimensions] = m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[0]]; // always use positive, vector dimension does not supposed to have notion of positive step size, so leaving as is // miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usable_dimensions],MI_POSITIVE); misize_t _sz; - miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usable_dimensions], &_sz); + miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_sz); numberOfComponents = _sz; - ++usable_dimensions; + ++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[usable_dimensions] = - m_MINCPImpl->m_MincFileDims[m_MINCPImpl->m_DimensionIndices[4]]; + 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[usable_dimensions], MI_POSITIVE); + miset_dimension_apparent_voxel_order(m_MINCPImpl->m_MincApparentDims[usableDimensions], MI_POSITIVE); misize_t _sz; - miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usable_dimensions], &_sz); + miget_dimension_size(m_MINCPImpl->m_MincApparentDims[usableDimensions], &_sz); numberOfComponents = _sz; - ++usable_dimensions; + ++usableDimensions; } // Set apparent dimension order to the MINC2 api - if (miset_apparent_dimension_order(m_MINCPImpl->m_Volume, usable_dimensions, m_MINCPImpl->m_MincApparentDims) < 0) + if (miset_apparent_dimension_order(m_MINCPImpl->m_Volume, usableDimensions, m_MINCPImpl->m_MincApparentDims) < 0) { itkExceptionMacro(" Can't set apparent dimension order!"); } - o_origin = dir_cos * origin; - - for (int i = 0; i < spatial_dimension_count; ++i) - { - this->SetOrigin(i, o_origin[i]); - } - miclass_t volume_data_class; if (miget_data_class(m_MINCPImpl->m_Volume, &volume_data_class) < 0) @@ -660,6 +678,9 @@ MINCImageIO::ReadImageInformation() const 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; @@ -944,22 +965,34 @@ MINCImageIO::WriteImageInformation() } // allocating dimensions - vnl_matrix dircosmatrix(nDims, nDims); - dircosmatrix.set_identity(); + vnl_matrix directionCosineMatrix(nDims, nDims); + directionCosineMatrix.set_identity(); vnl_vector origin(nDims); + // MINC stores direction cosines in PositiveCoordinateOrientation RAS + // need to convert to PositiveCoordinateOrientation LPS for ITK + vnl_matrix RAS_tofrom_LPS(nDims, nDims); + RAS_tofrom_LPS.set_identity(); + RAS_tofrom_LPS(0, 0) = -1.0; + RAS_tofrom_LPS(1, 1) = -1.0; + for (unsigned int i = 0; i < nDims; ++i) { for (unsigned int j = 0; j < nDims; ++j) { - dircosmatrix[i][j] = this->GetDirection(i)[j]; + directionCosineMatrix[i][j] = this->GetDirection(i)[j]; } origin[i] = this->GetOrigin(i); } - const vnl_matrix inverseDirectionCosines{ vnl_matrix_inverse(dircosmatrix).as_matrix() }; + const vnl_matrix inverseDirectionCosines{ vnl_matrix_inverse(directionCosineMatrix).as_matrix() }; origin *= inverseDirectionCosines; // transform to minc convention + + // Convert ITK direction cosines from PositiveCoordinateOrientation LPS to PositiveCoordinateOrientation RAS + if (this->m_RAStoLPS) + directionCosineMatrix *= RAS_tofrom_LPS; + for (unsigned int i = 0; i < nDims; ++i) { const unsigned int j = i + (nComp > 1 ? 1 : 0); @@ -968,7 +1001,7 @@ MINCImageIO::WriteImageInformation() { if (k < nDims) { - dir_cos[k] = dircosmatrix[i][k]; + dir_cos[k] = directionCosineMatrix[i][k]; } else { @@ -988,27 +1021,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; @@ -1016,10 +1043,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: @@ -1059,7 +1086,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 9c536baee36..be8f92310ac 100644 --- a/Modules/IO/MINC/test/itkMINCImageIOTest4.cxx +++ b/Modules/IO/MINC/test/itkMINCImageIOTest4.cxx @@ -23,11 +23,11 @@ #include "itkMINCImageIOFactory.h" #include "itkImageFileReader.h" #include "itkImageFileWriter.h" +#include "itksys/SystemTools.hxx" #include "itkImageMomentsCalculator.h" #include "itkStdStreamStateSave.h" #include "itkTestingMacros.h" - template int test_image_moments(const char * input_image, @@ -36,9 +36,17 @@ test_image_moments(const char * input_image, double mx, double my, double mz, - double epsilon) + double epsilon, + bool RAStofromLPS) { - // itk::MINCImageIO::Pointer mincIO1 = itk::MINCImageIO::New(); + if (RAStofromLPS) + { // need to flip expected moments to match flipped coordinate system + mx = -mx; + my = -my; + } + + itk::MINCImageIO::Pointer mincIO1 = itk::MINCImageIO::New(); + mincIO1->SetRAStoLPS(RAStofromLPS); using ReaderType = itk::ImageFileReader; @@ -50,7 +58,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 +99,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(RAStofromLPS); + if (itksys::SystemTools::StringEndsWith(output_image, ".mnc")) // HACK to enable use .mhd files + writer->SetImageIO(mincIO1); + writer->SetInput(reader->GetOutput()); writer->Update(); } @@ -105,16 +119,18 @@ itkMINCImageIOTest4(int argc, char * argv[]) // scope. const 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 RAStofromLPSTest = atoi(argv[3]); + bool RAStofromLPS = RAStofromLPSTest < 0 ? false : RAStofromLPSTest == 1; double total = 0.0; double mx = 0.0; @@ -123,14 +139,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 +164,14 @@ 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, RAStofromLPS) != + 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, RAStofromLPS) != + EXIT_SUCCESS) { ret = EXIT_FAILURE; } diff --git a/Modules/IO/TransformMINC/include/itkMINCTransformAdapter.h b/Modules/IO/TransformMINC/include/itkMINCTransformAdapter.h index cc0bfd110d0..8272203e181 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..69ea72267be 100644 --- a/Modules/IO/TransformMINC/src/itkMINCTransformIO.cxx +++ b/Modules/IO/TransformMINC/src/itkMINCTransformIO.cxx @@ -31,15 +31,18 @@ #include "itkVector.h" #include "itkDisplacementFieldTransform.h" #include "itkMetaDataObject.h" +#include "itkImageRegionIterator.h" +#include "vnl/vnl_vector_fixed.h" + namespace itk { template MINCTransformIOTemplate::MINCTransformIOTemplate() -{ - m_XFM_initialized = false; -} + : m_XFM_initialized(false) + , m_RAStoLPS(false) +{} template MINCTransformIOTemplate::~MINCTransformIOTemplate() @@ -93,13 +96,38 @@ MINCTransformIOTemplate::ReadOneTransform(VIO_General_tran ParametersType parameterArray; parameterArray.SetSize(12); + + Matrix affineMatrix; + + affineMatrix.SetIdentity(); + + // MINC stores transforms in PositiveCoordinateOrientation RAS + // need to convert to PositiveCoordinateOrientation LPS for ITK + Matrix RAStofromLPS; + RAStofromLPS.SetIdentity(); + RAStofromLPS(0, 0) = -1.0; + RAStofromLPS(1, 1) = -1.0; + + for (int j = 0; j < 3; ++j) + { + for (int i = 0; i < 3; ++i) + { + affineMatrix(i, j) = Transform_elem(*lin, j, i); + } + // shifts + affineMatrix(3, j) = Transform_elem(*lin, j, 3); + } + + if (this->m_RAStoLPS) // flip RAS PositiveCoordinateOrientation to LPS PositiveCoordinateOrientation + affineMatrix = RAStofromLPS * affineMatrix * RAStofromLPS; + for (int j = 0; j < 3; ++j) { for (int i = 0; i < 3; ++i) { - parameterArray.SetElement(i + j * 3, Transform_elem(*lin, j, i)); + parameterArray.SetElement(i + j * 3, affineMatrix(i, j)); } - parameterArray.SetElement(j + 9, Transform_elem(*lin, j, 3)); + parameterArray.SetElement(j + 9, affineMatrix(3, j)); } if (xfm->inverse_flag) @@ -138,14 +166,40 @@ MINCTransformIOTemplate::ReadOneTransform(VIO_General_tran using DisplacementFieldTransformType = DisplacementFieldTransform; using GridImageType = typename DisplacementFieldTransformType::DisplacementFieldType; using MincReaderType = ImageFileReader; + using OutputPixelType = vnl_vector_fixed; + const OutputPixelType RAStofromLPSVector = { -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, RAStofromLPSVector](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(), RAStofromLPSVector)); + oIt.Set(p); + } + }, + nullptr); + } TransformPointer transform; std::string transformTypeName = "DisplacementFieldTransform_"; @@ -155,11 +209,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 +284,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 RAStofromLPS; + Matrix affineMatrix; + + RAStofromLPS.SetIdentity(); + RAStofromLPS(0, 0) = -1.0; + RAStofromLPS(1, 1) = -1.0; + affineMatrix.SetIdentity(); + + + for (int j = 0; j < 3; ++j) + { + for (int i = 0; i < 3; ++i) + { + affineMatrix(j, i) = matrix(j, i); + } + affineMatrix(3, j) = offset[j]; + } + + if (this->m_RAStoLPS) + affineMatrix = RAStofromLPS * affineMatrix * RAStofromLPS; + 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) = affineMatrix(j, i); } - Transform_elem(lin, j, 3) = offset[j]; + Transform_elem(lin, j, 3) = affineMatrix(3, j); } // add 4th normalization row (not stored) Transform_elem(lin, 3, 3) = 1.0; @@ -245,29 +328,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 RAStofromLPSVector = { -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 rasGrid = GridImageType::New(); // flipped grid for RAS->LPS conversion + rasGrid->CopyInformation(grid); + rasGrid->SetRegions(grid->GetBufferedRegion()); + rasGrid->Allocate(true); + + itk::MultiThreaderBase::Pointer mt = itk::MultiThreaderBase::New(); + mt->ParallelizeImageRegion<3>( + grid->GetBufferedRegion(), + [grid, rasGrid, RAStofromLPSVector](const typename GridImageType::RegionType & region) { + itk::Vector p; + itk::ImageRegionConstIterator iIt(grid, region); + itk::ImageRegionIterator oIt(rasGrid, region); + for (; !iIt.IsAtEnd(); ++iIt, ++oIt) + { + p.SetVnlVector(element_product(iIt.Get().GetVnlVector(), RAStofromLPSVector)); + oIt.Set(p); + } + }, + nullptr); + + writer->SetInput(rasGrid); + } + 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 71f5e22ce8d..5aa2b16c746 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; constexpr 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; constexpr 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);