Skip to content

Commit

Permalink
BUG: Fix read/write of displacement field in NIFTI file format
Browse files Browse the repository at this point in the history 
NIFTI image containing displacement field (intent = 1006) was incorrectly read as scalar volume (all 3 vector components were set to the same value).

Fixed by reading NIFTI displacement vector files similarly to general-purpose vector images. For displacement vector files,
vector components are converted between NIFTI file's RAS coordinate system and ITK's LPS coordinate system,
unless this conversion is explicitly disabled by setting ConvertRASDisplacementVectors to false.

An ITK vector image can be written as NIFTI displacement vector file by setting "intent_code"="1006" in the metadata dictionary of the image.
  • Loading branch information
@lassoan
lassoan committed Jan 25, 2023
1 parent 17b0fde commit d52a682
Show file tree
Hide file tree
Showing 3 changed files with 173 additions and 9 deletions.
35 changes: 35 additions & 0 deletions Modules/IO/NIFTI/include/itkNiftiImageIO.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,6 +188,37 @@ class ITKIONIFTI_EXPORT NiftiImageIO : public ImageIOBase
itkSetMacro(LegacyAnalyze75Mode, NiftiImageIOEnums::Analyze75Flavor);
itkGetConstMacro(LegacyAnalyze75Mode, NiftiImageIOEnums::Analyze75Flavor);

/** Enable conversion of vector coordinates between RAS coordinate system (in NIFTI file) and
* LPS (ITK convention) when reading or writing a generic vector image file (intent code: 1007).
*
* This flag is disabled by default because vectors may store non-spatial information.
*
* By default an ITK vector image is written to NIFTI file as a generic vector image, unless
* "intent_code" field is set explicitly to set to "1006" (displacement vector).
*/
itkSetMacro(ConvertRASVectors, bool);
itkGetConstMacro(ConvertRASVectors, bool);
itkBooleanMacro(ConvertRASVectors);

/** Enable conversion of vector coordinates between RAS coordinate system (in NIFTI file) and
* LPS (ITK convention) when reading or writing a "displacement vector" file (intent code: 1006).
*
* This flag is enabled by default. Disabling it may be useful to avoid unnecessary conversions
* in applications that uses RAS coordinate system.
*
* To write a vector image as displacement vector: set "intent_code" to "1006" in the metadata
* dictionary of the input image.
*
\code
itk::MetaDataDictionary& dictionary = image->GetMetaDataDictionary();
itk::EncapsulateMetaData<std::string>(dictionary, "intent_code", "1006");
\endcode
*
*/
itkSetMacro(ConvertRASDisplacementVectors, bool);
itkGetConstMacro(ConvertRASDisplacementVectors, bool);
itkBooleanMacro(ConvertRASDisplacementVectors);

protected:
NiftiImageIO();
~NiftiImageIO() override;
Expand Down Expand Up @@ -238,6 +269,10 @@ class ITKIONIFTI_EXPORT NiftiImageIO : public ImageIOBase
double m_RescaleSlope{ 1.0 };
double m_RescaleIntercept{ 0.0 };

bool m_ConvertRAS{ false };
bool m_ConvertRASVectors{ false };
bool m_ConvertRASDisplacementVectors{ true };

IOComponentEnum m_OnDiskComponentType{ IOComponentEnum::UNKNOWNCOMPONENTTYPE };

NiftiImageIOEnums::Analyze75Flavor m_LegacyAnalyze75Mode{};
Expand Down
118 changes: 110 additions & 8 deletions Modules/IO/NIFTI/src/itkNiftiImageIO.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -478,6 +478,9 @@ NiftiImageIO::PrintSelf(std::ostream & os, Indent indent) const
os << indent << "RescaleIntercept: " << this->m_RescaleIntercept << std::endl;
os << indent << "OnDiskComponentType: " << this->m_OnDiskComponentType << std::endl;
os << indent << "LegacyAnalyze75Mode: " << this->m_LegacyAnalyze75Mode << std::endl;
os << indent << "ConvertRASVectors: " << m_ConvertRASVectors << std::endl;
os << indent << "ConvertRASDisplacementVectors: " << m_ConvertRASDisplacementVectors << std::endl;
os << indent << "ConvertRAS: " << m_ConvertRAS << std::endl;
}

