Skip to content

Commit

Permalink
Merge pull request #4864 from gdevenyi/MINC_IMAGE_RAS_LPS
Browse files Browse the repository at this point in the history
  • Loading branch information
@dzenanz
dzenanz authored Feb 22, 2025
2 parents a025ba0 + 18bd4b7 commit 33d1031
Show file tree
Hide file tree
Showing 13 changed files with 497 additions and 156 deletions.
6 changes: 6 additions & 0 deletions Modules/IO/MINC/include/itkMINCImageIO.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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

Expand Down
1 change: 1 addition & 0 deletions Modules/IO/MINC/src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@

set(ITKIOMINC_SRC itkMINCImageIO.cxx itkMINCImageIOFactory.cxx)

add_library(ITKIOMINC ${ITK_LIBRARY_BUILD_TYPE} ${ITKIOMINC_SRC})
Expand Down
121 changes: 77 additions & 44 deletions Modules/IO/MINC/src/itkMINCImageIO.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@

#include <memory> // For unique_ptr.


extern "C"
{
void
Expand Down Expand Up @@ -235,6 +236,7 @@ MINCImageIO::CloseVolume()

MINCImageIO::MINCImageIO()
: m_MINCPImpl(std::make_unique<MINCImageIOPImpl>())
, m_RAStoLPS(false)
{
for (auto & dimensionIndex : m_MINCPImpl->m_DimensionIndices)
{
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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<double, 3, 3>::GetIdentity();

Matrix<double, 3, 3> dir_cos{};
dir_cos.SetIdentity();
// Conversion matrix for MINC PositiveCoordinateOrientation RAS (LeftToRight, PosteriorToAnterior, InferiorToSuperior)
// to ITK PositiveCoordinateOrientation LPS (RightToLeft, AnteriorToPosterior, InferiorToSuperior)
auto RAStofromLPS = Matrix<double, 3, 3>::GetIdentity();
RAStofromLPS(0, 0) = -1.0;
RAStofromLPS(1, 1) = -1.0;
std::vector<double> dir_cos_temp(3);

Vector<double, 3> origin{};
Vector<double, 3> o_origin{};
Vector<double, 3> oOrigin{};

// minc api uses inverse order of dimensions , fastest varying are last
Vector<double, 3> sep;
Expand All @@ -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<double> _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)
{
Expand All @@ -486,52 +495,61 @@ MINCImageIO::ReadImageInformation()
sep[i - 1] = _sep;

this->SetDimensions(i - 1, static_cast<unsigned int>(_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)
Expand Down Expand Up @@ -660,6 +678,9 @@ MINCImageIO::ReadImageInformation()
const std::string classname(GetNameOfClass());
// EncapsulateMetaData<std::string>(thisDic,ITK_InputFilterName,
// classname);
// preserve information if the volume was PositiveCoordinateOrientation RAS to PositiveCoordinateOrientation LPS
// converted
EncapsulateMetaData<bool>(thisDic, "RAStoLPS", m_RAStoLPS);

// store history
size_t minc_history_length = 0;
Expand Down Expand Up @@ -944,22 +965,34 @@ MINCImageIO::WriteImageInformation()
}

// allocating dimensions
vnl_matrix<double> dircosmatrix(nDims, nDims);
dircosmatrix.set_identity();
vnl_matrix<double> directionCosineMatrix(nDims, nDims);
directionCosineMatrix.set_identity();
vnl_vector<double> origin(nDims);

// MINC stores direction cosines in PositiveCoordinateOrientation RAS
// need to convert to PositiveCoordinateOrientation LPS for ITK
vnl_matrix<double> 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<double> inverseDirectionCosines{ vnl_matrix_inverse<double>(dircosmatrix).as_matrix() };
const vnl_matrix<double> inverseDirectionCosines{ vnl_matrix_inverse<double>(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);
Expand All @@ -968,7 +1001,7 @@ MINCImageIO::WriteImageInformation()
{
if (k < nDims)
{
dir_cos[k] = dircosmatrix[i][k];
dir_cos[k] = directionCosineMatrix[i][k];
}
else
{
Expand All @@ -988,38 +1021,32 @@ 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;
// break;
// 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:
Expand Down Expand Up @@ -1059,7 +1086,6 @@ MINCImageIO::WriteImageInformation()
if (ExposeMetaData<std::string>(thisDic, "dimension_order", dimension_order))
{
// the format should be ((+|-)(X|Y|Z|V|T))*
// std::cout<<"Restoring original dimension order:"<<dimension_order<<std::endl;
if (dimension_order.length() == (minc_dimensions * 2))
{
dimorder_good = true;
Expand Down Expand Up @@ -1296,6 +1322,13 @@ MINCImageIO::WriteImageInformation()
// TODO: figure out what to do with it
}
}

// preserve information of MINC PositiveCoordinateOrientation RAS to ITK PositiveCoordinateOrientation LPS conversion
{
int tmp = (int)this->m_RAStoLPS;
miset_attr_values(m_MINCPImpl->m_Volume, MI_TYPE_INT, "itk", "RAStoLPS", 1, &tmp);
}

mifree_volume_props(hprops);
}

Expand Down
27 changes: 27 additions & 0 deletions Modules/IO/MINC/src/itkMINCImageIOConfigurePrivate.h.in
View file
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit 33d1031

Please sign in to comment.