From Benoit P, use clitkDicomRTStruct2Image with image with direction cosine

[clitk.git] / common / clitkDicomRTStruct2ImageFilter.cxx

diff --git a/common/clitkDicomRTStruct2ImageFilter.cxx b/common/clitkDicomRTStruct2ImageFilter.cxx
index 2e11fdc42a6ccd2a3243d8f6598b74c52cddd534..81c478e736c0c03677b1e472b555bed80effba31 100644 (file)
--- a/common/clitkDicomRTStruct2ImageFilter.cxx
+++ b/common/clitkDicomRTStruct2ImageFilter.cxx
@@ -23,6 +23,7 @@
 // clitk
 #include "clitkDicomRTStruct2ImageFilter.h"
 #include "clitkImageCommon.h"
+#include "vvImageWriter.h"
 
 // vtk
 #include <vtkVersion.h>
@@ -32,6 +33,7 @@
 #include <vtkLinearExtrusionFilter.h>
 #include <vtkMetaImageWriter.h>
 #include <vtkXMLPolyDataWriter.h>
+#include <vtkTransformPolyDataFilter.h>
 
 
 //--------------------------------------------------------------------
@@ -107,13 +109,19 @@ void clitk::DicomRTStruct2ImageFilter::SetImage(vvImage::Pointer image)
 {
   if (image->GetNumberOfDimensions() != 3) {
     std::cerr << "Error. Please provide a 3D image." << std::endl;
-    exit(0);
+    exit(EXIT_FAILURE);
   }
   mSpacing.resize(3);
   mOrigin.resize(3);
   mSize.resize(3);
   mDirection.resize(3);
-  mTransformMatrix = image->GetTransform()[0]->GetMatrix();
+  //mTransformMatrix = image->GetTransform()[0]->GetMatrix();
+  mTransformMatrix = vtkSmartPointer<vtkMatrix4x4>::New();
+  for(unsigned int i=0;i<4;i++) {
+      for(unsigned int j=0;j<4;j++) {
+          mTransformMatrix->SetElement(i,j,image->GetTransform()[0]->GetMatrix()->GetElement(i,j));
+      }
+  }
   for(unsigned int i=0; i<3; i++) {
     mSpacing[i] = image->GetSpacing()[i];
     mOrigin[i] = image->GetOrigin()[i];
@@ -132,7 +140,7 @@ void clitk::DicomRTStruct2ImageFilter::SetImageFilename(std::string f)
   itk::ImageIOBase::Pointer header = clitk::readImageHeader(f);
   if (header->GetNumberOfDimensions() < 3) {
     std::cerr << "Error. Please provide a 3D image instead of " << f << std::endl;
-    exit(0);
+    exit(EXIT_FAILURE);
   }
   if (header->GetNumberOfDimensions() > 3) {
     std::cerr << "Warning dimension > 3 are ignored" << std::endl;
@@ -149,6 +157,18 @@ void clitk::DicomRTStruct2ImageFilter::SetImageFilename(std::string f)
     for(unsigned int j=0; j<3; j++)
       mDirection[i][j] = header->GetDirection(i)[j];
   }
+  //cf. AddItkImage function in vvImage.txx
+  mTransformMatrix = vtkSmartPointer<vtkMatrix4x4>::New();
+  mTransformMatrix->Identity();
+  for(unsigned int i=0; i<3; i++) {
+      double tmp = 0;
+      for(unsigned int j=0; j<3; j++) {
+          mTransformMatrix->SetElement(i,j,mDirection[i][j]);
+          tmp -= mDirection[i][j] * mOrigin[j];
+      }
+      tmp += mOrigin[i];
+      mTransformMatrix->SetElement(i,3,tmp);
+  }
 }
 //--------------------------------------------------------------------
 
@@ -182,11 +202,11 @@ void clitk::DicomRTStruct2ImageFilter::Update()
 {
   if (!mROI) {
     std::cerr << "Error. No ROI set, please use SetROI." << std::endl;
-    exit(0);
+    exit(EXIT_FAILURE);
   }
   if (!ImageInfoIsSet()) {
     std::cerr << "Error. Please provide image info (spacing/origin) with SetImageFilename" << std::endl;
-    exit(0);
+    exit(EXIT_FAILURE);
   }
 
   // Get Mesh
@@ -206,7 +226,7 @@ void clitk::DicomRTStruct2ImageFilter::Update()
 
   // Get bounds
   double *bounds=mesh->GetBounds();
-
+  /*
   //Change mOrigin, mSize and mSpacing with respect to the directions
   // Spacing is influenced by input direction
   std::vector<double> tempSpacing;
@@ -238,7 +258,7 @@ void clitk::DicomRTStruct2ImageFilter::Update()
     }
     mSize[i] = lrint(tempSize[i]);
   }
-
+  */
   // Compute origin
   std::vector<double> origin;
   origin.resize(3);
@@ -259,7 +279,18 @@ void clitk::DicomRTStruct2ImageFilter::Update()
       extend[i] = mSize[i]-1;
     }
   }