bool
Expand Down Expand Up @@ -521,6 +524,36 @@ CastCopy(float * to, const void * from, size_t pixelcount)
}
}

// Internal function to convert vectors between RAS and LPS coordinate systems.
// Dimensions are CXYZT (ITK memory layout)
template <typename TBuffer>
void
ConvertRASToFromLPS_CXYZT(TBuffer * buffer, size_t size)
{
size_t numberOfVectors = size / 3;
for (size_t i = 0; i < numberOfVectors; ++i)
{
buffer[0] *= -1;
buffer[1] *= -1;
buffer += 3;
}
}

// Internal function to convert vectors between RAS and LPS coordinate systems.
// Dimensions are XYZTC (NIFTI memory layout)
template <typename TBuffer>
void
ConvertRASToFromLPS_XYZTC(TBuffer * buffer, size_t size)
{
// Flip the sign of the first two components (L<->R, P<->A)
// and keep the third component (S) unchanged.
size_t numberOfComponents = size / 3 * 2;
for (size_t i = 0; i < numberOfComponents; ++i)
{
buffer[i] *= -1;
}
}

void
NiftiImageIO::Read(void * buffer)
{
Expand Down Expand Up @@ -791,6 +824,27 @@ NiftiImageIO::Read(void * buffer)
}
}
}

if (this->m_ConvertRAS)
{
if (this->GetPixelType() != IOPixelEnum::VECTOR && this->GetPixelType() != IOPixelEnum::POINT)
{
itkExceptionMacro(<< "RAS conversion requires pixel to be 3-component vector or point. Current pixel type is "
<< numComponents << "-component " << this->GetPixelType() << ".");
}
switch (this->m_ComponentType)
{
case IOComponentEnum::FLOAT:
ConvertRASToFromLPS_CXYZT(static_cast<float *>(buffer), numElts * numComponents);
break;
case IOComponentEnum::DOUBLE:
ConvertRASToFromLPS_CXYZT(static_cast<double *>(buffer), numElts * numComponents);
break;
default:
itkExceptionMacro(<< "RAS conversion of datatype " << this->GetComponentTypeAsString(this->m_ComponentType)
<< " is not supported");
}
}
}

