Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: Fix MINC read/write to convert from/to RAS coordinate system #147

Closed
wants to merge 5 commits into from
Closed

WIP: Fix MINC read/write to convert from/to RAS coordinate system #147

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
59b72ce
BUG: Fix MINC read/write to convert from/to RAS coordinate system
gdevenyi Nov 14, 2018
983052b
BUG: Fix TransformMINC read/write from/to RAS coordinate system
gdevenyi Nov 14, 2018
ec3cff4
BUG: Fix grid TransformMINC read/write from/to RAS coordinate system
gdevenyi Nov 26, 2018
b8e93ae
PERF: Remove unndeeded matrix fills
gdevenyi Nov 30, 2018
666f747
BUG: Default to backwards compaitability
gdevenyi Dec 5, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions Modules/IO/MINC/include/itkMINCImageIO.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,11 @@ class ITKIOMINC_EXPORT MINCImageIO : public ImageIOBase
void SetCompressionLevel(int level);
int GetCompressionLevel() const;

/** Set the conversion between RAS and LPS coordinate systems.
Backwards compatibility feature for ezMINC tool */
itkSetMacro(ConvertCoordinatesToLPS, bool);
itkGetConstMacro(ConvertCoordinatesToLPS, bool);

/*-------- This part of the interface deals with reading data. ------ */