-
+  //Apply the transform to the mesh
+  vtkSmartPointer<vtkTransform> outputLabelmapGeometryTransform = vtkSmartPointer<vtkTransform>::New();
+  outputLabelmapGeometryTransform->SetMatrix(mTransformMatrix);
+  // Apparently the inverse is wrong...
+  //outputLabelmapGeometryTransform->Inverse();
+  vtkSmartPointer<vtkTransformPolyDataFilter> transformPolyDataFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
+#if VTK_MAJOR_VERSION <= 5
+  transformPolyDataFilter->SetInput(mesh);
+#else
+  transformPolyDataFilter->SetInputData(mesh);
+#endif
+  transformPolyDataFilter->SetTransform(outputLabelmapGeometryTransform);
   // Create new output image
   mBinaryImage = vtkSmartPointer<vtkImageData>::New();
 #if VTK_MAJOR_VERSION <= 5
@@ -281,11 +312,7 @@ void clitk::DicomRTStruct2ImageFilter::Update()
 
   // Extrude
   vtkSmartPointer<vtkLinearExtrusionFilter> extrude=vtkSmartPointer<vtkLinearExtrusionFilter>::New();
-#if VTK_MAJOR_VERSION <= 5
-  extrude->SetInput(mesh);
-#else
-  extrude->SetInputData(mesh);
-#endif
+  extrude->SetInputConnection(transformPolyDataFilter->GetOutputPort());
   ///We extrude in the -slice_spacing direction to respect the FOCAL convention (NEEDED !)
   extrude->SetVector(0, 0, -mSpacing[2]);
 
@@ -295,11 +322,7 @@ void clitk::DicomRTStruct2ImageFilter::Update()
   //http://www.nabble.com/Bug-in-vtkPolyDataToImageStencil--td23368312.html#a23370933
   sts->SetTolerance(0);
   sts->SetInformationInput(mBinaryImage);
-#if VTK_MAJOR_VERSION <= 5
-  sts->SetInput(extrude->GetOutput());
-#else
   sts->SetInputConnection(extrude->GetOutputPort(0));
-#endif
   //sts->SetInput(mesh);
 
   vtkSmartPointer<vtkImageStencil> stencil=vtkSmartPointer<vtkImageStencil>::New();
@@ -316,22 +339,26 @@ void clitk::DicomRTStruct2ImageFilter::Update()
   stencil->ReverseStencilOn();
   stencil->Update();
 
-  /*
-  vtkSmartPointer<vtkMetaImageWriter> w = vtkSmartPointer<vtkMetaImageWriter>::New();
-  w->SetInput(stencil->GetOutput());
-  w->SetFileName("binary2.mhd");
-  w->Write();
-  */
-
   mBinaryImage->ShallowCopy(stencil->GetOutput());
 
+  vvImage::Pointer vvBinaryImage = vvImage::New();
+  vtkSmartPointer<vtkTransform> vvBinaryImageT = vtkSmartPointer<vtkTransform>::New();
+  vvBinaryImageT->SetMatrix(mTransformMatrix);
+  vvBinaryImage->AddVtkImage(mBinaryImage, vvBinaryImageT);
+
   if (mWriteOutput) {
-    typedef itk::Image<unsigned char, 3> ImageType;
-    typedef itk::VTKImageToImageFilter<ImageType> ConnectorType;
-    ConnectorType::Pointer connector = ConnectorType::New();
-    connector->SetInput(GetOutput());
-    connector->Update();
-    clitk::writeImage<ImageType>(connector->GetOutput(), mOutputFilename);
+    //typedef itk::Image<unsigned char, 3> ImageType;
+    //typedef itk::VTKImageToImageFilter<ImageType> ConnectorType;
+    //ConnectorType::Pointer connector = ConnectorType::New();
+    //connector->SetInput(GetOutput());
+    //connector->Update();
+    //clitk::writeImage<ImageType>(connector->GetOutput(), mOutputFilename);
+    vvImageWriter::Pointer writer = vvImageWriter::New();
+    writer->SetInput(vvBinaryImage);
+    if (!vvBinaryImage->GetTransform().empty())
+        writer->SetSaveTransform(true);
+    writer->SetOutputFileName(mOutputFilename);
+    writer->Update();
   }
 }
 //--------------------------------------------------------------------