Add pixel type options for clitkMergeSequence (double and unsigned char)

[clitk.git] / common / vvImageReader.cxx

diff --git a/common/vvImageReader.cxx b/common/vvImageReader.cxx
index 0205355f6fafc40acd2dc82984ce11f1412dddc4..a2fc4e964ab8edeeef1e57cb0f6a4b742c371d44 100644 (file)
--- a/common/vvImageReader.cxx
+++ b/common/vvImageReader.cxx
@@ -19,9 +19,11 @@
 #define VVIMAGEREADER_CXX
 
 #include <itkImageFileReader.h>
+#include "gdcmImageHelper.h"
 #include "vvImageReader.h"
 #include "vvImageReader.txx"
 #include "clitkTransformUtilities.h"
+#include "clitkElastix.h"
 
 //------------------------------------------------------------------------------
 vvImageReader::vvImageReader()
@@ -31,6 +33,7 @@ vvImageReader::vvImageReader()
   mLastError = "";
   mType = UNDEFINEDIMAGETYPE;
   mSlice = 0;
+  mPatientCoordinateSystem = false;
 }
 //------------------------------------------------------------------------------
 
@@ -48,6 +51,14 @@ void vvImageReader::Update()
 //------------------------------------------------------------------------------
 
 
+//------------------------------------------------------------------------------
+void vvImageReader::SetPatientCoordinateSystem(bool patientCoordinateSystem)
+{
+  mPatientCoordinateSystem = patientCoordinateSystem;
+}
+//------------------------------------------------------------------------------
+
+
 //------------------------------------------------------------------------------
 void vvImageReader::Update(LoadedImageType type)
 {
@@ -56,9 +67,12 @@ 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);
   }
@@ -77,13 +91,13 @@ void vvImageReader::Update(int dim,std::string inputPixelType, LoadedImageType t
   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;
   }
@@ -110,37 +124,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;
-    }
+//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();
+//    }
+//  }
+//}
 //------------------------------------------------------------------------------
 
 
@@ -150,19 +164,48 @@ void vvImageReader::ReadMatImageTransform()
 {
   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++)
@@ -192,9 +235,22 @@ void vvImageReader::ReadMatImageTransform()
       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();
+      }
+    }
+
   }
 }
 //------------------------------------------------------------------------------