Authors belong to:
- University of LYON http://www.universite-lyon.fr/
- - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
+ - Léon Bérard cancer center http://www.centreleonberard.fr
- CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
This software is distributed WITHOUT ANY WARRANTY; without even
- BSD See included LICENSE.txt file
- CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
- ======================================================================-====*/
+ ===========================================================================**/
#ifndef VVIMAGEREADER_CXX
#define VVIMAGEREADER_CXX
#include <itkImageFileReader.h>
+#include "gdcmImageHelper.h"
#include "vvImageReader.h"
#include "vvImageReader.txx"
#include "clitkTransformUtilities.h"
+#include "clitkElastix.h"
//------------------------------------------------------------------------------
vvImageReader::vvImageReader()
mLastError = "";
mType = UNDEFINEDIMAGETYPE;
mSlice = 0;
+ mPatientCoordinateSystem = false;
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
+//------------------------------------------------------------------------------
+void vvImageReader::SetPatientCoordinateSystem(bool patientCoordinateSystem)
+{
+ mPatientCoordinateSystem = patientCoordinateSystem;
+}
+//------------------------------------------------------------------------------
+
+
//------------------------------------------------------------------------------
void vvImageReader::Update(LoadedImageType type)
{
mLastError="Unable to read file.";
} else {
reader->SetFileName(mInputFilenames[0]);
+ gdcm::ImageHelper::SetForcePixelSpacing(true);
reader->ReadImageInformation();
if (mInputFilenames.size() > 1)
Update(reader->GetNumberOfDimensions()+1,reader->GetComponentTypeAsString(reader->GetComponentType()),type);
+ else if (reader->GetNumberOfComponents() > 1 && type != VECTORFIELD && type != VECTORFIELDWITHTIME)
+ Update(reader->GetNumberOfDimensions()+1,reader->GetComponentTypeAsString(reader->GetComponentType()),VECTORPIXELIMAGE);
else
Update(reader->GetNumberOfDimensions(),reader->GetComponentTypeAsString(reader->GetComponentType()),type);
}
switch(mDim) {
case 2:
UpdateWithDim<2>(mInputPixelType);
- break;;
+ break;
case 3:
UpdateWithDim<3>(mInputPixelType);
- break;;
+ break;
case 4:
UpdateWithDim<4>(mInputPixelType);
- break;;
+ break;
default:
std::cerr << "dimension unknown in Update ! " << std::endl;
}
//------------------------------------------------------------------------------
//Read transformation in NKI format (Xdr, transposed, cm)
-void vvImageReader::ReadNkiImageTransform()
-{
- bool bRead=false;
- std::string filename = mInputFilenames[0]+".MACHINEORIENTATION";
- if(itksys::SystemTools::FileExists(filename.c_str())){
- typedef itk::ImageFileReader< itk::Image< double, 2 > > MatrixReaderType;
- MatrixReaderType::Pointer readerTransfo = MatrixReaderType::New();
- readerTransfo->SetFileName(filename);
- try {
- bRead = true;
- readerTransfo->Update();
- } catch( itk::ExceptionObject & err ) {
- bRead=false;
- std::cerr << "Cannot read " << filename << std::endl
- << "The error is: " << err << std::endl;
- }
+//void vvImageReader::ReadNkiImageTransform()
+//{
+// bool bRead=false;
+// std::string filename = mInputFilenames[0]+".MACHINEORIENTATION";
+// if(itksys::SystemTools::FileExists(filename.c_str())){
+// typedef itk::ImageFileReader< itk::Image< double, 2 > > MatrixReaderType;
+// MatrixReaderType::Pointer readerTransfo = MatrixReaderType::New();
+// readerTransfo->SetFileName(filename);
+// try {
+// bRead = true;
+// readerTransfo->Update();
+// } catch( itk::ExceptionObject & err ) {
+// bRead=false;
+// std::cerr << "Cannot read " << filename << std::endl
+// << "The error is: " << err << std::endl;
+// }
- if (bRead) {
- //Transpose matrix (NKI format)
- for(int j=0; j<4; j++)
- for(int i=0; i<4; i++)
- mImage->GetTransform()->GetMatrix()->SetElement(j,i,readerTransfo->GetOutput()->GetBufferPointer()[4*i+j]);
+// if (bRead) {
+// //Transpose matrix (NKI format)
+// for(int j=0; j<4; j++)
+// for(int i=0; i<4; i++)
+// mImage->GetTransform()->GetMatrix()->SetElement(j,i,readerTransfo->GetOutput()->GetBufferPointer()[4*i+j]);
- //From cm to mm
- for(int i=0; i<3; i++)
- mImage->GetTransform()->GetMatrix()->SetElement(i,3,10*mImage->GetTransform()->GetMatrix()->GetElement(i,3));
+// //From cm to mm
+// for(int i=0; i<3; i++)
+// mImage->GetTransform()->GetMatrix()->SetElement(i,3,10*mImage->GetTransform()->GetMatrix()->GetElement(i,3));
- mImage->GetTransform()->Inverse();
- }
- }
-}
+// mImage->GetTransform()->Inverse();
+// }
+// }
+//}
//------------------------------------------------------------------------------
//Read transformation in ASCII format
void vvImageReader::ReadMatImageTransform()
{
- std::string filename(mInputFilenames[0]+".mat");
+ std::string filename(mInputFilenames[0]);
+ std::string ext(itksys::SystemTools::GetFilenameLastExtension(filename));
+
+ // Try a ".mat" extension
+ if (ext.length() > 0) {
+ size_t pos = filename.rfind(ext);
+ filename.replace(pos, ext.length(), ".mat");
+ }
+ else
+ filename += ".mat";
std::ifstream f(filename.c_str());
+ itk::Matrix<double, 4, 4> itkMat;
+ bool itkMatRead = false;
if(f.is_open()) {
- f.close();
-
- itk::Matrix<double, 4, 4> itkMat = clitk::ReadMatrix3D(filename);
-
+ itkMatRead = true;
+
+ itkMat.SetIdentity();
+ try {
+ itkMat = clitk::ReadMatrix3D(filename);
+ }
+ catch (itk::ExceptionObject & err) {
+ itkWarningMacro(<< "Found " << filename
+ << " but this is not a 4x4 matrix so it is ignored.");
+ itkMatRead = false;
+ }
+ }
+ f.close();
+
+ // Try a ".elx" extension
+ filename = mInputFilenames[0];
+ if (ext.length() > 0) {
+ size_t pos = filename.rfind(ext);
+ filename.replace(pos, ext.length(), ".elx");
+ }
+ else
+ filename += ".elx";
+ f.open(filename.c_str());
+ if(!itkMatRead && f.is_open()) {
+ itkMatRead = true;
+ itkMat = clitk::createMatrixFromElastixFile<3>(filename, true);
+ }
+ f.close();
+
+ if(itkMatRead) {
vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
matrix->Identity();
for(int j=0; j<4; j++)
// file is inverted wrt the transformation given in the MAT file.
// We don't really know where the inversion takes place inside
// VTK or ITK. For what we could see in VV, the transformation
- // given in the MHD file seems "more correct" than that given in
+ // given in the MHD file seems "more correct" than that given in
// the MAT file. For this reason, we invert the matrix read from
// the MAT file before concatenating it to the current transformation.
// Still, it would be nice to find out what happens exactly between
// VTK and ITK...
//
- vtkSmartPointer<vtkMatrix4x4> inv_matrix = vtkSmartPointer<vtkMatrix4x4>::New();
+ // SR: 23/06/2011
+ // We actually use vtkImageReslice which takes the inverse transform
+ // as input. So it seems more correct to inverse the matrix given by
+ // itk::ImageBase< VImageDimension >::GetDirection() than the matrix
+ // in .mat which is indeed the inverse optimized by a typical
+ // affine registration algorithm.
+ // Validated using clitkAffineTransform --transform_grid
if (matrix->Determinant() == 0)
{
vtkGenericWarningMacro("Matrix in " << filename.c_str() << " cannot be inverted (determinant = 0)");
}
- else
- matrix->Invert(*matrix, *inv_matrix);
-
- mImage->GetTransform()->PostMultiply();
- mImage->GetTransform()->Concatenate(inv_matrix);
- mImage->GetTransform()->Update();
+
+ // TODO SR and BP: check on the list of transforms and not the first only
+ mImage->GetTransform()[0]->PreMultiply();
+ mImage->GetTransform()[0]->Concatenate(matrix);
+ mImage->GetTransform()[0]->Update();
+
+ //for image sequences, apply the transform to each images of the sequence
+ if (mImage->IsTimeSequence())
+ {
+ for (unsigned i = 1 ; i<mImage->GetTransform().size() ; i++)
+ {
+ mImage->GetTransform()[i]->PreMultiply();
+ mImage->GetTransform()[i]->Concatenate(matrix);
+ mImage->GetTransform()[i]->Update();
+ }
+ }
+
}
}
//------------------------------------------------------------------------------
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff