From 6938879c6f6b9993ee57bdadcd7d493ffa4f53cb Mon Sep 17 00:00:00 2001 From: Brian Anderson Date: Thu, 14 May 2026 12:37:41 -0700 Subject: [PATCH 1/2] BUG: Fix missing ROI slices in RT Struct rasterisation DicomRTStruct2ImageFilter used a global extrude+stencil with SetTolerance(0), which lost entire mask slices when the input contained degenerate sub-pixel contours (LUNG1-002 slice 88 with a 0.07x0.12 mm exclusion - 1/8 of a voxel). Replaced with per-slice rasterisation; the slab is centred on the destination voxel so the +0.5 mm Z hack is no longer needed. Also fixes a pre-existing copy-paste bug in UpdateDicomItem (Y was clobbered with Z, Z was never written), and softens the strict p[2] == mZ planarity assert to a tolerant warning that snaps off-plane points to the contour plane. Co-Authored-By: Claude Opus 4.7 (1M context) --- common/clitkDicomRTStruct2ImageFilter.cxx | 113 +++++++++++++++------- common/clitkDicomRT_Contour.cxx | 53 +++++----- common/clitkDicomRT_ROI.h | 1 + 3 files changed, 101 insertions(+), 66 deletions(-) diff --git a/common/clitkDicomRTStruct2ImageFilter.cxx b/common/clitkDicomRTStruct2ImageFilter.cxx index 2e11fdc4..898b057c 100644 --- a/common/clitkDicomRTStruct2ImageFilter.cxx +++ b/common/clitkDicomRTStruct2ImageFilter.cxx @@ -19,19 +19,26 @@ #include #include +#include +#include // clitk #include "clitkDicomRTStruct2ImageFilter.h" +#include "clitkDicomRT_Contour.h" +#include "clitkDicomRT_ROI.h" #include "clitkImageCommon.h" // vtk #include +#include #include #include #include #include #include #include +#include +#include //-------------------------------------------------------------------- @@ -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 extrude=vtkSmartPointer::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 sts=vtkSmartPointer::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 & contours = mROI->GetListOfContours(); + + // Group contour indices by destination mask slice index in Z. + std::map > sliceGroups; + for (size_t ci = 0; ci < contours.size(); ++ci) { + double cz = contours[ci]->GetZ(); + int sliceIdx = static_cast(lrint((cz - origin[2]) / mSpacing[2])); + if (sliceIdx < 0 || sliceIdx > static_cast(extend[2])) continue; + sliceGroups[sliceIdx].push_back(ci); + } + + for (std::map >::const_iterator git = sliceGroups.begin(); + git != sliceGroups.end(); ++git) { + int sliceIdx = git->first; + const std::vector & idxs = git->second; + if (idxs.empty()) continue; - vtkSmartPointer stencil=vtkSmartPointer::New(); + // Combine all contour polylines for this slice. + vtkSmartPointer append = vtkSmartPointer::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 tf = vtkSmartPointer::New(); + tf->Translate(0.0, 0.0, shiftZ); + vtkSmartPointer tfFilter = vtkSmartPointer::New(); + tfFilter->SetTransform(tf); + tfFilter->SetInputConnection(append->GetOutputPort(0)); + + vtkSmartPointer extrude = vtkSmartPointer::New(); + extrude->SetInputConnection(tfFilter->GetOutputPort(0)); + extrude->SetVector(0.0, 0.0, mSpacing[2]); + + vtkSmartPointer sts = vtkSmartPointer::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 stencil = vtkSmartPointer::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 w = vtkSmartPointer::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 ImageType; diff --git a/common/clitkDicomRT_Contour.cxx b/common/clitkDicomRT_Contour.cxx index 9a5042be..255199ba 100644 --- a/common/clitkDicomRT_Contour.cxx +++ b/common/clitkDicomRT_Contour.cxx @@ -19,6 +19,7 @@ #include "clitkDicomRT_Contour.h" #include +#include #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::New(); mData->SetDataTypeToDouble(); mData->SetNumberOfPoints(mNbOfPoints); + bool mZInitialized = false; for(unsigned int i=0; iSetPoint(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::New(); mData->SetDataTypeToDouble(); mData->SetNumberOfPoints(mNbOfPoints); + bool mZInitialized = false; for(unsigned int i=0; iSetPoint(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; diff --git a/common/clitkDicomRT_ROI.h b/common/clitkDicomRT_ROI.h index b4179332..491cecd4 100644 --- a/common/clitkDicomRT_ROI.h +++ b/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 & GetListOfContours() const { return mListOfContours; } void SetTransformMatrix(vtkMatrix4x4* matrix); From 5b97c84d2e8df40b5834debfd03cadca9ecd4866 Mon Sep 17 00:00:00 2001 From: Brian Anderson Date: Thu, 14 May 2026 12:42:25 -0700 Subject: [PATCH 2/2] ci: trigger build after enabling Actions on fork Co-Authored-By: Claude Opus 4.7 (1M context)