#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();
+// }
+// }
+//}
//------------------------------------------------------------------------------
{
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();
+ itkMatRead = true;
- itk::Matrix<double, 4, 4> itkMat = clitk::ReadMatrix3D(filename);
+ 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++)
vtkGenericWarningMacro("Matrix in " << filename.c_str() << " cannot be inverted (determinant = 0)");
}
- mImage->GetTransform()->PostMultiply();
- mImage->GetTransform()->Concatenate(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();
+ }
+ }
+
}
}
//------------------------------------------------------------------------------