X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=common%2FclitkElastix.h;h=bc333caef5d48d4c3f681a1cc0788ef5e394579b;hb=543b72e23ad001ac2a7743b9beacf48e2d0054ac;hp=ae5dc54a5485c07aa04a78d049adc0791f44c671;hpb=299afce0c4c4cc3524dccecf4ecffd46d7a66328;p=clitk.git

diff --git a/common/clitkElastix.h b/common/clitkElastix.h
index ae5dc54..bc333ca 100644
--- a/common/clitkElastix.h
+++ b/common/clitkElastix.h
@@ -19,6 +19,9 @@
 #ifndef clitkElastix_h
 #define clitkElastix_h
 
+#include <itkEuler3DTransform.h>
+#include <iterator>
+
 //--------------------------------------------------------------------
 namespace clitk {
 
@@ -63,80 +66,80 @@ GetValuesFromValue(const std::string & s,
 //-------------------------------------------------------------------
 template<unsigned int Dimension>
 typename itk::Matrix<double, Dimension+1, Dimension+1>
-createMatrixFromElastixFile(std::vector<std::string> & filename, bool verbose=true) {
+createMatrixFromElastixFile(std::string& filename, bool verbose=true) {
   if (Dimension != 3) {
     FATAL("Only 3D yet" << std::endl);
   }
-  typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
-
-  itk::CenteredEuler3DTransform<double>::Pointer mat = itk::CenteredEuler3DTransform<double>::New();
-  itk::CenteredEuler3DTransform<double>::Pointer previous;
-  for(uint j=0; j<filename.size(); j++) {
-
-    // Open file
-    if (verbose) std::cout << "Read elastix parameters in " << filename[j] << std::endl;
-    std::ifstream is;
-    clitk::openFileForReading(is, filename[j]);
-
-    // Check Transform
-    std::string s;
-    bool b = GetElastixValueFromTag(is, "Transform ", s);
-    if (!b) {
-      FATAL("Error must read 'Transform' in " << filename[j] << std::endl);
-    }
-    if (s != "EulerTransform") {
-      FATAL("Sorry only 'EulerTransform'" << std::endl);
-    }
+  typename itk::Matrix<double, Dimension+1, Dimension+1> matrix, init;
 
-    // FIXME check
-    //    (InitialTransformParametersFilename[j] "NoInitialTransform")
+  itk::Euler3DTransform<double>::Pointer mat = itk::Euler3DTransform<double>::New();
+  itk::Euler3DTransform<double>::Pointer previous;
 
-    // Get CenterOfRotationPoint
-    GetElastixValueFromTag(is, "CenterOfRotationPoint ", s); // space is needed
-    if (!b) {
-      FATAL("Error must read 'CenterOfRotationPoint' in " << filename[j] << std::endl);
-    }
-    std::vector<std::string> cor;
-    GetValuesFromValue(s, cor);
+  // Open file
+  if (verbose) std::cout << "Read elastix parameters in " << filename << std::endl;
+  std::ifstream is;
+  clitk::openFileForReading(is, filename);
 
-    // Get Transformparameters
-    GetElastixValueFromTag(is, "TransformParameters ", s); // space is needed
-    if (!b) {
-      FATAL("Error must read 'TransformParameters' in " << filename[j] << std::endl);
-    }
-    std::vector<std::string> results;
-    GetValuesFromValue(s, results);
-
-    // construct a stream from the string
-    itk::CenteredEuler3DTransform<double>::ParametersType p;
-    p.SetSize(9);
-    for(uint i=0; i<3; i++)
-      p[i] = atof(results[i].c_str()); // Rotation
-    for(uint i=0; i<3; i++)
-      p[i+3] = atof(cor[i].c_str()); // Centre of rotation
-    for(uint i=0; i<3; i++)
-      p[i+6] = atof(results[i+3].c_str()); // Translation
-    mat->SetParameters(p);
-
-    if (verbose) {
-      std::cout << "Rotation      (deg) : " << rad2deg(p[0]) << " " << rad2deg(p[1]) << " " << rad2deg(p[2]) << std::endl;
-      std::cout << "Center of rot (phy) : " << p[3] << " " << p[4] << " " << p[5] << std::endl;
-      std::cout << "Translation   (phy) : " << p[6] << " " << p[7] << " " << p[8] << std::endl;
-    }
+  // Check Transform
+  std::string s;
+  bool b = GetElastixValueFromTag(is, "Transform ", s);
+  if (!b) {
+    FATAL("Error must read 'Transform' in " << filename << std::endl);
+  }
+  if (s != "EulerTransform") {
+    FATAL("Sorry only 'EulerTransform'" << std::endl);
+  }
 
-    // Compose with previous if needed
-    if (j!=0) {
-      mat->Compose(previous);
-      if (verbose) {
-        std::cout << "Composed rotation      (deg) : " << rad2deg(mat->GetAngleX()) << " " << rad2deg(mat->GetAngleY()) << " " << rad2deg(mat->GetAngleZ()) << std::endl;
-        std::cout << "Composed center of rot (phy) : " << mat->GetCenter() << std::endl;
-        std::cout << "Compsoed translation   (phy) : " << mat->GetTranslation() << std::endl;
-      }
-    }
-    // previous = mat->Clone(); // ITK4
-    previous = itk::CenteredEuler3DTransform<double>::New();
-    previous->SetParameters(mat->GetParameters());
+  // Get previous
+  b = GetElastixValueFromTag(is, "InitialTransformParametersFileName ", s);
+  if(s == "NoInitialTransform")
+    init.SetIdentity();
+  else
+    init = createMatrixFromElastixFile<Dimension>(s, verbose);
+
+  // Get CenterOfRotationPoint
+  GetElastixValueFromTag(is, "CenterOfRotationPoint ", s); // space is needed
+  if (!b) {
+    FATAL("Error must read 'CenterOfRotationPoint' in " << filename << std::endl);
   }
+  std::vector<std::string> cor;
+  GetValuesFromValue(s, cor);
+  itk::Euler3DTransform<double>::CenterType c;
+  for(uint i=0; i<3; i++)
+    c[i] = atof(cor[i].c_str());
+  mat->SetCenter(c);
+
+  // Get Transformparameters
+  GetElastixValueFromTag(is, "ComputeZYX ", s); // space is needed
+  mat->SetComputeZYX( s==std::string("true") );
+
+  // Get Transformparameters
+  GetElastixValueFromTag(is, "TransformParameters ", s); // space is needed
+  if (!b) {
+    FATAL("Error must read 'TransformParameters' in " << filename << std::endl);
+  }
+  std::vector<std::string> results;
+  GetValuesFromValue(s, results);
+
+  // construct a stream from the string
+  itk::Euler3DTransform<double>::ParametersType p;
+  p.SetSize(6);
+  for(uint i=0; i<3; i++)
+    p[i] = atof(results[i].c_str()); // Rotation
+  for(uint i=0; i<3; i++)
+    p[i+3] = atof(results[i+3].c_str()); // Translation
+  mat->SetParameters(p);
+
+  if (verbose) {
+    std::cout << "Rotation      (deg) : " << rad2deg(p[0]) << " " << rad2deg(p[1]) << " " << rad2deg(p[2]) << std::endl;
+    std::cout << "Center of rot (phy) : " << c[0] << " " << c[1] << " " << c[2] << std::endl;
+    std::cout << "Translation   (phy) : " << p[3] << " " << p[4] << " " << p[5] << std::endl;
+  }
+
+  previous = itk::Euler3DTransform<double>::New();
+  previous->SetParameters(mat->GetParameters());
+  previous->SetCenter(c);
+  previous->SetComputeZYX(mat->GetComputeZYX());
 
   mat = previous;
   for(uint i=0; i<3; i++)
@@ -148,7 +151,7 @@ createMatrixFromElastixFile(std::vector<std::string> & filename, bool verbose=tr
   matrix[2][3] = mat->GetOffset()[2];
   matrix[3][3] = 1;
 
-  return matrix;
+  return matrix*init;
 }
 }
 //-------------------------------------------------------------------