NiftiImageIOEnums::NiftiFileEnum
Expand Down Expand Up @@ -1043,7 +1097,8 @@ NiftiImageIO::ReadImageInformation()
// Check the intent code, it is a vector image, or matrix image, then this is
// not true.
//
if (this->m_NiftiImage->intent_code == NIFTI_INTENT_VECTOR ||
if (this->m_NiftiImage->intent_code == NIFTI_INTENT_DISPVECT ||
this->m_NiftiImage->intent_code == NIFTI_INTENT_VECTOR ||
this->m_NiftiImage->intent_code == NIFTI_INTENT_SYMMATRIX)
{
if (this->m_NiftiImage->dim[4] > 1)
Expand Down Expand Up @@ -1085,7 +1140,8 @@ NiftiImageIO::ReadImageInformation()
this->SetNumberOfComponents(1);
}

if (this->m_NiftiImage->intent_code == NIFTI_INTENT_VECTOR ||
if (this->m_NiftiImage->intent_code == NIFTI_INTENT_DISPVECT ||
this->m_NiftiImage->intent_code == NIFTI_INTENT_VECTOR ||
this->m_NiftiImage->intent_code == NIFTI_INTENT_SYMMATRIX)
{
this->SetNumberOfComponents(this->m_NiftiImage->dim[5]);
Expand Down Expand Up @@ -1187,13 +1243,19 @@ NiftiImageIO::ReadImageInformation()

// there are a wide variety of intents we ignore
// but a few wee need to care about
this->m_ConvertRAS = false;
switch (this->m_NiftiImage->intent_code)
{
case NIFTI_INTENT_SYMMATRIX:
this->SetPixelType(IOPixelEnum::SYMMETRICSECONDRANKTENSOR);
break;
case NIFTI_INTENT_DISPVECT:
this->SetPixelType(IOPixelEnum::VECTOR);
this->m_ConvertRAS = m_ConvertRASDisplacementVectors;
break;
case NIFTI_INTENT_VECTOR:
this->SetPixelType(IOPixelEnum::VECTOR);
this->m_ConvertRAS = m_ConvertRASVectors;
break;
case NIFTI_INTENT_NONE:
case NIFTI_INTENT_CORREL:
Expand Down Expand Up @@ -1223,7 +1285,6 @@ NiftiImageIO::ReadImageInformation()
case NIFTI_INTENT_LABEL:
case NIFTI_INTENT_NEURONAME:
case NIFTI_INTENT_GENMATRIX:
case NIFTI_INTENT_DISPVECT:
case NIFTI_INTENT_POINTSET:
case NIFTI_INTENT_TRIANGLE:
case NIFTI_INTENT_QUATERNION:
Expand Down Expand Up @@ -1520,14 +1581,16 @@ NiftiImageIO::WriteImageInformation()
}
const unsigned int numComponents = this->GetNumberOfComponents();

MetaDataDictionary & thisDic = this->GetMetaDataDictionary();

// TODO: Also need to check for RGB images where numComponets=3
if (numComponents > 1 && !(this->GetPixelType() == IOPixelEnum::COMPLEX && numComponents == 2) &&
!(this->GetPixelType() == IOPixelEnum::RGB && numComponents == 3) &&
!(this->GetPixelType() == IOPixelEnum::RGBA && numComponents == 4))
{
this->m_NiftiImage->ndim = 5; // This must be 5 for NIFTI_INTENT_VECTOR
this->m_NiftiImage->ndim = 5; // This must be 5 for NIFTI_INTENT_VECTOR and NIFTI_INTENT_DISPVECT
// images.
this->m_NiftiImage->dim[0] = 5; // This must be 5 for NIFTI_INTENT_VECTOR
this->m_NiftiImage->dim[0] = 5; // This must be 5 for NIFTI_INTENT_VECTOR and NIFTI_INTENT_DISPVECT
// images.
if (this->GetNumberOfDimensions() > 4)
{
Expand All @@ -1543,7 +1606,21 @@ NiftiImageIO::WriteImageInformation()
}
else
{
this->m_NiftiImage->intent_code = NIFTI_INTENT_VECTOR;
// Each voxel is a vector. Set intent code to NIFTI_INTENT_VECTOR (default)
// or NIFTI_INTENT_DISPVECT.
int intentCode = NIFTI_INTENT_VECTOR;
std::string intentCodeStrInMetaData;
if (itk::ExposeMetaData<std::string>(thisDic, "intent_code", intentCodeStrInMetaData))
{
std::istringstream is(intentCodeStrInMetaData);
int intentCodeInMetaData = -1;
is >> intentCodeInMetaData;
if (intentCodeInMetaData == NIFTI_INTENT_DISPVECT)
{
intentCode = NIFTI_INTENT_DISPVECT;
}
}
this->m_NiftiImage->intent_code = intentCode;
}
this->m_NiftiImage->nu = this->m_NiftiImage->dim[5] = this->GetNumberOfComponents();
if (this->GetNumberOfDimensions() < 4)
Expand Down Expand Up @@ -1708,8 +1785,7 @@ NiftiImageIO::WriteImageInformation()
// TODO: Note both arguments are the same, no need to distinguish between them.
this->SetNIfTIOrientationFromImageIO(this->GetNumberOfDimensions(), this->GetNumberOfDimensions());

MetaDataDictionary & thisDic = this->GetMetaDataDictionary();
std::string temp;
std::string temp;
if (itk::ExposeMetaData<std::string>(thisDic, "aux_file", temp))
{
if (temp.length() > 23)
Expand All @@ -1721,6 +1797,10 @@ NiftiImageIO::WriteImageInformation()
strcpy(this->m_NiftiImage->aux_file, temp.c_str());
}
}

// Enable RAS conversion based on metadata and flags
this->m_ConvertRAS = (m_ConvertRASVectors && this->m_NiftiImage->intent_code == NIFTI_INTENT_VECTOR) ||
(m_ConvertRASDisplacementVectors && this->m_NiftiImage->intent_code == NIFTI_INTENT_DISPVECT);
}

namespace
Expand Down Expand Up @@ -2329,6 +2409,28 @@ NiftiImageIO::Write(const void * buffer)
}
}
}

