BUG: Fix missing ROI slices when RT Struct contains sub-pixel contours

common/clitkDicomRTStruct2ImageFilter.cxx

@@ -19,19 +19,26 @@
 
 #include <iterator>
 #include <algorithm>
+#include <map>
+#include <cmath>
 
 // clitk
 #include "clitkDicomRTStruct2ImageFilter.h"
+#include "clitkDicomRT_Contour.h"
+#include "clitkDicomRT_ROI.h"
 #include "clitkImageCommon.h"
 
 // vtk
 #include <vtkVersion.h>
+#include <vtkAppendPolyData.h>
 #include <vtkPolyDataToImageStencil.h>
 #include <vtkSmartPointer.h>
 #include <vtkImageStencil.h>
 #include <vtkLinearExtrusionFilter.h>
 #include <vtkMetaImageWriter.h>
 #include <vtkXMLPolyDataWriter.h>
+#include <vtkTransform.h>
+#include <vtkTransformPolyDataFilter.h>
 
 
 //--------------------------------------------------------------------
@@ -279,51 +286,83 @@ void clitk::DicomRTStruct2ImageFilter::Update()
   memset(mBinaryImage->GetScalarPointer(), 0,
          mBinaryImage->GetDimensions()[0]*mBinaryImage->GetDimensions()[1]*mBinaryImage->GetDimensions()[2]*sizeof(unsigned char));
 
-  // Extrude
-  vtkSmartPointer<vtkLinearExtrusionFilter> extrude=vtkSmartPointer<vtkLinearExtrusionFilter>::New();
-#if VTK_MAJOR_VERSION <= 5
-  extrude->SetInput(mesh);
-#else
-  extrude->SetInputData(mesh);
-#endif
-  ///We extrude in the -slice_spacing direction to respect the FOCAL convention (NEEDED !)
-  extrude->SetVector(0, 0, -mSpacing[2]);
-
-  // Binarization
-  vtkSmartPointer<vtkPolyDataToImageStencil> sts=vtkSmartPointer<vtkPolyDataToImageStencil>::New();
-  //The following line is extremely important
-  //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);
+  // Per-slice rasterisation.
+  //
+  // Each Z slice's closed contours are rasterised independently so that
+  // sub-voxel / degenerate contours on one slice cannot corrupt the parity
+  // count for neighbouring slices (the original implementation appended
+  // every contour into one polydata, extruded it globally, and ran the
+  // stencil with tolerance 0 -- which dropped entire slices when a tiny
+  // sub-pixel contour was present, e.g. LUNG1-002 slice 88 of NSCLC-Radiomics).
+  //
+  // For each slice we extrude the contours symmetrically around their Z
+  // (-spacing/2, +spacing/2), so the destination voxel centre sits exactly
+  // in the middle of the slab.  The even-odd rule inside vtkPolyDataToImageStencil
+  // still produces holes for nested contours on the same slice.
+  const std::vector<clitk::DicomRT_Contour::Pointer> & contours = mROI->GetListOfContours();
+
+  // Group contour indices by destination mask slice index in Z.
+  std::map<int, std::vector<size_t> > sliceGroups;
+  for (size_t ci = 0; ci < contours.size(); ++ci) {
+    double cz = contours[ci]->GetZ();
+    int sliceIdx = static_cast<int>(lrint((cz - origin[2]) / mSpacing[2]));
+    if (sliceIdx < 0 || sliceIdx > static_cast<int>(extend[2])) continue;
+    sliceGroups[sliceIdx].push_back(ci);
+  }
+
+  for (std::map<int, std::vector<size_t> >::const_iterator git = sliceGroups.begin();
+       git != sliceGroups.end(); ++git) {
+    int sliceIdx = git->first;
+    const std::vector<size_t> & idxs = git->second;
+    if (idxs.empty()) continue;
 
-  vtkSmartPointer<vtkImageStencil> stencil=vtkSmartPointer<vtkImageStencil>::New();
+    // Combine all contour polylines for this slice.
+    vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
+    for (size_t k = 0; k < idxs.size(); ++k) {
 #if VTK_MAJOR_VERSION <= 5
-  stencil->SetStencil(sts->GetOutput());
+      append->AddInput(contours[idxs[k]]->GetMesh());
 #else
-  stencil->SetStencilConnection(sts->GetOutputPort(0));
+      append->AddInputData(contours[idxs[k]]->GetMesh());
 #endif
+    }
+    append->Update();
+
+    // Centre the future extruded slab on the destination voxel centre.
+    // Shift = (voxelCenterZ - sourceZ) - mSpacing[2]/2, then extrude by +mSpacing[2].
+    double voxelCenterZ = origin[2] + sliceIdx * mSpacing[2];
+    double sourceZ = contours[idxs[0]]->GetZ();
+    double shiftZ = (voxelCenterZ - sourceZ) - mSpacing[2] / 2.0;
+
+    vtkSmartPointer<vtkTransform> tf = vtkSmartPointer<vtkTransform>::New();
+    tf->Translate(0.0, 0.0, shiftZ);
+    vtkSmartPointer<vtkTransformPolyDataFilter> tfFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
+    tfFilter->SetTransform(tf);
+    tfFilter->SetInputConnection(append->GetOutputPort(0));
+
+    vtkSmartPointer<vtkLinearExtrusionFilter> extrude = vtkSmartPointer<vtkLinearExtrusionFilter>::New();
+    extrude->SetInputConnection(tfFilter->GetOutputPort(0));
+    extrude->SetVector(0.0, 0.0, mSpacing[2]);
+
+    vtkSmartPointer<vtkPolyDataToImageStencil> sts = vtkSmartPointer<vtkPolyDataToImageStencil>::New();
+    // Keep VTK's default tolerance (~1e-3) so boundary voxels are decided
+    // deterministically.  The old SetTolerance(0) was the trigger for the
+    // missing-slice bug when combined with degenerate sub-voxel input.
+    sts->SetInformationInput(mBinaryImage);
+    sts->SetInputConnection(extrude->GetOutputPort(0));
+
+    vtkSmartPointer<vtkImageStencil> stencil = vtkSmartPointer<vtkImageStencil>::New();
+    stencil->SetStencilConnection(sts->GetOutputPort(0));
 #if VTK_MAJOR_VERSION <= 5
-  stencil->SetInput(mBinaryImage);
+    stencil->SetInput(mBinaryImage);
 #else
-  stencil->SetInputData(mBinaryImage);
+    stencil->SetInputData(mBinaryImage);
 #endif
-  stencil->ReverseStencilOn();
-  stencil->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());
+    // Accumulate into mBinaryImage so the next slice's pass sees previous fills.
+    mBinaryImage->ShallowCopy(stencil->GetOutput());
+  }
 
   if (mWriteOutput) {
     typedef itk::Image<unsigned char, 3> ImageType;

common/clitkDicomRT_Contour.cxx

@@ -19,6 +19,7 @@
 
 #include "clitkDicomRT_Contour.h"
 #include <vtkCellArray.h>
+#include <cmath>
 
 #if GDCM_MAJOR_VERSION >= 2
 #include "gdcmAttribute.h"
@@ -82,11 +83,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
+    points[i*3+2] = p[2];
   }
 
   // Get attribute
