X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=vv%2FvvImageReader.cxx;h=25e235a08edb3b713b347b5a5a1147aaed7fb061;hb=ff1d5b8af965e75ecdbc1b0cbff0dad6f3b3c511;hp=8bbdf251b57355a40b466065e913596e6ec0d7b9;hpb=91342e828163ad49dbed39426ed4521dbfc64629;p=clitk.git

diff --git a/vv/vvImageReader.cxx b/vv/vvImageReader.cxx
index 8bbdf25..25e235a 100644
--- a/vv/vvImageReader.cxx
+++ b/vv/vvImageReader.cxx
@@ -124,28 +124,33 @@ void vvImageReader::SetInputFilenames(const std::vector<std::string> & filenames
 //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();
+    }
   }
 }
 //------------------------------------------------------------------------------
@@ -161,10 +166,35 @@ void vvImageReader::ReadMatImageTransform()
     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();
   }
 }
 //------------------------------------------------------------------------------