if (this->m_ConvertRAS)
{
if (this->GetPixelType() != IOPixelEnum::VECTOR && this->GetPixelType() != IOPixelEnum::POINT)
{
itkExceptionMacro(<< "RAS conversion requires pixel to be 3-component vector or point. Current pixel type is "
<< numComponents << "-component " << this->GetPixelType() << ".");
}
switch (this->m_ComponentType)
{
case IOComponentEnum::FLOAT:
ConvertRASToFromLPS_XYZTC(reinterpret_cast<float *>(nifti_buf.get()), numComponents * seriesdist);
break;
case IOComponentEnum::DOUBLE:
ConvertRASToFromLPS_XYZTC(reinterpret_cast<double *>(nifti_buf.get()), numComponents * seriesdist);
break;
default:
itkExceptionMacro(<< "RAS conversion of datatype " << this->GetComponentTypeAsString(this->m_ComponentType)
<< " is not supported");
}
}

delete[] vecOrder;
dumpdata(buffer);
// Need a const cast here so that we don't have to copy the memory for
Expand Down
29 changes: 28 additions & 1 deletion Modules/IO/NIFTI/test/itkNiftiImageIOTest3.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ Decrement(ScalarType & value, std::enable_if_t<std::numeric_limits<ScalarType>::

template <typename ScalarType, unsigned int TVecLength, unsigned int TDimension>
int
TestImageOfVectors(const std::string & fname)
TestImageOfVectors(const std::string & fname, const std::string & intentCode = "")
{
constexpr int dimsize = 2;
/** Deformation field pixel type. */
Expand Down Expand Up @@ -143,6 +143,11 @@ TestImageOfVectors(const std::string & fname)
}
}
}
if (!intentCode.empty())
{
itk::MetaDataDictionary & dictionary = vi->GetMetaDataDictionary();
itk::EncapsulateMetaData<std::string>(dictionary, "intent_code", intentCode);
}
try
{
itk::IOTestHelper::WriteImage<VectorImageType, itk::NiftiImageIO>(vi, fname);
Expand Down Expand Up @@ -181,6 +186,24 @@ TestImageOfVectors(const std::string & fname)
return EXIT_FAILURE;
}
bool same = true;
if (!intentCode.empty())
{
const itk::MetaDataDictionary & dictionary = readback->GetMetaDataDictionary();
std::string readIntentCode;
if (itk::ExposeMetaData<std::string>(dictionary, "intent_code", readIntentCode))
{
if (readIntentCode != intentCode)
{
std::cout << "intent_code is different: " << readIntentCode << " != " << intentCode << std::endl;
same = false;
}
}
else
{
std::cout << "The read image should have an intent_code in its dictionary" << std::endl;
same = false;
}
}
if (readback->GetOrigin() != vi->GetOrigin())
{
std::cout << "Origin is different: " << readback->GetOrigin() << " != " << vi->GetOrigin() << std::endl;
Expand Down Expand Up @@ -295,5 +318,9 @@ itkNiftiImageIOTest3(int argc, char * argv[])
success |= TestImageOfVectors<float, 4, 4>(std::string("testVectorImage_float_4_4.nii.gz"));
success |= TestImageOfVectors<double, 3, 3>(std::string("testVectorImage_double_3_3.nii.gz"));

// Test reading/writing as displacement field (NIFTI intent code = 1006)
success |= TestImageOfVectors<double, 3, 1>(std::string("testDispacementImage_double.nii.gz"), std::string("1006"));
success |= TestImageOfVectors<float, 3, 1>(std::string("testDisplacementImage_float.nii.gz"), std::string("1006"));

return success;
}

0 comments on commit d52a682

Please sign in to comment.