//Read transformation in NKI format (Xdr, transposed, cm)
void vvImageReader::ReadNkiImageTransform()
{
- bool bRead=true;
- typedef itk::ImageFileReader< itk::Image< double, 2 > > MatrixReaderType;
- MatrixReaderType::Pointer readerTransfo = MatrixReaderType::New();
- readerTransfo->SetFileName(mInputFilenames[0]+".MACHINEORIENTATION");
- try {
- readerTransfo->Update();
- } catch( itk::ExceptionObject & err ) {
- bRead=false;
- }
-
- 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));
-
- mImage->GetTransform()->Inverse();
- mImage->UpdateReslice();
+ 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]);
+
+ //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();
+ }
}
}
//------------------------------------------------------------------------------
f.close();
itk::Matrix<double, 4, 4> itkMat = clitk::ReadMatrix3D(filename);
+
+ vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
+ matrix->Identity();
for(int j=0; j<4; j++)
for(int i=0; i<4; i++)
- mImage->GetTransform()->GetMatrix()->SetElement(j,i,itkMat[j][i]);
- mImage->UpdateReslice();
+ matrix->SetElement(j,i,itkMat[j][i]);
+
+ // RP: 14/03/2011
+ // For some reason, the transformation matrix given in the MHD
+ // 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
+ // 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();
+ 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();
}
}
//------------------------------------------------------------------------------