Bunch of changes

common/CMakeLists.txt

@@ -72,6 +72,25 @@ endif()
 #=========================================================
 
 target_link_libraries(clitkCommon ${VTK_LIBRARIES} ${ITK_LIBRARIES})
+# I have some linking problems on linux if VTK 7.1.1 (at least) is built statically - It asks for Qt5::X11Extras - I do not know why
+# First workaround - link clitkCommon with Qt5::X11Extras - Qt5::X11Extras is not necessary if VTK is built dynamically
+# Second workaround - oblige the user to recompile VTK as a dynamic library (better solution)
+if(UNIX AND NOT APPLE)
+    if(vv_QT_VERSION VERSION_EQUAL "5") #5
+        if(VTK_VERSION VERSION_EQUAL "7.1.1") #7.1.1
+            #message(${VTK_DIR})
+            list(GET VTK_LIBRARIES 0 FIRST_VTK_ELEMENT)
+            #message(${FIRST_VTK_ELEMENT})
+            file(GLOB FIRST_VTK_LIB ${VTK_DIR}/../../*${FIRST_VTK_ELEMENT}*.a)
+            #message(${FIRST_VTK_LIB})
+            if(EXISTS ${FIRST_VTK_LIB})
+                message(FATAL_ERROR "VTK is built as a statically library - you need to recompile VTK and/or ITK as dynamic libraries")
+                #find_package(Qt5X11Extras REQUIRED)
+                #target_link_libraries(clitkCommon Qt5::X11Extras)
+            endif()
+        endif()
+    endif()
+endif()
 
 add_library(clitkDicomRTStruct STATIC
   clitkDicomRT_Contour.cxx

common/clitkCommon.h

@@ -76,7 +76,7 @@ namespace clitk {
   //--------------------------------------------------------------------
   // when everything goes wrong
 #define WHEREAMI "[ " << __FILE__  << " ] line " << __LINE__
-#define FATAL(a) { std::cerr << "ERROR in " << WHEREAMI << ": " << a; exit(0); }
+#define FATAL(a) { std::cerr << "ERROR in " << WHEREAMI << ": " << a << std::endl; exit(0); }
 
   //--------------------------------------------------------------------
   // GGO with modified struct name

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();
   }
 }
 //--------------------------------------------------------------------
@@ -341,7 +368,10 @@ void clitk::DicomRTStruct2ImageFilter::Update()
 //--------------------------------------------------------------------
 vtkImageData * clitk::DicomRTStruct2ImageFilter::GetOutput()
 {
-  assert(mBinaryImage);
+  //assert(mBinaryImage);
+  if (mBinaryImage == NULL) {
+      FATAL("The binary RTStruct image is NULL");
+  }
   return mBinaryImage;
 }
 //--------------------------------------------------------------------

common/clitkDicomRTStruct2ImageFilter.h

@@ -72,17 +72,17 @@ namespace clitk {
 
 //--------------------------------------------------------------------
 
-template <int Dimension> 
-typename itk::Image<unsigned char,Dimension>::ConstPointer clitk::DicomRTStruct2ImageFilter::GetITKOutput()
-{
-  assert(mBinaryImage);
-  typedef itk::Image<unsigned char,Dimension> ConnectorImageType;
-  typedef itk::VTKImageToImageFilter <ConnectorImageType> ConnectorType;
-  typename ConnectorType::Pointer connector = ConnectorType::New();
-  connector->SetInput(mBinaryImage);
-  connector->Update();
-  return connector->GetOutput();
-}
+//template <int Dimension>
+//typename itk::Image<unsigned char,Dimension>::ConstPointer clitk::DicomRTStruct2ImageFilter::GetITKOutput()
+//{
+//  assert(mBinaryImage);
+//  typedef itk::Image<unsigned char,Dimension> ConnectorImageType;
+//  typedef itk::VTKImageToImageFilter <ConnectorImageType> ConnectorType;
+//  typename ConnectorType::Pointer connector = ConnectorType::New();
+//  connector->SetInput(mBinaryImage);
+//  connector->Update();
+//  return connector->GetOutput();
+//}
 //--------------------------------------------------------------------
 #endif // CLITKDICOMRT_TRUCT2IMAGEFILTER_H
 

common/clitkDicomRT_Contour.cxx

@@ -82,11 +82,7 @@ void clitk::DicomRT_Contour::UpdateDicomItem()
     double * p = mData->GetPoint(i);
     points[i*3] = p[0];
     points[i*3+1] = p[1];
-#if VTK_MAJOR_VERSION <= 5
-    points[i*3+1] = p[2];
-#else
     points[i*3+1] = p[2]-0.5;
-#endif
   }
 
   // Get attribute
@@ -96,7 +92,7 @@ void clitk::DicomRT_Contour::UpdateDicomItem()
   at.SetFromDataElement( contourdata );
 
   // Set attribute
-  at.SetValues(&points[0], points.size(), false);
+  at.SetValues(&points[0], points.size());
   DD(at.GetValues()[0]);
 
   DD("replace");
@@ -151,7 +147,10 @@ bool clitk::DicomRT_Contour::Read(gdcm::Item * item)
   at.SetFromDataElement( contourdata );
   const double* points = at.GetValues();
   //  unsigned int npts = at.GetNumberOfValues() / 3;
-  assert(at.GetNumberOfValues() == static_cast<unsigned int>(mNbOfPoints)*3);
+  //assert(at.GetNumberOfValues() == static_cast<unsigned int>(mNbOfPoints)*3);
+  if (at.GetNumberOfValues() != static_cast<unsigned int>(mNbOfPoints)*3) {
+      FATAL("The number of contour points is inconsistent with the number of triplets defining the contour");
+  }
 
   // Organize values
   mData = vtkSmartPointer<vtkPoints>::New();
@@ -161,19 +160,16 @@ bool clitk::DicomRT_Contour::Read(gdcm::Item * item)
     double p[3];
     p[0] = points[i*3];
     p[1] = points[i*3+1];
-#if VTK_MAJOR_VERSION <= 5
-    p[2] = points[i*3+2];
-#else
     p[2] = points[i*3+2]+0.5;
-#endif
     mData->SetPoint(i, p);
     if (mZ == -1) mZ = p[2];
-    if (p[2] != mZ) {
+    if (std::fabs(p[2] - mZ) > mTolerance) {
       DD(i);
       DD(p[2]);
       DD(mZ);
-      std::cout << "ERROR ! contour not in the same slice" << std::endl;
-      assert(p[2] == mZ);
+      //std::cout << "ERROR ! contour not in the same slice" << std::endl;
+      //assert(p[2] == mZ);
+      FATAL("ERROR ! contour not in the same slice");
     }
   }
 
@@ -202,7 +198,10 @@ bool clitk::DicomRT_Contour::Read(gdcm::SQItem * item)
 
   // Read values [Contour Data]
   std::vector<float> points = parse_string<float>(item->GetEntryValue(0x3006,0x0050),'\\');
-  assert(points.size() == static_cast<unsigned int>(mNbOfPoints)*3);
+  //assert(points.size() == static_cast<unsigned int>(mNbOfPoints)*3);
+  if (points.size() != static_cast<unsigned int>(mNbOfPoints)*3) {
+      FATAL("The number of contour points is inconsistent with the number of triplets defining the contour");
+  }
 
   // Organize values
   mData = vtkSmartPointer<vtkPoints>::New();
@@ -212,19 +211,16 @@ bool clitk::DicomRT_Contour::Read(gdcm::SQItem * item)
     double p[3];
     p[0] = points[i*3];
     p[1] = points[i*3+1];
-#if VTK_MAJOR_VERSION <= 5
-    p[2] = points[i*3+2];
-#else
     p[2] = points[i*3+2]+0.5;
-#endif
     mData->SetPoint(i, p);
     if (mZ == -1) mZ = p[2];
-    if (p[2] != mZ) {
+    if (std::fabs(p[2] - mZ) > mTolerance) {
       DD(i);
       DD(p[2]);
       DD(mZ);
-      std::cout << "ERROR ! contour not in the same slice" << std::endl;
-      assert(p[2] == mZ);
+      //std::cout << "ERROR ! contour not in the same slice" << std::endl;
+      //assert(p[2] == mZ);
+      FATAL("ERROR ! contour not in the same slice");
     }
   }
 
@@ -259,8 +255,17 @@ void clitk::DicomRT_Contour::SetTransformMatrix(vtkMatrix4x4* matrix)
   mTransformMatrix = matrix;
 }
 //--------------------------------------------------------------------
-
-
+//--------------------------------------------------------------------
+double clitk::DicomRT_Contour::GetTolerance()
+{
+  return mTolerance;
+}
+//--------------------------------------------------------------------
+void clitk::DicomRT_Contour::SetTolerance(double tol)
+{
+  mTolerance = tol;
+}
+//--------------------------------------------------------------------
 //--------------------------------------------------------------------
 void clitk::DicomRT_Contour::ComputeMeshFromDataPoints()
 {
@@ -271,10 +276,10 @@ void clitk::DicomRT_Contour::ComputeMeshFromDataPoints()
   mMesh->SetPoints(mPoints);
   vtkIdType ids[2];
   for (unsigned int idx=0 ; idx<mNbOfPoints ; idx++) {
-    double pointIn[4];
-    for (unsigned int j=0 ; j<3; ++j)
-      pointIn[j] = mData->GetPoint(idx)[j];
-    pointIn[3] = 1.0;
+    //double pointIn[4];
+    //for (unsigned int j=0 ; j<3; ++j)
+    //  pointIn[j] = mData->GetPoint(idx)[j];
+    //pointIn[3] = 1.0;
     /*double pointOut[4];
     mTransformMatrix->MultiplyPoint(pointIn, pointOut);
     std::cout << pointOut[0] << " " << pointOut[1] << " " << pointOut[2] << " " << pointOut[3] << std::endl;