X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=common%2FvvImageReader.cxx;h=c795cd9299697c7a34c784c1973541c7012e777d;hb=5c614f89dec012c7e8fee33c09fe7c1b955f5640;hp=4aad9b88099fed601a7a24a5aef4f2089824c45e;hpb=62c20109b437bef1b9df091dfe7728a0fdf98bbc;p=clitk.git

diff --git a/common/vvImageReader.cxx b/common/vvImageReader.cxx
index 4aad9b8..c795cd9 100644
--- a/common/vvImageReader.cxx
+++ b/common/vvImageReader.cxx
@@ -3,7 +3,7 @@
 
   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
@@ -14,7 +14,7 @@
 
   - 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
 
@@ -110,37 +110,37 @@ void vvImageReader::SetInputFilenames(const std::vector<std::string> & filenames
 
 //------------------------------------------------------------------------------
 //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;
-    }
-
-    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();
-    }
-  }
-}
+//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]);
+
+//      //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();
+//    }
+//  }
+//}
 //------------------------------------------------------------------------------
 
 
@@ -148,13 +148,21 @@ void vvImageReader::ReadNkiImageTransform()
 //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));
+  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());
   if(f.is_open()) {
     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++)
@@ -166,23 +174,28 @@ void vvImageReader::ReadMatImageTransform()
     // 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]->PostMultiply();
+    mImage->GetTransform()[0]->Concatenate(matrix);
+    mImage->GetTransform()[0]->Update();
   }
 }
 //------------------------------------------------------------------------------