/** Determine the file type. Returns true if this ImageIO can read the
Expand Down Expand Up @@ -139,6 +144,8 @@ class ITKIOMINC_EXPORT MINCImageIO : public ImageIOBase

MINCImageIOPImpl *m_MINCPImpl;

bool m_ConvertCoordinatesToLPS;

MatrixType m_DirectionCosines;

// complex type images, composed of complex numbers
Expand Down
30 changes: 30 additions & 0 deletions Modules/IO/MINC/src/itkMINCImageIO.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -252,6 +252,7 @@ MINCImageIO::MINCImageIO()
}

this->m_UseCompression = true;
this->m_ConvertCoordinatesToLPS = false;
this->m_MINCPImpl->m_CompressionLevel = 4; // Range 0-9; 0 = no file compression, 9 =
// maximum file compression
this->m_MINCPImpl->m_Volume_type = MI_TYPE_FLOAT;
Expand Down Expand Up @@ -456,6 +457,13 @@ void MINCImageIO::ReadImageInformation()
dir_cos.Fill(0.0);
dir_cos.SetIdentity();

//MINC stores direction cosines in RAS, need to convert to LPS for ITK
Matrix< double, 3,3 > RAS_tofrom_LPS;
RAS_tofrom_LPS.SetIdentity();
RAS_tofrom_LPS(0,0) = -1.0;
RAS_tofrom_LPS(1,1) = -1.0;
std::vector< double > dir_cos_temp(3);

Vector< double,3> origin,sep;
Vector< double,3> o_origin;
origin.Fill(0.0);
Expand Down Expand Up @@ -497,6 +505,9 @@ void MINCImageIO::ReadImageInformation()
}
}

if(this->m_ConvertCoordinatesToLPS)
dir_cos = RAS_tofrom_LPS*dir_cos;

if(this->m_MINCPImpl->m_DimensionIndices[0]!=-1) // have vector dimension
{
//micopy_dimension(this->m_MINCPImpl->m_MincFileDims[this->m_MINCPImpl->m_DimensionIndices[0]],&apparent_dimension_order[usable_dimensions]);
Expand Down Expand Up @@ -530,7 +541,17 @@ void MINCImageIO::ReadImageInformation()
o_origin=dir_cos*origin;

for(int i=0; i<spatial_dimension_count; i++)
{
this->SetOrigin(i,o_origin[i]);
if(this->m_ConvertCoordinatesToLPS)
{
for( unsigned int j=0; j<3; j++ )
{
dir_cos_temp[j] = dir_cos[j][i];
}
this->SetDirection(i,dir_cos_temp);
}
}

miclass_t volume_data_class;

Expand Down Expand Up @@ -924,6 +945,12 @@ void MINCImageIO::WriteImageInformation()
dircosmatrix.set_identity();
vnl_vector<double> origin(nDims);

//MINC stores direction cosines in RAS, need to convert to 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++ )
Expand All @@ -936,6 +963,9 @@ void MINCImageIO::WriteImageInformation()
vnl_matrix< double > inverseDirectionCosines = vnl_matrix_inverse< double >(dircosmatrix);
origin *= inverseDirectionCosines; //transform to minc convention

if(this->m_ConvertCoordinatesToLPS)
dircosmatrix *= RAS_tofrom_LPS;

for (unsigned int i = 0; i < nDims; i++ )
{
unsigned int j=i+(nComp>1 ? 1 : 0);
Expand Down
124 changes: 115 additions & 9 deletions Modules/IO/TransformMINC/include/itkMINCTransformIO.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,8 @@
#include "itkVector.h"
#include "itkDisplacementFieldTransform.h"
#include "itkMetaDataObject.h"
#include "itkImageRegionIterator.h"
#include "vnl/vnl_vector_fixed.h"

namespace itk
{
Expand Down Expand Up @@ -98,15 +100,44 @@ MINCTransformIOTemplate<TParametersValueType>
ParametersType parameterArray;
parameterArray.SetSize(12);


Matrix< double, 3,3 > RAS_affine_transform;
Matrix< double, 3,3 > LPS_affine_transform;
RAS_affine_transform.SetIdentity();

//MINC stores transforms in RAS, need to convert to LPS for ITK
Matrix< double, 3,3 > RAS_tofrom_LPS;
RAS_tofrom_LPS.SetIdentity();
RAS_tofrom_LPS(0,0) = -1.0;
RAS_tofrom_LPS(1,1) = -1.0;

for(int j = 0; j < 3; ++j)
{
for(int i = 0; i < 3; ++i)
{
parameterArray.SetElement(i+j*3, Transform_elem(*lin,j,i));
RAS_affine_transform(i,j) = Transform_elem(*lin,j,i);
}
parameterArray.SetElement(j+9, Transform_elem(*lin,j,3));
}

LPS_affine_transform = RAS_tofrom_LPS*RAS_affine_transform*RAS_tofrom_LPS;

for(int j = 0; j < 3; ++j)
{
for(int i = 0; i < 3; ++i)
{
parameterArray.SetElement(i+j*3, LPS_affine_transform(i,j));
}
//Here, again, RAS to LPS for the translations
if( (j == 2) )
{
parameterArray.SetElement(j+9, Transform_elem(*lin,j,3));
}
else
{
parameterArray.SetElement(j+9, -Transform_elem(*lin,j,3));
}
}

if(xfm->inverse_flag)
{
using AffineTransformType = AffineTransform<TParametersValueType, 3>;
Expand Down Expand Up @@ -143,6 +174,8 @@ MINCTransformIOTemplate<TParametersValueType>
using DisplacementFieldTransformType = DisplacementFieldTransform<TParametersValueType, 3>;
using GridImageType = typename DisplacementFieldTransformType::DisplacementFieldType;
using MincReaderType = ImageFileReader< GridImageType >;
using OutputPixelType = vnl_vector_fixed<TParametersValueType, 3>;
const OutputPixelType RAS_tofrom_LPS_vector = {-1, -1, 1};

MINCImageIO::Pointer mincIO = MINCImageIO::New();
typename MincReaderType::Pointer reader = MincReaderType::New();
Expand All @@ -152,6 +185,26 @@ MINCTransformIOTemplate<TParametersValueType>

typename GridImageType::Pointer grid = reader->GetOutput();

typename GridImageType::Pointer LPSgrid = GridImageType::New();
LPSgrid->CopyInformation(grid);
LPSgrid->SetRegions(grid->GetBufferedRegion());
LPSgrid->Allocate(true);

itk::MultiThreaderBase::Pointer mt = itk::MultiThreaderBase::New();
mt->ParallelizeImageRegion<3>(grid->GetBufferedRegion(),
[grid, LPSgrid, RAS_tofrom_LPS_vector](const typename GridImageType::RegionType & region)
{
itk::Vector<TParametersValueType, 3> p;
itk::ImageRegionConstIterator< GridImageType > iIt(grid, region);
itk::ImageRegionIterator< GridImageType > oIt(LPSgrid, region);
for (; !iIt.IsAtEnd(); ++iIt, ++oIt)
{
p.SetVnlVector(element_product(iIt.Get().GetVnlVector(),RAS_tofrom_LPS_vector));
oIt.Set(p);
}
},
nullptr);

TransformPointer transform;
std::string transformTypeName = "DisplacementFieldTransform_";
transformTypeName += typeNameString;
Expand All @@ -160,11 +213,11 @@ MINCTransformIOTemplate<TParametersValueType>
auto * gridTransform = static_cast< DisplacementFieldTransformType* >( transform.GetPointer());
if( xfm->inverse_flag ) //TODO: invert grid transform?
{
gridTransform->SetInverseDisplacementField( grid );
gridTransform->SetInverseDisplacementField( LPSgrid );
}
else
{
gridTransform->SetDisplacementField( grid );
gridTransform->SetDisplacementField( LPSgrid );
}

this->GetReadTransformList().push_back( transform );
Expand Down Expand Up @@ -231,14 +284,42 @@ MINCTransformIOTemplate<TParametersValueType>
MatrixType matrix = matrixOffsetTransform->GetMatrix();
OffsetType offset = matrixOffsetTransform->GetOffset();

//MINC stores everything in RAS, need to convert from LPS
Matrix< double, 3,3 > RAS_tofrom_LPS;
Matrix< double, 3,3 > RAS_affine_transform;
Matrix< double, 3,3 > LPS_affine_transform;
RAS_tofrom_LPS.SetIdentity();
RAS_tofrom_LPS(0,0) = -1.0;
RAS_tofrom_LPS(1,1) = -1.0;


for(int j=0; j < 3; ++j)
{
for(int i=0; i < 3; ++i)
{
Transform_elem(lin,j,i)=matrix(j,i);
LPS_affine_transform(j,i) = matrix(j,i);
}
Transform_elem(lin,j,3)=offset[j];
}

RAS_affine_transform = RAS_tofrom_LPS * LPS_affine_transform * RAS_tofrom_LPS;

for(int j=0; j < 3; ++j)
{
for(int i=0; i < 3; ++i)
{
Transform_elem(lin,j,i)=RAS_affine_transform(j,i);
}

//Here, again RAS to LPS for the translations
if( (j == 2) )
{
Transform_elem(lin,j,3)=offset[j];
}
else
{
Transform_elem(lin,j,3)=-offset[j];
}
}
//add 4th normalization row (not stored)
Transform_elem(lin,3,3)=1.0;

Expand All @@ -254,6 +335,10 @@ MINCTransformIOTemplate<TParametersValueType>
using DisplacementFieldTransformType = DisplacementFieldTransform<TParametersValueType, 3>;
using GridImageType = typename DisplacementFieldTransformType::DisplacementFieldType;
using MincWriterType = ImageFileWriter< GridImageType >;
typename GridImageType::Pointer RASgrid = GridImageType::New();
typename GridImageType::Pointer LPSgrid = GridImageType::New();
using OutputPixelType = vnl_vector_fixed<TParametersValueType, 3>;
const OutputPixelType RAS_tofrom_LPS_vector = {-1, -1, 1};
auto * _grid_transform = static_cast< DisplacementFieldTransformType* >( const_cast< TransformType* >( curTransform ));
char tmp[1024];
sprintf(tmp,"%s_grid_%d.mnc",xfm_file_base,serial);
Expand All @@ -266,17 +351,38 @@ MINCTransformIOTemplate<TParametersValueType>

if( _grid_transform->GetDisplacementField() )
{
writer->SetInput( _grid_transform->GetModifiableDisplacementField() );
LPSgrid = _grid_transform->GetModifiableDisplacementField();
}
else if( _grid_transform->GetInverseDisplacementField() )
{
writer->SetInput( _grid_transform->GetModifiableInverseDisplacementField() );
_inverse_grid=true;
LPSgrid = _grid_transform->GetModifiableInverseDisplacementField();
_inverse_grid=true;
}
else
{
itkExceptionMacro(<< "Trying to write-out displacement transform without displacement field");
}

RASgrid->CopyInformation(LPSgrid);
RASgrid->SetRegions(LPSgrid->GetBufferedRegion());
RASgrid->Allocate(true);

itk::MultiThreaderBase::Pointer mt = itk::MultiThreaderBase::New();
mt->ParallelizeImageRegion<3>(LPSgrid->GetBufferedRegion(),
[LPSgrid, RASgrid, RAS_tofrom_LPS_vector](const typename GridImageType::RegionType & region)
{
itk::Vector<TParametersValueType, 3> p;
itk::ImageRegionConstIterator< GridImageType > iIt(LPSgrid, region);
itk::ImageRegionIterator< GridImageType > oIt(RASgrid, region);
for (; !iIt.IsAtEnd(); ++iIt, ++oIt)
{
p.SetVnlVector(element_product(iIt.Get().GetVnlVector(),RAS_tofrom_LPS_vector));
oIt.Set(p);
}
},
nullptr);

writer->SetInput( RASgrid );
writer->Update();

xfm.push_back( VIO_General_transform() );
Expand Down