@@ -157,24 +154,23 @@ bool clitk::DicomRT_Contour::Read(gdcm::Item * item)
   mData = vtkSmartPointer<vtkPoints>::New();
   mData->SetDataTypeToDouble();
   mData->SetNumberOfPoints(mNbOfPoints);
+  bool mZInitialized = false;
   for(unsigned int i=0; i<mNbOfPoints; i++) {
     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) {
-      DD(i);
-      DD(p[2]);
-      DD(mZ);
-      std::cout << "ERROR ! contour not in the same slice" << std::endl;
-      assert(p[2] == mZ);
+    if (!mZInitialized) {
+      mZ = p[2];
+      mZInitialized = true;
+    }
+    else if (std::abs(p[2] - mZ) > 1e-3) {
+      std::cerr << "Warning: contour point " << i << " has Z=" << p[2]
+                << ", differs from contour plane Z=" << mZ
+                << " by " << (p[2]-mZ) << " mm. Snapping to plane." << std::endl;
+      p[2] = mZ;
     }
+    mData->SetPoint(i, p);
   }
 
   return true;
@@ -208,24 +204,23 @@ bool clitk::DicomRT_Contour::Read(gdcm::SQItem * item)
   mData = vtkSmartPointer<vtkPoints>::New();
   mData->SetDataTypeToDouble();
   mData->SetNumberOfPoints(mNbOfPoints);
+  bool mZInitialized = false;
   for(unsigned int i=0; i<mNbOfPoints; i++) {
     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) {
-      DD(i);
-      DD(p[2]);
-      DD(mZ);
-      std::cout << "ERROR ! contour not in the same slice" << std::endl;
-      assert(p[2] == mZ);
+    if (!mZInitialized) {
+      mZ = p[2];
+      mZInitialized = true;
+    }
+    else if (std::abs(p[2] - mZ) > 1e-3) {
+      std::cerr << "Warning: contour point " << i << " has Z=" << p[2]
+                << ", differs from contour plane Z=" << mZ
+                << " by " << (p[2]-mZ) << " mm. Snapping to plane." << std::endl;
+      p[2] = mZ;
     }
+    mData->SetPoint(i, p);
   }
 
   return true;

common/clitkDicomRT_ROI.h

@@ -66,6 +66,7 @@ class DicomRT_ROI : public itk::LightObject
 
   void SetImage(vvImage::Pointer im);
   DicomRT_Contour* GetContour(int n);
+  const std::vector<DicomRT_Contour::Pointer> & GetListOfContours() const { return mListOfContours; }
 
   void SetTransformMatrix(vtkMatrix4x4* matrix);