diff --git a/Modules/IO/MGHIO/CMakeLists.txt b/Modules/IO/MGHIO/CMakeLists.txt new file mode 100644 index 000000000000..e0254d3d4fb1 --- /dev/null +++ b/Modules/IO/MGHIO/CMakeLists.txt @@ -0,0 +1,3 @@ +project(MGHIO) +set(MGHIO_LIBRARIES MGHIO) +itk_module_impl() diff --git a/Modules/IO/MGHIO/README.md b/Modules/IO/MGHIO/README.md new file mode 100644 index 000000000000..23a2cf2761d2 --- /dev/null +++ b/Modules/IO/MGHIO/README.md @@ -0,0 +1,66 @@ +# MGHIO + +In-tree ITK module providing an `itk::MGHImageIO` plugin for the +`.mgh` / `.mgz` (MGH) volume format used by the FreeSurfer +neuroimaging toolchain. The factory `itk::MGHImageIOFactory` +auto-registers with `itk::ImageIOFactory` so any +`itk::ImageFileReader` / `itk::ImageFileWriter` transparently picks +up `.mgh` / `.mgz` files when the module is enabled. + +## Origin + +Ingested from the standalone remote module +[**InsightSoftwareConsortium/itkMGHImageIO**](https://github.com/InsightSoftwareConsortium/itkMGHImageIO) +on 2026-04-27, at upstream commit +[`af74507c`](https://github.com/InsightSoftwareConsortium/itkMGHImageIO/commit/af74507c1ddc82722637c37d1f5d169e3000553a). +The upstream repository will be archived read-only after this PR +merges; it remains reachable at the URL above. + +## What lives here + +Per the v3 ingestion strategy (see +`Utilities/Maintenance/RemoteModuleIngest/INGESTION_STRATEGY.md`), only +paths matching the narrow whitelist (code, headers, tests, wrapping, +module CMake) crossed the merge boundary: + +- `include/` — public C++ headers (`itkMGHImageIO.h`, `itkMGHImageIOFactory.h`). +- `src/` — non-template implementation. +- `test/` — CTest drivers and content-link stubs. +- `wrapping/` — Python wrapping descriptors. +- `CMakeLists.txt`, `itk-module.cmake` — build + module descriptors. + +Every surviving commit preserves original authorship; `git blame` +walks across the merge boundary to upstream authors going back to the +module's earliest history. + +## What was intentionally left upstream + +Everything outside the whitelist stays in the archived upstream repo. +If you need any of it, clone +. + +| Content in upstream | Why it did not ingest | +|---|---| +| `.github/`, `azure-pipelines.yml`, `Dockerfile` | Standalone-build CI scaffolding; not useful in-tree. | +| `CTestConfig.cmake`, `pyproject.toml`, `LICENSE`, `.clang-format` | Packaging / CI / style scaffolding superseded by ITK root. | +| Historical `Old/`, `paper/`, `docs/`, demo asset trees | Stripped by the whitelist to protect ITK's git pack from bloat. | + +## Long-form documentation + +- **Format reference** — see the FreeSurfer documentation: + . +- **Standalone build + history** — see the archived upstream at + . + +## Compliance level + +Compliance Level 2 (Alpha — niche execution environment dependence: +FreeSurfer-format I/O). Carried forward from the +`MODULE_COMPLIANCE_LEVEL 2` declaration in the previous +`Modules/Remote/MGHIO.remote.cmake`. + +## Content-link status + +The 3 input content-links under `test/MD5/` are already in `.cid` +(IPFS Content Identifier) form; no `.md5` → `.cid` normalization is +required for this module. diff --git a/Modules/IO/MGHIO/include/itkMGHImageIO.h b/Modules/IO/MGHIO/include/itkMGHImageIO.h new file mode 100644 index 000000000000..a3b3497b3eac --- /dev/null +++ b/Modules/IO/MGHIO/include/itkMGHImageIO.h @@ -0,0 +1,160 @@ +/*========================================================================= + * + * 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 itkMGHImageIO_h +#define itkMGHImageIO_h + +#include "itkMatrix.h" +#include "itkImageIOBase.h" +#include +#include "itk_zlib.h" + +#include "MGHIOExport.h" + +namespace itk +{ +/** \class MGHImageIO + * + * \author Hans J. Johnson + * \brief Class that defines how to read MGH file format. + * Originally developed as part of the Slicer software + * package under grants XXXX + * + * \ingroup IOFilters + * \ingroup MGHIO + */ +class MGHIO_EXPORT MGHImageIO : public ImageIOBase +{ +public: + /** Standard class type alias. */ + using Self = MGHImageIO; + using Superclass = ImageIOBase; + using Pointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkOverrideGetNameOfClassMacro(MGHImageIO); + + /*-------- This part of the interfaces deals with reading data. ----- */ + + /** Determine if the file can be read with this ImageIO implementation. + * \param FileNameToRead The name of the file to test for reading. + * \post Sets classes ImageIOBase::m_FileName variable to be FileNameToWrite + * \return Returns true if this ImageIO can read the file specified. + */ + bool + CanReadFile(const char * FileNameToRead) override; + + /** Set the spacing and dimension information for the set filename. */ + void + ReadImageInformation() override; + + /** Reads the data from disk into the memory buffer provided. */ + void + Read(void * pData) override; + + /*-------- This part of the interfaces deals with writing data. ----- */ + + /** Determine if the file can be written with this ImageIO implementation. + * \param FileNameToWrite The name of the file to test for writing. + * \post Sets classes ImageIOBase::m_FileName variable to be FileNameToWrite + * \return Returns true if this ImageIO can write the file specified. + */ + bool + CanWriteFile(const char * name) override; + + /** Set the spacing and dimension information for the set filename. */ + void + WriteImageInformation() override; + + /** Writes the data to disk from the memory buffer provided. Make sure + * that the IORegions has been set properly. */ + void + Write(const void * buffer) override; + +protected: + MGHImageIO(); + ~MGHImageIO() override; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + void + ReadVolumeHeader(); + +private: + static const int MRI_UCHAR = 0; + static const int MRI_INT = 1; + static const int MRI_FLOAT = 3; + static const int MRI_SHORT = 4; + static const int MRI_TENSOR = 6; + static const unsigned int FS_DIMENSION_HEADER_SIZE = sizeof(int) * 7; + static const unsigned int FS_RAS_HEADER_SIZE = (sizeof(float) * 15) + sizeof(short); + static const unsigned int FS_UNUSED_HEADER_SIZE = 256 - FS_RAS_HEADER_SIZE; + static const unsigned int FS_WHOLE_HEADER_SIZE = + FS_RAS_HEADER_SIZE + FS_DIMENSION_HEADER_SIZE + FS_UNUSED_HEADER_SIZE; + + /** check if a filename is for a compressed file */ + bool + IsCompressedFilename(const std::string fname); + /// processes the actual data buffer + void + SwapBytesIfNecessary(void * const buffer, const unsigned long numberOfPixels); + + /// examines the direction cosines and creates encapsulation data + // void MriDirCos(); + + void + WriteHeader(); + + void + WriteData(const void * buffer); + + void + PermuteFrameValues(const void * buffer, char * tempmemory); + + unsigned int + GetComponentSize() const override; + + std::string + GetOrientation(itk::Matrix directions); + + bool m_IsCompressed; + gzFile m_GZFile; + std::ofstream m_Output; + + // Utility function to assist with writing to disk in the + // proper format. TInType is static_cast type. + template + int + TWrite(const TInType inValue); + template + int + TRead(TOutType & outValue); + + int + TWrite(const char * buf, const unsigned long count); + void + OpenFile(); + void + CloseFile(); +}; +} // end namespace itk + +#endif // itkMGHImageIO_h diff --git a/Modules/IO/MGHIO/include/itkMGHImageIOFactory.h b/Modules/IO/MGHIO/include/itkMGHImageIOFactory.h new file mode 100644 index 000000000000..b9b3e28d208d --- /dev/null +++ b/Modules/IO/MGHIO/include/itkMGHImageIOFactory.h @@ -0,0 +1,74 @@ +/*========================================================================= + * + * 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 itkMGHImageIOFactory_h +#define itkMGHImageIOFactory_h + +#include "itkObjectFactoryBase.h" +#include "itkImageIOBase.h" + +#include "MGHIOExport.h" + +namespace itk +{ +/** \class MGHImageIOFactory + * \brief Create instances of MGHImageIO objects using an object factory. + * \ingroup MGHIO + */ +class MGHIO_EXPORT MGHImageIOFactory : public ObjectFactoryBase +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(MGHImageIOFactory); + + /** Standard class type alias */ + using Self = MGHImageIOFactory; + using Superclass = ObjectFactoryBase; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Class methods used to interface with the registered factories **/ + const char * + GetITKSourceVersion() const override; + + const char * + GetDescription() const override; + + /** Method for class instantiation **/ + itkFactorylessNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkOverrideGetNameOfClassMacro(MGHImageIOFactory); + + /** Register one factory of this type */ + static void + RegisterOneFactory() + { + MGHImageIOFactory::Pointer MGHFactory = MGHImageIOFactory::New(); + + ObjectFactoryBase::RegisterFactoryInternal(MGHFactory); + } + +protected: + MGHImageIOFactory(); + ~MGHImageIOFactory() override; + void + PrintSelf(std::ostream & os, Indent indent) const override; +}; +} // end namespace itk + +#endif /// itkMGHImageIOFactory_h diff --git a/Modules/IO/MGHIO/itk-module.cmake b/Modules/IO/MGHIO/itk-module.cmake new file mode 100644 index 000000000000..7e148c79b7c7 --- /dev/null +++ b/Modules/IO/MGHIO/itk-module.cmake @@ -0,0 +1,20 @@ +set( + DOCUMENTATION + "This modules contains an ImageIO class to read or write the + MGH file format that is an integral part of FreeSurfer based tools." +) + +itk_module( + MGHIO + ENABLE_SHARED + DEPENDS + ITKIOImageBase + ITKZLIB + TEST_DEPENDS + ITKTestKernel + ITKTransform + EXCLUDE_FROM_DEFAULT + FACTORY_NAMES + ImageIO::MGH + DESCRIPTION "${DOCUMENTATION}" +) diff --git a/Modules/IO/MGHIO/src/CMakeLists.txt b/Modules/IO/MGHIO/src/CMakeLists.txt new file mode 100644 index 000000000000..932468486b05 --- /dev/null +++ b/Modules/IO/MGHIO/src/CMakeLists.txt @@ -0,0 +1,7 @@ +set( + MGHIO_SRC + itkMGHImageIOFactory.cxx + itkMGHImageIO.cxx +) + +itk_module_add_library(MGHIO ${MGHIO_SRC}) diff --git a/Modules/IO/MGHIO/src/itkMGHImageIO.cxx b/Modules/IO/MGHIO/src/itkMGHImageIO.cxx new file mode 100644 index 000000000000..0180f474f659 --- /dev/null +++ b/Modules/IO/MGHIO/src/itkMGHImageIO.cxx @@ -0,0 +1,809 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "itkMGHImageIO.h" +#include "itkByteSwapper.h" +#include "itkMetaDataObject.h" +#include "itksys/SystemTools.hxx" + +namespace itk +{ + +// VALID file extensions +static const std::string __MGH_EXT(".mgh"); +static const std::string __MGZ_EXT(".mgz"); +static const std::string __GZ_EXT(".gz"); + +using MatrixType = itk::Matrix; +using VectorType = itk::Vector; + +static MatrixType +GetRAS2LPS() +{ + MatrixType RAS2LAS; + RAS2LAS.SetIdentity(); + RAS2LAS[0][0] = -1.0; + RAS2LAS[1][1] = -1.0; + RAS2LAS[2][2] = 1.0; + return RAS2LAS; +} + +// ------------------------------- +// +// Convert to BE +// +// ------------------------------- + +template +int +MGHImageIO ::TRead(TOutType & outValue) +{ + TDiskType onDiskValue; + const int result = ::gzread(this->m_GZFile, &onDiskValue, sizeof(TDiskType)); + itk::ByteSwapper::SwapFromSystemToBigEndian(&onDiskValue); + outValue = static_cast(onDiskValue); + return result; +} + +template +int +MGHImageIO ::TWrite(const TInType inValue) +{ + auto onDiskValue = static_cast(inValue); + itk::ByteSwapper::SwapFromSystemToBigEndian(&onDiskValue); + if (this->m_IsCompressed) + { + return ::gzwrite(this->m_GZFile, &onDiskValue, sizeof(TDiskType)); + } + else + { + this->m_Output.write(reinterpret_cast(&onDiskValue), sizeof(TDiskType)); + return this->m_Output.good() ? sizeof(TDiskType) : 0; + } +} + +int +MGHImageIO ::TWrite(const char * buf, const unsigned long count) +{ + if (this->m_IsCompressed) + { + const unsigned long int result = ::gzwrite(this->m_GZFile, buf, count); + if (result != count) + { + itkExceptionMacro(<< " Failed to write " << count << ", only wrote " << result); + } + return result; + } + else + { + this->m_Output.write(buf, count); + return this->m_Output.good() ? count : 0; + } +} +// -------------------------------------- +// +// MGHImageIO +// + + +MGHImageIO ::MGHImageIO() +{ + this->SetNumberOfDimensions(3); + std::fill(m_Dimensions.begin(), m_Dimensions.end(), 0U); + m_ByteOrder = (ByteSwapper::SystemIsBigEndian()) ? IOByteOrderEnum::BigEndian : IOByteOrderEnum::LittleEndian; +} + +MGHImageIO ::~MGHImageIO() +{ + // Nothing to do in destructor +} + +bool +MGHImageIO ::IsCompressedFilename(const std::string fname) +{ + // + // Originally MGHImageIO allowed ".mgh", ".mgz", + // "mgh.gz" and ".gz" + // + // The '.gz' extension is too ambiguous. It collides with NIfTI + // (.nii.gz) and given that there's no consistent 'magic number' in + // the header, the MGH library will blindly try and read any '.gz' + // file until it crashes or detects invalid data. + // + // So the bare '.gz' extension is no longer recognized. + const std::string lastExtension = itksys::SystemTools::GetFilenameLastExtension(fname.c_str()); + if (lastExtension == __MGZ_EXT) + { + return true; + } + if (lastExtension == __GZ_EXT) + { + const std::string fnameWithoutLastExtension = itksys::SystemTools::GetFilenameWithoutLastExtension(fname); + const std::string penultimateExtension = itksys::SystemTools::GetFilenameLastExtension(fnameWithoutLastExtension); + if (penultimateExtension == __MGH_EXT) + { + return true; + } + } + return false; +} + +bool +MGHImageIO ::CanReadFile(const char * FileNameToRead) +{ + const std::string filename(FileNameToRead); + + if (filename.empty()) + { + itkExceptionMacro(<< "A FileName must be specified."); + return false; + } + + // check if the correct extension is given by the user + const std::string extension = itksys::SystemTools::GetFilenameLastExtension(filename.c_str()); + if (extension == __MGH_EXT || this->IsCompressedFilename(filename)) + { + return true; + } + return false; +} + +void +MGHImageIO ::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + const std::string strSep = ", "; + + os << indent << "Data Dimensions: (" << m_Dimensions[0] << strSep << m_Dimensions[1] << strSep << m_Dimensions[2] + << ")\n" + << indent << "Data Spacing: (" << m_Spacing[0] << strSep << m_Spacing[1] << strSep << m_Spacing[2] << ")\n" + << indent << "Scalar Type: " << m_ComponentType << std::endl + << indent << "Number of Frames: " << m_NumberOfComponents << std::endl; + + os << indent << "RAS to IJK matrix: " << std::endl; +} + +void +MGHImageIO ::ReadImageInformation() +{ + this->m_GZFile = gzopen(m_FileName.c_str(), "rb"); + if (!this->m_GZFile) + { + itkExceptionMacro(<< "Can't find/open file: " << m_FileName); + return; + } + + ReadVolumeHeader(); + gzclose(this->m_GZFile); +} + +void +MGHImageIO ::ReadVolumeHeader() +{ + // check file reading + if (!this->m_GZFile) + { + itkExceptionMacro(<< "Can't find/open file: " << this->m_FileName); + return; + } + int version; + this->TRead(version); + this->TRead(m_Dimensions[0]); + this->TRead(m_Dimensions[1]); + this->TRead(m_Dimensions[2]); + // next is nframes + this->TRead(m_NumberOfComponents); + int type; + this->TRead(type); + int dof; + this->TRead(dof); // what does this do? + + // Convert type to an ITK type + switch (type) + { + case MRI_UCHAR: + { + m_ComponentType = IOComponentEnum::UCHAR; + } + break; + case MRI_INT: + { + m_ComponentType = IOComponentEnum::INT; + } + break; + case MRI_FLOAT: + { + m_ComponentType = IOComponentEnum::FLOAT; + } + break; + case MRI_SHORT: + { + m_ComponentType = IOComponentEnum::SHORT; + } + break; + case MRI_TENSOR: + { + m_ComponentType = IOComponentEnum::FLOAT; + m_NumberOfComponents = 9; + } + break; + default: + itkExceptionMacro(<< " Unknown data type " << type << " using float by default."); + m_ComponentType = IOComponentEnum::FLOAT; + } + + // Next short says whether RAS registration information is good. + // If so, read the voxel size and then the matrix + short RASgood; + this->TRead(RASgood); + if (RASgood) + { + for (size_t nSpacing = 0; nSpacing < 3; ++nSpacing) + { + this->TRead(m_Spacing[nSpacing]); + } + /* + From https://www.nmr.mgh.harvard.edu/~tosa/#coords: + To go from freesurfer voxel coordinates to RAS coordinates, they use: + translation: t_r, t_a, t_s is defined using c_r, c_a, c_s centre voxel position in RAS + rotation: direction cosines x_(r,a,s), y_(r,a,s), z_(r,a,s) + voxel size for scale: s_x, s_y, s_z + + [ x_r y_r z_r t_r][s_x 0 0 0] + [ x_a y_a z_a t_a][0 s_y 0 0] + [ x_s y_s z_s t_s][0 0 s_z 0] + [ 0 0 0 1 ][0 0 0 1] + Voxel center is a column matrix, multipled from the right + [v_x] + [v_y] + [v_z] + [ 1 ] + + In the MGH header, they hold: + x_r x_a x_s + y_r y_a y_s + z_r z_a z_s + c_r c_a c_s + */ + MatrixType MGHdirMatrix; + // reading in x_r x_a x_s y_r y_a y_s z_r z_a z_s and putting it into the + // matrix as: + // x_r y_r z_r + // x_a y_a z_a + // x_s y_s z_s + for (unsigned int c = 0; c < 3; ++c) + { + for (unsigned int r = 0; r < 3; ++r) // NOTE: Data stored row-major form, so traverse rows first + { + this->TRead(MGHdirMatrix[r][c]); + } + } + // getting orientation must be done before RAS->LAS conversion + // but not used const std::string orientation = GetOrientation( MGHdirMatrix ); + // convert the coordinates from RAS to LPS, as the ITK archetype assumes + // LPS volumes. volume orientation not related to scan order, always convert + MatrixType ITKdirMatrix = GetRAS2LPS() * MGHdirMatrix; // Pre-multipy by conversion + for (size_t c = 0; c < 3; ++c) + { + // now take x_r, x_a, x_s out of the matrix and set it to the direction + // vector 0, same for y_* and direction vector 1, z_* and vector 2 + std::vector vDir; + for (size_t r = 0; r < 3; ++r) + { + vDir.push_back(ITKdirMatrix[r][c]); + } + SetDirection(c, vDir); + } + + VectorType MGHcenterVoxel; // c_r c_a c_s + for (size_t ui = 0; ui < 3; ++ui) + { + this->TRead(MGHcenterVoxel[ui]); + } + // convert C to from RAS to LPS + VectorType ITKcenterVoxel = GetRAS2LPS() * MGHcenterVoxel; + + + // MriDirCos(); // convert direction cosines + // finally, store the origin of the image -> only works + // if the image is properly oriented in the sequel + // + // computed in CORONAL orientation = ITK_COORDINATE_ORIENTATION_LIA + VectorType fc; + fc[0] = m_Dimensions[0] * 0.5; + fc[1] = m_Dimensions[1] * 0.5; + fc[2] = m_Dimensions[2] * 0.5; + MatrixType spcing; + spcing[0][0] = m_Spacing[0]; + spcing[1][1] = m_Spacing[1]; + spcing[2][2] = m_Spacing[2]; + + VectorType ITKorigin = ITKcenterVoxel - (ITKdirMatrix * spcing * fc); + m_Origin[0] = ITKorigin[0]; + m_Origin[1] = ITKorigin[1]; + m_Origin[2] = ITKorigin[2]; + } + + // ================== + // read tags at the end of file + const size_t numValues = m_Dimensions[0] * m_Dimensions[1] * m_Dimensions[2]; + gzseek( + this->m_GZFile, FS_WHOLE_HEADER_SIZE + (m_NumberOfComponents * numValues * this->GetComponentSize()), SEEK_SET); + + float fBufTR; + // read TR, Flip, TE, FI, FOV + if (this->TRead(fBufTR)) + { + itk::MetaDataDictionary & thisDic = this->GetMetaDataDictionary(); + itk::EncapsulateMetaData(thisDic, std::string("TR"), fBufTR); + + // try to read flipAngle + float fBufFA; + if (this->TRead(fBufFA)) + { + itk::EncapsulateMetaData(thisDic, std::string("FlipAngle"), fBufFA); + // TE + float fBufTE; + if (this->TRead(fBufTE)) + { + itk::EncapsulateMetaData(thisDic, std::string("TE"), fBufTE); + // TI + float fBufTI; + if (this->TRead(fBufTI)) + { + itk::EncapsulateMetaData(thisDic, std::string("TI"), fBufTI); + // FOV + float fBufFOV; + if (this->TRead(fBufFOV)) + { + itk::EncapsulateMetaData(thisDic, std::string("FoV"), fBufFOV); + } + } + } + } + } +} + +void +MGHImageIO ::Read(void * pData) +{ + this->m_GZFile = gzopen(m_FileName.c_str(), "rb"); + if (!this->m_GZFile) + { + itkExceptionMacro(<< "Can't find/open file: " << m_FileName); + return; + } + + const unsigned long numPixels = m_Dimensions[0] * m_Dimensions[1] * m_Dimensions[2]; + + const unsigned int componentSize(this->GetComponentSize()); + + // check that the offset is actually computed wrt. the beginning + gzseek(this->m_GZFile, FS_WHOLE_HEADER_SIZE, SEEK_SET); + + const unsigned int frameSize = numPixels * componentSize; + + if (m_NumberOfComponents > 1) + { + auto * pBuffer = new char[frameSize]; + + const unsigned int pixelSize = componentSize * m_NumberOfComponents; + for (unsigned int frameIndex = 0; frameIndex < m_NumberOfComponents; ++frameIndex) + { + // read current frame + gzread(this->m_GZFile, pBuffer, frameSize); + // copy memory location in the final buffer + + auto * pSrc = (char *)pBuffer; + auto * pDst = (char *)pData; + + pDst += frameIndex * componentSize; + for (unsigned int ui = 0; ui < numPixels; ++ui, pSrc += componentSize, pDst += pixelSize) + { + for (unsigned int byteCount = 0; byteCount < componentSize; ++byteCount) + { + *(pDst + byteCount) = *(pSrc + byteCount); + } + } // next ui + } // next frameIndex + + // clear resources + delete[] pBuffer; + } + else + { + gzread(this->m_GZFile, pData, frameSize); + } + + gzclose(this->m_GZFile); + + SwapBytesIfNecessary(pData, numPixels * m_NumberOfComponents); + +} // end Read function + +void +MGHImageIO ::SwapBytesIfNecessary(void * const buffer, const unsigned long numberOfPixels) +{ + // NOTE: If machine order is little endian, and the data needs to be + // swapped, the SwapFromBigEndianToSystem is equivalent to + // SwapFromSystemToBigEndian. + + switch (m_ComponentType) + { + case IOComponentEnum::UCHAR: + { + ByteSwapper::SwapRangeFromSystemToBigEndian((unsigned char *)buffer, numberOfPixels); + } + break; + case IOComponentEnum::SHORT: + { + ByteSwapper::SwapRangeFromSystemToBigEndian((short *)buffer, numberOfPixels); + } + break; + case IOComponentEnum::INT: + { + ByteSwapper::SwapRangeFromSystemToBigEndian((int *)buffer, numberOfPixels); + } + break; + case IOComponentEnum::FLOAT: + { + ByteSwapper::SwapRangeFromSystemToBigEndian((float *)buffer, numberOfPixels); + } + break; + default: + ExceptionObject exception(__FILE__, __LINE__); + exception.SetDescription("Pixel Type Unknown"); + throw exception; + } // end switch +} + +bool +MGHImageIO ::CanWriteFile(const char * name) +{ + const std::string filename(name); + + if (filename.empty()) + { + itkExceptionMacro(<< "A FileName must be specified."); + return false; + } + + // TODO: Use ITK Extension extractor + const std::string extension = itksys::SystemTools::GetFilenameExtension(filename.c_str()); + if (extension != __MGH_EXT && !this->IsCompressedFilename(filename)) + { + return false; + } + return true; +} + +void +MGHImageIO ::OpenFile() +{ + if (this->m_IsCompressed) + { + this->m_GZFile = gzopen(m_FileName.c_str(), "wb"); + if (this->m_GZFile == nullptr) + { + itkExceptionMacro(<< " Failed to open gzFile for writing"); + itkExceptionMacro(<< " File cannot be written"); + } + } + else + { + this->m_Output.open(m_FileName.c_str(), std::ios::out | std::ios::binary); + if (this->m_Output.fail()) + { + itkExceptionMacro(<< " File cannot be written"); + } + } +} +void +MGHImageIO ::CloseFile() +{ + if (this->m_IsCompressed) + { + gzclose(this->m_GZFile); + } + else + { + this->m_Output.close(); + } +} +void +MGHImageIO ::WriteImageInformation() +{ + this->m_IsCompressed = this->IsCompressedFilename(this->m_FileName); + this->OpenFile(); + this->WriteHeader(); + this->CloseFile(); +} + +void +MGHImageIO ::Write(const void * buffer) +{ + this->m_IsCompressed = this->IsCompressedFilename(this->m_FileName); + this->OpenFile(); + this->WriteHeader(); + this->WriteData(buffer); + this->CloseFile(); +} + +void +MGHImageIO ::WriteHeader() +{ + // version + constexpr int mghVersion = 1; + this->TWrite(mghVersion); + // dimensions + for (size_t ui = 0; ui < 3; ++ui) + { + this->TWrite(m_Dimensions[ui]); + } + + // nframes + this->TWrite(m_NumberOfComponents); + + // type + switch (m_ComponentType) + { + case IOComponentEnum::UCHAR: + { + this->TWrite(MRI_UCHAR); + } + break; + case IOComponentEnum::INT: + { + this->TWrite(MRI_INT); + } + break; + case IOComponentEnum::FLOAT: + { + this->TWrite(MRI_FLOAT); + } + break; + case IOComponentEnum::SHORT: + { + this->TWrite(MRI_SHORT); + } + break; + default: + itkExceptionMacro(<< "MGHImageIO supports unsigned char, int, float and short"); + } + + // dof !?! -> default value = 1 + this->TWrite(1); + + // write RAS and voxel size info + // for now, RAS flag will be good + // in the future, check if the m_Directions matrix is a permutation matrix + this->TWrite(1); + // spacing + for (unsigned int ui = 0; ui < 3; ++ui) + { + this->TWrite(m_Spacing[ui]); + } + + // get directions matrix + MatrixType ITKdirMatrix; + for (unsigned int c = 0; c < 3; ++c) + { + const std::vector & dir_line = GetDirection(c); + for (unsigned int r = 0; r < 3; ++r) + { + ITKdirMatrix[r][c] = dir_line[r]; + } + } + MatrixType MGHdirMatrix = GetRAS2LPS() * ITKdirMatrix; // Pre-multipy by RAS2LPS + // writing in x_r x_a x_s y_r y_a y_s z_r z_a z_s and putting it into the + // matrix as: + // x_r y_r z_r + // x_a y_a z_a + // x_s y_s z_s + for (unsigned int c = 0; c < 3; ++c) + { + for (unsigned int r = 0; r < 3; ++r) // Write in row-major order + { + this->TWrite(MGHdirMatrix[r][c]); + } + } + + VectorType fc; + fc[0] = m_Dimensions[0] * 0.5; + fc[1] = m_Dimensions[1] * 0.5; + fc[2] = m_Dimensions[2] * 0.5; + VectorType ITKorigin; + ITKorigin[0] = m_Origin[0]; + ITKorigin[1] = m_Origin[1]; + ITKorigin[2] = m_Origin[2]; + MatrixType spcing; + spcing.SetIdentity(); + spcing[0][0] = m_Spacing[0]; + spcing[1][1] = m_Spacing[1]; + spcing[2][2] = m_Spacing[2]; + const VectorType ITKcenterVoxel = ITKorigin + (ITKdirMatrix * spcing * fc); + const VectorType MGHcenterVoxel = GetRAS2LPS() * ITKcenterVoxel; + for (size_t ui = 0; ui < 3; ++ui) + { + this->TWrite(MGHcenterVoxel[ui]); + } + // fill the rest of the buffer with zeros + const char zerobyte(0); + for (size_t i = 0; i < FS_UNUSED_HEADER_SIZE; ++i) + { + this->TWrite(zerobyte); + } +} + +void +MGHImageIO ::WriteData(const void * buffer) +{ + // swap bytes if necessary + const unsigned int numPixels = m_Dimensions[0] * m_Dimensions[1] * m_Dimensions[2]; + const unsigned long int numvalues = numPixels * m_NumberOfComponents; + const unsigned long int numbytes = this->GetComponentSize() * numvalues; + + auto * tempmemory = new char[numbytes]; + + // re-arrange data in frames + if (m_NumberOfComponents > 1) + { + PermuteFrameValues(buffer, tempmemory); + } + else + { + memcpy(tempmemory, buffer, numbytes); + } + + this->SwapBytesIfNecessary(tempmemory, numvalues); + + this->TWrite(tempmemory, this->GetImageSizeInBytes()); + delete[] tempmemory; + + // if present, the scan parameters are present at the end of the file, so now's the time to write them + itk::MetaDataDictionary & thisDic = this->GetMetaDataDictionary(); + + float fScanBuffer = 0.0F; + if (ExposeMetaData(thisDic, "TR", fScanBuffer)) + { + this->TWrite(fScanBuffer); + } // end TR + if (ExposeMetaData(thisDic, "FlipAngle", fScanBuffer)) + { + this->TWrite(fScanBuffer); + } // end FlipAngle + + if (ExposeMetaData(thisDic, "TE", fScanBuffer)) + { + this->TWrite(fScanBuffer); + } // end TE + + if (ExposeMetaData(thisDic, "TI", fScanBuffer)) + { + this->TWrite(fScanBuffer); + } // end TI + + if (ExposeMetaData(thisDic, "FoV", fScanBuffer)) + { + this->TWrite(fScanBuffer); + } // end FoV + + // no need to close the stream +} + +void +MGHImageIO ::PermuteFrameValues(const void * buffer, char * tempmemory) +{ + const unsigned int numPixels = m_Dimensions[0] * m_Dimensions[1] * m_Dimensions[2]; + const unsigned int valueSize(this->GetComponentSize()); + const unsigned int frameSize = numPixels * valueSize; + + const auto * pSrc = (const char *)buffer; + auto * pDst = (char *)tempmemory; + + for (unsigned int pixelIndex = 0; pixelIndex < numPixels; ++pixelIndex, pDst += valueSize) + { + for (unsigned int componentIndex = 0; componentIndex < m_NumberOfComponents; ++componentIndex, pSrc += valueSize) + { + std::copy(pSrc, pSrc + valueSize, pDst + frameSize * componentIndex); + } // next component index + } // next pixelIndex +} + +unsigned int +MGHImageIO ::GetComponentSize() const +{ + unsigned int returnValue; + switch (m_ComponentType) + { + case IOComponentEnum::UCHAR: + { + returnValue = sizeof(unsigned char); + } + break; + case IOComponentEnum::SHORT: + { + returnValue = sizeof(short); + } + break; + case IOComponentEnum::INT: + { + returnValue = sizeof(int); + } + break; + case IOComponentEnum::FLOAT: + { + returnValue = sizeof(float); + } + break; + default: + itkExceptionMacro(<< "MGHImageIO supports unsigned char, int, float and short"); + } + return returnValue; +} + +/** + * Examines the direction cosines and creates an Orientation String. + * The Orientation String is a three character string indicating the primary + * direction of each axis in the 3d matrix. The characters can be L,R,A,P,I,S. + **/ +std::string +MGHImageIO ::GetOrientation(itk::Matrix directions) +{ + std::string orientation(""); + for (int cAxes = 0; cAxes < 3; cAxes++) + { + const double sag = directions(0, cAxes); // LR axis + const double cor = directions(1, cAxes); // PA axis + const double ax = directions(2, cAxes); // IS axis + if (itk::Math::abs(sag) > itk::Math::abs(cor) && itk::Math::abs(sag) > itk::Math::abs(ax)) + { + if (sag > 0) + { + orientation += "R"; + } + else + { + orientation += "L"; + } + continue; + } + if (itk::Math::abs(cor) > itk::Math::abs(ax)) + { + if (cor > 0) + { + orientation += "A"; + } + else + { + orientation += "P"; + } + continue; + } + if (ax > 0) + { + orientation += "S"; + } + else + { + orientation += "I"; + } + } + return orientation; +} +} // end namespace itk diff --git a/Modules/IO/MGHIO/src/itkMGHImageIOFactory.cxx b/Modules/IO/MGHIO/src/itkMGHImageIOFactory.cxx new file mode 100644 index 000000000000..5d9aa3233668 --- /dev/null +++ b/Modules/IO/MGHIO/src/itkMGHImageIOFactory.cxx @@ -0,0 +1,64 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "itkMGHImageIOFactory.h" +#include "itkMGHImageIO.h" +#include "itkVersion.h" + + +namespace itk +{ +void +MGHImageIOFactory::PrintSelf(std::ostream &, Indent) const +{} + +MGHImageIOFactory::MGHImageIOFactory() +{ + this->RegisterOverride( + "itkImageIOBase", "itkMGHImageIO", "MGH Image IO", true, CreateObjectFunction::New()); +} + +MGHImageIOFactory::~MGHImageIOFactory() = default; + +const char * +MGHImageIOFactory::GetITKSourceVersion() const +{ + return ITK_SOURCE_VERSION; +} + +const char * +MGHImageIOFactory::GetDescription() const +{ + return "MGH ImageIO Factory, allows the loading of MGH/MGZ images into Insight"; +} +// Undocumented API used to register during static initialization. +// DO NOT CALL DIRECTLY. + +static bool MGHImageIOFactoryHasBeenRegistered = false; + +void MGHIO_EXPORT +MGHImageIOFactoryRegister__Private() +{ + if (!MGHImageIOFactoryHasBeenRegistered) + { + MGHImageIOFactoryHasBeenRegistered = true; + MGHImageIOFactory::RegisterOneFactory(); + } +} + +} // end namespace itk diff --git a/Modules/IO/MGHIO/test/CMakeLists.txt b/Modules/IO/MGHIO/test/CMakeLists.txt new file mode 100644 index 000000000000..bd38ae8a6623 --- /dev/null +++ b/Modules/IO/MGHIO/test/CMakeLists.txt @@ -0,0 +1,67 @@ +itk_module_test() +set(MGHIOTests itkMGHImageIOTest.cxx) + +set(MGH_DATA_ROOT ${CMAKE_CURRENT_SOURCE_DIR}/MD5) + +createtestdriver(MGHIO "${MGHIO-Test_LIBRARIES}" "${MGHIOTests}") + +itk_add_test( + NAME MGHFactoryCreationTest + COMMAND + MGHIOTestDriver + itkMGHImageIOTest + ${ITK_TEST_OUTPUT_DIR} + FactoryCreationTest +) + +itk_add_test( + NAME MGHReadImagesTest_mgz + COMMAND + MGHIOTestDriver + itkMGHImageIOTest + ${ITK_TEST_OUTPUT_DIR} + ReadImagesTest + DATA{${MGH_DATA_ROOT}/T1.mgz} + TEST.mgz +) + +itk_add_test( + NAME MGHReadImagesTest_mgh + COMMAND + MGHIOTestDriver + itkMGHImageIOTest + ${ITK_TEST_OUTPUT_DIR} + ReadImagesTest + DATA{${MGH_DATA_ROOT}/T1_uncompresed.mgh} + TEST.mgh +) + +itk_add_test( + NAME MGHReadImagesTest_mgh.gz + COMMAND + MGHIOTestDriver + itkMGHImageIOTest + ${ITK_TEST_OUTPUT_DIR} + ReadImagesTest + DATA{${MGH_DATA_ROOT}/T1_longname.mgh.gz} + TEST.mgh.gz +) + +itk_add_test( + NAME itkMGHIOInternalTests + COMMAND + MGHIOTestDriver + itkMGHImageIOTest + ${ITK_TEST_OUTPUT_DIR} + TestReadWriteOfSmallImageOfAllTypes +) +itk_add_test( + NAME itkMGHIOOriginTest + COMMAND + MGHIOTestDriver + itkMGHImageIOTest + ${ITK_TEST_OUTPUT_DIR} + TestOriginWriteTest + DATA{${MGH_DATA_ROOT}/T1_longname.mgh.gz} + OriginTEST.mgh.gz +) diff --git a/Modules/IO/MGHIO/test/MD5/T1.mgz.cid b/Modules/IO/MGHIO/test/MD5/T1.mgz.cid new file mode 100644 index 000000000000..a416ac4df70b --- /dev/null +++ b/Modules/IO/MGHIO/test/MD5/T1.mgz.cid @@ -0,0 +1 @@ +bafkreihdhvkddcmmtfgd5wugnrssr5cvzm4ab76lqj53e6f6tybyqrs3xy diff --git a/Modules/IO/MGHIO/test/MD5/T1_longname.mgh.gz.cid b/Modules/IO/MGHIO/test/MD5/T1_longname.mgh.gz.cid new file mode 100644 index 000000000000..a416ac4df70b --- /dev/null +++ b/Modules/IO/MGHIO/test/MD5/T1_longname.mgh.gz.cid @@ -0,0 +1 @@ +bafkreihdhvkddcmmtfgd5wugnrssr5cvzm4ab76lqj53e6f6tybyqrs3xy diff --git a/Modules/IO/MGHIO/test/MD5/T1_uncompresed.mgh.cid b/Modules/IO/MGHIO/test/MD5/T1_uncompresed.mgh.cid new file mode 100644 index 000000000000..f517f3270109 --- /dev/null +++ b/Modules/IO/MGHIO/test/MD5/T1_uncompresed.mgh.cid @@ -0,0 +1 @@ +bafkreib3qleljdy2t7czavhmkyzst5hlpobk7rt5pyix2n74hiaspyhfzi diff --git a/Modules/IO/MGHIO/test/itkMGHImageIOTest.cxx b/Modules/IO/MGHIO/test/itkMGHImageIOTest.cxx new file mode 100644 index 000000000000..3e711920bdc3 --- /dev/null +++ b/Modules/IO/MGHIO/test/itkMGHImageIOTest.cxx @@ -0,0 +1,154 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "itkMGHImageIOFactory.h" +#include "itkRandomImageSource.h" +#include "itkImageFileReader.h" +#include "itksys/SystemTools.hxx" +#include "itkMetaDataObject.h" +#include "itkIOCommon.h" + +#include "itkIOTestHelper.h" +#include "itkMGHImageIOTest.h" +#include + +int +itkMGHImageIOTest(int argc, char * argv[]) +{ + itk::ObjectFactoryBase::UnRegisterAllFactories(); + itk::MGHImageIOFactory::RegisterOneFactory(); + if (argc < 3) + { + std::cerr + << "ERROR: Incorrect number of arguments [ImageFileNameRead ImageFileNameWrite]" + << std::endl; + std::cerr << "Only " << argc << " arguments given." << std::endl; + return EXIT_FAILURE; + } + for (int i = 0; i < argc; ++i) + { + std::cout << i << " argv= " << argv[i] << std::endl; + } + // + // first argument is passing in the writable directory to do all testing + itksys::SystemTools::ChangeDirectory(argv[1]); + + static bool firstTime = true; + if (firstTime) + { + itk::ObjectFactoryBase::RegisterFactory(itk::MGHImageIOFactory::New()); + firstTime = false; + } + const std::string TestMode(argv[2]); + + bool returnSucceeded = true; + if (TestMode == std::string("FactoryCreationTest")) + // Tests added to increase code coverage. + { + itk::MGHImageIOFactory::Pointer MyFactoryTest = itk::MGHImageIOFactory::New(); + if (MyFactoryTest.IsNull()) + { + returnSucceeded &= false; + } + // This was made a protected function. MyFactoryTest->PrintSelf(std::cout,0); + } + else if (TestMode == std::string("TestReadWriteOfSmallImageOfAllTypes")) + { + std::string fn("test.mgz"); + // TODO: Need to test with images of non-identity direction cosigns, spacing, origin + returnSucceeded &= itkMGHImageIOTestReadWriteTest(fn, 3, "null", true); + returnSucceeded &= itkMGHImageIOTestReadWriteTest(fn, 3, "null", true); + returnSucceeded &= itkMGHImageIOTestReadWriteTest(fn, 3, "null", true); + returnSucceeded &= itkMGHImageIOTestReadWriteTest(fn, 3, "null", true); + returnSucceeded &= itkMGHImageIOTestReadWriteTest, 3>(fn, 3, "null", true); + } + else if (TestMode == std::string("ReadImagesTest")) // This is a mechanism for reading unsigned int images for + // testing. + { + using ImageType = itk::Image; + const std::string imageToBeRead(argv[3]); + const std::string imageToBeWritten(argv[4]); + try + { + std::cout << "Reading Image: " << imageToBeRead << std::endl; + ImageType::Pointer input = itk::IOTestHelper::ReadImage(imageToBeRead); + std::cout << input << std::endl; + itk::ImageFileWriter::Pointer testFactoryWriter = itk::ImageFileWriter::New(); + testFactoryWriter->SetFileName(imageToBeWritten); + testFactoryWriter->SetInput(input); + testFactoryWriter->Update(); + itk::ImageFileReader::Pointer testFactoryReader = itk::ImageFileReader::New(); + testFactoryReader->SetFileName(imageToBeWritten); + testFactoryReader->Update(); + ImageType::Pointer new_image = testFactoryReader->GetOutput(); + } + catch (itk::ExceptionObject & e) + { + e.Print(std::cerr); + returnSucceeded &= false; + } + } + else if (TestMode == "TestOriginWriteTest") + { + using ImageType = itk::Image; + ImageType::Pointer input; + const std::string imageToBeRead(argv[3]); + const std::string imageToBeWritten(argv[4]); + try + { + std::cout << "Reading Image: " << imageToBeRead << std::endl; + input = itk::IOTestHelper::ReadImage(imageToBeRead); + std::cout << input << std::endl; + + ImageType::PointType reference_origin; + reference_origin[0] = -123.4; + reference_origin[1] = 456.7; + reference_origin[2] = -890.0; + input->SetOrigin(reference_origin); + + itk::ImageFileWriter::Pointer testFactoryWriter = itk::ImageFileWriter::New(); + + testFactoryWriter->SetFileName(imageToBeWritten); + testFactoryWriter->SetInput(input); + testFactoryWriter->Update(); + itk::ImageFileReader::Pointer testFactoryReader = itk::ImageFileReader::New(); + testFactoryReader->SetFileName(imageToBeWritten); + testFactoryReader->Update(); + ImageType::Pointer new_image = testFactoryReader->GetOutput(); + const ImageType::PointType test_origin = new_image->GetOrigin(); + const double dist = reference_origin.EuclideanDistanceTo(test_origin); + if (dist > 1.0E-4) + { + std::cerr << std::setprecision(10) << "Origin written and origin read do not match: " + << "written: " << reference_origin << " read: " << test_origin << " distance: " << dist << std::endl; + returnSucceeded &= false; + } + } + catch (itk::ExceptionObject & e) + { + e.Print(std::cerr); + returnSucceeded &= false; + } + } + else + { + std::cerr << "Invalid TestMode : " << TestMode << std::endl; + returnSucceeded &= false; + } + return (returnSucceeded == true) ? EXIT_SUCCESS : EXIT_FAILURE; +} diff --git a/Modules/IO/MGHIO/test/itkMGHImageIOTest.h b/Modules/IO/MGHIO/test/itkMGHImageIOTest.h new file mode 100644 index 000000000000..66c3c756823d --- /dev/null +++ b/Modules/IO/MGHIO/test/itkMGHImageIOTest.h @@ -0,0 +1,314 @@ +/*========================================================================= + * + * 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 itkMGHImageIOTest_h +#define itkMGHImageIOTest_h + +#include +#include +#include +#include "itkImageRegionIterator.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkMGHImageIO.h" +#include "itkImage.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkRandomImageSource.h" +#include "itkDiffusionTensor3D.h" + +#include "itkEuler3DTransform.h" + + +template +typename itk::Image::Pointer +itkMGHImageIOTestGenerateRandomImage(const unsigned int size) +{ + using EulerTransformType = itk::Euler3DTransform; + EulerTransformType::Pointer eulerTransform = EulerTransformType::New(); + eulerTransform->SetIdentity(); + + // 15 degrees in radians + const double angleX = 15.0 * std::atan(1.0) / 45.0; + // 10 degrees in radians + const double angleY = 10.0 * std::atan(1.0) / 45.0; + // 5 degrees in radians + const double angleZ = 5.0 * std::atan(1.0) / 45.0; + eulerTransform->SetRotation(angleX, angleY, angleZ); + + using ImageType = itk::Image; + + typename ImageType::DirectionType permuteDirections; + const double dir1[3] = { -1., 0., 0. }; + const double dir2[3] = { 0., 0., 1. }; + const double dir3[3] = { 0., -1., 0. }; + permuteDirections.GetVnlMatrix().set_column(0, dir1); + permuteDirections.GetVnlMatrix().set_column(1, dir2); + permuteDirections.GetVnlMatrix().set_column(2, dir3); + + typename ImageType::DirectionType direction = eulerTransform->GetMatrix() * permuteDirections; + + for (unsigned int i = 0; i < VImageDimension; ++i) + { + for (unsigned int j = 0; j < VImageDimension; ++j) + { + direction[i][j] = static_cast(direction[i][j]); // Truncate for testing purposes + } + } + + typename ImageType::SizeType sz; + typename ImageType::SpacingType spacing; + typename ImageType::PointType origin; + + for (unsigned int i = 0; i < VImageDimension; ++i) + { + sz[i] = size; + spacing[i] = static_cast(i + 1.234567); + origin[i] = static_cast(1234.5); + } + + typename itk::RandomImageSource::Pointer source = itk::RandomImageSource::New(); + + source->SetDirection(direction); + source->SetSize(sz); + source->SetOrigin(origin); + source->SetSpacing(spacing); + + source->Update(); + typename ImageType::Pointer outImage = source->GetOutput(); + + { + itk::MetaDataDictionary & thisDic = outImage->GetMetaDataDictionary(); + // Add meta data to dictionary + // set TR, Flip, TE, FI, FOV //TODO: Add code that verifies these values + float fBufTR = 2.0F; + float fBufFA = 89.1F; + float fBufTE = 1.5F; + float fBufTI = 0.75F; + float fBufFOV = 321.0F; + itk::EncapsulateMetaData(thisDic, std::string("TR"), fBufTR); + // try to read flipAngle + itk::EncapsulateMetaData(thisDic, std::string("FlipAngle"), fBufFA); + // TE + itk::EncapsulateMetaData(thisDic, std::string("TE"), fBufTE); + // TI + itk::EncapsulateMetaData(thisDic, std::string("TI"), fBufTI); + // FOV + itk::EncapsulateMetaData(thisDic, std::string("FoV"), fBufFOV); + } + return outImage; +} + +// Template specialization for itkDiffusionTensor3D +template <> +itk::Image, 3>::Pointer +itkMGHImageIOTestGenerateRandomImage, 3>(unsigned int size) +{ + using TensorImageType = itk::Image, 3>; + using TensorImagePointer = TensorImageType::Pointer; + + const double dir1[3] = { 0, -1, 0 }; + const double dir2[3] = { 1, 0, 0 }; + const double dir3[3] = { 0, 0, -1 }; + + TensorImageType::DirectionType direction; + direction.GetVnlMatrix().set_column(0, dir1); + direction.GetVnlMatrix().set_column(1, dir2); + direction.GetVnlMatrix().set_column(2, dir3); + + TensorImageType::SizeType sz; + TensorImageType::SpacingType spacing; + TensorImageType::PointType origin; + for (unsigned int i = 0; i < 3; ++i) + { + sz[i] = size; + spacing[i] = static_cast(i * 3.21 + 1.0); + origin[i] = static_cast(1.23456); + } + + TensorImagePointer tensorImage = TensorImageType::New(); + tensorImage->SetRegions(sz); + tensorImage->SetSpacing(spacing); + tensorImage->SetOrigin(origin); + tensorImage->Allocate(); + tensorImage->SetDirection(direction); + + itk::DiffusionTensor3D pix; + using ScalarImageType = itk::Image; + using ScalarImagePointer = ScalarImageType::Pointer; + + std::vector scalarImageVec; + for (unsigned i = 0; i < pix.Size(); ++i) + { + ScalarImagePointer tempImage = itkMGHImageIOTestGenerateRandomImage(size); + scalarImageVec.push_back(tempImage); + } + itk::ImageRegionIteratorWithIndex tensorIt(tensorImage, tensorImage->GetLargestPossibleRegion()); + for (tensorIt.GoToBegin(); !tensorIt.IsAtEnd(); ++tensorIt) + { + for (unsigned int i = 0; i < pix.Size(); ++i) + { + const TensorImageType::IndexType & index = tensorIt.GetIndex(); + pix.SetNthComponent(i, scalarImageVec[i]->GetPixel(index)); + } + tensorIt.Set(pix); + } + return tensorImage; +} + + +template +bool +itkMGHImageIOTestReadWriteTest(std::string fn, unsigned int size, std::string inputFile, bool compression = false) +{ + using ImageType = itk::Image; + + typename itk::ImageFileReader::Pointer reader = itk::ImageFileReader::New(); + typename itk::ImageFileWriter::Pointer writer = itk::ImageFileWriter::New(); + + itk::MGHImageIO::Pointer io = itk::MGHImageIO::New(); + reader->SetImageIO(io); + writer->SetImageIO(io); + + typename ImageType::Pointer reference_image; + + if (inputFile != "null") + { + typename itk::ImageFileReader::Pointer tmpReader = itk::ImageFileReader::New(); + tmpReader->SetImageIO(io); + tmpReader->SetFileName(inputFile.c_str()); + try + { + tmpReader->Update(); + std::cout << "DONE READING INPUT IMAGE" << std::endl; + } + catch (itk::ExceptionObject & e) + { + std::cerr << e << std::endl; + return EXIT_FAILURE; + } + reference_image = tmpReader->GetOutput(); + } + else + { + // Generate a random reference_image. + reference_image = itkMGHImageIOTestGenerateRandomImage(size); + } + + // Write, then read the reference_image. + try + { + writer->SetFileName(fn.c_str()); + if (compression == true) + { + writer->UseCompressionOn(); + } + else + { + writer->UseCompressionOff(); + } + reader->SetFileName(fn.c_str()); + // writer->SetFileName("testDebug.mhd"); + // reader->SetFileName("testDebug.mhd"); + } + catch (itk::ExceptionObject & e) + { + std::cerr << e << std::endl; + return EXIT_FAILURE; + } + + writer->SetInput(reference_image); + reference_image->Print(std::cout); + std::cout << "----------" << std::endl; + + try + { + writer->Update(); + std::cout << "DONE WRITING TEST IMAGE" << std::endl; + reader->Update(); + std::cout << "DONE READING TEST IMAGE" << std::endl; + } + catch (itk::ExceptionObject & e) + { + std::cerr << "Exception in file reader or writer " << std::endl; + std::cerr << e.GetDescription() << std::endl; + std::cerr << e.GetLocation() << std::endl; + return EXIT_FAILURE; + } + + // Print the reference_image information. + + typename ImageType::Pointer test_image = reader->GetOutput(); + test_image->Print(std::cout); + std::cout << std::endl; + + bool isFailingPixelValues = false; + bool isFailingOrigin = false; + bool isFailingSpacing = false; + bool isFailingDirection = false; + // Compare input and output images. + itk::ImageRegionIterator a(reference_image, reference_image->GetRequestedRegion()); + itk::ImageRegionIterator b(test_image, test_image->GetRequestedRegion()); + for (a.GoToBegin(), b.GoToBegin(); !a.IsAtEnd(); ++a, ++b) + { + if (b.Get() != a.Get()) + { + std::cerr << "At index " << b.GetIndex() << " value " << b.Get() << " should be " << a.Get() << std::endl; + isFailingPixelValues = true; + } + } + for (int idx = 0; idx < 3; ++idx) + { + // itk::Math::FloatAlmostEqual( floatRepresentationfx1.asFloat, floatRepresentationfx2.asFloat, 4, 0.1f) + if (!itk::Math::FloatAlmostEqual(test_image->GetSpacing()[idx], reference_image->GetSpacing()[idx], 4, 1e-7)) + { + isFailingSpacing = true; + } + if (!itk::Math::FloatAlmostEqual(test_image->GetOrigin()[idx], reference_image->GetOrigin()[idx], 100000, 1e-1)) + { + isFailingOrigin = true; + } + for (int idx2 = 0; idx2 < 3; ++idx2) + { + if (!itk::Math::FloatAlmostEqual( + test_image->GetDirection()[idx][idx2], reference_image->GetDirection()[idx][idx2], 4, 1e-7)) + { + isFailingDirection = true; + } + } + } + std::cerr << std::fixed << std::setprecision(10) << std::endl; + // Report failure dianostics + if (isFailingSpacing) + { + std::cerr << "ERROR: Invalid Spacing: " << test_image->GetSpacing() << " ! = " << reference_image->GetSpacing() + << std::endl; + } + if (isFailingOrigin) + { + std::cerr << "ERROR: Invalid Origin: " << test_image->GetOrigin() << " ! = " << reference_image->GetOrigin() + << std::endl; + } + if (isFailingDirection) + { + std::cerr << "ERROR: Invalid Direction: \n" + << test_image->GetDirection() << " ! = \n" + << reference_image->GetDirection() << std::endl; + } + return !(isFailingPixelValues || isFailingOrigin || isFailingSpacing || isFailingDirection); +} + +#endif // itkMGHImageIOTest_h_ diff --git a/Modules/IO/MGHIO/wrapping/CMakeLists.txt b/Modules/IO/MGHIO/wrapping/CMakeLists.txt new file mode 100644 index 000000000000..c2bacc18dc8e --- /dev/null +++ b/Modules/IO/MGHIO/wrapping/CMakeLists.txt @@ -0,0 +1,3 @@ +itk_wrap_module(MGHIO) +itk_auto_load_submodules() +itk_end_wrap_module() diff --git a/Modules/IO/MGHIO/wrapping/itkMGHImageIO.wrap b/Modules/IO/MGHIO/wrapping/itkMGHImageIO.wrap new file mode 100644 index 000000000000..e43d19a9a403 --- /dev/null +++ b/Modules/IO/MGHIO/wrapping/itkMGHImageIO.wrap @@ -0,0 +1,2 @@ +itk_wrap_simple_class("itk::MGHImageIO" POINTER) +itk_wrap_simple_class("itk::MGHImageIOFactory" POINTER) diff --git a/Modules/Remote/MGHIO.remote.cmake b/Modules/Remote/MGHIO.remote.cmake deleted file mode 100644 index 9d91c36b707f..000000000000 --- a/Modules/Remote/MGHIO.remote.cmake +++ /dev/null @@ -1,50 +0,0 @@ -#-- # Grading Level Criteria Report -#-- EVALUATION DATE: 2020-03-01 -#-- EVALUATORS: [<>,<>] -#-- -#-- ## Compliance level 5 star (AKA ITK main modules, or remote modules that could become core modules) -#-- - [ ] Widespread community dependance -#-- - [ ] Above 90% code coverage -#-- - [ ] CI dashboards and testing monitored rigorously -#-- - [ ] Key API features are exposed in wrapping interface -#-- - [ ] All requirements of Levels 4,3,2,1 -#-- -#-- ## Compliance Level 4 star (Very high-quality code, perhaps small community dependance) -#-- - [ ] Meets all ITK code style standards -#-- - [ ] No external requirements beyond those needed by ITK proper -#-- - [ ] Builds and passes tests on all supported platforms within 1 month of each core tagged release -#-- - [ ] Windows Shared Library Build with Visual Studio -#-- - [ ] Mac with clang compiller -#-- - [ ] Linux with gcc compiler -#-- - [ ] Active developer community dedicated to maintaining code-base -#-- - [ ] 75% code coverage demonstrated for testing suite -#-- - [ ] Continuous integration testing performed -#-- - [ ] All requirements of Levels 3,2,1 -#-- -#-- ## Compliance Level 3 star (Quality beta code) -#-- - [ ] API | executable interface is considered mostly stable and feature complete -#-- - [ ] 10% C0-code coverage demonstrated for testing suite -#-- - [ ] Some tests exist and pass on at least some platform -#-- - [X] All requirements of Levels 2,1 -#-- -#-- ## Compliance Level 2 star (Alpha code feature API development or niche community/execution environment dependance ) -#-- - [X] Compiles for at least 1 niche set of execution envirionments, and perhaps others -#-- (may depend on specific external tools like a java environment, or specific external libraries to work ) -#-- - [X] All requirements of Levels 1 -#-- -#-- ## Compliance Level 1 star (Pre-alpha features under development and code of unknown quality) -#-- - [X] Code complies on at least 1 platform -#-- -#-- ## Compliance Level 0 star ( Code/Feature of known poor-quality or deprecated status ) -#-- - [ ] Code reviewed and explicitly identified as not recommended for use -#-- -#-- ### Please document here any justification for the criteria above -# Code style enforced by clang-format on 2020-02-19, and clang-tidy modernizations completed - -itk_fetch_module( - MGHIO - "MGHIO ImageIO plugin for ITK" - MODULE_COMPLIANCE_LEVEL 2 - GIT_REPOSITORY https://github.com/InsightSoftwareConsortium/itkMGHImageIO.git - GIT_TAG ef3756db3016378688c81ab7144126d9a7eb0904 - ) diff --git a/Utilities/Maintenance/RemoteModuleIngest/normalize-ingest-commits.py b/Utilities/Maintenance/RemoteModuleIngest/normalize-ingest-commits.py index b29288440d8e..3308386a379e 100755 --- a/Utilities/Maintenance/RemoteModuleIngest/normalize-ingest-commits.py +++ b/Utilities/Maintenance/RemoteModuleIngest/normalize-ingest-commits.py @@ -124,12 +124,28 @@ def normalize_message(subj: str, body: str) -> tuple[str, str, bool]: def replay(base: str, *, dry_run: bool, run_pre_commit: bool) -> int: - out = required(["git", "rev-list", "--reverse", f"{base}..HEAD"]) + # Skip merge commits. ``git cherry-pick`` rejects them without ``-m N``, + # and even with ``-m 1`` a merge that has been linearized by replaying + # its non-merge ancestors would re-apply that content. The non-merge + # commits we keep in this list collectively reproduce the final tree + # without the merge nodes; ``git blame`` after the merge boundary + # still walks back to original authors because authorship is preserved. + out = required(["git", "rev-list", "--reverse", "--no-merges", f"{base}..HEAD"]) commits = out.split() if not commits: print("Nothing to do — branch contains no commits beyond base.") return 0 + # For visibility: report how many merges were dropped from the replay. + all_count = len(required(["git", "rev-list", "--reverse", f"{base}..HEAD"]).split()) + merge_count = all_count - len(commits) + if merge_count: + print( + f"Skipping {merge_count} merge commit(s); " + f"their content arrives via the {len(commits)} non-merge commits.", + file=sys.stderr, + ) + print(f"Replaying {len(commits)} commits onto {base}...", file=sys.stderr) if dry_run: @@ -162,31 +178,40 @@ def replay(base: str, *, dry_run: bool, run_pre_commit: bool) -> int: for sha in commits: meta = required( - ["git", "show", "-s", "--format=%an%x00%ae%x00%aI%x00%s%x00%b", sha] + ["git", "show", "-s", "--format=%an%x00%ae%x00%aI%x00%P%x00%s%x00%b", sha] ) - an, ae, ad, subj, body = meta.split("\x00", 4) + an, ae, ad, parents, subj, body = meta.split("\x00", 5) + # Defense in depth: ``--no-merges`` should have stripped these + # already, but if a merge commit ever reaches here, treat it as + # mainline-relative so cherry-pick succeeds rather than crashing. + is_merge = len(parents.split()) > 1 # ``-X theirs`` biases toward the cherry-picked commit's content # whenever pre-commit's auto-fix on an earlier commit clashes # with the base state this commit expects. Pre-commit re-runs # below and re-normalizes the merged tree. - cp = run( - [ - "git", - "cherry-pick", - "--allow-empty", - "--no-commit", - "--strategy=recursive", - "-X", - "theirs", - sha, - ] - ) + cp_args = [ + "git", + "cherry-pick", + "--allow-empty", + "--no-commit", + "--strategy=recursive", + "-X", + "theirs", + ] + if is_merge: + cp_args += ["-m", "1"] + cp_args.append(sha) + cp = run(cp_args) if cp.returncode != 0: # Fall back to default 3-way merge — sometimes a true # content conflict needs human review. run(["git", "cherry-pick", "--abort"]) - cp2 = run(["git", "cherry-pick", "--allow-empty", "--no-commit", sha]) + cp2_args = ["git", "cherry-pick", "--allow-empty", "--no-commit"] + if is_merge: + cp2_args += ["-m", "1"] + cp2_args.append(sha) + cp2 = run(cp2_args) if cp2.returncode != 0: sys.stderr.write(f"cherry-pick failed for {sha[:10]}:\n{cp2.stderr}") run(["git", "cherry-pick", "--abort"]) diff --git a/pyproject.toml b/pyproject.toml index 94c367e0ffbb..75c5d8149738 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -55,6 +55,7 @@ cmd = '''cmake -DModule_AnisotropicDiffusionLBR:BOOL=ON -DModule_Montage:BOOL=ON -DModule_GenericLabelInterpolator:BOOL=ON + -DModule_MGHIO:BOOL=ON -DITK_COMPUTER_MEMORY_SIZE:STRING=11 -DModule_StructuralSimilarity:BOOL=ON''' description = "Configure ITK for CI (with ccache compiler launcher)"