From: David Sarrut Date: Mon, 3 Oct 2011 06:27:21 +0000 (+0200) Subject: Merge branch 'master' of /home/dsarrut/clitk3.server X-Git-Tag: v1.3.0~206 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=5b568893e14e5d955fa14f653bd176b54c6aee0c;hp=94cb14004ceab7bca9c892a435ace9f902e661b4;p=clitk.git Merge branch 'master' of /home/dsarrut/clitk3.server --- diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index 0263f3a..9f92d71 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -49,7 +49,7 @@ ADD_LIBRARY(clitkDicomRTStruct STATIC clitkDicomRT_Contour.cxx clitkDicomRT_ROI.cxx clitkDicomRT_StructureSet.cxx - clitkDicomRT_ROI_ConvertToImageFilter.cxx + clitkDicomRTStruct2ImageFilter.cxx ) TARGET_LINK_LIBRARIES(clitkDicomRTStruct vtkHybrid) diff --git a/common/clitkDicomRT_ROI_ConvertToImageFilter.cxx b/common/clitkDicomRTStruct2ImageFilter.cxx similarity index 87% rename from common/clitkDicomRT_ROI_ConvertToImageFilter.cxx rename to common/clitkDicomRTStruct2ImageFilter.cxx index e3b3472..1541cb2 100644 --- a/common/clitkDicomRT_ROI_ConvertToImageFilter.cxx +++ b/common/clitkDicomRTStruct2ImageFilter.cxx @@ -21,7 +21,7 @@ #include // clitk -#include "clitkDicomRT_ROI_ConvertToImageFilter.h" +#include "clitkDicomRTStruct2ImageFilter.h" #include "clitkImageCommon.h" // vtk @@ -33,7 +33,7 @@ //-------------------------------------------------------------------- -clitk::DicomRT_ROI_ConvertToImageFilter::DicomRT_ROI_ConvertToImageFilter() +clitk::DicomRTStruct2ImageFilter::DicomRTStruct2ImageFilter() { mROI = NULL; mWriteOutput = false; @@ -43,19 +43,19 @@ clitk::DicomRT_ROI_ConvertToImageFilter::DicomRT_ROI_ConvertToImageFilter() //-------------------------------------------------------------------- -clitk::DicomRT_ROI_ConvertToImageFilter::~DicomRT_ROI_ConvertToImageFilter() +clitk::DicomRTStruct2ImageFilter::~DicomRTStruct2ImageFilter() { } //-------------------------------------------------------------------- -bool clitk::DicomRT_ROI_ConvertToImageFilter::ImageInfoIsSet() const +bool clitk::DicomRTStruct2ImageFilter::ImageInfoIsSet() const { return mSize.size() && mSpacing.size() && mOrigin.size(); } //-------------------------------------------------------------------- -void clitk::DicomRT_ROI_ConvertToImageFilter::SetROI(clitk::DicomRT_ROI * roi) +void clitk::DicomRTStruct2ImageFilter::SetROI(clitk::DicomRT_ROI * roi) { mROI = roi; } @@ -63,7 +63,7 @@ void clitk::DicomRT_ROI_ConvertToImageFilter::SetROI(clitk::DicomRT_ROI * roi) //-------------------------------------------------------------------- -void clitk::DicomRT_ROI_ConvertToImageFilter::SetCropMaskEnabled(bool b) +void clitk::DicomRTStruct2ImageFilter::SetCropMaskEnabled(bool b) { mCropMask = b; } @@ -71,7 +71,7 @@ void clitk::DicomRT_ROI_ConvertToImageFilter::SetCropMaskEnabled(bool b) //-------------------------------------------------------------------- -void clitk::DicomRT_ROI_ConvertToImageFilter::SetOutputImageFilename(std::string s) +void clitk::DicomRTStruct2ImageFilter::SetOutputImageFilename(std::string s) { mOutputFilename = s; mWriteOutput = true; @@ -80,7 +80,7 @@ void clitk::DicomRT_ROI_ConvertToImageFilter::SetOutputImageFilename(std::string //-------------------------------------------------------------------- -void clitk::DicomRT_ROI_ConvertToImageFilter::SetImageFilename(std::string f) +void clitk::DicomRTStruct2ImageFilter::SetImageFilename(std::string f) { itk::ImageIOBase::Pointer header = clitk::readImageHeader(f); if (header->GetNumberOfDimensions() < 3) { @@ -101,23 +101,23 @@ void clitk::DicomRT_ROI_ConvertToImageFilter::SetImageFilename(std::string f) } //-------------------------------------------------------------------- -void clitk::DicomRT_ROI_ConvertToImageFilter::SetOutputOrigin(const double* origin) +void clitk::DicomRTStruct2ImageFilter::SetOutputOrigin(const double* origin) { std::copy(origin,origin+3,std::back_inserter(mOrigin)); } //-------------------------------------------------------------------- -void clitk::DicomRT_ROI_ConvertToImageFilter::SetOutputSpacing(const double* spacing) +void clitk::DicomRTStruct2ImageFilter::SetOutputSpacing(const double* spacing) { std::copy(spacing,spacing+3,std::back_inserter(mSpacing)); } //-------------------------------------------------------------------- -void clitk::DicomRT_ROI_ConvertToImageFilter::SetOutputSize(const unsigned long* size) +void clitk::DicomRTStruct2ImageFilter::SetOutputSize(const unsigned long* size) { std::copy(size,size+3,std::back_inserter(mSize)); } //-------------------------------------------------------------------- -void clitk::DicomRT_ROI_ConvertToImageFilter::Update() +void clitk::DicomRTStruct2ImageFilter::Update() { if (!mROI) { std::cerr << "Error. No ROI set, please use SetROI." << std::endl; @@ -221,7 +221,7 @@ void clitk::DicomRT_ROI_ConvertToImageFilter::Update() //-------------------------------------------------------------------- -vtkImageData * clitk::DicomRT_ROI_ConvertToImageFilter::GetOutput() +vtkImageData * clitk::DicomRTStruct2ImageFilter::GetOutput() { assert(mBinaryImage); return mBinaryImage; diff --git a/common/clitkDicomRT_ROI_ConvertToImageFilter.h b/common/clitkDicomRTStruct2ImageFilter.h similarity index 89% rename from common/clitkDicomRT_ROI_ConvertToImageFilter.h rename to common/clitkDicomRTStruct2ImageFilter.h index cb06933..1add60e 100644 --- a/common/clitkDicomRT_ROI_ConvertToImageFilter.h +++ b/common/clitkDicomRTStruct2ImageFilter.h @@ -17,8 +17,8 @@ =========================================================================*/ -#ifndef CLITKDICOMRT_ROI_CONVERTTOIMAGEFILTER_H -#define CLITKDICOMRT_ROI_CONVERTTOIMAGEFILTER_H +#ifndef CLITKDICOMRTSTRUCT2IMAGEFILTER_H +#define CLITKDICOMRT_TRUCT2IMAGEFILTER_H #include "clitkDicomRT_ROI.h" #include "clitkImageCommon.h" @@ -29,11 +29,11 @@ namespace clitk { //-------------------------------------------------------------------- - class DicomRT_ROI_ConvertToImageFilter { + class DicomRTStruct2ImageFilter { public: - DicomRT_ROI_ConvertToImageFilter(); - ~DicomRT_ROI_ConvertToImageFilter(); + DicomRTStruct2ImageFilter(); + ~DicomRTStruct2ImageFilter(); void SetROI(clitk::DicomRT_ROI * roi); ///This is used to create a mask with the same characteristics as an input image @@ -66,7 +66,7 @@ namespace clitk { //-------------------------------------------------------------------- template -typename itk::Image::ConstPointer clitk::DicomRT_ROI_ConvertToImageFilter::GetITKOutput() +typename itk::Image::ConstPointer clitk::DicomRTStruct2ImageFilter::GetITKOutput() { assert(mBinaryImage); typedef itk::Image ConnectorImageType; @@ -77,5 +77,5 @@ typename itk::Image::ConstPointer clitk::DicomRT_ROI_Co return connector->GetOutput(); } //-------------------------------------------------------------------- -#endif // CLITKDICOMRT_ROI_CONVERTTOIMAGEFILTER_H +#endif // CLITKDICOMRT_TRUCT2IMAGEFILTER_H diff --git a/common/clitkDicomRT_Contour.cxx b/common/clitkDicomRT_Contour.cxx index 0cc62b0..9874563 100644 --- a/common/clitkDicomRT_Contour.cxx +++ b/common/clitkDicomRT_Contour.cxx @@ -52,11 +52,73 @@ void clitk::DicomRT_Contour::Print(std::ostream & os) const //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +void clitk::DicomRT_Contour::UpdateDicomItem() +{ + DD("DicomRT_Contour::UpdateDicomItem"); + + gdcm::DataSet & nestedds = mItem->GetNestedDataSet(); + + //NON CONSIDER CONTOUR ITEM NOT SEQ ITEM ? + + // Contour type [Contour Geometric Type] + gdcm::Attribute<0x3006,0x0042> contgeotype; + contgeotype.SetFromDataSet( nestedds ); + + // Number of points [Number of Contour Points] + gdcm::Attribute<0x3006,0x0046> numcontpoints; + numcontpoints.SetFromDataSet( nestedds ); + DD(mNbOfPoints); + mNbOfPoints = numcontpoints.GetValue(); + DD(mNbOfPoints); + + // Contour dicom values from DataPoints + //ComputeDataPointsFromMesh(); + uint nb = mData->GetNumberOfPoints(); + std::vector points; + points.resize(mNbOfPoints*3); + for(unsigned int i=0; iGetPoint(i); + points[i*3] = p[0]; + points[i*3+1] = p[1]; + points[i*3+2] = p[2]; + } + + // Get attribute + gdcm::Attribute<0x3006,0x0050> at; + gdcm::Tag tcontourdata(0x3006,0x0050); + gdcm::DataElement contourdata = nestedds.GetDataElement( tcontourdata ); + at.SetFromDataElement( contourdata ); + + // Set attribute + at.SetValues(&points[0], points.size(), false); + DD(at.GetValues()[0]); + + DD("replace"); + nestedds.Replace(at.GetAsDataElement()); + DD("done"); + + // Change Number of points [Number of Contour Points] + numcontpoints.SetValue(nb); + nestedds.Replace(numcontpoints.GetAsDataElement()); + + // Test + gdcm::DataElement aa = nestedds.GetDataElement( tcontourdata ); + at.SetFromDataElement( aa ); + const double* bb = at.GetValues(); + DD(bb[0]); + +} +//-------------------------------------------------------------------- + + //-------------------------------------------------------------------- #if GDCM_MAJOR_VERSION == 2 -bool clitk::DicomRT_Contour::Read(gdcm::Item const & item) +bool clitk::DicomRT_Contour::Read(gdcm::Item * item) { - const gdcm::DataSet& nestedds2 = item.GetNestedDataSet(); + mItem = item; + + const gdcm::DataSet& nestedds2 = item->GetNestedDataSet(); // Contour type [Contour Geometric Type] gdcm::Attribute<0x3006,0x0042> contgeotype; @@ -82,7 +144,7 @@ bool clitk::DicomRT_Contour::Read(gdcm::Item const & item) const gdcm::DataElement & contourdata = nestedds2.GetDataElement( tcontourdata ); at.SetFromDataElement( contourdata ); const double* points = at.GetValues(); - + // unsigned int npts = at.GetNumberOfValues() / 3; assert(at.GetNumberOfValues() == static_cast(mNbOfPoints)*3); // Organize values @@ -162,7 +224,7 @@ bool clitk::DicomRT_Contour::Read(gdcm::SQItem * item) vtkPolyData * clitk::DicomRT_Contour::GetMesh() { if (!mMeshIsUpToDate) { - ComputeMesh(); + ComputeMeshFromDataPoints(); } return mMesh; } @@ -170,7 +232,15 @@ vtkPolyData * clitk::DicomRT_Contour::GetMesh() //-------------------------------------------------------------------- -void clitk::DicomRT_Contour::ComputeMesh() +void clitk::DicomRT_Contour::SetMesh(vtkPolyData * mesh) +{ + mMesh = mesh; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +void clitk::DicomRT_Contour::ComputeMeshFromDataPoints() { // DD("ComputeMesh Contour"); mMesh = vtkSmartPointer::New(); @@ -184,11 +254,35 @@ void clitk::DicomRT_Contour::ComputeMesh() mData->GetPoint(idx)[2]); ids[0]=idx; ids[1]=(ids[0]+1) % mNbOfPoints; //0-1,1-2,...,n-1-0 - // DD(ids[0]); -// DD(ids[1]); mMesh->GetLines()->InsertNextCell(2,ids); } - // DD(mMesh->GetNumberOfCells()); mMeshIsUpToDate = true; } //-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +void clitk::DicomRT_Contour::ComputeDataPointsFromMesh() +{ + DD("ComputeDataPointsFromMesh"); + + + /*todo + + mMesh = vtkSmartPointer::New(); + mMesh->Allocate(); //for cell structures + mPoints = vtkSmartPointer::New(); + mMesh->SetPoints(mPoints); + vtkIdType ids[2]; + for (unsigned int idx=0 ; idxGetPoints()->InsertNextPoint(mData->GetPoint(idx)[0], + mData->GetPoint(idx)[1], + mData->GetPoint(idx)[2]); + ids[0]=idx; + ids[1]=(ids[0]+1) % mNbOfPoints; //0-1,1-2,...,n-1-0 + mMesh->GetLines()->InsertNextCell(2,ids); + } + mMeshIsUpToDate = true; +*/ +} +//-------------------------------------------------------------------- diff --git a/common/clitkDicomRT_Contour.h b/common/clitkDicomRT_Contour.h index 0484da6..91d75d6 100644 --- a/common/clitkDicomRT_Contour.h +++ b/common/clitkDicomRT_Contour.h @@ -42,17 +42,23 @@ public: itkNewMacro(Self); void Print(std::ostream & os = std::cout) const; + #if GDCM_MAJOR_VERSION == 2 - bool Read(gdcm::Item const & item); + bool Read(gdcm::Item * item); + void UpdateDicomItem(); #else bool Read(gdcm::SQItem * item); #endif + vtkPolyData * GetMesh(); + void SetMesh(vtkPolyData * mesh); vtkPoints * GetPoints() {return mData;} double GetZ() const {return mZ;} + protected: - void ComputeMesh(); + void ComputeMeshFromDataPoints(); + void ComputeDataPointsFromMesh(); unsigned int mNbOfPoints; std::string mType; vtkSmartPointer mData; @@ -61,6 +67,8 @@ protected: bool mMeshIsUpToDate; ///Z location of the contour double mZ; + + gdcm::Item * mItem; private: DicomRT_Contour(); diff --git a/common/clitkDicomRT_ROI.cxx b/common/clitkDicomRT_ROI.cxx index 88bab72..75fc8e9 100644 --- a/common/clitkDicomRT_ROI.cxx +++ b/common/clitkDicomRT_ROI.cxx @@ -20,6 +20,8 @@ #include "clitkDicomRT_ROI.h" #include #include +#include +#include #if GDCM_MAJOR_VERSION == 2 #include "gdcmAttribute.h" @@ -31,11 +33,13 @@ clitk::DicomRT_ROI::DicomRT_ROI() { mName = "NoName"; mNumber = -1; + mImage = NULL; mColor.resize(3); mColor[0] = mColor[1] = mColor[2] = 0; mMeshIsUpToDate = false; mBackgroundValue = 0; mForegroundValue = 1; + SetDicomUptodateFlag(false); } //-------------------------------------------------------------------- @@ -134,21 +138,38 @@ double clitk::DicomRT_ROI::GetForegroundValueLabelImage() const //-------------------------------------------------------------------- #if GDCM_MAJOR_VERSION == 2 -void clitk::DicomRT_ROI::Read(std::map & rois, gdcm::Item const & item) +bool clitk::DicomRT_ROI::Read(gdcm::Item * itemInfo, gdcm::Item * itemContour) { - const gdcm::DataSet& nestedds = item.GetNestedDataSet(); - + // Keep dicom item + mItemInfo = itemInfo; + mItemContour = itemContour; + // DD(mItemInfo); + + // ROI number [Referenced ROI Number] + const gdcm::DataSet & nesteddsInfo = mItemInfo->GetNestedDataSet(); + gdcm::Attribute<0x3006,0x0022> roinumber; + roinumber.SetFromDataSet( nesteddsInfo ); + int nb1 = roinumber.GetValue(); + + // Check this is the same with the other item + const gdcm::DataSet & nestedds = mItemContour->GetNestedDataSet(); gdcm::Attribute<0x3006,0x0084> referencedroinumber; referencedroinumber.SetFromDataSet( nestedds ); - // Change number if needed - - // TODO + int nb2 = referencedroinumber.GetValue(); + + // Must never be different + if (nb1 != nb2) { + DD(nb2); + DD(nb1); + FATAL("nb1 must equal nb2" << std::endl); + } + mNumber = nb1; - // ROI number [Referenced ROI Number] - mNumber = referencedroinumber.GetValue(); - - // Retrieve ROI Name - mName = rois[mNumber]; + // Retrieve ROI Name (in the info item) + gdcm::Attribute<0x3006,0x26> roiname; + roiname.SetFromDataSet( nesteddsInfo ); + mName = roiname.GetValue(); + // DD(mName); // ROI Color [ROI Display Color] gdcm::Attribute<0x3006,0x002a> color = {}; @@ -163,10 +184,12 @@ void clitk::DicomRT_ROI::Read(std::map & rois, gdcm::Item cons if( !nestedds.FindDataElement( tcsq ) ) { std::cerr << "Warning. Could not read contour for structure <" << mName << ">, number" << mNumber << " ? I ignore it" << std::endl; - return; + SetDicomUptodateFlag(true); + return false; } const gdcm::DataElement& csq = nestedds.GetDataElement( tcsq ); - gdcm::SmartPointer sqi2 = csq.GetValueAsSQ(); + mContoursSequenceOfItems = csq.GetValueAsSQ(); + gdcm::SmartPointer & sqi2 = mContoursSequenceOfItems; if( !sqi2 || !sqi2->GetNumberOfItems() ) { } @@ -174,22 +197,19 @@ void clitk::DicomRT_ROI::Read(std::map & rois, gdcm::Item cons for(unsigned int i = 0; i < nitems; ++i) { - const gdcm::Item & j = sqi2->GetItem(i+1); // Item start at #1 + gdcm::Item & j = sqi2->GetItem(i+1); // Item start at #1 DicomRT_Contour::Pointer c = DicomRT_Contour::New(); - bool b = c->Read(j); + bool b = c->Read(&j); if (b) { mListOfContours.push_back(c); } } + SetDicomUptodateFlag(true); + return true; } #else void clitk::DicomRT_ROI::Read(std::map & rois, gdcm::SQItem * item) { - - // Change number if needed - - // TODO - // ROI number [Referenced ROI Number] mNumber = atoi(item->GetEntryValue(0x3006,0x0084).c_str()); @@ -215,6 +235,7 @@ void clitk::DicomRT_ROI::Read(std::map & rois, gdcm::SQItem * else { std::cerr << "Warning. Could not read contour for structure <" << mName << ">, number" << mNumber << " ? I ignore it" << std::endl; } + SetDicomUptodateFlag(true); } #endif //-------------------------------------------------------------------- @@ -232,18 +253,23 @@ void clitk::DicomRT_ROI::SetImage(vvImage * image) vtkPolyData * clitk::DicomRT_ROI::GetMesh() { if (!mMeshIsUpToDate) { - ComputeMesh(); + ComputeMeshFromContour(); } return mMesh; } +//-------------------------------------------------------------------- + + //-------------------------------------------------------------------- clitk::DicomRT_Contour * clitk::DicomRT_ROI::GetContour(int n) { return mListOfContours[n]; } +//-------------------------------------------------------------------- + //-------------------------------------------------------------------- -void clitk::DicomRT_ROI::ComputeMesh() +void clitk::DicomRT_ROI::ComputeMeshFromContour() { vtkSmartPointer append = vtkSmartPointer::New(); for(unsigned int i=0; i + gdcm::Attribute<0x3006,0x26> roiname; + roiname.SetValue(GetName()); + gdcm::DataElement de = roiname.GetAsDataElement(); + gdcm::DataSet & ds = mItemInfo->GetNestedDataSet(); + ds.Replace(de); + + // From MESH to CONTOURS + ComputeContoursFromImage(); + + // Update contours + DD(mListOfContours.size()); + for(uint i=0; iUpdateDicomItem();//mItemContour); + } + + // Nb of contours + unsigned int nitems = mContoursSequenceOfItems->GetNumberOfItems(); + DD(nitems); + + // Write [Contour Sequence] = 0x3006,0x0040) + gdcm::DataSet & dsc = mItemContour->GetNestedDataSet(); + gdcm::Tag tcsq(0x3006,0x0040); + const gdcm::DataElement& csq = dsc.GetDataElement( tcsq ); + gdcm::DataElement dec(csq); + dec.SetValue(*mContoursSequenceOfItems); + dsc.Replace(dec); + + gdcm::DataSet & a = mContoursSequenceOfItems->GetItem(1).GetNestedDataSet(); + gdcm::Attribute<0x3006,0x0050> at; + gdcm::Tag tcontourdata(0x3006,0x0050); + gdcm::DataElement contourdata = a.GetDataElement( tcontourdata ); + at.SetFromDataElement( contourdata ); + const double* points = at.GetValues(); + DD(points[0]); + +} +//-------------------------------------------------------------------- +#endif + //-------------------------------------------------------------------- void clitk::DicomRT_ROI::SetFromBinaryImage(vvImage * image, int n, std::string name, @@ -290,3 +369,49 @@ vvImage * clitk::DicomRT_ROI::GetImage() const return mImage; } //-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +void clitk::DicomRT_ROI::ComputeContoursFromImage() +{ + DD("ComputeMeshFromImage"); + + // Check that an image is loaded + if (!mImage) return; + + // Only consider 3D here + if (mImage->GetNumberOfDimensions() != 3) { + FATAL("DicomRT_ROI::ComputeMeshFromImage only work with 3D images"); + } + + // Get the VTK image + vtkImageData * image = mImage->GetVTKImages()[0]; + + // Get initial extend for the clipping + vtkSmartPointer clipper = vtkSmartPointer::New(); + clipper->SetInput(image); + int* extent = image->GetExtent(); + DDV(extent, 6); + // std::vector extend; + + // Prepare the marching squares + vtkSmartPointer squares = vtkSmartPointer::New(); + squares->SetInput(clipper->GetOutput()); + + // Loop on slice + uint n = image->GetDimensions()[2]; + DD(n); + for(uint i=0; iGetOrigin()[2]+i*image->GetSpacing()[2]; + DDV(extent, 6); + clipper->SetOutputWholeExtent(extent[0],extent[1],extent[2], + extent[3],extent[4],extent[5]); + + + squares->Update(); + DD(squares->GetNumberOfContours()); + mListOfContours[i]->SetMesh(squares->GetOutput()); + } +} +//-------------------------------------------------------------------- diff --git a/common/clitkDicomRT_ROI.h b/common/clitkDicomRT_ROI.h index 0fb3076..f0be23d 100644 --- a/common/clitkDicomRT_ROI.h +++ b/common/clitkDicomRT_ROI.h @@ -35,11 +35,6 @@ public: itkNewMacro(Self); void Print(std::ostream & os = std::cout) const; -#if GDCM_MAJOR_VERSION == 2 - void Read(std::map & rois, gdcm::Item const & item); -#else - void Read(std::map & rois, gdcm::SQItem * item); -#endif void SetFromBinaryImage(vvImage * image, int n, std::string name, std::vector color, @@ -64,10 +59,24 @@ public: void SetImage(vvImage * im); DicomRT_Contour* GetContour(int n); - // double GetContourSpacing() const {return mZDelta;} + // Compute a vtk mesh from the dicom contours + void ComputeMeshFromContour(); + void ComputeContoursFromImage(); + // Indicate if the mesh is uptodate according to the dicom + void SetDicomUptodateFlag(bool b) { m_DicomUptodateFlag = b; } + bool GetDicomUptoDateFlag() const { return m_DicomUptodateFlag; } + void SetName(std::string n) { mName = n; } + + // Read from DICOM RT STRUCT +#if GDCM_MAJOR_VERSION == 2 + bool Read(gdcm::Item * itemInfo, gdcm::Item * itemContour); + void UpdateDicomItem(); +#else + void Read(std::map & rois, gdcm::SQItem * item); +#endif + protected: - void ComputeMesh(); std::string mName; std::string mFilename; int mNumber; @@ -78,8 +87,13 @@ protected: vvImage::Pointer mImage; double mBackgroundValue; double mForegroundValue; - ///Spacing between two contours - // double mZDelta; + bool m_DicomUptodateFlag; + +#if GDCM_MAJOR_VERSION == 2 + gdcm::Item * mItemInfo; + gdcm::Item * mItemContour; + gdcm::SmartPointer mContoursSequenceOfItems; +#endif private: DicomRT_ROI(); diff --git a/common/clitkDicomRT_StructureSet.cxx b/common/clitkDicomRT_StructureSet.cxx index 8543fb7..4ac4705 100644 --- a/common/clitkDicomRT_StructureSet.cxx +++ b/common/clitkDicomRT_StructureSet.cxx @@ -20,10 +20,6 @@ #include "clitkDicomRT_StructureSet.h" #include #include "gdcmFile.h" -#if GDCM_MAJOR_VERSION == 2 -#include "gdcmReader.h" -#include "gdcmAttribute.h" -#endif //-------------------------------------------------------------------- clitk::DicomRT_StructureSet::DicomRT_StructureSet() @@ -35,6 +31,7 @@ clitk::DicomRT_StructureSet::DicomRT_StructureSet() mName = "NoName"; mDate = "NoDate"; mTime = "NoTime"; + mFile = NULL; } //-------------------------------------------------------------------- @@ -103,23 +100,21 @@ const std::string & clitk::DicomRT_StructureSet::GetTime() const //-------------------------------------------------------------------- -const std::vector & clitk::DicomRT_StructureSet::GetListOfROI() const -{ - return mListOfROI; -} +// const std::vector & clitk::DicomRT_StructureSet::GetListOfROI() const +// { +// return mListOfROI; +// } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -clitk::DicomRT_ROI* clitk::DicomRT_StructureSet::GetROI(int n) +clitk::DicomRT_ROI* clitk::DicomRT_StructureSet::GetROIFromROINumber(int n) { - if (mMapOfROIIndex.find(n) == mMapOfROIIndex.end()) { + if (mROIs.find(n) == mROIs.end()) { std::cerr << "No ROI number " << n << std::endl; return NULL; } - // DD(mListOfROI[mMapOfROIIndex[n]]->GetName()); - //DD(mListOfROI[mMapOfROIIndex[n]]->GetROINumber()); - return mListOfROI[mMapOfROIIndex[n]]; + return mROIs[n]; } //-------------------------------------------------------------------- @@ -133,30 +128,106 @@ void clitk::DicomRT_StructureSet::Print(std::ostream & os) const << "Struct Label = " << mLabel << std::endl << "Struct Name = " << mName << std::endl << "Struct Time = " << mTime << std::endl - << "Number of ROI = " << mListOfROI.size() << std::endl; - for(unsigned int i=0; iPrint(os); + << "Number of ROI = " << mROIs.size() << std::endl; + for(ROIConstIteratorType iter = mROIs.begin(); iter != mROIs.end(); iter++) { + iter->second->Print(os); } } //-------------------------------------------------------------------- +#if GDCM_MAJOR_VERSION == 2 //-------------------------------------------------------------------- -void clitk::DicomRT_StructureSet::Read(const std::string & filename) +int clitk::DicomRT_StructureSet::ReadROINumber(const gdcm::Item & item) +{ + // 0x3006,0x0022 = [ROI Number] + const gdcm::DataSet & nestedds = item.GetNestedDataSet(); + gdcm::Attribute<0x3006,0x0022> roinumber; + roinumber.SetFromDataSet( nestedds ); + return roinumber.GetValue(); +} +//-------------------------------------------------------------------- +#endif + +//-------------------------------------------------------------------- +void clitk::DicomRT_StructureSet::Write(const std::string & filename) { - // Open DICOM #if GDCM_MAJOR_VERSION == 2 - gdcm::Reader reader; - reader.SetFileName(filename.c_str()); - reader.Read(); + DD("DCM RT Writer"); + + // Assert that the gdcm file is still open (we can write only if it was readed) + if (mFile == NULL) { + //assert(mFile != NULL); + FATAL("Sorry, I can write DICOM only if it was read first from a file with 'Read' function"); + } - const gdcm::File & file = reader.GetFile(); - const gdcm::DataSet & ds = file.GetDataSet(); + // Loop and update each ROI + for(ROIIteratorType iter = mROIs.begin(); iter != mROIs.end(); iter++) { + iter->second->UpdateDicomItem(); + } + // Write [ Structure Set ROI Sequence ] = 0x3006,0x0020 + gdcm::DataSet & ds = mFile->GetDataSet(); + gdcm::Tag tssroisq(0x3006,0x0020); + const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq ); + gdcm::DataElement de(ssroisq); + de.SetValue(*mROIInfoSequenceOfItems); + ds.Replace(de); + + // Write [ ROI Contour Sequence ] = 0x3006,0x0039 + DD("ici"); + gdcm::Tag troicsq(0x3006,0x0039); + const gdcm::DataElement &roicsq = ds.GetDataElement( troicsq ); + gdcm::DataElement de2(roicsq); + de2.SetValue(*mROIContoursSequenceOfItems); + ds.Replace(de2); + + //DEBUG + gdcm::DataSet & a = mROIContoursSequenceOfItems->GetItem(1).GetNestedDataSet(); + gdcm::Tag tcsq(0x3006,0x0040); + const gdcm::DataElement& csq = a.GetDataElement( tcsq ); + gdcm::SmartPointer sqi2 = csq.GetValueAsSQ(); + gdcm::Item & j = sqi2->GetItem(1); + gdcm::DataSet & b = j.GetNestedDataSet(); + gdcm::Attribute<0x3006,0x0050> at; + gdcm::Tag tcontourdata(0x3006,0x0050); + gdcm::DataElement contourdata = b.GetDataElement( tcontourdata ); + at.SetFromDataElement( contourdata ); + const double* points = at.GetValues(); + DD(points[0]); + + + // Write dicom + gdcm::Writer writer; + //writer.CheckFileMetaInformationOff(); + writer.SetFileName(filename.c_str()); + writer.SetFile(*mFile); + DD("before write"); + writer.Write(); + DD("End write"); +#else + FATAL("Sorry not compatible with GDCM1, use GDCM2"); +#endif +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +void clitk::DicomRT_StructureSet::Read(const std::string & filename) +{ + // Open DICOM +#if GDCM_MAJOR_VERSION == 2 + // Read gdcm file + mReader = new gdcm::Reader; + mReader->SetFileName(filename.c_str()); + mReader->Read(); + mFile = &(mReader->GetFile()); + const gdcm::DataSet & ds = mFile->GetDataSet(); + // Check file type //Verify if the file is a RT-Structure-Set dicom file gdcm::MediaStorage ms; - ms.SetFromFile(file); + ms.SetFromFile(*mFile); if( ms != gdcm::MediaStorage::RTStructureSetStorage ) { std::cerr << "Error. the file " << filename @@ -196,12 +267,17 @@ void clitk::DicomRT_StructureSet::Read(const std::string & filename) mName = atname.GetValue(); mTime = time.GetValue(); + // Temporary store the list of items + std::map mMapOfROIInfo; + std::map mMapOfROIContours; + //---------------------------------- // Read all ROI Names and number // 0x3006,0x0020 = [ Structure Set ROI Sequence ] gdcm::Tag tssroisq(0x3006,0x0020); const gdcm::DataElement &ssroisq = ds.GetDataElement( tssroisq ); - gdcm::SmartPointer roi_seq = ssroisq.GetValueAsSQ(); + mROIInfoSequenceOfItems = ssroisq.GetValueAsSQ(); + gdcm::SmartPointer & roi_seq = mROIInfoSequenceOfItems; assert(roi_seq); // TODO error message for(unsigned int ridx = 0; ridx < roi_seq->GetNumberOfItems(); ++ridx) { @@ -210,13 +286,13 @@ void clitk::DicomRT_StructureSet::Read(const std::string & filename) gdcm::Attribute<0x3006,0x26> roiname; roiname.SetFromDataSet( nestedds ); - std::string name = roiname.GetValue(); // 0x3006,0x0026 = [ROI Name] - gdcm::Attribute<0x3006,0x0022> roinumber; - roinumber.SetFromDataSet( nestedds ); - int nb = roinumber.GetValue(); // 0x3006,0x0022 = [ROI Number] - // Change number if needed + std::string name = roiname.GetValue(); // 0x3006,0x0026 = [ROI Name] - //TODO + // 0x3006,0x0022 = [ROI Number] + int nb = ReadROINumber(item); + + // Store the item + mMapOfROIInfo[nb] = &item; // Check if such a number already exist if (mMapOfROIName.find(nb) != mMapOfROIName.end()) { @@ -226,43 +302,57 @@ void clitk::DicomRT_StructureSet::Read(const std::string & filename) // Add in map mMapOfROIName[nb] = name; } - // DD(mMapOfROIName.size()); //---------------------------------- - // Read all ROI + // Read all ROI item // 0x3006,0x0039 = [ ROI Contour Sequence ] gdcm::Tag troicsq(0x3006,0x0039); const gdcm::DataElement &roicsq = ds.GetDataElement( troicsq ); gdcm::SmartPointer roi_contour_seq = roicsq.GetValueAsSQ(); + mROIContoursSequenceOfItems = roi_contour_seq; assert(roi_contour_seq); // TODO error message - int n=0; - for(unsigned int ridx = 0; ridx < roi_contour_seq->GetNumberOfItems(); ++ridx) - { + for(unsigned int ridx = 0; ridx < roi_contour_seq->GetNumberOfItems(); ++ridx) { gdcm::Item & item = roi_contour_seq->GetItem( ridx + 1); // Item starts at 1 + // ROI number [Referenced ROI Number] + const gdcm::DataSet& nestedds = item.GetNestedDataSet(); + gdcm::Attribute<0x3006,0x0084> referencedroinumber; + referencedroinumber.SetFromDataSet( nestedds ); + int nb = referencedroinumber.GetValue(); + // Store the item + mMapOfROIContours[nb] = &item; + } + //---------------------------------- + // Create the ROIs + for(std::map::iterator i = mMapOfROIInfo.begin(); i != mMapOfROIInfo.end(); i++) { + int nb = i->first;//ReadROINumber(i);//mROIIndex[i]; + // Create the roi DicomRT_ROI::Pointer roi = DicomRT_ROI::New(); - roi->Read(mMapOfROIName, item); - mListOfROI.push_back(roi); - mMapOfROIIndex[roi->GetROINumber()] = n; - n++; - } + roi->Read(mMapOfROIInfo[nb], mMapOfROIContours[nb]); + // mListOfROI.push_back(roi); + // mMapOfROIIndex[nb] = i; + mROIs[nb] = roi; + } + //---------------------------------------------------------------------------------------- + //---------------------------------------------------------------------------------------- + //---------------------------------------------------------------------------------------- #else - gdcm::File reader; - reader.SetFileName(filename.c_str()); - reader.SetMaxSizeLoadEntry(16384); // Needed ... - reader.SetLoadMode(gdcm::LD_NOSHADOW); // don't load shadow tags (in order to save memory) - reader.Load(); - + mFile = new gdcm::File; + mFile->SetFileName(filename.c_str()); + mFile->SetMaxSizeLoadEntry(16384); // Needed ... + mFile->SetLoadMode(gdcm::LD_NOSHADOW); // don't load shadow tags (in order to save memory) + mFile->Load(); + // Check file type //Verify if the file is a RT-Structure-Set dicom file - if (!gdcm::Util::DicomStringEqual(reader.GetEntryValue(0x0008,0x0016),"1.2.840.10008.5.1.4.1.1.481.3")) { //SOP clas UID + if (!gdcm::Util::DicomStringEqual(mFile->GetEntryValue(0x0008,0x0016),"1.2.840.10008.5.1.4.1.1.481.3")) { //SOP clas UID std::cerr << "Error. the file " << filename << " is not a Dicom Struct ? (must have a SOP Class UID [0008|0016] = 1.2.840.10008.5.1.4.1.1.481.3 ==> [RT Structure Set Storage])" << std::endl; exit(0); } - if (!gdcm::Util::DicomStringEqual(reader.GetEntryValue(0x0008,0x0060),"RTSTRUCT")) { //SOP clas UID + if (!gdcm::Util::DicomStringEqual(mFile->GetEntryValue(0x0008,0x0060),"RTSTRUCT")) { //SOP clas UID std::cerr << "Error. the file " << filename << " is not a Dicom Struct ? (must have 0x0008,0x0060 = RTSTRUCT [RT Structure Set Storage])" << std::endl; @@ -270,30 +360,26 @@ void clitk::DicomRT_StructureSet::Read(const std::string & filename) } // Read global info - mStudyID = reader.GetValEntry(0x0020,0x0010)->GetValue(); - mStudyTime = reader.GetValEntry(0x008,0x0020)->GetValue(); - mStudyDate = reader.GetValEntry(0x008,0x0030)->GetValue(); - mLabel = reader.GetValEntry(0x3006,0x002)->GetValue(); - if (!reader.GetValEntry(0x3006,0x004)) { + mStudyID = mFile->GetValEntry(0x0020,0x0010)->GetValue(); + mStudyTime = mFile->GetValEntry(0x008,0x0020)->GetValue(); + mStudyDate = mFile->GetValEntry(0x008,0x0030)->GetValue(); + mLabel = mFile->GetValEntry(0x3006,0x002)->GetValue(); + if (!mFile->GetValEntry(0x3006,0x004)) { mName = "Anonymous"; } else { - mName = reader.GetValEntry(0x3006,0x004)->GetValue(); + mName = mFile->GetValEntry(0x3006,0x004)->GetValue(); } - mTime = reader.GetValEntry(0x3006,0x009)->GetValue(); + mTime = mFile->GetValEntry(0x3006,0x009)->GetValue(); //---------------------------------- // Read all ROI Names and number // 0x3006,0x0020 = [ Structure Set ROI Sequence ] - gdcm::SeqEntry * roi_seq=reader.GetSeqEntry(0x3006,0x0020); + gdcm::SeqEntry * roi_seq=mFile->GetSeqEntry(0x3006,0x0020); assert(roi_seq); // TODO error message for (gdcm::SQItem* r=roi_seq->GetFirstSQItem(); r!=0; r=roi_seq->GetNextSQItem()) { std::string name = r->GetEntryValue(0x3006,0x0026); // 0x3006,0x0026 = [ROI Name] int nb = atoi(r->GetEntryValue(0x3006,0x0022).c_str()); // 0x3006,0x0022 = [ROI Number] - // Change number if needed - - //TODO - // Check if such a number already exist if (mMapOfROIName.find(nb) != mMapOfROIName.end()) { std::cerr << "WARNING. A Roi already exist with the number " @@ -302,19 +388,17 @@ void clitk::DicomRT_StructureSet::Read(const std::string & filename) // Add in map mMapOfROIName[nb] = name; } - // DD(mMapOfROIName.size()); //---------------------------------- // Read all ROI // 0x3006,0x0039 = [ ROI Contour Sequence ] - gdcm::SeqEntry * roi_contour_seq=reader.GetSeqEntry(0x3006,0x0039); + gdcm::SeqEntry * roi_contour_seq=mFile->GetSeqEntry(0x3006,0x0039); assert(roi_contour_seq); // TODO error message int n=0; for (gdcm::SQItem* r=roi_contour_seq->GetFirstSQItem(); r!=0; r=roi_contour_seq->GetNextSQItem()) { DicomRT_ROI::Pointer roi = DicomRT_ROI::New(); roi->Read(mMapOfROIName, r); - mListOfROI.push_back(roi); - mMapOfROIIndex[roi->GetROINumber()] = n; + mROIs[roi->GetROINumber()] = roi; n++; } @@ -326,22 +410,19 @@ void clitk::DicomRT_StructureSet::Read(const std::string & filename) //-------------------------------------------------------------------- int clitk::DicomRT_StructureSet::AddBinaryImageAsNewROI(vvImage * im, std::string n) { - //DD("AddBinaryImageAsNewROI"); // Search max ROI number int max = -1; - for(unsigned int i=0; iGetROINumber() > max) - max = mListOfROI[i]->GetROINumber(); + for(ROIConstIteratorType iter = mROIs.begin(); iter != mROIs.end(); iter++) { + // for(unsigned int i=0; isecond; + if (roi->GetROINumber() > max) + max = roi->GetROINumber(); } - // DD(max); ++max; - //DD(max); // Compute name std::ostringstream oss; oss << vtksys::SystemTools::GetFilenameName(vtksys::SystemTools::GetFilenameWithoutLastExtension(n)); - // << "_roi_" << max << vtksys::SystemTools::GetFilenameLastExtension(n); - //DD(oss.str()); mMapOfROIName[max] = oss.str(); // Set color @@ -353,9 +434,7 @@ int clitk::DicomRT_StructureSet::AddBinaryImageAsNewROI(vvImage * im, std::strin // Create ROI DicomRT_ROI::Pointer roi = DicomRT_ROI::New(); roi->SetFromBinaryImage(im, max, oss.str(), color, n); - mListOfROI.push_back(roi); - mMapOfROIIndex[mListOfROI.size()-1] = max; - //DD(mMapOfROIIndex[mListOfROI.size()-1]); + mROIs[max] = roi; return max; } //-------------------------------------------------------------------- diff --git a/common/clitkDicomRT_StructureSet.h b/common/clitkDicomRT_StructureSet.h index 3dbb8a6..5e4cf97 100644 --- a/common/clitkDicomRT_StructureSet.h +++ b/common/clitkDicomRT_StructureSet.h @@ -20,10 +20,20 @@ #ifndef CLITKDICOMRT_STRUCTURESET_H #define CLITKDICOMRT_STRUCTURESET_H +// clitk #include "clitkCommon.h" #include "clitkDicomRT_ROI.h" + +// vv #include "vvImage.h" +// gdcm +#if GDCM_MAJOR_VERSION == 2 +#include "gdcmReader.h" +#include "gdcmWriter.h" +#include "gdcmAttribute.h" +#endif + namespace clitk { //-------------------------------------------------------------------- @@ -34,11 +44,15 @@ public: typedef itk::SmartPointer Pointer; itkNewMacro(Self); + typedef typename std::map::iterator ROIIteratorType; + typedef typename std::map::const_iterator ROIConstIteratorType; + void Print(std::ostream & os = std::cout) const; void Read(const std::string & filename); + void Write(const std::string & filename); - const std::vector & GetListOfROI() const; - clitk::DicomRT_ROI * GetROI(int n); + clitk::DicomRT_ROI * GetROIFromROINumber(int n); + std::map & GetROIs() { return mROIs; } const std::string & GetStudyID() const; const std::string & GetStudyTime() const; const std::string & GetStudyDate() const; @@ -49,6 +63,11 @@ public: int AddBinaryImageAsNewROI(vvImage * i, std::string name); +#if GDCM_MAJOR_VERSION == 2 + // Static + static int ReadROINumber(const gdcm::Item & item); +#endif + protected: std::string mStudyID; std::string mStudyTime; @@ -57,9 +76,15 @@ protected: std::string mName; std::string mDate; std::string mTime; + + std::map mROIs; std::map mMapOfROIName; - std::map mMapOfROIIndex; - std::vector mListOfROI; +#if GDCM_MAJOR_VERSION == 2 + gdcm::Reader * mReader; + gdcm::SmartPointer mROIInfoSequenceOfItems; + gdcm::SmartPointer mROIContoursSequenceOfItems; +#endif + gdcm::File * mFile; private: DicomRT_StructureSet(); diff --git a/common/clitkFilterBase.cxx b/common/clitkFilterBase.cxx index 482a82c..32dd9b4 100644 --- a/common/clitkFilterBase.cxx +++ b/common/clitkFilterBase.cxx @@ -29,6 +29,7 @@ clitk::FilterBase::FilterBase() VerboseWarningFlagOn(); VerboseWarningFlagOff(); VerboseMemoryFlagOff(); + VerboseImageSizeFlagOff(); SetWarning(""); VerboseWarningFlagOn(); m_IsCancelled = false; diff --git a/common/clitkFilterBase.h b/common/clitkFilterBase.h index 25cfc4c..29caaf0 100644 --- a/common/clitkFilterBase.h +++ b/common/clitkFilterBase.h @@ -69,6 +69,11 @@ namespace clitk { itkGetConstMacro(VerboseMemoryFlag, bool); itkBooleanMacro(VerboseMemoryFlag); + // Verbose ImageSize + itkSetMacro(VerboseImageSizeFlag, bool); + itkGetConstMacro(VerboseImageSizeFlag, bool); + itkBooleanMacro(VerboseImageSizeFlag); + // Steps management itkSetMacro(NumberOfSteps, int); itkGetConstMacro(NumberOfSteps, int); @@ -115,13 +120,14 @@ namespace clitk { virtual ~FilterBase() {} void StartNewStep(std::string s); template - void StopCurrentStep(typename TInternalImageType::Pointer p); + void StopCurrentStep(typename TInternalImageType::Pointer p, std::string txt=""); void StopCurrentStep(); bool m_VerboseFlag; bool m_VerboseOptionFlag; bool m_VerboseStepFlag; bool m_VerboseMemoryFlag; + bool m_VerboseImageSizeFlag; bool m_WriteStepFlag; int m_CurrentStepNumber; int m_NumberOfSteps; diff --git a/common/clitkFilterBase.txx b/common/clitkFilterBase.txx index 81568ed7..d6c028e 100644 --- a/common/clitkFilterBase.txx +++ b/common/clitkFilterBase.txx @@ -59,7 +59,7 @@ void clitk::FilterBase::VerboseOptionV(std::string name, int nb, OptionType * va //-------------------------------------------------------------------- template -void clitk::FilterBase::StopCurrentStep(typename TInternalImageType::Pointer p) +void clitk::FilterBase::StopCurrentStep(typename TInternalImageType::Pointer p, std::string txt) { StopCurrentStep(); if (m_WriteStepFlag) { @@ -68,6 +68,19 @@ void clitk::FilterBase::StopCurrentStep(typename TInternalImageType::Pointer p) clitk::writeImage(p, name.str()); } clitk::PrintMemory(GetVerboseMemoryFlag(), "End of step"); + if (GetVerboseImageSizeFlag()) { + std::ostream & os = std::cout; + int dim = p->GetImageDimension(); + int nb = 1; + os << txt << " size = "; + for(unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i] << "x"; + nb *= p->GetLargestPossibleRegion().GetSize()[i]; + } + os << p->GetLargestPossibleRegion().GetSize()[dim-1] << " "; + nb *= p->GetLargestPossibleRegion().GetSize()[dim-1]; + os << " pixels = " << nb << std::endl; + } } //-------------------------------------------------------------------- diff --git a/common/clitkImage2DicomRTStructFilter.h b/common/clitkImage2DicomRTStructFilter.h new file mode 100644 index 0000000..749a0be --- /dev/null +++ b/common/clitkImage2DicomRTStructFilter.h @@ -0,0 +1,64 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + Main authors : XX XX XX + + Authors belongs to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + - BSD http://www.opensource.org/licenses/bsd-license.php + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + + =========================================================================*/ + +#ifndef CLITKIMAGE2DICOMRTSTRUCTFILTER_H +#define CLITKIMAGE2DICOMRTSTRUCTFILTER_H + +// clitk +#include "clitkDicomRT_ROI.h" +#include "clitkImageCommon.h" +#include "clitkFilterBase.h" +#include "clitkDicomRT_StructureSet.h" + +namespace clitk { + + //-------------------------------------------------------------------- + template + class Image2DicomRTStructFilter: public clitk::FilterBase { + + public: + Image2DicomRTStructFilter(); + ~Image2DicomRTStructFilter(); + + typedef itk::Image ImageType; + typedef typename ImageType::Pointer ImagePointer; + typedef typename clitk::DicomRT_StructureSet::Pointer DicomRTStructPointer; + + // Set inputs + itkSetMacro(Input, ImagePointer); + itkGetConstMacro(Input, ImagePointer); + + // Run filter + void Update(); + + // Get output + itkGetConstMacro(DicomRTStruct, DicomRTStructPointer); + + protected: + ImagePointer m_Input; + DicomRTStructPointer m_DicomRTStruct; + }; + //-------------------------------------------------------------------- + +} // end namespace clitk + +#include "clitkImage2DicomRTStructFilter.txx" + +#endif // CLITKIMAGE2DICOMRTSTRUCTFILTER_H + diff --git a/common/clitkImage2DicomRTStructFilter.txx b/common/clitkImage2DicomRTStructFilter.txx new file mode 100644 index 0000000..6ab72e2 --- /dev/null +++ b/common/clitkImage2DicomRTStructFilter.txx @@ -0,0 +1,106 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + Main authors : XX XX XX + + Authors belongs to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + - BSD http://www.opensource.org/licenses/bsd-license.php + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + + =========================================================================*/ + +// std +#include +#include + +// clitk +#include "clitkImage2DicomRTStructFilter.h" +#include "clitkImageCommon.h" +#include "vvFromITK.h" + +// vtk +#include +#include +#include +#include +#include +#include + +// itk +#include +#include + +//-------------------------------------------------------------------- +template +clitk::Image2DicomRTStructFilter::Image2DicomRTStructFilter() +{ + DD("Image2DicomRTStructFilter Const"); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +clitk::Image2DicomRTStructFilter::~Image2DicomRTStructFilter() +{ + DD("Image2DicomRTStructFilter Destructor"); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void clitk::Image2DicomRTStructFilter::Update() +{ + DD("Image2DicomRTStructFilter::GenerateData"); + + // Read DicomRTStruct + std::string filename = "RS.zzQAnotmt_french01_.dcm"; + clitk::DicomRT_StructureSet::Pointer structset = clitk::DicomRT_StructureSet::New(); + structset->Read(filename); + + DD(structset->GetName()); + clitk::DicomRT_ROI * roi = structset->GetROIFromROINumber(1); // Aorta + DD(roi->GetName()); + DD(roi->GetROINumber()); + + // Add an image to the roi + vvImage::Pointer im = vvImageFromITK<3, PixelType>(m_Input); + roi->SetImage(im); + + // Get one contour + DD("Compute Mesh"); + roi->ComputeMeshFromImage(); + vtkSmartPointer mesh = roi->GetMesh(); + DD("done"); + + // Change the mesh (shift by 10); + // const vtkSmartPointer & points = mesh->GetPoints(); + // for(uint i=0; iGetNumberOfVerts (); i++) { + // DD(i); + // double * p = points->GetPoint(i); + // p[0] += 30; + // points->SetPoint(i, p); + // } + roi->SetName("TOTO"); + roi->SetDicomUptodateFlag(false); // indicate that dicom info must be updated from the mesh. + + // Convert to dicom ? + DD("TODO"); + + // Write + structset->Write("toto.dcm"); +} +//-------------------------------------------------------------------- + + + + diff --git a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h index 7eb531e..fa5e95b 100644 --- a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h +++ b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.h @@ -78,7 +78,7 @@ namespace clitk { typedef itk::Image FloatImageType; /** Orientation types */ - typedef enum { AtRightTo = 0, AtLeftTo = 1, + typedef enum { RightTo = 0, LeftTo = 1, AntTo = 2, PostTo = 3, InfTo = 4, SupTo = 5, Angle = 6 } OrientationTypeEnumeration; @@ -128,6 +128,12 @@ namespace clitk { itkSetMacro(CombineWithOrFlag, bool); itkBooleanMacro(CombineWithOrFlag); + // I dont want to verify inputs information + virtual void VerifyInputInformation() { } + + // For debug + void PrintOptions(); + protected: AddRelativePositionConstraintToLabelImageFilter(); virtual ~AddRelativePositionConstraintToLabelImageFilter() {} diff --git a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx index 54f4311..694910c 100644 --- a/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx +++ b/itk/clitkAddRelativePositionConstraintToLabelImageFilter.txx @@ -121,22 +121,22 @@ clitk::AddRelativePositionConstraintToLabelImageFilter:: AddOrientationTypeString(std::string t) { m_OrientationTypeString.push_back(t); - switch (t[0]) { - case 'L' : AddOrientationType(AtLeftTo); break; - case 'R' : AddOrientationType(AtRightTo);break; - case 'A' : AddOrientationType(AntTo);break; - case 'P' : AddOrientationType(PostTo);break; - case 'S' : AddOrientationType(SupTo);break; - case 'I' : AddOrientationType(InfTo);break; - case 'N': - if (t == "NotLeftTo") { AddOrientationType(AtLeftTo); InverseOrientationFlagOn(); break; } - if (t == "NotRightTo") { AddOrientationType(AtRightTo); InverseOrientationFlagOn(); break; } - if (t == "NotAntTo") { AddOrientationType(AntTo); InverseOrientationFlagOn(); break; } - if (t == "NotPostTo") { AddOrientationType(PostTo); InverseOrientationFlagOn(); break; } - if (t == "NotSupTo") { AddOrientationType(SupTo); InverseOrientationFlagOn(); break; } - if (t == "NotInfTo") { AddOrientationType(InfTo); InverseOrientationFlagOn(); break; } - default: clitkExceptionMacro("Error, you must provide L,R or A,P or S,I (or NotLeftTo, NotRightTo etc)"); - } + + if (t == "LeftTo") { AddOrientationType(LeftTo); return; } + if (t == "RightTo") { AddOrientationType(RightTo); return; } + if (t == "AntTo") { AddOrientationType(AntTo); return; } + if (t == "PostTo") { AddOrientationType(PostTo); return; } + if (t == "SupTo") { AddOrientationType(SupTo); return; } + if (t == "InfTo") { AddOrientationType(InfTo); return; } + + if (t == "NotLeftTo") { AddOrientationType(LeftTo); InverseOrientationFlagOn(); return; } + if (t == "NotRightTo") { AddOrientationType(RightTo); InverseOrientationFlagOn(); return; } + if (t == "NotAntTo") { AddOrientationType(AntTo); InverseOrientationFlagOn(); return; } + if (t == "NotPostTo") { AddOrientationType(PostTo); InverseOrientationFlagOn(); return; } + if (t == "NotSupTo") { AddOrientationType(SupTo); InverseOrientationFlagOn(); return; } + if (t == "NotInfTo") { AddOrientationType(InfTo); InverseOrientationFlagOn(); return; } + + clitkExceptionMacro("Error, you must provide LeftTo,RightTo or AntTo,PostTo or SupTo,InfTo (or NotLeftTo, NotRightTo etc) but you give " << t); } //-------------------------------------------------------------------- @@ -192,11 +192,11 @@ AddOrientationType(OrientationTypeEnumeration orientation) { m_OrientationType.push_back(orientation); switch (orientation) { - case AtRightTo: + case RightTo: m_Angle1.push_back(clitk::deg2rad(0)); m_Angle2.push_back(clitk::deg2rad(0)); break; - case AtLeftTo: + case LeftTo: m_Angle1.push_back(clitk::deg2rad(180)); m_Angle2.push_back(clitk::deg2rad(0)); break; @@ -230,6 +230,26 @@ AddOrientationType(OrientationTypeEnumeration orientation) //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::AddRelativePositionConstraintToLabelImageFilter:: +PrintOptions() +{ + DD((int)this->GetBackgroundValue()); + DD((int)this->GetObjectBackgroundValue()); + DDV(this->GetOrientationTypeString(), (uint)this->GetNumberOfAngles()); + DD(this->GetIntermediateSpacingFlag()); + DD(this->GetIntermediateSpacing()); + DD(this->GetFuzzyThreshold()); + DD(this->GetAutoCropFlag()); + DD(this->GetInverseOrientationFlag()); + DD(this->GetRemoveObjectFlag()); + DD(this->GetCombineWithOrFlag()); +} +//-------------------------------------------------------------------- + + //-------------------------------------------------------------------- template void @@ -240,6 +260,18 @@ GenerateData() clitkExceptionMacro("Add at least one orientation type"); } + if (GetVerboseOptionFlag()) { + for(int i=0; i(itk::ProcessObject::GetInput(0)); object = dynamic_cast(itk::ProcessObject::GetInput(1)); diff --git a/itk/clitkBooleanOperatorLabelImageFilter.h b/itk/clitkBooleanOperatorLabelImageFilter.h index 71eccda..adb7116 100644 --- a/itk/clitkBooleanOperatorLabelImageFilter.h +++ b/itk/clitkBooleanOperatorLabelImageFilter.h @@ -97,6 +97,9 @@ namespace clitk { itkStaticConstMacro(InputImage1Dimension, unsigned int, TInputImage1::ImageDimension); itkStaticConstMacro(InputImage2Dimension, unsigned int, TInputImage2::ImageDimension); itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension); + + // I dont want to verify inputs information + virtual void VerifyInputInformation() { } protected: BooleanOperatorLabelImageFilter(); diff --git a/itk/clitkCropLikeImageFilter.h b/itk/clitkCropLikeImageFilter.h index d3290ac..9281fca 100644 --- a/itk/clitkCropLikeImageFilter.h +++ b/itk/clitkCropLikeImageFilter.h @@ -65,6 +65,9 @@ namespace clitk { /** ImageDimension constants */ itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension); + // I dont want to verify inputs information + virtual void VerifyInputInformation() { } + protected: CropLikeImageFilter(); virtual ~CropLikeImageFilter() {} diff --git a/itk/clitkCropLikeImageFilter.txx b/itk/clitkCropLikeImageFilter.txx index 2877239..73dd7a7 100644 --- a/itk/clitkCropLikeImageFilter.txx +++ b/itk/clitkCropLikeImageFilter.txx @@ -21,9 +21,10 @@ // clitk #include "clitkCommon.h" +#include "clitkPasteImageFilter.h" // itk -#include "itkPasteImageFilter.h" +//#include "itkPasteImageFilter.h" //-------------------------------------------------------------------- template @@ -228,7 +229,7 @@ GenerateData() { output->FillBuffer(GetBackgroundValue()); // Paste image inside - typedef itk::PasteImageFilter PasteFilterType; + typedef clitk::PasteImageFilter PasteFilterType; typename PasteFilterType::Pointer pasteFilter = PasteFilterType::New(); //pasteFilter->ReleaseDataFlagOn(); // change nothing ? // pasteFilter->InPlaceOn(); // makt it seg fault diff --git a/itk/clitkPasteImageFilter.h b/itk/clitkPasteImageFilter.h new file mode 100644 index 0000000..72550cf --- /dev/null +++ b/itk/clitkPasteImageFilter.h @@ -0,0 +1,49 @@ +/*========================================================================= + * + * COPY OF itkPasteImageFilter to remove VerifyInputInformation + * + *=========================================================================*/ + +#ifndef __clitkPasteImageFilter_h +#define __clitkPasteImageFilter_h + +#include "itkPasteImageFilter.h" + +namespace clitk +{ + using namespace itk; + + template< class TInputImage, class TSourceImage = TInputImage, class TOutputImage = TInputImage > + class ITK_EXPORT PasteImageFilter: + public itk::PasteImageFilter< TInputImage, TSourceImage, TOutputImage > + { + public: + virtual void VerifyInputInformation() { } + + /** Standard class typedefs. */ + typedef PasteImageFilter Self; + typedef InPlaceImageFilter< TInputImage, TOutputImage > Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(PasteImageFilter, InPlaceImageFilter); + + protected: + PasteImageFilter(); + ~PasteImageFilter() {} + + private: + PasteImageFilter(const Self &); //purposely not implemented + void operator=(const Self &); //purposely not implemented + }; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkPasteImageFilter.hxx" +#endif + +#endif diff --git a/itk/clitkPasteImageFilter.hxx b/itk/clitkPasteImageFilter.hxx new file mode 100644 index 0000000..348d858 --- /dev/null +++ b/itk/clitkPasteImageFilter.hxx @@ -0,0 +1,26 @@ +/*========================================================================= + * + * COPY OF itkPasteImageFilter to remove VerifyInputInformation + * + *=========================================================================*/ + +#ifndef __clitkPasteImageFilter_hxx +#define __clitkPasteImageFilter_hxx + +#include "itkPasteImageFilter.h" + +namespace clitk +{ + template< class TInputImage, class TSourceImage, class TOutputImage > + PasteImageFilter< TInputImage, TSourceImage, TOutputImage > + ::PasteImageFilter() + { + this->ProcessObject::SetNumberOfRequiredInputs(2); + + this->InPlaceOff(); + this->m_DestinationIndex.Fill(0); + } + +} // end namespace clitk + +#endif diff --git a/itk/clitkSegmentationUtils.h b/itk/clitkSegmentationUtils.h index 6204c9c..f58b7c4 100644 --- a/itk/clitkSegmentationUtils.h +++ b/itk/clitkSegmentationUtils.h @@ -363,6 +363,10 @@ namespace clitk { void And(ImageType * input, const ImageType * object, typename ImageType::PixelType BG=0); + template + void Or(ImageType * input, + const ImageType * object, + typename ImageType::PixelType BG=0); //-------------------------------------------------------------------- diff --git a/itk/clitkSegmentationUtils.txx b/itk/clitkSegmentationUtils.txx index afcefcf..edb208a 100644 --- a/itk/clitkSegmentationUtils.txx +++ b/itk/clitkSegmentationUtils.txx @@ -353,9 +353,9 @@ namespace clitk { sliceRelPosFilter->AddOrientationTypeString(orientation); sliceRelPosFilter->SetIntermediateSpacingFlag((spacing != -1)); sliceRelPosFilter->SetIntermediateSpacing(spacing); - sliceRelPosFilter->SetUniqueConnectedComponentBySlice(uniqueConnectedComponent); - sliceRelPosFilter->CCLSelectionFlagOff(); - sliceRelPosFilter->SetUseASingleObjectConnectedComponentBySliceFlag(singleObjectCCL); + sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(uniqueConnectedComponent); + sliceRelPosFilter->ObjectCCLSelectionFlagOff(); + sliceRelPosFilter->SetUseTheLargestObjectCCLFlag(singleObjectCCL); // sliceRelPosFilter->SetInverseOrientationFlag(inverseflag); sliceRelPosFilter->SetAutoCropFlag(autocropFlag); sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); @@ -364,6 +364,7 @@ namespace clitk { } //-------------------------------------------------------------------- + //-------------------------------------------------------------------- template bool @@ -378,6 +379,7 @@ namespace clitk { } //-------------------------------------------------------------------- + //-------------------------------------------------------------------- template bool @@ -818,6 +820,7 @@ namespace clitk { int mainDirection, double offsetToKeep) { + assert((mainDirection==0) || (mainDirection==1)); typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; SliceIteratorType siter = SliceIteratorType(input, input->GetLargestPossibleRegion()); @@ -929,6 +932,34 @@ namespace clitk { //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + template + void + Or(ImageType * input, + const ImageType * object, + typename ImageType::PixelType BG) + { + typename ImageType::Pointer o; + bool resized=false; + if (!clitk::HaveSameSizeAndSpacing(input, object)) { + o = clitk::ResizeImageLike(object, input, BG); + resized = true; + } + + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(input); + if (resized) boolFilter->SetInput2(o); + else boolFilter->SetInput2(object); + boolFilter->SetBackgroundValue1(BG); + boolFilter->SetBackgroundValue2(BG); + boolFilter->SetOperationType(BoolFilterType::Or); + boolFilter->Update(); + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- template typename ImageType::Pointer diff --git a/itk/clitkSliceBySliceRelativePositionFilter.h b/itk/clitkSliceBySliceRelativePositionFilter.h index 53c5a87..f860c03 100644 --- a/itk/clitkSliceBySliceRelativePositionFilter.h +++ b/itk/clitkSliceBySliceRelativePositionFilter.h @@ -74,28 +74,28 @@ namespace clitk { itkGetConstMacro(Direction, int); itkSetMacro(Direction, int); - itkGetConstMacro(UniqueConnectedComponentBySlice, bool); - itkSetMacro(UniqueConnectedComponentBySlice, bool); - itkBooleanMacro(UniqueConnectedComponentBySlice); + itkGetConstMacro(UniqueConnectedComponentBySliceFlag, bool); + itkSetMacro(UniqueConnectedComponentBySliceFlag, bool); + itkBooleanMacro(UniqueConnectedComponentBySliceFlag); itkGetConstMacro(IgnoreEmptySliceObjectFlag, bool); itkSetMacro(IgnoreEmptySliceObjectFlag, bool); itkBooleanMacro(IgnoreEmptySliceObjectFlag); - itkGetConstMacro(UseASingleObjectConnectedComponentBySliceFlag, bool); - itkSetMacro(UseASingleObjectConnectedComponentBySliceFlag, bool); - itkBooleanMacro(UseASingleObjectConnectedComponentBySliceFlag); - - itkGetConstMacro(CCLSelectionFlag, bool); - itkSetMacro(CCLSelectionFlag, bool); - itkBooleanMacro(CCLSelectionFlag); - itkGetConstMacro(CCLSelectionDimension, int); - itkSetMacro(CCLSelectionDimension, int); - itkGetConstMacro(CCLSelectionDirection, int); - itkSetMacro(CCLSelectionDirection, int); - itkGetConstMacro(CCLSelectionIgnoreSingleCCLFlag, bool); - itkSetMacro(CCLSelectionIgnoreSingleCCLFlag, bool); - itkBooleanMacro(CCLSelectionIgnoreSingleCCLFlag); + itkGetConstMacro(UseTheLargestObjectCCLFlag, bool); + itkSetMacro(UseTheLargestObjectCCLFlag, bool); + itkBooleanMacro(UseTheLargestObjectCCLFlag); + + itkGetConstMacro(ObjectCCLSelectionFlag, bool); + itkSetMacro(ObjectCCLSelectionFlag, bool); + itkBooleanMacro(ObjectCCLSelectionFlag); + itkGetConstMacro(ObjectCCLSelectionDimension, int); + itkSetMacro(ObjectCCLSelectionDimension, int); + itkGetConstMacro(ObjectCCLSelectionDirection, int); + itkSetMacro(ObjectCCLSelectionDirection, int); + itkGetConstMacro(ObjectCCLSelectionIgnoreSingleCCLFlag, bool); + itkSetMacro(ObjectCCLSelectionIgnoreSingleCCLFlag, bool); + itkBooleanMacro(ObjectCCLSelectionIgnoreSingleCCLFlag); protected: SliceBySliceRelativePositionFilter(); @@ -109,14 +109,14 @@ namespace clitk { ImagePointer object; ImagePointer m_working_input; ImagePointer m_working_object; - bool m_UniqueConnectedComponentBySlice; + bool m_UniqueConnectedComponentBySliceFlag; int m_Direction; bool m_IgnoreEmptySliceObjectFlag; - bool m_UseASingleObjectConnectedComponentBySliceFlag; - bool m_CCLSelectionFlag; - int m_CCLSelectionDimension; - int m_CCLSelectionDirection; - bool m_CCLSelectionIgnoreSingleCCLFlag; + bool m_UseTheLargestObjectCCLFlag; + bool m_ObjectCCLSelectionFlag; + int m_ObjectCCLSelectionDimension; + int m_ObjectCCLSelectionDirection; + bool m_ObjectCCLSelectionIgnoreSingleCCLFlag; private: SliceBySliceRelativePositionFilter(const Self&); //purposely not implemented diff --git a/itk/clitkSliceBySliceRelativePositionFilter.txx b/itk/clitkSliceBySliceRelativePositionFilter.txx index f2452d8..d75cc88 100644 --- a/itk/clitkSliceBySliceRelativePositionFilter.txx +++ b/itk/clitkSliceBySliceRelativePositionFilter.txx @@ -31,16 +31,16 @@ SliceBySliceRelativePositionFilter(): clitk::AddRelativePositionConstraintToLabelImageFilter() { SetDirection(2); - UniqueConnectedComponentBySliceOff(); + UniqueConnectedComponentBySliceFlagOff(); SetIgnoreEmptySliceObjectFlag(false); - UseASingleObjectConnectedComponentBySliceFlagOn(); + UseTheLargestObjectCCLFlagOff(); this->VerboseStepFlagOff(); this->WriteStepFlagOff(); this->SetCombineWithOrFlag(false); - CCLSelectionFlagOff(); - SetCCLSelectionDimension(0); - SetCCLSelectionDirection(1); - CCLSelectionIgnoreSingleCCLFlagOff(); + ObjectCCLSelectionFlagOff(); + SetObjectCCLSelectionDimension(0); + SetObjectCCLSelectionDirection(1); + ObjectCCLSelectionIgnoreSingleCCLFlagOff(); } //-------------------------------------------------------------------- @@ -81,11 +81,15 @@ PrintOptions() DD(this->GetIntermediateSpacingFlag()); DD(this->GetIntermediateSpacing()); DD(this->GetFuzzyThreshold()); - DD(this->GetUniqueConnectedComponentBySlice()); + DD(this->GetUniqueConnectedComponentBySliceFlag()); DD(this->GetAutoCropFlag()); DD(this->GetInverseOrientationFlag()); DD(this->GetRemoveObjectFlag()); DD(this->GetCombineWithOrFlag()); + DD(this->GetUseTheLargestObjectCCLFlag()); + DD(this->GetObjectCCLSelectionFlag()); + DD(this->GetObjectCCLSelectionDimension()); + DD(this->GetObjectCCLSelectionIgnoreSingleCCLFlag()); } //-------------------------------------------------------------------- @@ -187,23 +191,23 @@ GenerateOutputInformation() if ((!GetIgnoreEmptySliceObjectFlag()) || (nb!=0)) { // Select or not a single CCL ? - if (GetUseASingleObjectConnectedComponentBySliceFlag()) { + if (GetUseTheLargestObjectCCLFlag()) { mObjectSlices[i] = KeepLabels(mObjectSlices[i], 0, 1, 1, 1, true); } // Select a single according to a position if more than one CCL - if (GetCCLSelectionFlag()) { + if (GetObjectCCLSelectionFlag()) { // if several CCL, choose the most extrema according a direction, // if not -> should we consider this slice ? if (nb<2) { - if (GetCCLSelectionIgnoreSingleCCLFlag()) { + if (GetObjectCCLSelectionIgnoreSingleCCLFlag()) { mObjectSlices[i] = SetBackground(mObjectSlices[i], mObjectSlices[i], 1, this->GetBackgroundValue(), true); } } - int dim = GetCCLSelectionDimension(); - int direction = GetCCLSelectionDirection(); + int dim = GetObjectCCLSelectionDimension(); + int direction = GetObjectCCLSelectionDirection(); std::vector centroids; ComputeCentroids(mObjectSlices[i], this->GetBackgroundValue(), centroids); uint index=1; @@ -222,7 +226,7 @@ GenerateOutputInformation() true); } } - } // end GetCCLSelectionFlag = true + } // end GetbjectCCLSelectionFlag = true // Relative position typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; @@ -234,10 +238,14 @@ GenerateOutputInformation() relPosFilter->SetInput(mInputSlices[i]); relPosFilter->SetInputObject(mObjectSlices[i]); relPosFilter->SetRemoveObjectFlag(this->GetRemoveObjectFlag()); + // This flag (InverseOrientation) *must* be set before + // AddOrientation because AddOrientation can change it. + relPosFilter->SetInverseOrientationFlag(this->GetInverseOrientationFlag()); for(int j=0; jGetNumberOfAngles(); j++) { relPosFilter->AddOrientationTypeString(this->GetOrientationTypeString(j)); + //DD(this->GetOrientationTypeString(j)); } - relPosFilter->SetInverseOrientationFlag(this->GetInverseOrientationFlag()); + //DD(this->GetInverseOrientationFlag()); //relPosFilter->SetOrientationType(this->GetOrientationType()); relPosFilter->SetIntermediateSpacing(this->GetIntermediateSpacing()); relPosFilter->SetIntermediateSpacingFlag(this->GetIntermediateSpacingFlag()); @@ -248,7 +256,7 @@ GenerateOutputInformation() mInputSlices[i] = relPosFilter->GetOutput(); // Select main CC if needed - if (GetUniqueConnectedComponentBySlice()) { + if (GetUniqueConnectedComponentBySliceFlag()) { mInputSlices[i] = Labelize(mInputSlices[i], 0, true, 1); mInputSlices[i] = KeepLabels(mInputSlices[i], 0, 1, 1, 1, true); } diff --git a/itk/itkImageToVTKImageFilter.h b/itk/itkImageToVTKImageFilter.h index 4edf4ae..a50e5b7 100644 --- a/itk/itkImageToVTKImageFilter.h +++ b/itk/itkImageToVTKImageFilter.h @@ -15,89 +15,87 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html ===========================================================================**/ -#ifndef __itkImageToVTKImageFilter_h -#define __itkImageToVTKImageFilter_h -#include "itkVTKImageExport.h" -#include "vtkImageImport.h" -#include "vtkImageData.h" - -namespace itk -{ - -/** \class ImageToVTKImageFilter - * \brief Converts an ITK image into a VTK image and plugs a - * itk data pipeline to a VTK datapipeline. - * - * This class puts together an itkVTKImageExporter and a vtkImageImporter. - * It takes care of the details related to the connection of ITK and VTK - * pipelines. The User will perceive this filter as an adaptor to which - * an itk::Image can be plugged as input and a vtkImage is produced as - * output. - * - * \ingroup ImageFilters - */ -template -class ITK_EXPORT ImageToVTKImageFilter : public ProcessObject -{ -public: - /** Standard class typedefs. */ - typedef ImageToVTKImageFilter Self; - typedef ProcessObject Superclass; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Run-time type information (and related methods). */ - itkTypeMacro(ImageToVTKImageFilter, ProcessObject); - - /** Some typedefs. */ - typedef TInputImage InputImageType; - typedef typename InputImageType::ConstPointer InputImagePointer; - typedef VTKImageExport< InputImageType> ExporterFilterType; - typedef typename ExporterFilterType::Pointer ExporterFilterPointer; - - /** Get the output in the form of a vtkImage. - This call is delegated to the internal vtkImageImporter filter */ - vtkImageData * GetOutput() const; - - /** Set the input in the form of an itk::Image */ - void SetInput( const InputImageType * ); - - /** Return the internal VTK image importer filter. - This is intended to facilitate users the access - to methods in the importer */ - vtkImageImport * GetImporter() const; - - /** Return the internal ITK image exporter filter. - This is intended to facilitate users the access - to methods in the exporter */ - ExporterFilterType * GetExporter() const; - - /** This call delegate the update to the importer */ - void Update(); - -protected: - ImageToVTKImageFilter(); - virtual ~ImageToVTKImageFilter(); - -private: - ImageToVTKImageFilter(const Self&); //purposely not implemented - void operator=(const Self&); //purposely not implemented - - ExporterFilterPointer m_Exporter; - vtkImageImport * m_Importer; - -}; - -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "itkImageToVTKImageFilter.txx" -#endif - -#endif - - - +#ifndef __itkImageToVTKImageFilter_h +#define __itkImageToVTKImageFilter_h +#include "itkVTKImageExport.h" +#include "vtkImageImport.h" +#include "vtkImageData.h" + +namespace itk +{ + +/** \class ImageToVTKImageFilter + * \brief Converts an ITK image into a VTK image and plugs a + * itk data pipeline to a VTK datapipeline. + * + * This class puts together an itkVTKImageExporter and a vtkImageImporter. + * It takes care of the details related to the connection of ITK and VTK + * pipelines. The User will perceive this filter as an adaptor to which + * an itk::Image can be plugged as input and a vtkImage is produced as + * output. + * + * \ingroup ImageFilters + */ +template +class ITK_EXPORT ImageToVTKImageFilter : public ProcessObject +{ +public: + /** Standard class typedefs. */ + typedef ImageToVTKImageFilter Self; + typedef ProcessObject Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(ImageToVTKImageFilter, ProcessObject); + + /** Some typedefs. */ + typedef TInputImage InputImageType; + typedef typename InputImageType::ConstPointer InputImagePointer; + typedef VTKImageExport< InputImageType> ExporterFilterType; + typedef typename ExporterFilterType::Pointer ExporterFilterPointer; + + /** Get the output in the form of a vtkImage. + This call is delegated to the internal vtkImageImporter filter */ + vtkImageData * GetOutput() const; + + /** Set the input in the form of an itk::Image */ + void SetInput( const InputImageType * ); + + /** Return the internal VTK image importer filter. + This is intended to facilitate users the access + to methods in the importer */ + vtkImageImport * GetImporter() const; + + /** Return the internal ITK image exporter filter. + This is intended to facilitate users the access + to methods in the exporter */ + ExporterFilterType * GetExporter() const; + + /** This call delegate the update to the importer */ + void Update(); + +protected: + ImageToVTKImageFilter(); + virtual ~ImageToVTKImageFilter(); + +private: + ImageToVTKImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + ExporterFilterPointer m_Exporter; + vtkImageImport * m_Importer; + +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkImageToVTKImageFilter.txx" +#endif + +#endif + diff --git a/segmentation/CMakeLists.txt b/segmentation/CMakeLists.txt index d6812e0..0cc7e2b 100644 --- a/segmentation/CMakeLists.txt +++ b/segmentation/CMakeLists.txt @@ -50,13 +50,15 @@ IF(CLITK_BUILD_SEGMENTATION) TARGET_LINK_LIBRARIES(clitkBool clitkCommon ${ITK_LIBRARIES} clitkSegmentationGgoLib) SET(SEGMENTATION_INSTALL ${SEGMENTATION_INSTALL} clitkBool) + WRAP_GGO(clitkRelativePosition_GGO_C ../tools/clitkRelativePosition.ggo) + WRAP_GGO(clitkExtractMediastinum_GGO_C clitkExtractMediastinum.ggo) - ADD_EXECUTABLE(clitkExtractMediastinum clitkExtractMediastinum.cxx ${clitkExtractMediastinum_GGO_C}) + ADD_EXECUTABLE(clitkExtractMediastinum clitkExtractMediastinum.cxx ${clitkExtractMediastinum_GGO_C} ${clitkRelativePosition_GGO_C}) TARGET_LINK_LIBRARIES(clitkExtractMediastinum clitkCommon clitkSegmentationGgoLib ${ITK_LIBRARIES}) SET(SEGMENTATION_INSTALL ${SEGMENTATION_INSTALL} clitkExtractMediastinum) WRAP_GGO(clitkExtractLymphStations_GGO_C clitkExtractLymphStations.ggo) - ADD_EXECUTABLE(clitkExtractLymphStations clitkExtractLymphStations.cxx clitkFilterWithAnatomicalFeatureDatabaseManagement.cxx ${clitkExtractLymphStations_GGO_C}) + ADD_EXECUTABLE(clitkExtractLymphStations clitkExtractLymphStations.cxx clitkFilterWithAnatomicalFeatureDatabaseManagement.cxx ${clitkExtractLymphStations_GGO_C} ${clitkRelativePosition_GGO_C}) TARGET_LINK_LIBRARIES(clitkExtractLymphStations clitkSegmentationGgoLib clitkCommon vtkHybrid) SET(SEGMENTATION_INSTALL ${SEGMENTATION_INSTALL} clitkExtractLymphStations) diff --git a/segmentation/clitkConnectedComponentLabeling.cxx b/segmentation/clitkConnectedComponentLabeling.cxx index 47ea003..c601e9f 100644 --- a/segmentation/clitkConnectedComponentLabeling.cxx +++ b/segmentation/clitkConnectedComponentLabeling.cxx @@ -24,26 +24,26 @@ //-------------------------------------------------------------------- int main(int argc, char * argv[]) { - clitk::PrintMemory(true, "start"); + // clitk::PrintMemory(true, "start"); // Init command line GGO(clitkConnectedComponentLabeling, args_info); CLITK_INIT; // Filter - clitk::PrintMemory(true, "before filter"); + //clitk::PrintMemory(true, "before filter"); typedef clitk::ConnectedComponentLabelingGenericFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetArgsInfo(args_info); - clitk::PrintMemory(true, "before update"); + //clitk::PrintMemory(true, "before update"); try { filter->Update(); } catch(std::runtime_error e) { std::cout << e.what() << std::endl; } - clitk::PrintMemory(true, "after filter"); + //clitk::PrintMemory(true, "after filter"); return EXIT_SUCCESS; } // This is the end, my friend diff --git a/segmentation/clitkConnectedComponentLabelingGenericFilter.txx b/segmentation/clitkConnectedComponentLabelingGenericFilter.txx index 3c57841..89154a5 100644 --- a/segmentation/clitkConnectedComponentLabelingGenericFilter.txx +++ b/segmentation/clitkConnectedComponentLabelingGenericFilter.txx @@ -32,7 +32,7 @@ template clitk::ConnectedComponentLabelingGenericFilter::ConnectedComponentLabelingGenericFilter(): ImageToImageGenericFilter("ConnectedComponentLabeling") { - // InitializeImageType<2>(); + InitializeImageType<2>(); InitializeImageType<3>(); //InitializeImageType<4>(); } @@ -47,7 +47,7 @@ void clitk::ConnectedComponentLabelingGenericFilter::InitializeIma ADD_IMAGE_TYPE(Dim, uchar); ADD_IMAGE_TYPE(Dim, short); // ADD_IMAGE_TYPE(Dim, int); - // ADD_IMAGE_TYPE(Dim, float); + ADD_IMAGE_TYPE(Dim, float); } //-------------------------------------------------------------------- @@ -72,14 +72,12 @@ template template void clitk::ConnectedComponentLabelingGenericFilter::UpdateWithInputImageType() { - DD("UpdateWithInputImageType"); - // Reading input typename ImageType::Pointer input = this->template GetInput(0); // Output image type typedef itk::Image OutputImageType; - PrintMemory(true, "initial"); + //PrintMemory(true, "initial"); typename OutputImageType::Pointer output; { @@ -95,11 +93,11 @@ void clitk::ConnectedComponentLabelingGenericFilter::UpdateWithInp // connectFilter->SetNumberOfThreads(8); connectFilter->Update(); temp = connectFilter->GetOutput(); - PrintMemory(true, "after udpate"); + // PrintMemory(true, "after udpate"); } - PrintMemory(true, "after CCL block"); - DD(input->GetReferenceCount()); - DD(temp->GetReferenceCount()); + // PrintMemory(true, "after CCL block"); + // DD(input->GetReferenceCount()); + // DD(temp->GetReferenceCount()); // Sort by size and remove too small area. typedef itk::RelabelComponentImageFilter RelabelFilterType; @@ -109,9 +107,9 @@ void clitk::ConnectedComponentLabelingGenericFilter::UpdateWithInp relabelFilter->SetMinimumObjectSize(mArgsInfo.minSize_arg); relabelFilter->Update(); - DD(mArgsInfo.inputBG_arg); - DD(mArgsInfo.full_flag); - DD(mArgsInfo.minSize_arg); + // DD(mArgsInfo.inputBG_arg); + // DD(mArgsInfo.full_flag); + // DD(mArgsInfo.minSize_arg); // Set information const std::vector & a @@ -120,16 +118,16 @@ void clitk::ConnectedComponentLabelingGenericFilter::UpdateWithInp for(unsigned int i=0; iGetSizeOfObjectsInPhysicalUnits(); m_OriginalNumberOfObjects = relabelFilter->GetOriginalNumberOfObjects(); - DD(m_OriginalNumberOfObjects); - DD(m_SizeOfObjectsInPhysicalUnits.size()); + // DD(m_OriginalNumberOfObjects); + // DD(m_SizeOfObjectsInPhysicalUnits.size()); output = relabelFilter->GetOutput(); } - PrintMemory(true, "after block"); + // PrintMemory(true, "after block"); // Write/Save results this->template SetNextOutput(output); - PrintMemory(true, "end filter "); + // PrintMemory(true, "end filter "); } //-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractLymphStation_1RL.txx b/segmentation/clitkExtractLymphStation_1RL.txx new file mode 100644 index 0000000..608fdc9 --- /dev/null +++ b/segmentation/clitkExtractLymphStation_1RL.txx @@ -0,0 +1,185 @@ + + +// vtk +#include +#include +#include + +// clitk +#include "clitkMeshToBinaryImageFilter.h" + +// itk +#include + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_1RL_SetDefaultValues() +{ + +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_1RL() +{ + if ((!CheckForStation("1R")) && (!CheckForStation("1L"))) return; + + StartNewStep("Stations 1RL"); + StartSubStep(); + + // Get the current support + StartNewStep("[Station 1RL] Get the current 1RL suppport"); + m_ListOfStations["1R"] = m_ListOfSupports["S1R"]; + m_ListOfStations["1L"] = m_ListOfSupports["S1L"]; + StopCurrentStep(m_ListOfStations["1R"]); + + // Specific processes + ExtractStation_1RL_Ant_Limits(); + ExtractStation_1RL_Post_Limits(); + m_Working_Support = m_ListOfStations["1R"]; + Remove_Structures(" 1R", "ScaleneMuscleAnt"); + Remove_Structures(" 1R", "CommonCarotidArteryRight"); + m_Working_Support = m_ListOfStations["1L"]; + Remove_Structures(" 1L", "ScaleneMuscleAnt"); + Remove_Structures(" 1L", "CommonCarotidArteryLeft"); + + // Generic RelativePosition processes + m_ListOfStations["1R"] = this->ApplyRelativePositionList("Station_1R", m_ListOfStations["1R"]); + m_ListOfStations["1L"] = this->ApplyRelativePositionList("Station_1L", m_ListOfStations["1L"]); + + // Store image filenames into AFDB + writeImage(m_ListOfStations["1R"], "seg/Station1R.mhd"); + writeImage(m_ListOfStations["1L"], "seg/Station1L.mhd"); + GetAFDB()->SetImageFilename("Station1R", "seg/Station1R.mhd"); + GetAFDB()->SetImageFilename("Station1L", "seg/Station1L.mhd"); + WriteAFDB(); + StopSubStep(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_1RL_Ant_Limits() +{ + // ----------------------------------------------------- + StartNewStep("[Station 1RL] anterior limits with Trachea and Thyroid"); + + /* + The idea here it to consider the most anterior points int the + Trachea or the Thyroid and cut all parts anterior to it. This is + an heuristic, not written explicitely in the articles. + */ + + MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); + MaskImagePointer Thyroid = GetAFDB()->template GetImage("Thyroid"); + MaskImagePointer S1R = m_ListOfStations["1R"]; + MaskImagePointer S1L = m_ListOfStations["1L"]; + MaskImagePointer support = m_ListOfSupports["S1RL"]; + + // Resize like S1R. Warning: for Thyroid, only Right part is thus + // taken into account + Trachea = clitk::ResizeImageLike(Trachea, S1R, GetBackgroundValue()); + Thyroid = clitk::ResizeImageLike(Thyroid, S1R, GetBackgroundValue()); + + // Search for most Ant point, slice by slice, between Trachea and Thyroid + std::vector Trachea_slices; + clitk::ExtractSlices(Trachea, 2, Trachea_slices); + std::vector Thyroid_slices; + clitk::ExtractSlices(Thyroid, 2, Thyroid_slices); + std::vector A; + std::vector B; + for(uint i=0; i(Trachea_slices[i], + GetBackgroundValue(), + 1, true, p); + FindExtremaPointInAGivenDirection(Thyroid_slices[i], + GetBackgroundValue(), + 1, true, q); + if (q[1] < p[1]) p = q; // Now p is the most ant. + // Add a little margin, 3mm + p[1] -= 3; + // Convert in 3D + ImagePointType x; //dummy + A.push_back(x); + B.push_back(x); + clitk::PointsUtils::Convert2DTo3D(p, Trachea, i, A[i]); + p[0] += 10; + clitk::PointsUtils::Convert2DTo3D(p, Trachea, i, B[i]); + } + + // Remove anterior to this line (keep +10 offset) + clitk::SliceBySliceSetBackgroundFromLineSeparation(S1R, A, B, + GetBackgroundValue(), 1, +10); + m_ListOfStations["1R"] = S1R; + + clitk::SliceBySliceSetBackgroundFromLineSeparation(S1L, A, B, + GetBackgroundValue(), 1, +10); + StopCurrentStep(S1L); + m_ListOfStations["1L"] = S1L; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_1RL_Post_Limits() +{ + // ----------------------------------------------------- + StartNewStep("[Station 1RL] Posterior limits with VertebralArtery"); + + MaskImagePointer VertebralArtery = GetAFDB()->template GetImage("VertebralArtery"); + MaskImagePointer S1R = m_ListOfStations["1R"]; + MaskImagePointer S1L = m_ListOfStations["1L"]; + + // Resize like S1R. + VertebralArtery = clitk::ResizeImageLike(VertebralArtery, S1R, GetBackgroundValue()); + + // Search for most Ant point, slice by slice, between Trachea and Thyroid + std::vector VertebralArtery_slices; + clitk::ExtractSlices(VertebralArtery, 2, VertebralArtery_slices); + std::vector A; + std::vector B; + for(uint i=0; i(VertebralArtery_slices[i], + GetBackgroundValue(), + 1, false, p); + // Add a little margin ? No. + //p[1] += 0; + //DD(p); + // Convert in 3D + ImagePointType x; //dummy + A.push_back(x); + B.push_back(x); + clitk::PointsUtils::Convert2DTo3D(p, VertebralArtery, i, A[i]); + p[0] += 10; + clitk::PointsUtils::Convert2DTo3D(p, VertebralArtery, i, B[i]); + } + + // Remove anterior to this line (keep -10 offset) + clitk::SliceBySliceSetBackgroundFromLineSeparation(S1R, A, B, + GetBackgroundValue(), 1, -10); + m_ListOfStations["1R"] = S1R; + + // Remove anterior to this line (keep -10 offset) + clitk::SliceBySliceSetBackgroundFromLineSeparation(S1L, A, B, + GetBackgroundValue(), 1, -10); + StopCurrentStep(S1L); + m_ListOfStations["1L"] = S1L; +} +//-------------------------------------------------------------------- + + diff --git a/segmentation/clitkExtractLymphStation_2RL.txx b/segmentation/clitkExtractLymphStation_2RL.txx index 71298eb..e31f2d5 100644 --- a/segmentation/clitkExtractLymphStation_2RL.txx +++ b/segmentation/clitkExtractLymphStation_2RL.txx @@ -11,287 +11,58 @@ // itk #include -//-------------------------------------------------------------------- -template -class comparePointsX { -public: - bool operator() (PointType i, PointType j) { return (i[0] -class comparePointsWithAngle { -public: - bool operator() (PairType i, PairType j) { return (i.second < j.second); } -}; -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void HypercubeCorners(std::vector > & out) { - std::vector > previous; - HypercubeCorners(previous); - out.resize(previous.size()*2); - for(unsigned int i=0; i p; - if (i -void HypercubeCorners<1>(std::vector > & out) { - out.resize(2); - out[0][0] = 0; - out[1][0] = 1; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image, - std::vector & bounds) -{ - // Get image max/min coordinates - const unsigned int dim=ImageType::ImageDimension; - typedef typename ImageType::PointType PointType; - typedef typename ImageType::IndexType IndexType; - PointType min_c, max_c; - IndexType min_i, max_i; - min_i = image->GetLargestPossibleRegion().GetIndex(); - for(unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i] + min_i[i]; - image->TransformIndexToPhysicalPoint(min_i, min_c); - image->TransformIndexToPhysicalPoint(max_i, max_c); - - // Get corners coordinates - HypercubeCorners(bounds); - for(unsigned int i=0; i void clitk::ExtractLymphStationsFilter:: ExtractStation_2RL_SetDefaultValues() { - SetFuzzyThreshold("2RL", "CommonCarotidArtery", 0.7); - SetFuzzyThreshold("2RL", "BrachioCephalicArtery", 0.7); - SetFuzzyThreshold("2RL", "BrachioCephalicVein", 0.3); - SetFuzzyThreshold("2RL", "Aorta", 0.7); - SetFuzzyThreshold("2RL", "SubclavianArteryRight", 0.5); - SetFuzzyThreshold("2RL", "SubclavianArteryLeft", 0.8); } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template +template void -clitk::ExtractLymphStationsFilter:: -ExtractStation_2RL_SI_Limits() +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL() { - // Apex of the chest or Sternum & Carina. - StartNewStep("[Station 2RL] Inf/Sup limits with Sternum and TopOfAorticArch/CaudalMarginOfLeftBrachiocephalicVein"); - - /* Rod says: "For the inferior border, unlike in Atlas – UM, there - is now a difference between 2R and 2L. 2R stops at the - intersection of the caudal margin of the innominate vein with the - trachea. 2L extends less inferiorly to the superior border of the - aortic arch." */ - - /* Get TopOfAorticArch and CaudalMarginOfLeftBrachiocephalicVein - - TopOfAorticArch -> can be obtain from Aorta -> most sup part. - - CaudalMarginOfLeftBrachiocephalicVein -> must inf part of BrachicephalicVein - */ - MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); - MaskImagePointType p = Aorta->GetOrigin(); // initialise to avoid warning - clitk::FindExtremaPointInAGivenDirection(Aorta, GetBackgroundValue(), 2, false, p); - double TopOfAorticArchZ=p[2]; - GetAFDB()->SetDouble("TopOfAorticArchZ", TopOfAorticArchZ); - - MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); - clitk::FindExtremaPointInAGivenDirection(BrachioCephalicVein, GetBackgroundValue(), 2, true, p); - double CaudalMarginOfLeftBrachiocephalicVeinZ=p[2]; - GetAFDB()->SetDouble("CaudalMarginOfLeftBrachiocephalicVeinZ", CaudalMarginOfLeftBrachiocephalicVeinZ); - - // First, cut on the most inferior part. Add one slice because this - // inf slice should not be included. - double inf = std::min(CaudalMarginOfLeftBrachiocephalicVeinZ, TopOfAorticArchZ) + m_Working_Support->GetSpacing()[2]; + if ((!CheckForStation("2R")) && (!CheckForStation("2L"))) return; - // Get Sternum and search for the upper position - MaskImagePointer Sternum = GetAFDB()->template GetImage("Sternum"); + StartNewStep("Stations 2RL"); + StartSubStep(); - // Search most sup point - clitk::FindExtremaPointInAGivenDirection(Sternum, GetBackgroundValue(), 2, false, p); - double m_SternumZ = p[2]; - // +Sternum->GetSpacing()[2]; // One more slice, because it is below this point // NOT HERE ?? WHY DIFFERENT FROM 3A ?? - - //* Crop support : - m_Working_Support = - clitk::CropImageAlongOneAxis(m_Working_Support, 2, - inf, m_SternumZ, true, - GetBackgroundValue()); - - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_2RL_Ant_Limits() -{ - // ----------------------------------------------------- - /* Rod says: "The anterior border, as with the Atlas – UM, is - posterior to the vessels (right subclavian vein, left - brachiocephalic vein, right brachiocephalic vein, left subclavian - artery, left common carotid artery and brachiocephalic trunk). - These vessels are not included in the nodal station. The anterior - border is drawn to the midpoint of the vessel and an imaginary - line joins the middle of these vessels. Between the vessels, - station 2 is in contact with station 3a." */ - - // ----------------------------------------------------- - StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery"); - - // Get CommonCarotidArtery - MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage("CommonCarotidArtery"); - - // Remove Ant to CommonCarotidArtery - typedef SliceBySliceRelativePositionFilter SliceRelPosFilterType; - typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); - sliceRelPosFilter->VerboseStepFlagOff(); - sliceRelPosFilter->WriteStepFlagOff(); - sliceRelPosFilter->SetInput(m_Working_Support); - sliceRelPosFilter->SetInputObject(CommonCarotidArtery); - sliceRelPosFilter->SetDirection(2); - sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "CommonCarotidArtery")); - sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); - sliceRelPosFilter->IntermediateSpacingFlagOn(); - sliceRelPosFilter->SetIntermediateSpacing(2); - sliceRelPosFilter->UniqueConnectedComponentBySliceOff(); - sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff(); - sliceRelPosFilter->AutoCropFlagOn(); - sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); - sliceRelPosFilter->RemoveObjectFlagOff(); - sliceRelPosFilter->Update(); - m_Working_Support = sliceRelPosFilter->GetOutput(); - - // End - StopCurrentStep(m_Working_Support); + // Get the current support + StartNewStep("[Station 2RL] Get the current 2RL suppport"); + m_ListOfStations["2R"] = m_ListOfSupports["S2R"]; + m_ListOfStations["2L"] = m_ListOfSupports["S2L"]; + StopCurrentStep(m_ListOfStations["2R"]); + + // Do the same limits for R & L + m_Working_Support = m_ListOfStations["2R"]; + ExtractStation_2RL_Ant_Limits("2R"); + ExtractStation_2RL_Remove_Structures(" 2R"); m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; - - // ----------------------------------------------------- - // Remove Ant to H line from the Ant most part of the - // CommonCarotidArtery until we reach the first slice of - // BrachioCephalicArtery - StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery, H line"); - - // First, find the first slice of BrachioCephalicArtery - MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage("BrachioCephalicArtery"); - MaskImagePointType p = BrachioCephalicArtery->GetOrigin(); // initialise to avoid warning - clitk::FindExtremaPointInAGivenDirection(BrachioCephalicArtery, GetBackgroundValue(), 2, false, p); - double TopOfBrachioCephalicArteryZ=p[2] + BrachioCephalicArtery->GetSpacing()[2]; // Add one slice - - // Remove CommonCarotidArtery below this point - CommonCarotidArtery = clitk::CropImageRemoveLowerThan(CommonCarotidArtery, 2, TopOfBrachioCephalicArteryZ, true, GetBackgroundValue()); - // Find most Ant points - std::vector ccaAntPositionA; - std::vector ccaAntPositionB; - clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(CommonCarotidArtery, - GetBackgroundValue(), 2, - 1, true, // Ant - 0, // Horizontal line - -3, // margin - ccaAntPositionA, - ccaAntPositionB); - // Cut ant to this line - clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, - ccaAntPositionA, - ccaAntPositionB, - GetBackgroundValue(), 1, 10); - - // End - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; + m_Working_Support = m_ListOfStations["2L"]; + ExtractStation_2RL_Ant_Limits("2L"); + ExtractStation_2RL_Remove_Structures(" 2L"); m_ListOfStations["2L"] = m_Working_Support; - // ----------------------------------------------------- - // Ant limit with the BrachioCephalicArtery - StartNewStep("[Station 2RL] Ant limits with BrachioCephalicArtery line"); - - // Remove Ant to BrachioCephalicArtery - m_Working_Support = - clitk::SliceBySliceRelativePosition(m_Working_Support, BrachioCephalicArtery, 2, - GetFuzzyThreshold("2RL", "BrachioCephalicArtery"), "NotAntTo", false, 2, true, false); - // End - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; + // Remove superior part to BrachioCephalicVein (used then by RelPos) + ExtractStation_2RL_Cut_BrachioCephalicVein_superiorly_when_it_split(); - // ----------------------------------------------------- - // Ant limit with the BrachioCephalicArtery H line - StartNewStep("[Station 2RL] Ant limits with BrachioCephalicArtery, Horizontal line"); + // Generic RelativePosition processes + m_ListOfStations["2R"] = this->ApplyRelativePositionList("Station_2R", m_ListOfStations["2R"]); + m_ListOfStations["2L"] = this->ApplyRelativePositionList("Station_2L", m_ListOfStations["2L"]); - // Find most Ant points - std::vector bctAntPositionA; - std::vector bctAntPositionB; - clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(BrachioCephalicArtery, - GetBackgroundValue(), 2, - 1, true, // Ant - 0, // Horizontal line - -1, // margin - bctAntPositionA, - bctAntPositionB); - // Cut ant to this line - clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, - bctAntPositionA, - bctAntPositionB, - GetBackgroundValue(), 1, 10); - // End - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; - - // ----------------------------------------------------- - // Ant limit with the BrachioCephalicVein - StartNewStep("[Station 2RL] Ant limits with BrachioCephalicVein"); - - // Remove Ant to BrachioCephalicVein - MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); - m_Working_Support = - clitk::SliceBySliceRelativePosition(m_Working_Support, BrachioCephalicVein, 2, - GetFuzzyThreshold("2RL", "BrachioCephalicVein"), "NotAntTo", false, 2, true, false); - // End - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; + // Store image filenames into AFDB + writeImage(m_ListOfStations["2R"], "seg/Station2R.mhd"); + writeImage(m_ListOfStations["2L"], "seg/Station2L.mhd"); + GetAFDB()->SetImageFilename("Station2R", "seg/Station2R.mhd"); + GetAFDB()->SetImageFilename("Station2L", "seg/Station2L.mhd"); + WriteAFDB(); + StopSubStep(); } //-------------------------------------------------------------------- @@ -300,299 +71,22 @@ ExtractStation_2RL_Ant_Limits() template void clitk::ExtractLymphStationsFilter:: -ExtractStation_2RL_Ant_Limits2() +ExtractStation_2RL_Ant_Limits(std::string s) { // ----------------------------------------------------- - /* Rod says: "The anterior border, as with the Atlas – UM, is - posterior to the vessels (right subclavian vein, left - brachiocephalic vein, right brachiocephalic vein, left subclavian - artery, left common carotid artery and brachiocephalic trunk). - These vessels are not included in the nodal station. The anterior - border is drawn to the midpoint of the vessel and an imaginary - line joins the middle of these vessels. Between the vessels, - station 2 is in contact with station 3a." */ - - // ----------------------------------------------------- - StartNewStep("[Station 2RL] Ant limits with vessels centroids"); - - /* Here, we consider the vessels form a kind of anterior barrier. We - link all vessels centroids and remove what is post to it. - - select the list of structure - vessel1 = BrachioCephalicArtery - vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right) - vessel3 = CommonCarotidArtery - vessel4 = SubclavianArtery - other = Thyroid - other = Aorta - - crop images as needed - - slice by slice, choose the CCL and extract centroids - - slice by slice, sort according to polar angle wrt Trachea centroid. - - slice by slice, link centroids and close contour - - remove outside this contour - - merge with support - */ - - // Read structures - MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage("BrachioCephalicArtery"); - MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); - MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage("CommonCarotidArtery"); - MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage("SubclavianArtery"); - MaskImagePointer Thyroid = GetAFDB()->template GetImage("Thyroid"); - MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); - MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); + StartNewStep("[Station "+s+"] Ant limits with vessels centroids"); - // Resize all structures like support - BrachioCephalicArtery = - clitk::ResizeImageLike(BrachioCephalicArtery, m_Working_Support, GetBackgroundValue()); - CommonCarotidArtery = - clitk::ResizeImageLike(CommonCarotidArtery, m_Working_Support, GetBackgroundValue()); - SubclavianArtery = - clitk::ResizeImageLike(SubclavianArtery, m_Working_Support, GetBackgroundValue()); - Thyroid = - clitk::ResizeImageLike(Thyroid, m_Working_Support, GetBackgroundValue()); - Aorta = - clitk::ResizeImageLike(Aorta, m_Working_Support, GetBackgroundValue()); - BrachioCephalicVein = - clitk::ResizeImageLike(BrachioCephalicVein, m_Working_Support, GetBackgroundValue()); - Trachea = - clitk::ResizeImageLike(Trachea, m_Working_Support, GetBackgroundValue()); - - // Extract slices - std::vector slices_BrachioCephalicArtery; - clitk::ExtractSlices(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery); - std::vector slices_BrachioCephalicVein; - clitk::ExtractSlices(BrachioCephalicVein, 2, slices_BrachioCephalicVein); - std::vector slices_CommonCarotidArtery; - clitk::ExtractSlices(CommonCarotidArtery, 2, slices_CommonCarotidArtery); - std::vector slices_SubclavianArtery; - clitk::ExtractSlices(SubclavianArtery, 2, slices_SubclavianArtery); - std::vector slices_Thyroid; - clitk::ExtractSlices(Thyroid, 2, slices_Thyroid); - std::vector slices_Aorta; - clitk::ExtractSlices(Aorta, 2, slices_Aorta); - std::vector slices_Trachea; - clitk::ExtractSlices(Trachea, 2, slices_Trachea); - unsigned int n= slices_BrachioCephalicArtery.size(); - - // Get the boundaries of one slice - std::vector bounds; - ComputeImageBoundariesCoordinates(slices_BrachioCephalicArtery[0], bounds); - - // For all slices, for all structures, find the centroid and build the contour - // List of 3D points (for debug) - std::vector p3D; - - vtkSmartPointer append = vtkSmartPointer::New(); - for(unsigned int i=0; i(slices_CommonCarotidArtery[i], - GetBackgroundValue(), true, 1); - slices_SubclavianArtery[i] = Labelize(slices_SubclavianArtery[i], - GetBackgroundValue(), true, 1); - slices_BrachioCephalicArtery[i] = Labelize(slices_BrachioCephalicArtery[i], - GetBackgroundValue(), true, 1); - slices_BrachioCephalicVein[i] = Labelize(slices_BrachioCephalicVein[i], - GetBackgroundValue(), true, 1); - slices_Thyroid[i] = Labelize(slices_Thyroid[i], - GetBackgroundValue(), true, 1); - slices_Aorta[i] = Labelize(slices_Aorta[i], - GetBackgroundValue(), true, 1); - - // Search centroids - std::vector points2D; - std::vector centroids1; - std::vector centroids2; - std::vector centroids3; - std::vector centroids4; - std::vector centroids5; - std::vector centroids6; - ComputeCentroids(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1); - ComputeCentroids(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2); - ComputeCentroids(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3); - ComputeCentroids(slices_Thyroid[i], GetBackgroundValue(), centroids4); - ComputeCentroids(slices_Aorta[i], GetBackgroundValue(), centroids5); - ComputeCentroids(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6); - - // BrachioCephalicVein -> when it is separated into two CCL, we - // only consider the most at Right one - if (centroids6.size() > 2) { - if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2); - else centroids6.erase(centroids6.begin()+1); - } - - // BrachioCephalicVein -> when SubclavianArtery has 2 CCL - // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein - if ((centroids3.size() ==1) && (centroids2.size() > 2)) { - centroids6.clear(); - } - - for(unsigned int j=1; j centroids_trachea; - ComputeCentroids(slices_Trachea[i], GetBackgroundValue(), centroids_trachea); - typedef std::pair PointAngleType; - std::vector angles; - for(unsigned int j=0; j0) angle = atan(y/x); - if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI; - if ((x<0) && (y<0)) angle = atan(y/x)-M_PI; - if (x==0) { - if (y>0) angle = M_PI/2.0; - if (y<0) angle = -M_PI/2.0; - if (y==0) angle = 0; - } - angle = clitk::rad2deg(angle); - // Angle is [-180;180] wrt the X axis. We change the X axis to - // be the vertical line, because we want to sort from Right to - // Left from Post to Ant. - if (angle>0) - angle = (270-angle); - if (angle<0) { - angle = -angle-90; - if (angle<0) angle = 360-angle; - } - PointAngleType a; - a.first = points2D[j]; - a.second = angle; - angles.push_back(a); - } - - // Do nothing if less than 2 points --> n - if (points2D.size() < 3) { //continue; - continue; - } - - // Sort points2D according to polar angles - std::sort(angles.begin(), angles.end(), comparePointsWithAngle()); - for(unsigned int j=0; j toadd; - unsigned int index = 0; - double dmax = 5; - while (indexdmax) { - - MaskSlicePointType b; - b[0] = a[0]+(c[0]-a[0])/2.0; - b[1] = a[1]+(c[1]-a[1])/2.0; - - // Compute distance to trachea centroid - MaskSlicePointType m = centroids_trachea[1]; - double da = m.EuclideanDistanceTo(a); - double db = m.EuclideanDistanceTo(b); - //double dc = m.EuclideanDistanceTo(c); - - // Mean distance, find point on the line from trachea centroid - // to b - double alpha = (da+db)/2.0; - MaskSlicePointType v; - double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2)); - v[0] = (b[0]-m[0])/n; - v[1] = (b[1]-m[1])/n; - MaskSlicePointType r; - r[0] = m[0]+alpha*v[0]; - r[1] = m[1]+alpha*v[1]; - points2D.insert(points2D.begin()+index+1, r); - } - else { - index++; - } - } - // DDV(points2D, points2D.size()); - - // Add some points to close the contour - // - H line towards Right - MaskSlicePointType p = points2D[0]; - p[0] = bounds[0][0]; - points2D.insert(points2D.begin(), p); - // - corner Right/Post - p = bounds[0]; - points2D.insert(points2D.begin(), p); - // - H line towards Left - p = points2D.back(); - p[0] = bounds[2][0]; - points2D.push_back(p); - // - corner Right/Post - p = bounds[2]; - points2D.push_back(p); - // Close contour with the first point - points2D.push_back(points2D[0]); - // DDV(points2D, points2D.size()); - - // Build 3D points from the 2D points - std::vector points3D; - clitk::PointsUtils::Convert2DListTo3DList(points2D, i, m_Working_Support, points3D); - for(unsigned int x=0; x mesh = Build3DMeshFrom2DContour(points3D); - append->AddInput(mesh); - } - - // DEBUG: write points3D - clitk::WriteListOfLandmarks(p3D, "vessels-centroids.txt"); - - // Build the final 3D mesh form the list 2D mesh - append->Update(); - vtkSmartPointer mesh = append->GetOutput(); - - // Debug, write the mesh - /* - vtkSmartPointer w = vtkSmartPointer::New(); - w->SetInput(mesh); - w->SetFileName("bidon.vtk"); - w->Write(); - */ - - // Compute a single binary 3D image from the list of contours - clitk::MeshToBinaryImageFilter::Pointer filter = - clitk::MeshToBinaryImageFilter::New(); - filter->SetMesh(mesh); - filter->SetLikeImage(m_Working_Support); - filter->Update(); - MaskImagePointer binarizedContour = filter->GetOutput(); - - // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse - ImagePointType p = p3D[2]; // This is the first centroid of the first slice - p[1] += 50; // 50 mm Post from this point must be kept - ImageIndexType index; - binarizedContour->TransformPhysicalPointToIndex(p, index); - bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue()); - + // WARNING, as I used "And" after, empty slice in binarizedContour + // lead to remove part of the support, although we want to keep + // unchanged. So we decide to ResizeImageLike but pad with + // ForegroundValue instead of BG + + // Get or compute the binary mask that separate Ant/Post part + // according to vessels + MaskImagePointer binarizedContour = FindAntPostVessels2(); + binarizedContour = clitk::ResizeImageLike(binarizedContour, + m_Working_Support, + GetForegroundValue()); // remove from support typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); @@ -601,17 +95,12 @@ ExtractStation_2RL_Ant_Limits2() boolFilter->SetInput2(binarizedContour); boolFilter->SetBackgroundValue1(GetBackgroundValue()); boolFilter->SetBackgroundValue2(GetBackgroundValue()); - if (isInside) - boolFilter->SetOperationType(BoolFilterType::And); - else - boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->SetOperationType(BoolFilterType::And); boolFilter->Update(); m_Working_Support = boolFilter->GetOutput(); - + // End StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; } //-------------------------------------------------------------------- @@ -620,35 +109,43 @@ ExtractStation_2RL_Ant_Limits2() template void clitk::ExtractLymphStationsFilter:: -ExtractStation_2RL_Post_Limits() +ExtractStation_2RL_Cut_BrachioCephalicVein_superiorly_when_it_split() { - StartNewStep("[Station 2RL] Post limits with post wall of Trachea"); - - // Get Trachea - MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); + // ----------------------------------------------------- + StartNewStep("[Station 2RL] Cut BrachioCephalicVein superiorly (when it split)"); - // Resize like the current support (to have the same number of slices) - Trachea = clitk::ResizeImageLike(Trachea, m_Working_Support, GetBackgroundValue()); + // Get BrachioCephalicVein + MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); + + // Remove the part superior to the slice where BrachioCephalicVein + // divide in two CCL + std::vector BCV_slices; + clitk::ExtractSlices(BrachioCephalicVein, 2, BCV_slices); + bool stop = false; + uint i=0; + while (!stop) { + // Count the number of CCL + int nb; + clitk::LabelizeAndCountNumberOfObjects(BCV_slices[i], GetBackgroundValue(), true, 1, nb); + if (nb>1) stop = true; + i++; + } + // Convert slice into coordinate + MaskImagePointType p; + MaskImageIndexType index; + index[0] = index[1] = 0; + index[2] = i; + BrachioCephalicVein->TransformIndexToPhysicalPoint(index, p); + BrachioCephalicVein = + clitk::CropImageRemoveGreaterThan(BrachioCephalicVein, 2, + p[2], true, + GetBackgroundValue()); - // Find extrema post positions - std::vector tracheaPostPositionsA; - std::vector tracheaPostPositionsB; - clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(Trachea, - GetBackgroundValue(), 2, - 1, false, // Post - 0, // Horizontal line - 1, - tracheaPostPositionsA, - tracheaPostPositionsB); - // Cut post to this line - clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, - tracheaPostPositionsA, - tracheaPostPositionsB, - GetBackgroundValue(), 1, -10); - // END - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; + // Now, insert this image in the AFDB (but do not store on disk) + GetAFDB()->template SetImage("BrachioCephalicVein_ForS2RL", "bidon", + BrachioCephalicVein, false); + // End + StopCurrentStep(BrachioCephalicVein); } //-------------------------------------------------------------------- @@ -682,173 +179,18 @@ Build3DMeshFrom2DContour(const std::vector & points) template void clitk::ExtractLymphStationsFilter:: -ExtractStation_2RL_LR_Limits() +ExtractStation_2RL_Remove_Structures(std::string s) { - // --------------------------------------------------------------------------- - StartNewStep("[Station 2RL] Left/Right limits with Aorta"); - MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); - // DD(GetFuzzyThreshold("2RL", "BrachioCephalicVein")); - m_Working_Support = - clitk::SliceBySliceRelativePosition(m_Working_Support, Aorta, 2, - GetFuzzyThreshold("2RL", "Aorta"), - "RightTo", false, 2, true, false); - // END - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; - - // --------------------------------------------------------------------------- - StartNewStep("[Station 2RL] Left/Right limits with SubclavianArtery (Right)"); - - // SliceBySliceRelativePosition + select CCL most at Right - MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage("SubclavianArtery"); - typedef SliceBySliceRelativePositionFilter SliceRelPosFilterType; - typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); - sliceRelPosFilter->VerboseStepFlagOff(); - sliceRelPosFilter->WriteStepFlagOff(); - sliceRelPosFilter->SetInput(m_Working_Support); - sliceRelPosFilter->SetInputObject(SubclavianArtery); - sliceRelPosFilter->SetDirection(2); - sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "SubclavianArteryRight")); - sliceRelPosFilter->AddOrientationTypeString("NotRightTo"); - sliceRelPosFilter->IntermediateSpacingFlagOn(); - sliceRelPosFilter->SetIntermediateSpacing(2); - sliceRelPosFilter->UniqueConnectedComponentBySliceOff(); - sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff(); - - sliceRelPosFilter->CCLSelectionFlagOn(); // select one CCL by slice - sliceRelPosFilter->SetCCLSelectionDimension(0); // select according to X (0) axis - sliceRelPosFilter->SetCCLSelectionDirection(-1); // select most at Right - sliceRelPosFilter->CCLSelectionIgnoreSingleCCLFlagOn(); // ignore if only one CCL - - sliceRelPosFilter->AutoCropFlagOn(); - sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); - sliceRelPosFilter->RemoveObjectFlagOff(); - sliceRelPosFilter->Update(); - m_Working_Support = sliceRelPosFilter->GetOutput(); - - // END - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; - - - // --------------------------------------------------------------------------- - StartNewStep("[Station 2RL] Left/Right limits with SubclavianArtery (Left)"); - - // SliceBySliceRelativePosition + select CCL most at Right - sliceRelPosFilter = SliceRelPosFilterType::New(); - sliceRelPosFilter->VerboseStepFlagOff(); - sliceRelPosFilter->WriteStepFlagOff(); - sliceRelPosFilter->SetInput(m_Working_Support); - sliceRelPosFilter->SetInputObject(SubclavianArtery); - sliceRelPosFilter->SetDirection(2); - sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "SubclavianArteryLeft")); - sliceRelPosFilter->AddOrientationTypeString("NotLeftTo"); - sliceRelPosFilter->IntermediateSpacingFlagOn(); - sliceRelPosFilter->SetIntermediateSpacing(2); - sliceRelPosFilter->UniqueConnectedComponentBySliceOff(); - sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff(); - - sliceRelPosFilter->CCLSelectionFlagOn(); // select one CCL by slice - sliceRelPosFilter->SetCCLSelectionDimension(0); // select according to X (0) axis - sliceRelPosFilter->SetCCLSelectionDirection(+1); // select most at Left - sliceRelPosFilter->CCLSelectionIgnoreSingleCCLFlagOff(); // do not ignore if only one CCL - - sliceRelPosFilter->AutoCropFlagOn(); - sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); - sliceRelPosFilter->RemoveObjectFlagOff(); - sliceRelPosFilter->Update(); - m_Working_Support = sliceRelPosFilter->GetOutput(); - - // END - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_2RL_Remove_Structures() -{ - Remove_Structures("2RL", "BrachioCephalicVein"); - Remove_Structures("2RL", "CommonCarotidArtery"); - Remove_Structures("2RL", "SubclavianArtery"); - Remove_Structures("2RL", "Thyroid"); - Remove_Structures("2RL", "Aorta"); - - // END - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support; + // m_Working_Support must be set + Remove_Structures(s, "BrachioCephalicVein"); + Remove_Structures(s, "BrachioCephalicArtery"); + Remove_Structures(s, "CommonCarotidArteryLeft"); + Remove_Structures(s, "CommonCarotidArteryRight"); + Remove_Structures(s, "SubclavianArteryLeft"); + Remove_Structures(s, "SubclavianArteryRight"); + Remove_Structures(s, "Thyroid"); + Remove_Structures(s, "Aorta"); } //-------------------------------------------------------------------- -//-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_2RL_SeparateRL() -{ - // --------------------------------------------------------------------------- - StartNewStep("[Station 2RL] Separate 2R/2L according to Trachea"); - - /*Rod says: - - "For station 2 there is a shift, dividing 2R from 2L, from midline - to the left paratracheal border." - - Algo: - - Consider Trachea SliceBySlice - - find extrema at Left - - add margins towards Right - - remove what is at Left of this line - */ - - // Get Trachea - MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); - - // Resize like the current support (to have the same number of slices) - Trachea = clitk::ResizeImageLike(Trachea, m_Working_Support, GetBackgroundValue()); - - // Find extrema post positions - std::vector tracheaLeftPositionsA; - std::vector tracheaLeftPositionsB; - clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(Trachea, - GetBackgroundValue(), 2, - 0, false, // Left - 1, // Vertical line - 1, // margins - tracheaLeftPositionsA, - tracheaLeftPositionsB); - // Copy support for R and L - typedef itk::ImageDuplicator DuplicatorType; - DuplicatorType::Pointer duplicator = DuplicatorType::New(); - duplicator->SetInputImage(m_Working_Support); - duplicator->Update(); - MaskImageType::Pointer m_Working_Support2 = duplicator->GetOutput(); - - // Cut post to this line for Right part - clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, - tracheaLeftPositionsA, - tracheaLeftPositionsB, - GetBackgroundValue(), 0, -10); - writeImage(m_Working_Support, "R.mhd"); - - // Cut post to this line for Left part - clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support2, - tracheaLeftPositionsA, - tracheaLeftPositionsB, - GetBackgroundValue(), 0, +10); - writeImage(m_Working_Support2, "L.mhd"); - - // END - StopCurrentStep(m_Working_Support); - m_ListOfStations["2R"] = m_Working_Support; - m_ListOfStations["2L"] = m_Working_Support2; -} -//-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractLymphStation_3A.txx b/segmentation/clitkExtractLymphStation_3A.txx index a55071c..ba2fa62 100644 --- a/segmentation/clitkExtractLymphStation_3A.txx +++ b/segmentation/clitkExtractLymphStation_3A.txx @@ -5,38 +5,155 @@ void clitk::ExtractLymphStationsFilter:: ExtractStation_3A_SetDefaultValues() { - SetFuzzyThreshold("3A", "Sternum", 0.5); - SetFuzzyThreshold("3A", "SubclavianArtery", 0.5); } //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_3A() +{ + if (!CheckForStation("3A")) return; + + StartNewStep("Station 3A"); + StartSubStep(); + + // Get the current support + StartNewStep("[Station 3A] Get the current 3A suppport"); + m_Working_Support = m_ListOfSupports["S3A"]; + m_ListOfStations["3A"] = m_Working_Support; + StopCurrentStep(m_Working_Support); + + ExtractStation_3A_Post_Left_Limits_With_Aorta_S5_Support(); + ExtractStation_3A_Post_Limits_With_Dilated_Aorta_S6_Support(); + ExtractStation_3A_AntPost_Superiorly(); + ExtractStation_3A_Remove_Structures(); + + // Generic RelativePosition processes + m_ListOfStations["3A"] = this->ApplyRelativePositionList("Station_3A", m_ListOfStations["3A"]); + + // Keep a single CCL + m_ListOfStations["3A"] = + clitk::SliceBySliceKeepMainCCL(m_ListOfStations["3A"], + GetBackgroundValue(), + GetForegroundValue()); + + // Store image filenames into AFDB + writeImage(m_ListOfStations["3A"], "seg/Station3A.mhd"); + GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd"); + WriteAFDB(); + StopSubStep(); +} +//-------------------------------------------------------------------- + + + //-------------------------------------------------------------------- template void clitk::ExtractLymphStationsFilter:: -ExtractStation_3A_SI_Limits() +ExtractStation_3A_Post_Left_Limits_With_Aorta_S5_Support() { - // Apex of the chest or Sternum & Carina. - StartNewStep("[Station 3A] Inf/Sup limits with Sternum and Carina"); + StartNewStep("[Station 3A] Post limits in S5 support according to Aorta"); - // Get Carina position (has been determined in Station8) - m_CarinaZ = GetAFDB()->GetDouble("CarinaZ"); + // Consider Aorta, remove Left/Post part ; only around S5 + // Get S5 support and Aorta + MaskImagePointer S5 = m_ListOfSupports["S5"]; + MaskImagePointer Aorta = GetAFDB()->template GetImage ("Aorta"); + Aorta = clitk::ResizeImageLike(Aorta, S5, GetBackgroundValue()); + + // Inferiorly, Aorta has two CCL that merge into a single one when + // S6 appears. Loop on Aorta slices, select the most ant one, detect + // the most ant point. + std::vector slices; + clitk::ExtractSlices(Aorta, 2, slices); + std::vector points; + for(uint i=0; i(slices[i], GetBackgroundValue(), false, 1); + std::vector c; + clitk::ComputeCentroids(slices[i], GetBackgroundValue(), c); + assert(c.size() == 3); // only 2 CCL + typename MaskSliceType::PixelType l; + if (c[1][1] > c[2][1]) { // We will remove the label=1 + l = 1; + } + else { + l = 2;// We will remove the label=2 + } + slices[i] = clitk::SetBackground(slices[i], slices[i], l, + GetBackgroundValue(), true); + + // Detect the most ant point + MaskSlicePointType p; + MaskImagePointType pA; + clitk::FindExtremaPointInAGivenDirection(slices[i], GetBackgroundValue(), 1, true, p); + // Set the X coordinate to the X coordinate of the centroid + if (l==1) p[0] = c[2][0]; + else p[0] = c[1][0]; + + // Convert in 3D and store + clitk::PointsUtils::Convert2DTo3D(p, Aorta, i, pA); + points.push_back(pA); + } - // Get Sternum and search for the upper position - MaskImagePointer Sternum = GetAFDB()->template GetImage("Sternum"); + // DEBUG + // MaskImagePointer o = clitk::JoinSlices(slices, Aorta, 2); + // writeImage(o, "o.mhd"); + // clitk::WriteListOfLandmarks(points, "Ant-Aorta.txt"); + + // Remove Post/Left to this point + m_Working_Support = + clitk::SliceBySliceSetBackgroundFromPoints(m_Working_Support, + GetBackgroundValue(), 2, + points, + true, // Set BG if X greater than point[x], and + true); // if Y greater than point[y] + + StopCurrentStep(m_Working_Support); + m_ListOfStations["3A"] = m_Working_Support; +} +//-------------------------------------------------------------------- - // Search most sup point - MaskImagePointType ps = Sternum->GetOrigin(); // initialise to avoid warning - clitk::FindExtremaPointInAGivenDirection(Sternum, GetBackgroundValue(), 2, false, ps); - double m_SternumZ = ps[2]+Sternum->GetSpacing()[2]; // One more slice, because it is below this point - //* Crop support : +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_3A_Post_Limits_With_Dilated_Aorta_S6_Support() +{ + StartNewStep("[Station 3A] Post limits with dilated Aorta"); + + // Consider Aorta + MaskImagePointer Aorta = GetAFDB()->template GetImage ("Aorta"); + + // Limits the support to S6 + MaskImagePointer S6 = m_ListOfSupports["S6"]; + Aorta = clitk::ResizeImageLike(Aorta, S6, GetBackgroundValue()); + + // Extend 1cm anteriorly + MaskImagePointType radius; // in mm + radius[0] = 10; + radius[1] = 10; + radius[2] = 0; // required + Aorta = clitk::Dilate(Aorta, radius, GetBackgroundValue(), GetForegroundValue(), false); + + // Now, insert this image in the AFDB (but do not store on disk) + GetAFDB()->template SetImage("Aorta_Dilated_Anteriorly", "bidon", + Aorta, false); + /* + // Not Post to m_Working_Support = - clitk::CropImageAlongOneAxis(m_Working_Support, 2, - m_CarinaZ, m_SternumZ, true, - GetBackgroundValue()); + clitk::SliceBySliceRelativePosition(m_Working_Support, Aorta, 2, + GetFuzzyThreshold("3A", "Aorta"), + "NotPostTo", true, + Aorta->GetSpacing()[0], false, false); + */ + + StopCurrentStep(m_Working_Support); m_ListOfStations["3A"] = m_Working_Support; } @@ -47,16 +164,29 @@ ExtractStation_3A_SI_Limits() template void clitk::ExtractLymphStationsFilter:: -ExtractStation_3A_Ant_Limits() +ExtractStation_3A_AntPost_Superiorly() { - StartNewStep("[Station 3A] Ant limits with Sternum"); + StartNewStep("[Station 3A] Post limits superiorly"); + + // Get or compute the binary mask that separate Ant/Post part + // according to vessels + MaskImagePointer binarizedContour = FindAntPostVessels2(); + binarizedContour = clitk::ResizeImageLike(binarizedContour, + m_Working_Support, + GetBackgroundValue()); - // Get Sternum, keep posterior part. - MaskImagePointer Sternum = GetAFDB()->template GetImage("Sternum"); - m_Working_Support = - clitk::SliceBySliceRelativePosition(m_Working_Support, Sternum, 2, - GetFuzzyThreshold("3A", "Sternum"), "PostTo", - false, 3, true, false); + // remove from support + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(m_Working_Support); + boolFilter->SetInput2(binarizedContour); + boolFilter->SetBackgroundValue1(GetBackgroundValue()); + boolFilter->SetBackgroundValue2(GetBackgroundValue()); + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + m_Working_Support = boolFilter->GetOutput(); + StopCurrentStep(m_Working_Support); m_ListOfStations["3A"] = m_Working_Support; } @@ -67,16 +197,64 @@ ExtractStation_3A_Ant_Limits() template void clitk::ExtractLymphStationsFilter:: -ExtractStation_3A_Post_Limits() +ExtractStation_3A_Remove_Structures() { - StartNewStep("[Station 3A] Post limits with SubclavianArtery"); + Remove_Structures(" 3A", "Aorta"); + Remove_Structures(" 3A", "SubclavianArteryLeft"); + Remove_Structures(" 3A", "SubclavianArteryRight"); + Remove_Structures(" 3A", "Thyroid"); + Remove_Structures(" 3A", "CommonCarotidArteryLeft"); + Remove_Structures(" 3A", "CommonCarotidArteryRight"); + Remove_Structures(" 3A", "BrachioCephalicArtery"); + // Remove_Structures("3A", "Bones"); --> should be in extractmediastinum + // Remove_Structures("3A", "BrachioCephalicVein"); ? + + StartNewStep("[Station 3A] Remove part of BrachioCephalicVein"); + // resize like support, extract slices + // while single CCL -> remove + // when two remove only the most post + MaskImagePointer BrachioCephalicVein = + GetAFDB()->template GetImage ("BrachioCephalicVein"); + BrachioCephalicVein = clitk::ResizeImageLike(BrachioCephalicVein, + m_Working_Support, + GetBackgroundValue()); + std::vector slices; + std::vector slices_BCV; + clitk::ExtractSlices(m_Working_Support, 2, slices); + clitk::ExtractSlices(BrachioCephalicVein, 2, slices_BCV); + for(uint i=0; i(slices_BCV[i], 0, true, 1); + + // Compute centroids + std::vector centroids; + ComputeCentroids(slices_BCV[i], GetBackgroundValue(), centroids); + + // If several centroid, select the one most anterior + if (centroids.size() > 2) { + // Only keep the one most post + typename MaskSliceType::PixelType label; + if (centroids[1][1] > centroids[2][1]) { + label = 2; + } + else { + label = 1; + } + // "remove" the CCL + slices_BCV[i] = clitk::SetBackground(slices_BCV[i], + slices_BCV[i], + label, + GetBackgroundValue(), + true); + } + + // Remove from the support + clitk::AndNot(slices[i], slices_BCV[i], GetBackgroundValue()); + } + + // Joint + m_Working_Support = clitk::JoinSlices(slices, m_Working_Support, 2); - // Get Sternum, keep posterior part. - MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage("SubclavianArtery"); - m_Working_Support = - clitk::SliceBySliceRelativePosition(m_Working_Support, SubclavianArtery, 2, - GetFuzzyThreshold("3A", "SubclavianArtery"), "AntTo", - false, 3, true, false); StopCurrentStep(m_Working_Support); m_ListOfStations["3A"] = m_Working_Support; } diff --git a/segmentation/clitkExtractLymphStation_3P.txx b/segmentation/clitkExtractLymphStation_3P.txx index 1c9fdc0..7120498 100644 --- a/segmentation/clitkExtractLymphStation_3P.txx +++ b/segmentation/clitkExtractLymphStation_3P.txx @@ -10,40 +10,41 @@ ExtractStation_3P_SetDefaultValues() //-------------------------------------------------------------------- -template +template void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3P_SI_Limits() +clitk::ExtractLymphStationsFilter:: +ExtractStation_3P() { - /* - Apex of the chest & Carina. - */ - StartNewStep("[Station 3P] Inf/Sup limits with apex of the chest and carina"); + if (!CheckForStation("3P")) return; - // Get Carina position (has been determined in Station8) - m_CarinaZ = GetAFDB()->GetDouble("CarinaZ"); - - // Get Apex of the Chest. The "lungs" structure is autocroped so we - // consider the most superior point. - MaskImagePointer Lungs = GetAFDB()->template GetImage("Lungs"); - MaskImageIndexType index = Lungs->GetLargestPossibleRegion().GetIndex(); - index += Lungs->GetLargestPossibleRegion().GetSize(); - MaskImagePointType p; - Lungs->TransformIndexToPhysicalPoint(index, p); - p[2] -= 5; // 5 mm below because the autocrop is slightly above the apex - double m_ApexOfTheChest = p[2]; - - /* Crop support : - Superior limit = carina - Inferior limit = Apex of the chest */ - m_Working_Support = - clitk::CropImageAlongOneAxis(m_Working_Support, 2, - m_CarinaZ, - m_ApexOfTheChest, true, - GetBackgroundValue()); + StartNewStep("Station 3P"); + StartSubStep(); + // Get the current support + StartNewStep("[Station 3P] Get the current 3P suppport"); + m_Working_Support = m_ListOfSupports["S3P"]; + m_ListOfStations["3P"] = m_Working_Support; StopCurrentStep(m_Working_Support); + + // ExtractStation_3P_LR_inf_Limits(); // <-- done in RelPosList + + // Generic RelativePosition processes + m_ListOfStations["3P"] = this->ApplyRelativePositionList("Station_3P", m_ListOfStations["3P"]); + m_Working_Support = m_ListOfStations["3P"]; + ExtractStation_8_Single_CCL_Limits(); // YES 8 ! + ExtractStation_3P_Remove_Structures(); // after CCL m_ListOfStations["3P"] = m_Working_Support; + + // Old stuff + // LR limits superiorly => not here for the moment because not clear in the def + // ExtractStation_3P_LR_sup_Limits_2(); //TODO + // ExtractStation_3P_LR_sup_Limits(); // old version to change + + // Store image filenames into AFDB + writeImage(m_ListOfStations["3P"], "seg/Station3P.mhd"); + GetAFDB()->SetImageFilename("Station3P", "seg/Station3P.mhd"); + WriteAFDB(); + StopSubStep(); } //-------------------------------------------------------------------- @@ -54,172 +55,10 @@ void clitk::ExtractLymphStationsFilter:: ExtractStation_3P_Remove_Structures() { - /* - 3 - (sup) remove Aorta, VB, subcl - not LR subcl ! -> a séparer LR ? - (inf) remove Eso, Aorta, Azygousvein - */ - StartNewStep("[Station 3P] Remove some structures."); - - Remove_Structures("3P", "Aorta"); - Remove_Structures("3P", "VertebralBody"); - Remove_Structures("3P", "SubclavianArtery"); Remove_Structures("3P", "Esophagus"); - Remove_Structures("3P", "Azygousvein"); Remove_Structures("3P", "Thyroid"); - Remove_Structures("3P", "VertebralArtery"); - - StopCurrentStep(m_Working_Support); - m_ListOfStations["3P"] = m_Working_Support; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3P_Ant_Limits() -{ - /* - Ant Post limit : - - Anteriorly, it is in contact with the posterior aspect of Stations - 1–2 superiorly (Fig. 2A–C) and with Stations 4R and 4L inferiorly - (Fig. 2D–I and 3A–C). The anterior limit of Station 3P is kept - posterior to the trachea, which is defined by an imaginary - horizontal line running along the posterior wall of the trachea - (Fig. 2B,E, red line). Posteriorly, it is delineated along the - anterior and lateral borders of the vertebral body until an - imaginary horizontal line running 1 cm posteriorly to the - anterior border of the vertebral body (Fig. 2D). - - 1 - post to the trachea : define an imaginary line just post the - Trachea and remove what is anterior. - - */ - StartNewStep("[Station 3P] Ant limits with Trachea "); - - // Load Trachea - MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); - - // Crop like current support (need by SliceBySliceSetBackgroundFromLineSeparation after) - Trachea = - clitk::ResizeImageLike(Trachea, m_Working_Support, GetBackgroundValue()); - - // Slice by slice, determine the most post point of the trachea (A) - // and consider a second point on the same horizontal line (B) - std::vector p_A; - std::vector p_B; - std::vector slices; - clitk::ExtractSlices(Trachea, 2, slices); - MaskImagePointType p; - typename MaskSliceType::PointType sp; - for(uint i=0; i(slices[i], 0, true, 100); - slices[i] = KeepLabels(slices[i], GetBackgroundValue(), - GetForegroundValue(), 1, 1, true); - // Find most posterior point - clitk::FindExtremaPointInAGivenDirection(slices[i], GetBackgroundValue(), - 1, false, sp); - clitk::PointsUtils::Convert2DTo3D(sp, Trachea, i, p); - p_A.push_back(p); - p[0] -= 20; - p_B.push_back(p); - } - clitk::WriteListOfLandmarks(p_A, "trachea-post.txt"); - - // Remove Ant part above those lines - clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, - p_A, p_B, - GetBackgroundValue(), - 1, 10); - // End, release memory - GetAFDB()->template ReleaseImage("Trachea"); - StopCurrentStep(m_Working_Support); - m_ListOfStations["3P"] = m_Working_Support; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3P_Post_Limits() -{ - /* - Ant Post limit : - - Anteriorly, it is in contact with the posterior aspect of Stations - 1–2 superiorly (Fig. 2A–C) and with Stations 4R and 4L inferiorly - (Fig. 2D–I and 3A–C). The anterior limit of Station 3P is kept - posterior to the trachea, which is defined by an imaginary - horizontal line running along the posterior wall of the trachea - (Fig. 2B,E, red line). Posteriorly, it is delineated along the - anterior and lateral borders of the vertebral body until an - imaginary horizontal line running 1 cm posteriorly to the - anterior border of the vertebral body (Fig. 2D). - - 2 - post to the trachea : define an imaginary line just post the - Trachea and remove what is anterior. - - */ - StartNewStep("[Station 3P] Post limits with VertebralBody "); - - // Load VertebralBody - MaskImagePointer VertebralBody = GetAFDB()->template GetImage("VertebralBody"); - - // Crop like current support (need by SliceBySliceSetBackgroundFromLineSeparation after) - VertebralBody = clitk::ResizeImageLike(VertebralBody, m_Working_Support, GetBackgroundValue()); - - writeImage(VertebralBody, "vb.mhd"); - - // Compute VertebralBody most Ant position (again because slices - // changes). Slice by slice, determine the most post point of the - // trachea (A) and consider a second point on the same horizontal - // line (B) - std::vector p_A; - std::vector p_B; - std::vector slices; - clitk::ExtractSlices(VertebralBody, 2, slices); - MaskImagePointType p; - typename MaskSliceType::PointType sp; - for(uint i=0; i(slices[i], GetBackgroundValue(), - 1, true, sp); - - // If the VertebralBody stop superiorly before the end of - // m_Working_Support, we consider the same previous point. - if (!found) { - p = p_A.back(); - p[2] += VertebralBody->GetSpacing()[2]; - p_A.push_back(p); - p = p_B.back(); - p[2] += VertebralBody->GetSpacing()[2]; - p_B.push_back(p); - } - else { - clitk::PointsUtils::Convert2DTo3D(sp, VertebralBody, i, p); - p[1] += 10; // Consider 10 mm more post - p_A.push_back(p); - p[0] -= 20; - p_B.push_back(p); - } - } - clitk::WriteListOfLandmarks(p_A, "vb-ant.txt"); - - - // Remove Ant part above those lines - clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, - p_A, p_B, - GetBackgroundValue(), - 1, -10); - writeImage(m_Working_Support, "a.mhd"); - + Remove_Structures("3P", "VertebralArtery"); // (inside the station) StopCurrentStep(m_Working_Support); m_ListOfStations["3P"] = m_Working_Support; } @@ -327,7 +166,7 @@ ExtractStation_3P_LR_sup_Limits() relPosFilter->SetBackgroundValue(GetBackgroundValue()); relPosFilter->SetInput(slices_support[i]); relPosFilter->SetInputObject(object); - relPosFilter->AddOrientationTypeString("R"); + relPosFilter->AddOrientationTypeString("RightTo"); relPosFilter->SetInverseOrientationFlag(true); // relPosFilter->SetIntermediateSpacing(3); relPosFilter->SetIntermediateSpacingFlag(false); @@ -343,7 +182,7 @@ ExtractStation_3P_LR_sup_Limits() relPosFilter->SetBackgroundValue(GetBackgroundValue()); relPosFilter->SetInput(slices_support[i]); relPosFilter->SetInputObject(object); - relPosFilter->AddOrientationTypeString("A"); + relPosFilter->AddOrientationTypeString("AntTo"); relPosFilter->SetInverseOrientationFlag(true); // relPosFilter->SetIntermediateSpacing(3); relPosFilter->SetIntermediateSpacingFlag(false); @@ -422,29 +261,6 @@ template void clitk::ExtractLymphStationsFilter:: ExtractStation_3P_LR_sup_Limits_2() -{ - /* - Use VertebralArtery to limit. - - - StartNewStep("[Station 3P] Left/Right limits with VertebralArtery"); - - // Load structures - MaskImagePointer VertebralArtery = GetAFDB()->template GetImage("VertebralArtery"); - - clitk::AndNot(m_Working_Support, VertebralArtery); - - StopCurrentStep(m_Working_Support); - m_ListOfStations["3P"] = m_Working_Support; - */ -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3P_LR_inf_Limits() { /* "On the right side, the limit is defined by the air–soft-tissue @@ -452,44 +268,16 @@ ExtractStation_3P_LR_inf_Limits() interface superiorly (Fig. 2A–C) and the aorta inferiorly (Figs. 2D–I and 3A–C)." - inf : not Right to Azygousvein + sup : + Resizelike support : Trachea, SubclavianArtery + Trachea, slice by slice, get centroid trachea + SubclavianArtery, slice by slice, CCL + prendre la 1ère à L et R, not at Left + */ - StartNewStep("[Station 3P] Left/Right limits (inferior part) with Azygousvein and Aorta"); - - // Load structures - MaskImagePointer AzygousVein = GetAFDB()->template GetImage("AzygousVein"); - MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); - - typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; - typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); - relPosFilter->VerboseStepFlagOff(); - relPosFilter->WriteStepFlagOff(); - relPosFilter->SetBackgroundValue(GetBackgroundValue()); - relPosFilter->SetInput(m_Working_Support); - relPosFilter->SetInputObject(AzygousVein); - relPosFilter->AddOrientationTypeString("R"); - relPosFilter->SetInverseOrientationFlag(true); - // relPosFilter->SetIntermediateSpacing(3); - relPosFilter->SetIntermediateSpacingFlag(false); - relPosFilter->SetFuzzyThreshold(0.7); - relPosFilter->AutoCropFlagOn(); - relPosFilter->Update(); - m_Working_Support = relPosFilter->GetOutput(); - - writeImage(m_Working_Support, "before-L-aorta.mhd"); - relPosFilter->SetInput(m_Working_Support); - relPosFilter->SetInputObject(Aorta); - relPosFilter->AddOrientationTypeString("L"); - relPosFilter->SetInverseOrientationFlag(true); - // relPosFilter->SetIntermediateSpacing(3); - relPosFilter->SetIntermediateSpacingFlag(false); - relPosFilter->SetFuzzyThreshold(0.7); - relPosFilter->AutoCropFlagOn(); - relPosFilter->Update(); - m_Working_Support = relPosFilter->GetOutput(); - writeImage(m_Working_Support, "after-L-aorta.mhd"); + // StartNewStep("[Station 3P] Left/Right limits (superior part) "); + - StopCurrentStep(m_Working_Support); - m_ListOfStations["3P"] = m_Working_Support; } //-------------------------------------------------------------------- + diff --git a/segmentation/clitkExtractLymphStation_7.txx b/segmentation/clitkExtractLymphStation_7.txx index b12fa52..d6df428 100644 --- a/segmentation/clitkExtractLymphStation_7.txx +++ b/segmentation/clitkExtractLymphStation_7.txx @@ -65,7 +65,7 @@ ExtractStation_7_SI_Limits() A, B, 2, 0, false); // Get the CarinaZ position - m_CarinaZ = FindCarinaSlicePosition(); + double m_CarinaZ = FindCarina(); // Crop support m_Working_Support = @@ -354,7 +354,7 @@ ExtractStation_7_RL_Limits_OLD() sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); sliceRelPosFilter->SetIntermediateSpacingFlag(true); sliceRelPosFilter->SetIntermediateSpacing(3); - sliceRelPosFilter->SetUniqueConnectedComponentBySlice(false); + sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(false); sliceRelPosFilter->SetAutoCropFlag(false); sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); sliceRelPosFilter->Update(); @@ -380,7 +380,7 @@ ExtractStation_7_RL_Limits_OLD() sliceRelPosFilter->AddOrientationTypeString("NotPostTo"); sliceRelPosFilter->SetIntermediateSpacingFlag(true); sliceRelPosFilter->SetIntermediateSpacing(3); - sliceRelPosFilter->SetUniqueConnectedComponentBySlice(false); + sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(false); sliceRelPosFilter->SetAutoCropFlag(false); sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); sliceRelPosFilter->Update(); @@ -403,7 +403,7 @@ ExtractStation_7_RL_Limits_OLD() sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); sliceRelPosFilter->SetIntermediateSpacingFlag(true); sliceRelPosFilter->SetIntermediateSpacing(3); - sliceRelPosFilter->SetUniqueConnectedComponentBySlice(false); + sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(false); sliceRelPosFilter->SetAutoCropFlag(false); sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); sliceRelPosFilter->Update(); @@ -421,7 +421,7 @@ ExtractStation_7_RL_Limits_OLD() sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); sliceRelPosFilter->SetIntermediateSpacingFlag(true); sliceRelPosFilter->SetIntermediateSpacing(3); - sliceRelPosFilter->SetUniqueConnectedComponentBySlice(false); + sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(false); sliceRelPosFilter->SetAutoCropFlag(false); sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); sliceRelPosFilter->Update(); @@ -439,7 +439,7 @@ ExtractStation_7_RL_Limits_OLD() sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); sliceRelPosFilter->SetIntermediateSpacingFlag(true); sliceRelPosFilter->SetIntermediateSpacing(3); - sliceRelPosFilter->SetUniqueConnectedComponentBySlice(false); + sliceRelPosFilter->SetUniqueConnectedComponentBySliceFlag(false); sliceRelPosFilter->SetAutoCropFlag(true); sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); sliceRelPosFilter->Update(); diff --git a/segmentation/clitkExtractLymphStation_8.txx b/segmentation/clitkExtractLymphStation_8.txx index 7402a62..d4bbfc3 100644 --- a/segmentation/clitkExtractLymphStation_8.txx +++ b/segmentation/clitkExtractLymphStation_8.txx @@ -24,19 +24,22 @@ void clitk::ExtractLymphStationsFilter:: ExtractStation_8() { - if (CheckForStation("8")) { - ExtractStation_8_SI_Limits(); // OK, validated - ExtractStation_8_Ant_Limits(); // OK, validated - ExtractStation_8_Left_Sup_Limits(); // OK, validated - ExtractStation_8_Left_Inf_Limits(); // OK, validated - ExtractStation_8_Single_CCL_Limits(); // OK, validated - ExtractStation_8_Remove_Structures(); // OK, validated - - // Store image filenames into AFDB - writeImage(m_ListOfStations["8"], "seg/Station8.mhd"); - GetAFDB()->SetImageFilename("Station8", "seg/Station8.mhd"); - WriteAFDB(); - } + if (!CheckForStation("8")) return; + + StartNewStep("Station 8"); + StartSubStep(); + ExtractStation_8_SI_Limits(); // OK, validated + ExtractStation_8_Ant_Limits(); // OK, validated + ExtractStation_8_Left_Sup_Limits(); // OK, validated + ExtractStation_8_Left_Inf_Limits(); // OK, validated + ExtractStation_8_Single_CCL_Limits(); // OK, validated + ExtractStation_8_Remove_Structures(); // OK, validated + + // Store image filenames into AFDB + writeImage(m_ListOfStations["8"], "seg/Station8.mhd"); + GetAFDB()->SetImageFilename("Station8", "seg/Station8.mhd"); + WriteAFDB(); + StopSubStep(); } //-------------------------------------------------------------------- @@ -154,7 +157,7 @@ ExtractStation_8_Ant_Limits() relPosFilter->AddOrientationTypeString("PostTo"); // relPosFilter->InverseOrientationFlagOff(); relPosFilter->SetDirection(2); // Z axis - relPosFilter->UniqueConnectedComponentBySliceOff(); + relPosFilter->UniqueConnectedComponentBySliceFlagOff(); relPosFilter->SetIntermediateSpacing(3); relPosFilter->IntermediateSpacingFlagOn(); relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("8", "Esophagus")); diff --git a/segmentation/clitkExtractLymphStation_Supports.txx b/segmentation/clitkExtractLymphStation_Supports.txx index 903948d..ad2e957 100644 --- a/segmentation/clitkExtractLymphStation_Supports.txx +++ b/segmentation/clitkExtractLymphStation_Supports.txx @@ -11,64 +11,574 @@ ExtractStationSupports() DD("ExtractStationSupports"); // Get initial Mediastinum - m_Working_Support = m_Mediastinum = GetAFDB()->template GetImage("Mediastinum", true); + m_Working_Support = m_Mediastinum = this->GetAFDB()->template GetImage("Mediastinum", true); - // Superior limits: CricoidCartilag - // Inferior limits: lung - // (the Mediastinum support already stop at this limit) - - // Consider above Carina - m_CarinaZ = FindCarinaSlicePosition(); + // Consider sup/inf to Carina + double m_CarinaZ = FindCarina(); MaskImagePointer m_Support_Superior_to_Carina = clitk::CropImageRemoveLowerThan(m_Working_Support, 2, m_CarinaZ, true, GetBackgroundValue()); MaskImagePointer m_Support_Inferior_to_Carina = clitk::CropImageRemoveGreaterThan(m_Working_Support, 2, m_CarinaZ, true, GetBackgroundValue()); + m_ListOfSupports["Support_Superior_to_Carina"] = m_Support_Superior_to_Carina; + m_ListOfSupports["Support_Inferior_to_Carina"] = m_Support_Inferior_to_Carina; + writeImage(m_Support_Inferior_to_Carina, "seg/Support_Inf_Carina.mhd"); + this->GetAFDB()->SetImageFilename("Support_Inf_Carina", "seg/Support_Inf_Carina.mhd"); + writeImage(m_Support_Superior_to_Carina, "seg/Support_Sup_Carina.mhd"); + this->GetAFDB()->SetImageFilename("Support_Sup_Carina", "seg/Support_Sup_Carina.mhd"); + + // S1RL + Support_SupInf_S1RL(); + Support_LeftRight_S1R_S1L(); + + // S2RL + Support_SupInf_S2R_S2L(); + Support_LeftRight_S2R_S2L(); + + // S4RL + Support_SupInf_S4R_S4L(); + Support_LeftRight_S4R_S4L(); + + // Post limits of S1,S2,S4 + Support_Post_S1S2S4(); + + // S3AP + Support_S3P(); + Support_S3A(); + + // S5, S6 + Support_S5(); + Support_S6(); + + // Below Carina S7,8,9,10 + m_ListOfSupports["S7"] = clitk::Clone(m_Support_Inferior_to_Carina); + m_ListOfSupports["S8"] = clitk::Clone(m_Support_Inferior_to_Carina); + m_ListOfSupports["S9"] = clitk::Clone(m_Support_Inferior_to_Carina); + m_ListOfSupports["S10"] = clitk::Clone(m_Support_Inferior_to_Carina); + m_ListOfSupports["S11"] = clitk::Clone(m_Support_Inferior_to_Carina); + + // Store image filenames into AFDB + writeImage(m_ListOfSupports["S1R"], "seg/Support_S1R.mhd"); + this->GetAFDB()->SetImageFilename("Support_S1R", "seg/Support_S1R.mhd"); + writeImage(m_ListOfSupports["S1L"], "seg/Support_S1L.mhd"); + this->GetAFDB()->SetImageFilename("Support_S1L", "seg/Support_S1L.mhd"); + + writeImage(m_ListOfSupports["S2L"], "seg/Support_S2L.mhd"); + this->GetAFDB()->SetImageFilename("Support_S2L", "seg/Support_S2L.mhd"); + writeImage(m_ListOfSupports["S2R"], "seg/Support_S2R.mhd"); + this->GetAFDB()->SetImageFilename("Support_S2R", "seg/Support_S2R.mhd"); + + writeImage(m_ListOfSupports["S3P"], "seg/Support_S3P.mhd"); + this->GetAFDB()->SetImageFilename("Support_S3P", "seg/Support_S3P.mhd"); + writeImage(m_ListOfSupports["S3A"], "seg/Support_S3A.mhd"); + this->GetAFDB()->SetImageFilename("Support_S3A", "seg/Support_S3A.mhd"); + + writeImage(m_ListOfSupports["S4L"], "seg/Support_S4L.mhd"); + this->GetAFDB()->SetImageFilename("Support_S4L", "seg/Support_S4L.mhd"); + writeImage(m_ListOfSupports["S4R"], "seg/Support_S4R.mhd"); + this->GetAFDB()->SetImageFilename("Support_S4R", "seg/Support_S4R.mhd"); + + writeImage(m_ListOfSupports["S5"], "seg/Support_S5.mhd"); + this->GetAFDB()->SetImageFilename("Support_S5", "seg/Support_S5.mhd"); + writeImage(m_ListOfSupports["S6"], "seg/Support_S6.mhd"); + this->GetAFDB()->SetImageFilename("Support_S6", "seg/Support_S6.mhd"); + + writeImage(m_ListOfSupports["S7"], "seg/Support_S7.mhd"); + this->GetAFDB()->SetImageFilename("Support_S7", "seg/Support_S7.mhd"); + + writeImage(m_ListOfSupports["S8"], "seg/Support_S8.mhd"); + this->GetAFDB()->SetImageFilename("Support_S8", "seg/Support_S8.mhd"); + + writeImage(m_ListOfSupports["S9"], "seg/Support_S9.mhd"); + this->GetAFDB()->SetImageFilename("Support_S9", "seg/Support_S9.mhd"); + + writeImage(m_ListOfSupports["S10"], "seg/Support_S10.mhd"); + this->GetAFDB()->SetImageFilename("Support_S10", "seg/Support_S10.mhd"); + + writeImage(m_ListOfSupports["S11"], "seg/Support_S11.mhd"); + this->GetAFDB()->SetImageFilename("Support_S11", "seg/Support_S11.mhd"); + WriteAFDB(); +} +//-------------------------------------------------------------------- - // Consider only Superior to Carina - m_Working_Support = m_Support_Superior_to_Carina; +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_SupInf_S1RL() +{ // Step : S1RL - StartNewStep("[Support] sup-inf S1RL"); + StartNewStep("[Support] Sup-Inf S1RL"); /* - Lower border: clavicles bilaterally and, in the midline, the upper - border of the manubrium, 1R designates right-sided nodes, 1L, - left-sided nodes in this region - 2R: Upper border: apex of the right lung and pleural space, and in the midline, the upper border of the manubrium 2L: Upper border: apex of the left lung and pleural space, and in the midline, the upper border of the manubrium + + => apex / manubrium = up Sternum + */ + m_Working_Support = m_ListOfSupports["Support_Superior_to_Carina"]; + MaskImagePointer Sternum = this->GetAFDB()->template GetImage ("Sternum"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(Sternum, GetBackgroundValue(), 2, false, p); + MaskImagePointer S1RL = + clitk::CropImageRemoveLowerThan(m_Working_Support, 2, + p[2], true, GetBackgroundValue()); + m_Working_Support = + clitk::CropImageRemoveGreaterThan(m_Working_Support, 2, + p[2], true, GetBackgroundValue()); + m_ListOfSupports["S1RL"] = S1RL; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_LeftRight_S1R_S1L() +{ + // Step S1RL : Left-Right + StartNewStep("[Support] Left-Right S1R S1L"); + std::vector A; + std::vector B; + // Search for centroid positions of trachea + MaskImagePointer Trachea = this->GetAFDB()->template GetImage ("Trachea"); + MaskImagePointer S1RL = m_ListOfSupports["S1RL"]; + Trachea = clitk::ResizeImageLike(Trachea, S1RL, GetBackgroundValue()); + std::vector slices; + clitk::ExtractSlices(Trachea, 2, slices); + for(uint i=0; i(slices[i], 0, false, 10); + std::vector c; + clitk::ComputeCentroids(slices[i], GetBackgroundValue(), c); + ImagePointType a,b; + clitk::PointsUtils::Convert2DTo3D(c[1], Trachea, i, a); + A.push_back(a); + b = a; + b[1] += 50; + B.push_back(b); + } + clitk::WriteListOfLandmarks(A, "S1LR_A.txt"); + clitk::WriteListOfLandmarks(B, "S1LR_B.txt"); + + // Clone support + MaskImagePointer S1R = clitk::Clone(S1RL); + MaskImagePointer S1L = clitk::Clone(S1RL); + + // Right part + clitk::SliceBySliceSetBackgroundFromLineSeparation(S1R, A, B, + GetBackgroundValue(), 0, -10); + S1R = clitk::AutoCrop(S1R, GetBackgroundValue()); + m_ListOfSupports["S1R"] = S1R; + + // Left part + clitk::SliceBySliceSetBackgroundFromLineSeparation(S1L, A, B, + GetBackgroundValue(), 0, 10); + S1L = clitk::AutoCrop(S1L, GetBackgroundValue()); + m_ListOfSupports["S1L"] = S1L; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_SupInf_S2R_S2L() +{ + // Step : S2RL Sup-Inf limits + /* + 2R Lower border: intersection of caudal margin of innominate vein with + the trachea + 2L Lower border: superior border of the aortic arch */ + StartNewStep("[Support] Sup-Inf S2RL"); + m_Working_Support = m_ListOfSupports["Support_Superior_to_Carina"]; + + // S2R Caudal Margin Of Left BrachiocephalicVein + MaskImagePointer BrachioCephalicVein = this->GetAFDB()->template GetImage("BrachioCephalicVein"); + MaskImagePointType p; + clitk::FindExtremaPointInAGivenDirection(BrachioCephalicVein, GetBackgroundValue(), 2, true, p); + // I add slightly more than a slice + double CaudalMarginOfLeftBrachiocephalicVeinZ=p[2]+ 1.1*m_Working_Support->GetSpacing()[2]; + this->GetAFDB()->SetDouble("CaudalMarginOfLeftBrachiocephalicVeinZ", CaudalMarginOfLeftBrachiocephalicVeinZ); + MaskImagePointer S2R = + clitk::CropImageRemoveLowerThan(m_Working_Support, 2, + CaudalMarginOfLeftBrachiocephalicVeinZ, true, + GetBackgroundValue()); + // S2L : Top Of Aortic Arch + MaskImagePointer Aorta = this->GetAFDB()->template GetImage("Aorta"); + clitk::FindExtremaPointInAGivenDirection(Aorta, GetBackgroundValue(), 2, false, p); + // Save the TopOfAorticArchZ + this->GetAFDB()->SetDouble("TopOfAorticArchZ", p[2]); + // I substract slightly more than a slice to respect delineation + double TopOfAorticArchZ=p[2]- 1.1*m_Working_Support->GetSpacing()[2]; + // this->GetAFDB()->SetDouble("TopOfAorticArchZ", TopOfAorticArchZ); + MaskImagePointer S2L = + clitk::CropImageRemoveLowerThan(m_Working_Support, 2, + TopOfAorticArchZ, true, + GetBackgroundValue()); + /* + // S2RL: Superior support, I use inferior part of S1RL + MaskImagePointer S1L = m_ListOfSupports["S1L"]; + clitk::FindExtremaPointInAGivenDirection(S1L, GetBackgroundValue(), 2, true, p); + DD(p); + S2L = + clitk::CropImageRemoveGreaterThan(S2L, 2, + p[2], true, + GetBackgroundValue()); + MaskImagePointer S1R = m_ListOfSupports["S1R"]; + clitk::FindExtremaPointInAGivenDirection(S1R, GetBackgroundValue(), 2, true, p); + DD(p); + S2R = + clitk::CropImageRemoveGreaterThan(S2R, 2, + p[2], true, + GetBackgroundValue()); + */ + // Superior limits, use Sternum (but not strictly inf to S1RL + MaskImagePointer Sternum = this->GetAFDB()->template GetImage ("Sternum"); + clitk::FindExtremaPointInAGivenDirection(Sternum, GetBackgroundValue(), 2, false, p); + // Add one slice + p[2] = p[2] + m_Working_Support->GetSpacing()[2]; + S2L = + clitk::CropImageRemoveGreaterThan(S2L, 2, + p[2], true, GetBackgroundValue()); + S2R = + clitk::CropImageRemoveGreaterThan(S2R, 2, + p[2], true, GetBackgroundValue()); + // The is the end + m_ListOfSupports["S2L"] = S2L; + m_ListOfSupports["S2R"] = S2R; +} +//-------------------------------------------------------------------- - // // LeftRight cut along trachea - // MaskImagePointer Trachea = GetAFDB()->GetImage("Trachea"); - // // build a ant-post line for each slice - // MaskImagePointer m_Support_SupRight = - // clitk::CropImageRemoveLowerThan(m_Working_Support, 2, - // m_CarinaZ, true, GetBackgroundValue()); +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_LeftRight_S2R_S2L() +{ + // --------------------------------------------------------------------------- + /* Step : S2RL LeftRight + As for lymph node station 4R, 2R includes nodes extending to the + left lateral border of the trachea + Rod says: "For station 2 there is a shift, dividing 2R from 2L, from midline + to the left paratracheal border." + */ + StartNewStep("[Support] Separate 2R/2L according to Trachea"); + MaskImagePointer S2R = m_ListOfSupports["S2R"]; + MaskImagePointer S2L = m_ListOfSupports["S2L"]; + S2R = LimitsWithTrachea(S2R, 0, 1, -10); + S2L = LimitsWithTrachea(S2L, 0, 1, 10); + m_ListOfSupports["S2R"] = S2R; + m_ListOfSupports["S2L"] = S2L; + this->GetAFDB()->template ReleaseImage("Trachea"); +} +//-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_SupInf_S4R_S4L() +{ + // --------------------------------------------------------------------------- + /* Step : S4RL Sup-Inf + - start at the end of 2R and 2L + - stop ? + - 4R + Rod says : "The inferior border is at the lower border of the azygous vein." + Rod says : difficulties + (was : "ends at the upper lobe bronchus or where the right pulmonary artery + crosses the midline of the mediastinum ") + - 4L + Rod says : "The lower border is to upper margin of the left main pulmonary artery." + (was LLL bronchus) + */ + StartNewStep("[Support] Sup-Inf limits of 4R/4L"); + // Start from the support, remove 2R and 2L + MaskImagePointer S4RL = clitk::Clone(m_Working_Support); + MaskImagePointer S2R = m_ListOfSupports["S2R"]; + MaskImagePointer S2L = m_ListOfSupports["S2L"]; + clitk::AndNot(S4RL, S2R, GetBackgroundValue()); + clitk::AndNot(S4RL, S2L, GetBackgroundValue()); + S4RL = clitk::AutoCrop(S4RL, GetBackgroundValue()); - // Store image filenames into AFDB - m_ListOfSupports["S1R"] = m_Working_Support; - writeImage(m_ListOfSupports["S1R"], "seg/Support_S1R.mhd"); - GetAFDB()->SetImageFilename("Support_S1R", "seg/Support_S1R.mhd"); - WriteAFDB(); + // Copy, stop 4R at AzygousVein and 4L at LeftPulmonaryArtery + MaskImagePointer S4R = clitk::Clone(S4RL); + MaskImagePointer S4L = clitk::Clone(S4RL); + + // Get AzygousVein and limit according to LowerBorderAzygousVein + MaskImagePointer LowerBorderAzygousVein + = this->GetAFDB()->template GetImage("LowerBorderAzygousVein"); + std::vector c; + clitk::ComputeCentroids(LowerBorderAzygousVein, GetBackgroundValue(), c); + S4R = clitk::CropImageRemoveLowerThan(S4RL, 2, + c[1][2], true, GetBackgroundValue()); + S4R = clitk::AutoCrop(S4R, GetBackgroundValue()); + m_ListOfSupports["S4R"] = S4R; + + + // Limit according to LeftPulmonaryArtery + MaskImagePointer LeftPulmonaryArtery + = this->GetAFDB()->template GetImage("LeftPulmonaryArtery"); + MaskImagePointType p; + clitk::FindExtremaPointInAGivenDirection(LeftPulmonaryArtery, GetBackgroundValue(), + 2, false, p); + S4L = clitk::CropImageRemoveLowerThan(S4RL, 2, + p[2], true, GetBackgroundValue()); + S4L = clitk::AutoCrop(S4L, GetBackgroundValue()); + m_ListOfSupports["S4L"] = S4L; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_LeftRight_S4R_S4L() +{ + // --------------------------------------------------------------------------- + /* Step : S4RL LeftRight + + - 4R: includes right paratracheal nodes, and pretracheal nodes + extending to the left lateral border of trachea + + - 4L: includes nodes to the left of the left lateral border of + the trachea, medial to the ligamentum arteriosum + + => same than 2RL + */ + StartNewStep("[Support] Left Right separation of 4R/4L"); + + MaskImagePointer S4R = m_ListOfSupports["S4R"]; + MaskImagePointer S4L = m_ListOfSupports["S4L"]; + S4R = LimitsWithTrachea(S4R, 0, 1, -10); + S4L = LimitsWithTrachea(S4L, 0, 1, 10); + m_ListOfSupports["S4R"] = S4R; + m_ListOfSupports["S4L"] = S4L; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +typename clitk::ExtractLymphStationsFilter::MaskImagePointer +clitk::ExtractLymphStationsFilter:: +LimitsWithTrachea(MaskImageType * input, int extremaDirection, int lineDirection, + double offset) +{ + MaskImagePointType min, max; + GetMinMaxBoundary(input, min, max); + return LimitsWithTrachea(input, extremaDirection, lineDirection, offset, max[2]); +} +template +typename clitk::ExtractLymphStationsFilter::MaskImagePointer +clitk::ExtractLymphStationsFilter:: +LimitsWithTrachea(MaskImageType * input, int extremaDirection, int lineDirection, + double offset, double maxSupPosition) +{ + /* + Take the input mask, consider the trachea and limit according to + Left border of the trachea. Keep at Left or at Right according to + the offset + */ + // Read the trachea + MaskImagePointer Trachea = this->GetAFDB()->template GetImage("Trachea"); + + // Find extrema post positions + std::vector tracheaLeftPositionsA; + std::vector tracheaLeftPositionsB; + + // Crop Trachea only on the Sup-Inf axes, without autocrop + // Trachea = clitk::ResizeImageLike(Trachea, input, GetBackgroundValue()); + MaskImagePointType min, max; + GetMinMaxBoundary(input, min, max); + Trachea = clitk::CropImageAlongOneAxis(Trachea, 2, min[2], max[2], + false, GetBackgroundValue()); + + // Select the main CCL (because of bronchus) + Trachea = SliceBySliceKeepMainCCL(Trachea, GetBackgroundValue(), GetForegroundValue()); + + // Slice by slice, build the separation line + clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(Trachea, + GetBackgroundValue(), 2, + extremaDirection, false, // Left + lineDirection, // Vertical line + 1, // margins + tracheaLeftPositionsA, + tracheaLeftPositionsB); + // Do not consider trachea above the limit + int indexMax=tracheaLeftPositionsA.size(); + for(uint i=0; i maxSupPosition) { + indexMax = i; + i = tracheaLeftPositionsA.size(); // stop loop + } + } + tracheaLeftPositionsA.erase(tracheaLeftPositionsA.begin()+indexMax, tracheaLeftPositionsA.end()); + tracheaLeftPositionsB.erase(tracheaLeftPositionsB.begin()+indexMax, tracheaLeftPositionsB.end()); + + // Cut post to this line + clitk::SliceBySliceSetBackgroundFromLineSeparation(input, + tracheaLeftPositionsA, + tracheaLeftPositionsB, + GetBackgroundValue(), + extremaDirection, offset); + MaskImagePointer output = clitk::AutoCrop(input, GetBackgroundValue()); + return output; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_Post_S1S2S4() +{ + StartNewStep("[Support] Post limits of S1RL, S2RL, S4RL"); + + double m_ApexOfTheChest = FindApexOfTheChest(); + + // Post limits with Trachea + MaskImagePointer S1R = m_ListOfSupports["S1R"]; + MaskImagePointer S1L = m_ListOfSupports["S1L"]; + MaskImagePointer S2R = m_ListOfSupports["S2R"]; + MaskImagePointer S2L = m_ListOfSupports["S2L"]; + MaskImagePointer S4R = m_ListOfSupports["S4R"]; + MaskImagePointer S4L = m_ListOfSupports["S4L"]; + S1L = LimitsWithTrachea(S1L, 1, 0, -10, m_ApexOfTheChest); + S1R = LimitsWithTrachea(S1R, 1, 0, -10, m_ApexOfTheChest); + S2R = LimitsWithTrachea(S2R, 1, 0, -10, m_ApexOfTheChest); + S2L = LimitsWithTrachea(S2L, 1, 0, -10, m_ApexOfTheChest); + S4R = LimitsWithTrachea(S4R, 1, 0, -10, m_ApexOfTheChest); + S4L = LimitsWithTrachea(S4L, 1, 0, -10, m_ApexOfTheChest); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_S3P() +{ + StartNewStep("[Support] Ant limits of S3P and Post limits of S1RL, S2RL, S4RL"); + + // Initial S3P support + MaskImagePointer S3P = clitk::Clone(m_ListOfSupports["Support_Superior_to_Carina"]); + + // Stop at Lung Apex + double m_ApexOfTheChest = FindApexOfTheChest(); + S3P = + clitk::CropImageRemoveGreaterThan(S3P, 2, + m_ApexOfTheChest, true, + GetBackgroundValue()); + // Ant limits with Trachea + S3P = LimitsWithTrachea(S3P, 1, 0, 10); + m_ListOfSupports["S3P"] = S3P; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_S3A() +{ + StartNewStep("[Support] Sup-Inf and Post limits for S3A"); + + // Initial S3A support + MaskImagePointer S3A = clitk::Clone(m_ListOfSupports["Support_Superior_to_Carina"]); + + // Stop at Lung Apex or like S2/S1 (upper border Sternum - manubrium) ? + + //double m_ApexOfTheChest = FindApexOfTheChest(); + + MaskImagePointer Sternum = this->GetAFDB()->template GetImage ("Sternum"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(Sternum, GetBackgroundValue(), 2, false, p); + + S3A = + clitk::CropImageRemoveGreaterThan(S3A, 2, + //m_ApexOfTheChest + p[2], true, + GetBackgroundValue()); + // Ant limits with Trachea + S3A = LimitsWithTrachea(S3A, 1, 0, -10); + m_ListOfSupports["S3A"] = S3A; } //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_S5() +{ + StartNewStep("[Support] Sup-Inf limits S5 with aorta"); + + // Initial S5 support + MaskImagePointer S5 = clitk::Clone(this->GetAFDB()->template GetImage("Mediastinum", true)); + + // Sup limits with Aorta + double sup = FindInferiorBorderOfAorticArch(); + + // Inf limits with "upper rim of the left main pulmonary artery" + // For the moment only, it will change. + MaskImagePointer PulmonaryTrunk = this->GetAFDB()->template GetImage("PulmonaryTrunk"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(PulmonaryTrunk, GetBackgroundValue(), 2, false, p); + + // Cut Sup/Inf + S5 = clitk::CropImageAlongOneAxis(S5, 2, p[2], sup, true, GetBackgroundValue()); + + m_ListOfSupports["S5"] = S5; +} +//-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +Support_S6() +{ + StartNewStep("[Support] Sup-Inf limits S6 with aorta"); + + // Initial S6 support like S3A + MaskImagePointer S6 = clitk::Clone(m_ListOfSupports["S3A"]); + + // Inf Sup limits with Aorta + double sup = FindSuperiorBorderOfAorticArch(); + double inf = FindInferiorBorderOfAorticArch(); + + // Cut Sup/Inf + S6 = clitk::CropImageAlongOneAxis(S6, 2, inf, sup, true, GetBackgroundValue()); + + m_ListOfSupports["S6"] = S6; +} +//-------------------------------------------------------------------- + diff --git a/segmentation/clitkExtractLymphStations.ggo b/segmentation/clitkExtractLymphStations.ggo index 96ac994..82ed33b 100644 --- a/segmentation/clitkExtractLymphStations.ggo +++ b/segmentation/clitkExtractLymphStations.ggo @@ -3,44 +3,41 @@ package "clitkExtractLymphStations" version "1.0" purpose "Extract LymphStations with help of TODO" -option "config" - "Config file" string no -option "imagetypes" - "Display allowed image types" flag off -option "verbose" v "Verbose" flag off -option "verboseStep" - "Verbose each step" flag off -option "writeStep" w "Write image at each step" flag off -option "verboseOption" - "Display options values" flag off -option "verboseWarningOff" - "Do not display warning" flag off -option "verboseMemory" - "Display memory usage" flag off +option "config" - "Config file" string no +option "imagetypes" - "Display allowed image types" flag off +option "verbose" v "Verbose" flag off +option "verboseStep" - "Verbose each step" flag off +option "writeStep" w "Write image at each step" flag off +option "verboseOption" - "Display options values" flag off +option "verboseWarningOff" - "Do not display warning" flag off +option "verboseMemory" - "Display memory usage" flag off section "I/O" -option "afdb" a "Input Anatomical Feature DB" string no -option "input" i "Input filename" string no -option "output" o "Output lungs mask filename" string no -option "station" - "Force to compute station even if already exist in the DB" string no multiple +option "afdb" a "Input Anatomical Feature DB" string no +option "input" i "Input filename" string no +option "output" o "Output lungs mask filename" string no + +section "Options for all stations" +option "station" - "Force to compute station even if already exist in the DB" string no multiple +option "nosupport" - "Do not recompute the station supports if already available" flag off +option "relpos" - "List of filenames for relativeposition operations" string no multiple + +# section "Options for Station 3A" +# section "Options for Station 2RL" +# section "Options for Station 1RL" section "Options for Station 8" -option "S8_esophagusDilatationForAnt" - "Dilatation of esophagus, in mm, for 'anterior' limits (default=15,2,1)" double no multiple +option "S8_esophagusDilatationForAnt" - "Dilatation of esophagus, in mm, for 'anterior' limits (default=15,2,1)" double no multiple option "S8_esophagusDilatationForRight" - "Dilatation of esophagus, in mm, for 'right' limits (default=5,10,1)" double no multiple -option "S8_ft_Esophagus" - "Fuzzy Threshold for 'Post' to dilated Esophagus" double default="0.5" no +option "S8_ft_Esophagus" - "Fuzzy Threshold for 'Post' to dilated Esophagus" double default="0.5" no section "Options for Station 7" -option "S7_ft_Bronchi" - "Fuzzy Threshold for Left/Right bronchi" double default="0.1" no -option "S7_ft_LeftSuperiorPulmonaryVein" - "Fuzzy Threshold for LeftSuperiorPulmonaryVein" double default="0.3" no +option "S7_ft_Bronchi" - "Fuzzy Threshold for Left/Right bronchi" double default="0.1" no +option "S7_ft_LeftSuperiorPulmonaryVein" - "Fuzzy Threshold for LeftSuperiorPulmonaryVein" double default="0.3" no option "S7_ft_RightSuperiorPulmonaryVein" - "Fuzzy Threshold for RightSuperiorPulmonaryVein" double default="0.2" no -option "S7_ft_RightPulmonaryArtery" - "Fuzzy Threshold for RightPulmonaryArtery" double default="0.3" no -option "S7_ft_LeftPulmonaryArtery" - "Fuzzy Threshold for LeftPulmonaryArtery (NOT USEFUL YET)" double default="0.5" no -option "S7_ft_SVC" - "Fuzzy Threshold for SVC" double default="0.2" no -option "S7_UseMostInferiorPartOnly" - "Inferior separation with S8." flag off - -section "Options for Station 3A" -option "S3A_ft_Sternum" - "Fuzzy Threshold for Post to Sternum" double default="0.5" no -option "S3A_ft_SubclavianArtery" - "Fuzzy Threshold for Ant to SubclavianArtery" double default="0.2" no - -section "Options for Station 2RL" -option "S2RL_ft_CommonCarotidArtery" - "Threshold for NotAntTo to CommonCarotidArtery" double default="0.7" no -option "S2RL_ft_BrachioCephalicTrunk" - "Threshold for NotAntTo to BrachioCephalicTrunk" double default="0.7" no -option "S2RL_ft_BrachioCephalicVein" - "Threshold for NotAntTo to BrachioCephalicVein" double default="0.3" no -option "S2RL_ft_Aorta" - "Threshold for RightTo to Aorta" double default="0.7" no -option "S2RL_ft_SubclavianArteryRight" - "Threshold for NotLeft to SubclavianArteryRight" double default="0.5" no -option "S2RL_ft_SubclavianArteryLeft" - "Threshold for NotRight to SubclavianArteryLeft" double default="0.8" no +option "S7_ft_RightPulmonaryArtery" - "Fuzzy Threshold for RightPulmonaryArtery" double default="0.3" no +option "S7_ft_LeftPulmonaryArtery" - "Fuzzy Threshold for LeftPulmonaryArtery (NOT USEFUL YET)" double default="0.5" no +option "S7_ft_SVC" - "Fuzzy Threshold for SVC" double default="0.2" no +option "S7_UseMostInferiorPartOnly" - "Inferior separation with S8." flag off + diff --git a/segmentation/clitkExtractLymphStationsFilter.h b/segmentation/clitkExtractLymphStationsFilter.h index dd99c5f..60a73e0 100644 --- a/segmentation/clitkExtractLymphStationsFilter.h +++ b/segmentation/clitkExtractLymphStationsFilter.h @@ -20,8 +20,7 @@ #define CLITKEXTRACTLYMPHSTATIONSFILTER_H // clitk -#include "clitkFilterBase.h" -#include "clitkFilterWithAnatomicalFeatureDatabaseManagement.h" +#include "clitkStructuresExtractionFilter.h" // vtk #include @@ -37,14 +36,12 @@ namespace clitk { template class ITK_EXPORT ExtractLymphStationsFilter: - public virtual clitk::FilterBase, - public clitk::FilterWithAnatomicalFeatureDatabaseManagement, - public itk::ImageToImageFilter > + public clitk::StructuresExtractionFilter { public: /** Standard class typedefs. */ - typedef itk::ImageToImageFilter > Superclass; + typedef clitk::StructuresExtractionFilter Superclass; typedef ExtractLymphStationsFilter Self; typedef itk::SmartPointer Pointer; typedef itk::SmartPointer ConstPointer; @@ -83,15 +80,13 @@ namespace clitk { /** ImageDimension constants */ itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension); FILTERBASE_INIT; - + itkGetConstMacro(BackgroundValue, MaskImagePixelType); itkGetConstMacro(ForegroundValue, MaskImagePixelType); itkSetMacro(BackgroundValue, MaskImagePixelType); itkSetMacro(ForegroundValue, MaskImagePixelType); // Station 8 - // itkSetMacro(DistanceMaxToAnteriorPartOfTheSpine, double); - //itkGetConstMacro(DistanceMaxToAnteriorPartOfTheSpine, double); itkSetMacro(EsophagusDiltationForAnt, MaskImagePointType); itkGetConstMacro(EsophagusDiltationForAnt, MaskImagePointType); itkSetMacro(EsophagusDiltationForRight, MaskImagePointType); @@ -109,6 +104,11 @@ namespace clitk { void AddComputeStation(std::string station) ; void SetFuzzyThreshold(std::string station, std::string tag, double value); double GetFuzzyThreshold(std::string station, std::string tag); + void SetThreshold(std::string station, std::string tag, double value); + double GetThreshold(std::string station, std::string tag); + itkGetConstMacro(ComputeStationsSupportsFlag, bool); + itkSetMacro(ComputeStationsSupportsFlag, bool); + itkBooleanMacro(ComputeStationsSupportsFlag); protected: ExtractLymphStationsFilter(); @@ -118,6 +118,17 @@ namespace clitk { virtual void GenerateInputRequestedRegion(); virtual void GenerateData(); + // To avoid repeat "this->" + AnatomicalFeatureDatabase * GetAFDB() { return clitk::FilterWithAnatomicalFeatureDatabaseManagement::GetAFDB(); } + void WriteAFDB() { clitk::FilterWithAnatomicalFeatureDatabaseManagement::WriteAFDB(); } + void LoadAFDB() { clitk::FilterWithAnatomicalFeatureDatabaseManagement::LoadAFDB(); } + void StartNewStep(std::string s) { clitk::FilterBase::StartNewStep(s); } + void StartSubStep() { clitk::FilterBase::StartSubStep(); } + template + void StopCurrentStep(typename TInternalImageType::Pointer p, std::string txt="") { clitk::FilterBase::StopCurrentStep(p, txt); } + void StopCurrentStep() {clitk::FilterBase::StopCurrentStep(); } + void StopSubStep() {clitk::FilterBase::StopSubStep(); } + ImageConstPointer m_Input; MaskImagePointer m_Mediastinum; MaskImagePointer m_Working_Support; @@ -131,21 +142,43 @@ namespace clitk { void Remove_Structures(std::string station, std::string s); // Functions common to several stations - void FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B); - double FindCarinaSlicePosition(); + double FindCarina(); + double FindApexOfTheChest(); + double FindSuperiorBorderOfAorticArch(); + double FindInferiorBorderOfAorticArch(); void FindLeftAndRightBronchi(); + void FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B); + MaskImagePointer FindAntPostVesselsOLD(); + MaskImagePointer FindAntPostVessels2(); // Global parameters typedef std::map FuzzyThresholdByStructureType; std::map m_FuzzyThreshold; + typedef std::map ThresholdByStructureType; + std::map m_Threshold; // Station's supports void ExtractStationSupports(); - + void Support_SupInf_S1RL(); + void Support_LeftRight_S1R_S1L(); + void Support_SupInf_S2R_S2L(); + void Support_LeftRight_S2R_S2L(); + void Support_SupInf_S4R_S4L(); + void Support_LeftRight_S4R_S4L(); + void Support_Post_S1S2S4(); + void Support_S3P(); + void Support_S3A(); + void Support_S5(); + void Support_S6(); + + MaskImagePointer LimitsWithTrachea(MaskImageType * input, + int extremaDirection, int lineDirection, + double offset, double maxSupPosition); + MaskImagePointer LimitsWithTrachea(MaskImageType * input, + int extremaDirection, int lineDirection, + double offset); // Station 8 - // double m_DistanceMaxToAnteriorPartOfTheSpine; double m_DiaphragmInferiorLimit; - double m_CarinaZ; double m_OriginOfRightMiddleLobeBronchusZ; double m_InjectedThresholdForS8; MaskImagePointer m_Esophagus; @@ -164,32 +197,43 @@ namespace clitk { // Station 3P void ExtractStation_3P(); void ExtractStation_3P_SetDefaultValues(); - void ExtractStation_3P_SI_Limits(); + void ExtractStation_3P_LR_inf_Limits(); + void ExtractStation_3P_LR_sup_Limits_2(); void ExtractStation_3P_Remove_Structures(); - void ExtractStation_3P_Ant_Limits(); - void ExtractStation_3P_Post_Limits(); void ExtractStation_3P_LR_sup_Limits(); - void ExtractStation_3P_LR_sup_Limits_2(); - void ExtractStation_3P_LR_inf_Limits(); + + // Station 3A + void ExtractStation_3A(); + void ExtractStation_3A_SetDefaultValues(); + void ExtractStation_3A_Post_Left_Limits_With_Aorta_S5_Support(); + void ExtractStation_3A_Post_Limits_With_Dilated_Aorta_S6_Support(); + void ExtractStation_3A_AntPost_Superiorly(); + void ExtractStation_3A_Remove_Structures(); // Station 2RL void ExtractStation_2RL(); void ExtractStation_2RL_SetDefaultValues(); - void ExtractStation_2RL_SI_Limits(); - void ExtractStation_2RL_Ant_Limits(); - void ExtractStation_2RL_Ant_Limits2(); - void ExtractStation_2RL_Post_Limits(); - void ExtractStation_2RL_LR_Limits(); - void ExtractStation_2RL_Remove_Structures(); - void ExtractStation_2RL_SeparateRL(); + void ExtractStation_2RL_Ant_Limits(std::string s); + void ExtractStation_2RL_Remove_Structures(std::string s); + void ExtractStation_2RL_Cut_BrachioCephalicVein_superiorly_when_it_split(); vtkSmartPointer Build3DMeshFrom2DContour(const std::vector & points); - // Station 3A - void ExtractStation_3A(); - void ExtractStation_3A_SetDefaultValues(); - void ExtractStation_3A_SI_Limits(); - void ExtractStation_3A_Ant_Limits(); - void ExtractStation_3A_Post_Limits(); + // Station 1RL + void ExtractStation_1RL(); + void ExtractStation_1RL_SetDefaultValues(); + void ExtractStation_1RL_Ant_Limits(); + void ExtractStation_1RL_Post_Limits(); + + // Station 4RL + void ExtractStation_4RL_SetDefaultValues(); + void ExtractStation_4RL(); + void ExtractStation_4RL_SI_Limits(); + void ExtractStation_4RL_LR_Limits(); + void ExtractStation_4RL_AP_Limits(); + MaskImagePointer m_RightSupport; + MaskImagePointer m_LeftSupport; + + // Station 7 void ExtractStation_7(); @@ -197,11 +241,11 @@ namespace clitk { void ExtractStation_7_SI_Limits(); void ExtractStation_7_RL_Interior_Limits(); - void ExtractStation_7_RL_Limits_OLD(); void ExtractStation_7_Posterior_Limits(); void ExtractStation_7_Remove_Structures(); bool m_S7_UseMostInferiorPartOnlyFlag; + bool m_ComputeStationsSupportsFlag; MaskImagePointer m_Working_Trachea; MaskImagePointer m_LeftBronchus; MaskImagePointer m_RightBronchus; @@ -219,14 +263,6 @@ namespace clitk { ListOfPointsType & LR, ListOfPointsType & Ant, ListOfPointsType & Post); - // Station 4RL - void ExtractStation_4RL(); - void ExtractStation_4RL_SI_Limits(); - void ExtractStation_4RL_LR_Limits(); - void ExtractStation_4RL_AP_Limits(); - MaskImagePointer m_RightSupport; - MaskImagePointer m_LeftSupport; - private: ExtractLymphStationsFilter(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented @@ -240,12 +276,13 @@ namespace clitk { #ifndef ITK_MANUAL_INSTANTIATION #include "clitkExtractLymphStationsFilter.txx" #include "clitkExtractLymphStation_Supports.txx" -#include "clitkExtractLymphStation_8.txx" #include "clitkExtractLymphStation_3P.txx" #include "clitkExtractLymphStation_2RL.txx" #include "clitkExtractLymphStation_3A.txx" -#include "clitkExtractLymphStation_7.txx" #include "clitkExtractLymphStation_4RL.txx" +#include "clitkExtractLymphStation_1RL.txx" +#include "clitkExtractLymphStation_8.txx" +#include "clitkExtractLymphStation_7.txx" #endif #endif diff --git a/segmentation/clitkExtractLymphStationsFilter.txx b/segmentation/clitkExtractLymphStationsFilter.txx index 58f2991..f12d3bc 100644 --- a/segmentation/clitkExtractLymphStationsFilter.txx +++ b/segmentation/clitkExtractLymphStationsFilter.txx @@ -27,6 +27,7 @@ #include "clitkAutoCropFilter.h" #include "clitkSegmentationUtils.h" #include "clitkSliceBySliceRelativePositionFilter.h" +#include "clitkMeshToBinaryImageFilter.h" // itk #include @@ -40,24 +41,31 @@ // itk ENST #include "RelativePositionPropImageFilter.h" +// vtk +#include +#include +#include + //-------------------------------------------------------------------- template clitk::ExtractLymphStationsFilter:: ExtractLymphStationsFilter(): - clitk::FilterBase(), - clitk::FilterWithAnatomicalFeatureDatabaseManagement(), - itk::ImageToImageFilter() + clitk::StructuresExtractionFilter() { this->SetNumberOfRequiredInputs(1); SetBackgroundValue(0); SetForegroundValue(1); + ComputeStationsSupportsFlagOn(); // Default values - ExtractStation_8_SetDefaultValues(); ExtractStation_3P_SetDefaultValues(); ExtractStation_2RL_SetDefaultValues(); ExtractStation_3A_SetDefaultValues(); + ExtractStation_1RL_SetDefaultValues(); + ExtractStation_4RL_SetDefaultValues(); + ExtractStation_7_SetDefaultValues(); + ExtractStation_8_SetDefaultValues(); } //-------------------------------------------------------------------- @@ -68,55 +76,78 @@ void clitk::ExtractLymphStationsFilter:: GenerateOutputInformation() { // Get inputs - LoadAFDB(); + this->LoadAFDB(); m_Input = dynamic_cast(itk::ProcessObject::GetInput(0)); - m_Mediastinum = GetAFDB()->template GetImage ("Mediastinum"); + m_Mediastinum = this->GetAFDB()->template GetImage ("Mediastinum"); + + // Clean some computer landmarks to force the recomputation + this->GetAFDB()->RemoveTag("AntPostVesselsSeparation"); // Global supports for stations - StartNewStep("Supports for stations"); - StartSubStep(); - ExtractStationSupports(); - StopSubStep(); + bool supportsExist = true; + try { + m_ListOfSupports["S1R"] = this->GetAFDB()->template GetImage("Support_S1R"); + m_ListOfSupports["S1L"] = this->GetAFDB()->template GetImage("Support_S1L"); + m_ListOfSupports["S2R"] = this->GetAFDB()->template GetImage("Support_S2R"); + m_ListOfSupports["S2L"] = this->GetAFDB()->template GetImage("Support_S2L"); + m_ListOfSupports["S4R"] = this->GetAFDB()->template GetImage("Support_S4R"); + m_ListOfSupports["S4L"] = this->GetAFDB()->template GetImage("Support_S4L"); + + m_ListOfSupports["S3A"] = this->GetAFDB()->template GetImage("Support_S3A"); + m_ListOfSupports["S3P"] = this->GetAFDB()->template GetImage("Support_S3P"); + m_ListOfSupports["S5"] = this->GetAFDB()->template GetImage("Support_S5"); + m_ListOfSupports["S6"] = this->GetAFDB()->template GetImage("Support_S6"); + m_ListOfSupports["S7"] = this->GetAFDB()->template GetImage("Support_S7"); + m_ListOfSupports["S8"] = this->GetAFDB()->template GetImage("Support_S8"); + m_ListOfSupports["S9"] = this->GetAFDB()->template GetImage("Support_S9"); + m_ListOfSupports["S10"] = this->GetAFDB()->template GetImage("Support_S10"); + m_ListOfSupports["S11"] = this->GetAFDB()->template GetImage("Support_S11"); + } catch(clitk::ExceptionObject o) { + supportsExist = false; + } - // Extract Station8 - StartNewStep("Station 8"); - StartSubStep(); - ExtractStation_8(); - StopSubStep(); + if (!supportsExist || GetComputeStationsSupportsFlag()) { + this->StartNewStep("Supports for stations"); + this->StartSubStep(); + this->GetAFDB()->RemoveTag("CarinaZ"); + this->GetAFDB()->RemoveTag("ApexOfTheChestZ"); + this->GetAFDB()->RemoveTag("ApexOfTheChest"); + this->GetAFDB()->RemoveTag("RightBronchus"); + this->GetAFDB()->RemoveTag("LeftBronchus"); + this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArchZ"); + this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArch"); + this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArchZ"); + this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArch"); + ExtractStationSupports(); + this->StopSubStep(); + } // Extract Station3P - StartNewStep("Station 3P"); - StartSubStep(); ExtractStation_3P(); - StopSubStep(); - - // Extract Station2RL - /* - StartNewStep("Station 2RL"); - StartSubStep(); - ExtractStation_2RL(); - StopSubStep(); - */ // Extract Station3A - StartNewStep("Station 3A"); - StartSubStep(); ExtractStation_3A(); - StopSubStep(); + + // Extract Station2RL + ExtractStation_2RL(); + + // Extract Station1RL + ExtractStation_1RL(); + + // Extract Station1RL + ExtractStation_4RL(); + + // ---------- TODO ----------------------- + + // Extract Station8 + ExtractStation_8(); // Extract Station7 - StartNewStep("Station 7"); - StartSubStep(); + this->StartNewStep("Station 7"); + this->StartSubStep(); ExtractStation_7(); - StopSubStep(); - - if (0) { // temporary suppress - // Extract Station4RL - StartNewStep("Station 4RL"); - StartSubStep(); - //ExtractStation_4RL(); - StopSubStep(); - } + this->StopSubStep(); + } //-------------------------------------------------------------------- @@ -154,8 +185,8 @@ CheckForStation(std::string station) // Check if station already exist in DB bool found = false; - if (GetAFDB()->TagExist(s)) { - m_ListOfStations[station] = GetAFDB()->template GetImage(s); + if (this->GetAFDB()->TagExist(s)) { + m_ListOfStations[station] = this->GetAFDB()->template GetImage(s); found = true; } @@ -164,7 +195,7 @@ CheckForStation(std::string station) std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl; } if (!found || GetComputeStation(station)) { - m_Working_Support = m_Mediastinum = GetAFDB()->template GetImage("Mediastinum", true); + m_Working_Support = m_Mediastinum = this->GetAFDB()->template GetImage("Mediastinum", true); return true; } else { @@ -197,88 +228,80 @@ AddComputeStation(std::string station) //-------------------------------------------------------------------- + //-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3P() -{ - if (CheckForStation("3P")) { - ExtractStation_3P_SI_Limits(); - ExtractStation_3P_Ant_Limits(); - ExtractStation_3P_Post_Limits(); - ExtractStation_3P_LR_sup_Limits(); - // ExtractStation_3P_LR_sup_Limits_2(); - ExtractStation_3P_LR_inf_Limits(); - ExtractStation_8_Single_CCL_Limits(); // YES 8 ! - ExtractStation_3P_Remove_Structures(); // after CCL - // Store image filenames into AFDB - writeImage(m_ListOfStations["3P"], "seg/Station3P.mhd"); - GetAFDB()->SetImageFilename("Station3P", "seg/Station3P.mhd"); - WriteAFDB(); - } -} +template +class comparePointsX { +public: + bool operator() (PointType i, PointType j) { return (i[0] -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_3A() -{ - if (CheckForStation("3A")) { - ExtractStation_3A_SI_Limits(); - ExtractStation_3A_Ant_Limits(); - ExtractStation_3A_Post_Limits(); - // Store image filenames into AFDB - writeImage(m_ListOfStations["3A"], "seg/Station3A.mhd"); - GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd"); - WriteAFDB(); +template +class comparePointsWithAngle { +public: + bool operator() (PairType i, PairType j) { return (i.second < j.second); } +}; +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void HypercubeCorners(std::vector > & out) { + std::vector > previous; + HypercubeCorners(previous); + out.resize(previous.size()*2); + for(unsigned int i=0; i p; + if (i -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_2RL() -{ - if (CheckForStation("2RL")) { - ExtractStation_2RL_SI_Limits(); - ExtractStation_2RL_Post_Limits(); - ExtractStation_2RL_Ant_Limits2(); - ExtractStation_2RL_Ant_Limits(); - ExtractStation_2RL_LR_Limits(); - ExtractStation_2RL_Remove_Structures(); - ExtractStation_2RL_SeparateRL(); - - // Store image filenames into AFDB - writeImage(m_ListOfStations["2R"], "seg/Station2R.mhd"); - writeImage(m_ListOfStations["2L"], "seg/Station2L.mhd"); - GetAFDB()->SetImageFilename("Station2R", "seg/Station2R.mhd"); - GetAFDB()->SetImageFilename("Station2L", "seg/Station2L.mhd"); - WriteAFDB(); - } +template<> +void HypercubeCorners<1>(std::vector > & out) { + out.resize(2); + out[0][0] = 0; + out[1][0] = 1; } //-------------------------------------------------------------------- //-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -ExtractStation_4RL() { - DD("TODO"); - exit(0); - /* - WARNING ONLY 4R FIRST !!! (not same inf limits) - */ - ExtractStation_4RL_SI_Limits(); - ExtractStation_4RL_LR_Limits(); - ExtractStation_4RL_AP_Limits(); +template +void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image, + std::vector & bounds) +{ + // Get image max/min coordinates + const unsigned int dim=ImageType::ImageDimension; + typedef typename ImageType::PointType PointType; + typedef typename ImageType::IndexType IndexType; + PointType min_c, max_c; + IndexType min_i, max_i; + min_i = image->GetLargestPossibleRegion().GetIndex(); + for(unsigned int i=0; iGetLargestPossibleRegion().GetSize()[i] + min_i[i]; + image->TransformIndexToPhysicalPoint(min_i, min_c); + image->TransformIndexToPhysicalPoint(max_i, max_c); + + // Get corners coordinates + HypercubeCorners(bounds); + for(unsigned int i=0; i:: Remove_Structures(std::string station, std::string s) { try { - StartNewStep("[Station"+station+"] Remove "+s); - MaskImagePointer Structure = GetAFDB()->template GetImage(s); + this->StartNewStep("[Station"+station+"] Remove "+s); + MaskImagePointer Structure = this->GetAFDB()->template GetImage(s); clitk::AndNot(m_Working_Support, Structure, GetBackgroundValue()); } catch(clitk::ExceptionObject e) { @@ -312,6 +335,17 @@ SetFuzzyThreshold(std::string station, std::string tag, double value) //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +SetThreshold(std::string station, std::string tag, double value) +{ + m_Threshold[station][tag] = value; +} +//-------------------------------------------------------------------- + + //-------------------------------------------------------------------- template double @@ -333,6 +367,27 @@ GetFuzzyThreshold(std::string station, std::string tag) //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +double +clitk::ExtractLymphStationsFilter:: +GetThreshold(std::string station, std::string tag) +{ + if (m_Threshold.find(station) == m_Threshold.end()) { + clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+")."); + return 0.0; + } + + if (m_Threshold[station].find(tag) == m_Threshold[station].end()) { + clitkExceptionMacro("Could not find options "+tag+" in the list of Threshold for station "+station+"."); + return 0.0; + } + + return m_Threshold[station][tag]; +} +//-------------------------------------------------------------------- + + //-------------------------------------------------------------------- template void @@ -344,22 +399,22 @@ FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B) // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus try { - GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A); - GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B); + this->GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A); + this->GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B); } catch(clitk::ExceptionObject & o) { DD("FindLineForS7S8Separation"); // Load LeftLowerLobeBronchus and get centroid point MaskImagePointer LeftLowerLobeBronchus = - GetAFDB()->template GetImage ("LeftLowerLobeBronchus"); + this->GetAFDB()->template GetImage ("LeftLowerLobeBronchus"); std::vector c; clitk::ComputeCentroids(LeftLowerLobeBronchus, GetBackgroundValue(), c); A = c[1]; // Load RightMiddleLobeBronchus and get superior point (not centroid here) MaskImagePointer RightMiddleLobeBronchus = - GetAFDB()->template GetImage ("RightMiddleLobeBronchus"); + this->GetAFDB()->template GetImage ("RightMiddleLobeBronchus"); bool b = FindExtremaPointInAGivenDirection(RightMiddleLobeBronchus, GetBackgroundValue(), 2, false, B); @@ -368,8 +423,8 @@ FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B) } // Insert into the DB - GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A); - GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B); + this->GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A); + this->GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B); } } //-------------------------------------------------------------------- @@ -379,28 +434,28 @@ FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B) template double clitk::ExtractLymphStationsFilter:: -FindCarinaSlicePosition() +FindCarina() { double z; try { - z = GetAFDB()->GetDouble("CarinaZ"); + z = this->GetAFDB()->GetDouble("CarinaZ"); } catch(clitk::ExceptionObject e) { DD("FindCarinaSlicePosition"); // Get Carina - MaskImagePointer Carina = GetAFDB()->template GetImage("Carina"); + MaskImagePointer Carina = this->GetAFDB()->template GetImage("Carina"); // Get Centroid and Z value std::vector centroids; clitk::ComputeCentroids(Carina, GetBackgroundValue(), centroids); // We dont need Carina structure from now - Carina->Delete(); + this->GetAFDB()->template ReleaseImage("Carina"); // Put inside the AFDB - GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]); - GetAFDB()->SetDouble("CarinaZ", centroids[1][2]); - WriteAFDB(); + this->GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]); + this->GetAFDB()->SetDouble("CarinaZ", centroids[1][2]); + this->WriteAFDB(); z = centroids[1][2]; } return z; @@ -408,6 +463,38 @@ FindCarinaSlicePosition() //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +double +clitk::ExtractLymphStationsFilter:: +FindApexOfTheChest() +{ + double z; + try { + z = this->GetAFDB()->GetDouble("ApexOfTheChestZ"); + } + catch(clitk::ExceptionObject e) { + DD("FindApexOfTheChestPosition"); + MaskImagePointer Lungs = this->GetAFDB()->template GetImage("Lungs"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(Lungs, GetBackgroundValue(), 2, false, p); + + // We dont need Lungs structure from now + this->GetAFDB()->template ReleaseImage("Lungs"); + + // Put inside the AFDB + this->GetAFDB()->SetPoint3D("ApexOfTheChest", p); + p[2] -= 5; // We consider 5 mm lower + this->GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]); + this->WriteAFDB(); + z = p[2]; + } + return z; +} +//-------------------------------------------------------------------- + + //-------------------------------------------------------------------- template void @@ -415,8 +502,8 @@ clitk::ExtractLymphStationsFilter:: FindLeftAndRightBronchi() { try { - m_RightBronchus = GetAFDB()->template GetImage ("RightBronchus"); - m_LeftBronchus = GetAFDB()->template GetImage ("LeftBronchus"); + m_RightBronchus = this->GetAFDB()->template GetImage ("RightBronchus"); + m_LeftBronchus = this->GetAFDB()->template GetImage ("LeftBronchus"); } catch(clitk::ExceptionObject & o) { @@ -425,10 +512,10 @@ FindLeftAndRightBronchi() // a Left and Right bronchus. // Get the trachea - MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); + MaskImagePointer Trachea = this->GetAFDB()->template GetImage("Trachea"); // Get the Carina position - m_CarinaZ = FindCarinaSlicePosition(); + double m_CarinaZ = FindCarina(); // Consider only inferiorly to the Carina MaskImagePointer m_Working_Trachea = @@ -492,12 +579,813 @@ FindLeftAndRightBronchi() LeftBronchus = clitk::AutoCrop(LeftBronchus, GetBackgroundValue()); // Insert int AFDB if need after - GetAFDB()->template SetImage ("RightBronchus", "seg/rightBronchus.mhd", + this->GetAFDB()->template SetImage ("RightBronchus", "seg/rightBronchus.mhd", RightBronchus, true); - GetAFDB()->template SetImage ("LeftBronchus", "seg/leftBronchus.mhd", + this->GetAFDB()->template SetImage ("LeftBronchus", "seg/leftBronchus.mhd", LeftBronchus, true); } } //-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +double +clitk::ExtractLymphStationsFilter:: +FindSuperiorBorderOfAorticArch() +{ + double z; + try { + z = this->GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ"); + } + catch(clitk::ExceptionObject e) { + DD("FindSuperiorBorderOfAorticArch"); + MaskImagePointer Aorta = this->GetAFDB()->template GetImage("Aorta"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(Aorta, GetBackgroundValue(), 2, false, p); + p[2] += Aorta->GetSpacing()[2]; // the slice above + + // We dont need Lungs structure from now + this->GetAFDB()->template ReleaseImage("Aorta"); + + // Put inside the AFDB + this->GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p); + this->GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]); + this->WriteAFDB(); + z = p[2]; + } + return z; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +double +clitk::ExtractLymphStationsFilter:: +FindInferiorBorderOfAorticArch() +{ + double z; + try { + z = this->GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ"); + } + catch(clitk::ExceptionObject e) { + DD("FindInferiorBorderOfAorticArch"); + MaskImagePointer Aorta = this->GetAFDB()->template GetImage("Aorta"); + std::vector slices; + clitk::ExtractSlices(Aorta, 2, slices); + bool found=false; + uint i = slices.size()-1; + while (!found) { + slices[i] = Labelize(slices[i], 0, false, 10); + std::vector c; + clitk::ComputeCentroids(slices[i], GetBackgroundValue(), c); + if (c.size()>2) { + found = true; + } + else { + i--; + } + } + MaskImageIndexType index; + index[0] = index[1] = 0.0; + index[2] = i+1; + MaskImagePointType lower; + Aorta->TransformIndexToPhysicalPoint(index, lower); + + // We dont need Lungs structure from now + this->GetAFDB()->template ReleaseImage("Aorta"); + + // Put inside the AFDB + this->GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower); + this->GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]); + this->WriteAFDB(); + z = lower[2]; + } + return z; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +typename clitk::ExtractLymphStationsFilter::MaskImagePointer +clitk::ExtractLymphStationsFilter:: +FindAntPostVesselsOLD() +{ + // ----------------------------------------------------- + /* Rod says: "The anterior border, as with the Atlas – UM, is + posterior to the vessels (right subclavian vein, left + brachiocephalic vein, right brachiocephalic vein, left subclavian + artery, left common carotid artery and brachiocephalic trunk). + These vessels are not included in the nodal station. The anterior + border is drawn to the midpoint of the vessel and an imaginary + line joins the middle of these vessels. Between the vessels, + station 2 is in contact with station 3a." */ + + // Check if no already done + bool found = true; + MaskImagePointer binarizedContour; + try { + binarizedContour = this->GetAFDB()->template GetImage ("AntPostVesselsSeparation"); + } + catch(clitk::ExceptionObject e) { + found = false; + } + if (found) { + return binarizedContour; + } + + /* Here, we consider the vessels form a kind of anterior barrier. We + link all vessels centroids and remove what is post to it. + - select the list of structure + vessel1 = BrachioCephalicArtery + vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right) + vessel3 = CommonCarotidArtery + vessel4 = SubclavianArtery + other = Thyroid + other = Aorta + - crop images as needed + - slice by slice, choose the CCL and extract centroids + - slice by slice, sort according to polar angle wrt Trachea centroid. + - slice by slice, link centroids and close contour + - remove outside this contour + - merge with support + */ + + // Read structures + MaskImagePointer BrachioCephalicArtery = this->GetAFDB()->template GetImage("BrachioCephalicArtery"); + MaskImagePointer BrachioCephalicVein = this->GetAFDB()->template GetImage("BrachioCephalicVein"); + MaskImagePointer CommonCarotidArtery = this->GetAFDB()->template GetImage("CommonCarotidArtery"); + MaskImagePointer SubclavianArtery = this->GetAFDB()->template GetImage("SubclavianArtery"); + MaskImagePointer Thyroid = this->GetAFDB()->template GetImage("Thyroid"); + MaskImagePointer Aorta = this->GetAFDB()->template GetImage("Aorta"); + MaskImagePointer Trachea = this->GetAFDB()->template GetImage("Trachea"); + + // Create a temporay support + // From first slice of BrachioCephalicVein to end of 3A + MaskImagePointer support = this->GetAFDB()->template GetImage("Support_Sup_Carina"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(BrachioCephalicVein, GetBackgroundValue(), 2, true, p); + double inf = p [2]; + clitk::FindExtremaPointInAGivenDirection(this->GetAFDB()->template GetImage("Support_S3A"), + GetBackgroundValue(), 2, false, p); + double sup = p [2]; + support = clitk::CropImageAlongOneAxis(support, 2, inf, sup, + false, GetBackgroundValue()); + + // Resize all structures like support + BrachioCephalicArtery = + clitk::ResizeImageLike(BrachioCephalicArtery, support, GetBackgroundValue()); + CommonCarotidArtery = + clitk::ResizeImageLike(CommonCarotidArtery, support, GetBackgroundValue()); + SubclavianArtery = + clitk::ResizeImageLike(SubclavianArtery, support, GetBackgroundValue()); + Thyroid = + clitk::ResizeImageLike(Thyroid, support, GetBackgroundValue()); + Aorta = + clitk::ResizeImageLike(Aorta, support, GetBackgroundValue()); + BrachioCephalicVein = + clitk::ResizeImageLike(BrachioCephalicVein, support, GetBackgroundValue()); + Trachea = + clitk::ResizeImageLike(Trachea, support, GetBackgroundValue()); + + // Extract slices + std::vector slices_BrachioCephalicArtery; + clitk::ExtractSlices(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery); + std::vector slices_BrachioCephalicVein; + clitk::ExtractSlices(BrachioCephalicVein, 2, slices_BrachioCephalicVein); + std::vector slices_CommonCarotidArtery; + clitk::ExtractSlices(CommonCarotidArtery, 2, slices_CommonCarotidArtery); + std::vector slices_SubclavianArtery; + clitk::ExtractSlices(SubclavianArtery, 2, slices_SubclavianArtery); + std::vector slices_Thyroid; + clitk::ExtractSlices(Thyroid, 2, slices_Thyroid); + std::vector slices_Aorta; + clitk::ExtractSlices(Aorta, 2, slices_Aorta); + std::vector slices_Trachea; + clitk::ExtractSlices(Trachea, 2, slices_Trachea); + unsigned int n= slices_BrachioCephalicArtery.size(); + + // Get the boundaries of one slice + std::vector bounds; + ComputeImageBoundariesCoordinates(slices_BrachioCephalicArtery[0], bounds); + + // For all slices, for all structures, find the centroid and build the contour + // List of 3D points (for debug) + std::vector p3D; + + vtkSmartPointer append = vtkSmartPointer::New(); + for(unsigned int i=0; i(slices_CommonCarotidArtery[i], + GetBackgroundValue(), true, 1); + slices_SubclavianArtery[i] = Labelize(slices_SubclavianArtery[i], + GetBackgroundValue(), true, 1); + slices_BrachioCephalicArtery[i] = Labelize(slices_BrachioCephalicArtery[i], + GetBackgroundValue(), true, 1); + slices_BrachioCephalicVein[i] = Labelize(slices_BrachioCephalicVein[i], + GetBackgroundValue(), true, 1); + slices_Thyroid[i] = Labelize(slices_Thyroid[i], + GetBackgroundValue(), true, 1); + slices_Aorta[i] = Labelize(slices_Aorta[i], + GetBackgroundValue(), true, 1); + + // Search centroids + std::vector points2D; + std::vector centroids1; + std::vector centroids2; + std::vector centroids3; + std::vector centroids4; + std::vector centroids5; + std::vector centroids6; + ComputeCentroids(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1); + ComputeCentroids(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2); + ComputeCentroids(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3); + ComputeCentroids(slices_Thyroid[i], GetBackgroundValue(), centroids4); + ComputeCentroids(slices_Aorta[i], GetBackgroundValue(), centroids5); + ComputeCentroids(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6); + + // BrachioCephalicVein -> when it is separated into two CCL, we + // only consider the most at Right one + if (centroids6.size() > 2) { + if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2); + else centroids6.erase(centroids6.begin()+1); + } + + // BrachioCephalicVein -> when SubclavianArtery has 2 CCL + // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein + if ((centroids3.size() ==1) && (centroids2.size() > 2)) { + centroids6.clear(); + } + + for(unsigned int j=1; j centroids_trachea; + ComputeCentroids(slices_Trachea[i], GetBackgroundValue(), centroids_trachea); + typedef std::pair PointAngleType; + std::vector angles; + for(unsigned int j=0; j0) angle = atan(y/x); + if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI; + if ((x<0) && (y<0)) angle = atan(y/x)-M_PI; + if (x==0) { + if (y>0) angle = M_PI/2.0; + if (y<0) angle = -M_PI/2.0; + if (y==0) angle = 0; + } + angle = clitk::rad2deg(angle); + // Angle is [-180;180] wrt the X axis. We change the X axis to + // be the vertical line, because we want to sort from Right to + // Left from Post to Ant. + if (angle>0) + angle = (270-angle); + if (angle<0) { + angle = -angle-90; + if (angle<0) angle = 360-angle; + } + PointAngleType a; + a.first = points2D[j]; + a.second = angle; + angles.push_back(a); + } + + // Do nothing if less than 2 points --> n + if (points2D.size() < 3) { //continue; + continue; + } + + // Sort points2D according to polar angles + std::sort(angles.begin(), angles.end(), comparePointsWithAngle()); + for(unsigned int j=0; j toadd; + unsigned int index = 0; + double dmax = 5; + while (indexdmax) { + + MaskSlicePointType b; + b[0] = a[0]+(c[0]-a[0])/2.0; + b[1] = a[1]+(c[1]-a[1])/2.0; + + // Compute distance to trachea centroid + MaskSlicePointType m = centroids_trachea[1]; + double da = m.EuclideanDistanceTo(a); + double db = m.EuclideanDistanceTo(b); + //double dc = m.EuclideanDistanceTo(c); + + // Mean distance, find point on the line from trachea centroid + // to b + double alpha = (da+db)/2.0; + MaskSlicePointType v; + double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2)); + v[0] = (b[0]-m[0])/n; + v[1] = (b[1]-m[1])/n; + MaskSlicePointType r; + r[0] = m[0]+alpha*v[0]; + r[1] = m[1]+alpha*v[1]; + points2D.insert(points2D.begin()+index+1, r); + } + else { + index++; + } + } + // DDV(points2D, points2D.size()); + + // Add some points to close the contour + // - H line towards Right + MaskSlicePointType p = points2D[0]; + p[0] = bounds[0][0]; + points2D.insert(points2D.begin(), p); + // - corner Right/Post + p = bounds[0]; + points2D.insert(points2D.begin(), p); + // - H line towards Left + p = points2D.back(); + p[0] = bounds[2][0]; + points2D.push_back(p); + // - corner Right/Post + p = bounds[2]; + points2D.push_back(p); + // Close contour with the first point + points2D.push_back(points2D[0]); + // DDV(points2D, points2D.size()); + + // Build 3D points from the 2D points + std::vector points3D; + clitk::PointsUtils::Convert2DListTo3DList(points2D, i, support, points3D); + for(unsigned int x=0; x mesh = Build3DMeshFrom2DContour(points3D); + append->AddInput(mesh); + } + + // DEBUG: write points3D + clitk::WriteListOfLandmarks(p3D, "vessels-centroids.txt"); + + // Build the final 3D mesh form the list 2D mesh + append->Update(); + vtkSmartPointer mesh = append->GetOutput(); + + // Debug, write the mesh + /* + vtkSmartPointer w = vtkSmartPointer::New(); + w->SetInput(mesh); + w->SetFileName("bidon.vtk"); + w->Write(); + */ + + // Compute a single binary 3D image from the list of contours + clitk::MeshToBinaryImageFilter::Pointer filter = + clitk::MeshToBinaryImageFilter::New(); + filter->SetMesh(mesh); + filter->SetLikeImage(support); + filter->Update(); + binarizedContour = filter->GetOutput(); + + // Crop + clitk::FindExtremaPointInAGivenDirection(binarizedContour, GetForegroundValue(), 2, true, p); + inf = p[2]; + DD(p); + clitk::FindExtremaPointInAGivenDirection(binarizedContour, GetForegroundValue(), 2, false, p); + sup = p[2]; + DD(p); + binarizedContour = clitk::CropImageAlongOneAxis(binarizedContour, 2, inf, sup, + false, GetBackgroundValue()); + // Update the AFDB + writeImage(binarizedContour, "seg/AntPostVesselsSeparation.mhd"); + this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd"); + this->WriteAFDB(); + return binarizedContour; + + /* + // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse + ImagePointType p = p3D[2]; // This is the first centroid of the first slice + p[1] += 50; // 50 mm Post from this point must be kept + ImageIndexType index; + binarizedContour->TransformPhysicalPointToIndex(p, index); + bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue()); + + // remove from support + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(m_Working_Support); + boolFilter->SetInput2(binarizedContour); + boolFilter->SetBackgroundValue1(GetBackgroundValue()); + boolFilter->SetBackgroundValue2(GetBackgroundValue()); + if (isInside) + boolFilter->SetOperationType(BoolFilterType::And); + else + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + m_Working_Support = boolFilter->GetOutput(); + */ + + // End + //StopCurrentStep(m_Working_Support); + //m_ListOfStations["2R"] = m_Working_Support; + //m_ListOfStations["2L"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +typename clitk::ExtractLymphStationsFilter::MaskImagePointer +clitk::ExtractLymphStationsFilter:: +FindAntPostVessels2() +{ + // ----------------------------------------------------- + /* Rod says: "The anterior border, as with the Atlas – UM, is + posterior to the vessels (right subclavian vein, left + brachiocephalic vein, right brachiocephalic vein, left subclavian + artery, left common carotid artery and brachiocephalic trunk). + These vessels are not included in the nodal station. The anterior + border is drawn to the midpoint of the vessel and an imaginary + line joins the middle of these vessels. Between the vessels, + station 2 is in contact with station 3a." */ + + // Check if no already done + bool found = true; + MaskImagePointer binarizedContour; + try { + binarizedContour = this->GetAFDB()->template GetImage ("AntPostVesselsSeparation"); + } + catch(clitk::ExceptionObject e) { + found = false; + } + if (found) { + return binarizedContour; + } + + /* Here, we consider the vessels form a kind of anterior barrier. We + link all vessels centroids and remove what is post to it. + - select the list of structure + vessel1 = BrachioCephalicArtery + vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right) + vessel3 = CommonCarotidArtery + vessel4 = SubclavianArtery + other = Thyroid + other = Aorta + - crop images as needed + - slice by slice, choose the CCL and extract centroids + - slice by slice, sort according to polar angle wrt Trachea centroid. + - slice by slice, link centroids and close contour + - remove outside this contour + - merge with support + */ + + // Read structures + std::map MapOfStructures; + typedef std::map::iterator MapIter; + MapOfStructures["BrachioCephalicArtery"] = this->GetAFDB()->template GetImage("BrachioCephalicArtery"); + MapOfStructures["BrachioCephalicVein"] = this->GetAFDB()->template GetImage("BrachioCephalicVein"); + MapOfStructures["CommonCarotidArteryLeft"] = this->GetAFDB()->template GetImage("CommonCarotidArteryLeft"); + MapOfStructures["CommonCarotidArteryRight"] = this->GetAFDB()->template GetImage("CommonCarotidArteryRight"); + MapOfStructures["SubclavianArteryLeft"] = this->GetAFDB()->template GetImage("SubclavianArteryLeft"); + MapOfStructures["SubclavianArteryRight"] = this->GetAFDB()->template GetImage("SubclavianArteryRight"); + MapOfStructures["Thyroid"] = this->GetAFDB()->template GetImage("Thyroid"); + MapOfStructures["Aorta"] = this->GetAFDB()->template GetImage("Aorta"); + MapOfStructures["Trachea"] = this->GetAFDB()->template GetImage("Trachea"); + + std::vector ListOfStructuresNames; + + // Create a temporay support + // From first slice of BrachioCephalicVein to end of 3A or end of 2RL + MaskImagePointer support = this->GetAFDB()->template GetImage("Support_Sup_Carina"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(MapOfStructures["BrachioCephalicVein"], + GetBackgroundValue(), 2, true, p); + double inf = p[2]; + clitk::FindExtremaPointInAGivenDirection(this->GetAFDB()->template GetImage("Support_S3A"), + GetBackgroundValue(), 2, false, p); + MaskImagePointType p2; + clitk::FindExtremaPointInAGivenDirection(this->GetAFDB()->template GetImage("Support_S2L"), + GetBackgroundValue(), 2, false, p2); + if (p2[2] > p[2]) p = p2; + + double sup = p[2]+support->GetSpacing()[2];//one slice more ? + support = clitk::CropImageAlongOneAxis(support, 2, inf, sup, + false, GetBackgroundValue()); + + // Resize all structures like support + for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) { + it->second = clitk::ResizeImageLike(it->second, support, GetBackgroundValue()); + } + + // Extract slices + typedef std::vector SlicesType; + std::map MapOfSlices; + for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) { + SlicesType s; + clitk::ExtractSlices(it->second, 2, s); + MapOfSlices[it->first] = s; + } + + unsigned int n= MapOfSlices["Trachea"].size(); + + // Get the boundaries of one slice + std::vector bounds; + ComputeImageBoundariesCoordinates(MapOfSlices["Trachea"][0], bounds); + + // LOOP ON SLICES + // For all slices, for all structures, find the centroid and build the contour + // List of 3D points (for debug) + std::vector p3D; + vtkSmartPointer append = vtkSmartPointer::New(); + for(unsigned int i=0; ifirst][i]; + s = clitk::Labelize(s, GetBackgroundValue(), true, 1); + } + + // Search centroids + std::vector points2D; + typedef std::vector CentroidsType; + std::map MapOfCentroids; + for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) { + std::string structure = it->first; + MaskSlicePointer & s = MapOfSlices[structure][i]; + CentroidsType c; + clitk::ComputeCentroids(s, GetBackgroundValue(), c); + MapOfCentroids[structure] = c; + } + + + // BrachioCephalicVein -> when it is separated into two CCL, we + // only consider the most at Right one + CentroidsType & c = MapOfCentroids["BrachioCephalicVein"]; + if (c.size() > 2) { + if (c[1][0] < c[2][0]) c.erase(c.begin()+2); + else c.erase(c.begin()+1); + } + + /* + // BrachioCephalicVein -> when SubclavianArtery has 2 CCL + // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein + if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1) + && (MapOfCentroids["SubclavianArtery"].size() > 2)) { + MapOfCentroids["BrachioCephalicVein"].clear(); + } + */ + + // Add all 2D points + for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) { + std::string structure = it->first; + if (structure != "Trachea") { // do not add centroid of trachea + CentroidsType & c = MapOfCentroids[structure]; + for(unsigned int j=1; j centroids_trachea; + //ComputeCentroids(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea); + CentroidsType centroids_trachea = MapOfCentroids["Trachea"]; + typedef std::pair PointAngleType; + std::vector angles; + for(unsigned int j=0; j0) angle = atan(y/x); + if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI; + if ((x<0) && (y<0)) angle = atan(y/x)-M_PI; + if (x==0) { + if (y>0) angle = M_PI/2.0; + if (y<0) angle = -M_PI/2.0; + if (y==0) angle = 0; + } + angle = clitk::rad2deg(angle); + // Angle is [-180;180] wrt the X axis. We change the X axis to + // be the vertical line, because we want to sort from Right to + // Left from Post to Ant. + if (angle>0) + angle = (270-angle); + if (angle<0) { + angle = -angle-90; + if (angle<0) angle = 360-angle; + } + PointAngleType a; + a.first = points2D[j]; + a.second = angle; + angles.push_back(a); + } + + // Do nothing if less than 2 points --> n + if (points2D.size() < 3) { //continue; + continue; + } + + // Sort points2D according to polar angles + std::sort(angles.begin(), angles.end(), comparePointsWithAngle()); + for(unsigned int j=0; j toadd; + unsigned int index = 0; + double dmax = 5; + while (indexdmax) { + + MaskSlicePointType b; + b[0] = a[0]+(c[0]-a[0])/2.0; + b[1] = a[1]+(c[1]-a[1])/2.0; + + // Compute distance to trachea centroid + MaskSlicePointType m = centroids_trachea[1]; + double da = m.EuclideanDistanceTo(a); + double db = m.EuclideanDistanceTo(b); + //double dc = m.EuclideanDistanceTo(c); + + // Mean distance, find point on the line from trachea centroid + // to b + double alpha = (da+db)/2.0; + MaskSlicePointType v; + double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2)); + v[0] = (b[0]-m[0])/n; + v[1] = (b[1]-m[1])/n; + MaskSlicePointType r; + r[0] = m[0]+alpha*v[0]; + r[1] = m[1]+alpha*v[1]; + points2D.insert(points2D.begin()+index+1, r); + } + else { + index++; + } + } + // DDV(points2D, points2D.size()); + + // Add some points to close the contour + // - H line towards Right + MaskSlicePointType p = points2D[0]; + p[0] = bounds[0][0]; + points2D.insert(points2D.begin(), p); + // - corner Right/Post + p = bounds[0]; + points2D.insert(points2D.begin(), p); + // - H line towards Left + p = points2D.back(); + p[0] = bounds[2][0]; + points2D.push_back(p); + // - corner Right/Post + p = bounds[2]; + points2D.push_back(p); + // Close contour with the first point + points2D.push_back(points2D[0]); + // DDV(points2D, points2D.size()); + + // Build 3D points from the 2D points + std::vector points3D; + clitk::PointsUtils::Convert2DListTo3DList(points2D, i, support, points3D); + for(unsigned int x=0; x mesh = Build3DMeshFrom2DContour(points3D); + append->AddInput(mesh); + // if (i ==n-1) { // last slice + // clitk::PointsUtils::Convert2DListTo3DList(points2D, i+1, support, points3D); + // vtkSmartPointer mesh = Build3DMeshFrom2DContour(points3D); + // append->AddInput(mesh); + // } + } + + // DEBUG: write points3D + clitk::WriteListOfLandmarks(p3D, "vessels-centroids.txt"); + + // Build the final 3D mesh form the list 2D mesh + append->Update(); + vtkSmartPointer mesh = append->GetOutput(); + + // Debug, write the mesh + /* + vtkSmartPointer w = vtkSmartPointer::New(); + w->SetInput(mesh); + w->SetFileName("bidon.vtk"); + w->Write(); + */ + + // Compute a single binary 3D image from the list of contours + clitk::MeshToBinaryImageFilter::Pointer filter = + clitk::MeshToBinaryImageFilter::New(); + filter->SetMesh(mesh); + filter->SetLikeImage(support); + filter->Update(); + binarizedContour = filter->GetOutput(); + + // Crop + clitk::FindExtremaPointInAGivenDirection(binarizedContour, + GetForegroundValue(), 2, true, p); + inf = p[2]; + clitk::FindExtremaPointInAGivenDirection(binarizedContour, + GetForegroundValue(), 2, false, p); + sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice + binarizedContour = clitk::CropImageAlongOneAxis(binarizedContour, 2, inf, sup, + false, GetBackgroundValue()); + + // Update the AFDB + writeImage(binarizedContour, "seg/AntPostVesselsSeparation.mhd"); + this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd"); + this->WriteAFDB(); + return binarizedContour; + + /* + // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse + ImagePointType p = p3D[2]; // This is the first centroid of the first slice + p[1] += 50; // 50 mm Post from this point must be kept + ImageIndexType index; + binarizedContour->TransformPhysicalPointToIndex(p, index); + bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue()); + + // remove from support + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(m_Working_Support); + boolFilter->SetInput2(binarizedContour); + boolFilter->SetBackgroundValue1(GetBackgroundValue()); + boolFilter->SetBackgroundValue2(GetBackgroundValue()); + if (isInside) + boolFilter->SetOperationType(BoolFilterType::And); + else + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + m_Working_Support = boolFilter->GetOutput(); + */ + + // End + //StopCurrentStep(m_Working_Support); + //m_ListOfStations["2R"] = m_Working_Support; + //m_ListOfStations["2L"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + + #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX diff --git a/segmentation/clitkExtractLymphStationsGenericFilter.txx b/segmentation/clitkExtractLymphStationsGenericFilter.txx index 8b6167e..317c78a 100644 --- a/segmentation/clitkExtractLymphStationsGenericFilter.txx +++ b/segmentation/clitkExtractLymphStationsGenericFilter.txx @@ -69,6 +69,8 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) f->SetVerboseMemoryFlag(mArgsInfo.verboseMemory_flag); f->SetAFDBFilename(mArgsInfo.afdb_arg); + f->SetComputeStationsSupportsFlag(!mArgsInfo.nosupport_flag); + // Station 8 //f->SetDistanceMaxToAnteriorPartOfTheSpine(mArgsInfo.S8_maxAntSpine_arg); f->SetFuzzyThreshold("8", "Esophagus", mArgsInfo.S8_ft_Esophagus_arg); @@ -111,6 +113,8 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) for(uint i=0; iAddComputeStation(mArgsInfo.station_arg[i]); + // Station 3A + // Station 7 f->SetFuzzyThreshold("7", "Bronchi", mArgsInfo.S7_ft_Bronchi_arg); f->SetFuzzyThreshold("7", "LeftSuperiorPulmonaryVein", mArgsInfo.S7_ft_LeftSuperiorPulmonaryVein_arg); @@ -120,17 +124,13 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) f->SetFuzzyThreshold("7", "SVC", mArgsInfo.S7_ft_SVC_arg); f->SetS7_UseMostInferiorPartOnlyFlag(mArgsInfo.S7_UseMostInferiorPartOnly_flag); - // Station 3A - f->SetFuzzyThreshold("3A", "Sternum", mArgsInfo.S3A_ft_Sternum_arg); - f->SetFuzzyThreshold("3A", "SubclavianArtery", mArgsInfo.S3A_ft_SubclavianArtery_arg); - // Station 2RL - f->SetFuzzyThreshold("2RL", "CommonCarotidArtery", mArgsInfo.S2RL_ft_CommonCarotidArtery_arg); - f->SetFuzzyThreshold("2RL", "BrachioCephalicTrunk", mArgsInfo.S2RL_ft_BrachioCephalicTrunk_arg); - f->SetFuzzyThreshold("2RL", "BrachioCephalicVein", mArgsInfo.S2RL_ft_BrachioCephalicVein_arg); - f->SetFuzzyThreshold("2RL", "Aorta", mArgsInfo.S2RL_ft_Aorta_arg); - f->SetFuzzyThreshold("2RL", "SubclavianArteryLeft", mArgsInfo.S2RL_ft_SubclavianArteryLeft_arg); - f->SetFuzzyThreshold("2RL", "SubclavianArteryRight", mArgsInfo.S2RL_ft_SubclavianArteryRight_arg); + + // Station 1RL + + // Set RelativePositionList filenames + for(uint i=0; iAddRelativePositionListFilename(mArgsInfo.relpos_arg[i]); } //-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractMediastinum.ggo b/segmentation/clitkExtractMediastinum.ggo index bbb473d..bd2bbe4 100644 --- a/segmentation/clitkExtractMediastinum.ggo +++ b/segmentation/clitkExtractMediastinum.ggo @@ -17,22 +17,8 @@ section "I/O" option "input" i "Input CT filename" string yes option "output" o "Output lungs mask filename" string yes option "afdb" a "Input Anatomical Feature DB (needs Patient, Lungs, Trachea, VertebralBody)" string no +option "relpos" - "List of filenames for relativeposition operations" string no multiple + option "useBones" - "If set : do use bones mask (when image is not injected)" flag off option "maxAntSpine" - "Distance max to anterior part of the VertebralBody (spine) in mm" double no default="10" -option "spacing" - "Intermediate resampling spacing" double no default="6" -option "ft_LR_lungs" - "RelPos LR limits lungs" double no default="0.3" -option "ft_bones" - "RelPos AP limits bones" double no default="0.6" -option "ft_inf_lungs" - "RelPos inf limits bones" double no default="0.05" -option "ft_ant_sternum" - "RelPos ant limits sternum" double no default="0.5" - -#section "Step 1 : Left/Right limits with lungs" - -section "Step 2 : Ant/Post limits with bones" -#option "antSpine" - "Distance max to anterior part of the spine in mm" double no default="10" - -#section "Step 3 : Inf limits with Lung" - -# section "Step x : threshold for removing bones and injected parts" -# option "upper" - "Upper threshold" double no default="150" -# option "lower" - "Lower threshold" double no default="-200" \ No newline at end of file diff --git a/segmentation/clitkExtractMediastinumFilter.h b/segmentation/clitkExtractMediastinumFilter.h index 61abdee..3144c95 100644 --- a/segmentation/clitkExtractMediastinumFilter.h +++ b/segmentation/clitkExtractMediastinumFilter.h @@ -19,8 +19,7 @@ #ifndef CLITKEXTRACTMEDIASTINUMFILTER_H #define CLITKEXTRACTMEDIASTINUMFILTER_H -#include "clitkFilterBase.h" -#include "clitkFilterWithAnatomicalFeatureDatabaseManagement.h" +#include "clitkStructuresExtractionFilter.h" namespace clitk { @@ -38,9 +37,7 @@ namespace clitk { template class ITK_EXPORT ExtractMediastinumFilter: - public virtual clitk::FilterBase, - public clitk::FilterWithAnatomicalFeatureDatabaseManagement, - public itk::ImageToImageFilter > + public clitk::StructuresExtractionFilter { public: @@ -69,7 +66,8 @@ namespace clitk { typedef typename MaskSliceType::PointType MaskSlicePointType; /** Standard class typedefs. */ - typedef itk::ImageToImageFilter Superclass; + // typedef itk::ImageToImageFilter Superclass; + typedef clitk::StructuresExtractionFilter Superclass; typedef ExtractMediastinumFilter Self; typedef itk::SmartPointer Pointer; typedef itk::SmartPointer ConstPointer; @@ -88,8 +86,8 @@ namespace clitk { MaskImagePixelType fgLeftLung=1, MaskImagePixelType fgRightLung=2); void SetInputBonesLabelImage(const MaskImageType * image, MaskImagePixelType bg=0); void SetInputTracheaLabelImage(const MaskImageType * image, MaskImagePixelType bg=0); - - // Output filename (for AFBD) + + // Output filename (for AFBD) itkSetMacro(OutputMediastinumFilename, std::string); itkGetConstMacro(OutputMediastinumFilename, std::string); @@ -106,9 +104,6 @@ namespace clitk { itkSetMacro(BackgroundValueBones, MaskImagePixelType); itkGetConstMacro(BackgroundValueBones, MaskImagePixelType); - itkGetConstMacro(BackgroundValue, MaskImagePixelType); - itkGetConstMacro(ForegroundValue, MaskImagePixelType); - itkSetMacro(ForegroundValueLeftLung, MaskImagePixelType); itkGetConstMacro(ForegroundValueLeftLung, MaskImagePixelType); @@ -138,10 +133,7 @@ namespace clitk { virtual void GenerateOutputInformation(); virtual void GenerateInputRequestedRegion(); virtual void GenerateData(); - - itkSetMacro(BackgroundValue, MaskImagePixelType); - itkSetMacro(ForegroundValue, MaskImagePixelType); - + MaskImagePixelType m_BackgroundValuePatient; MaskImagePixelType m_BackgroundValueLung; MaskImagePixelType m_BackgroundValueBones; @@ -149,9 +141,6 @@ namespace clitk { MaskImagePixelType m_ForegroundValueLeftLung; MaskImagePixelType m_ForegroundValueRightLung; - MaskImagePixelType m_BackgroundValue; - MaskImagePixelType m_ForegroundValue; - MaskImagePointer output; MaskImagePointer patient; MaskImagePointer lung; diff --git a/segmentation/clitkExtractMediastinumFilter.txx b/segmentation/clitkExtractMediastinumFilter.txx index 67bc305..312fd14 100644 --- a/segmentation/clitkExtractMediastinumFilter.txx +++ b/segmentation/clitkExtractMediastinumFilter.txx @@ -45,9 +45,7 @@ template clitk::ExtractMediastinumFilter:: ExtractMediastinumFilter(): - clitk::FilterBase(), - clitk::FilterWithAnatomicalFeatureDatabaseManagement(), - itk::ImageToImageFilter() + clitk::StructuresExtractionFilter() { this->SetNumberOfRequiredInputs(1); SetBackgroundValuePatient(0); @@ -55,12 +53,6 @@ ExtractMediastinumFilter(): SetBackgroundValueBones(0); SetForegroundValueLeftLung(1); SetForegroundValueRightLung(2); - SetBackgroundValue(0); - SetForegroundValue(1); - SetIntermediateSpacing(6); - SetFuzzyThreshold("LR_lungs", 0.3); - SetFuzzyThreshold("bones", 0.6); - SetFuzzyThreshold("inf_lungs", 0.05); SetDistanceMaxToAnteriorPartOfTheVertebralBody(10); SetOutputMediastinumFilename("mediastinum.mhd"); UseBonesOff(); @@ -156,7 +148,6 @@ GenerateInputRequestedRegion() // Superclass::GenerateInputRequestedRegion(); // DD("End GenerateInputRequestedRegion"); } - //-------------------------------------------------------------------- @@ -171,6 +162,7 @@ SetInput(const ImageType * image) //-------------------------------------------------------------------- + //-------------------------------------------------------------------- template void @@ -181,19 +173,21 @@ GenerateOutputInformation() { //-------------------------------------------------------------------- // Get input pointers - LoadAFDB(); + this->LoadAFDB(); ImageConstPointer input = dynamic_cast(itk::ProcessObject::GetInput(0)); MaskImagePointer patient, lung, bones, trachea; - patient = GetAFDB()->template GetImage ("Patient"); - lung = GetAFDB()->template GetImage ("Lungs"); + patient = this->GetAFDB()->template GetImage ("Patient"); + lung = this->GetAFDB()->template GetImage ("Lungs"); if (GetUseBones()) { - bones = GetAFDB()->template GetImage ("Bones"); + bones = this->GetAFDB()->template GetImage ("Bones"); } - trachea = GetAFDB()->template GetImage ("Trachea"); - + trachea = this->GetAFDB()->template GetImage ("Trachea"); + + //this->VerboseImageSizeFlagOn(); + //-------------------------------------------------------------------- // Step : Crop support (patient) to lung extend in RL - StartNewStep("Crop support like lungs along LR"); + this->StartNewStep("[Mediastinum] Crop support like lungs along LR"); typedef clitk::CropLikeImageFilter CropFilterType; typename CropFilterType::Pointer cropFilter = CropFilterType::New(); cropFilter->SetInput(patient); @@ -205,7 +199,7 @@ GenerateOutputInformation() { //-------------------------------------------------------------------- // Step : Crop support (previous) to bones extend in AP if (GetUseBones()) { - StartNewStep("Crop support like bones along AP"); + this->StartNewStep("[Mediastinum] Crop support like bones along AP"); cropFilter = CropFilterType::New(); cropFilter->SetInput(output); cropFilter->SetCropLikeImage(bones, 1);// Indicate that we only crop in Y (Ant-Post) axe @@ -216,7 +210,7 @@ GenerateOutputInformation() { //-------------------------------------------------------------------- // Step : patient minus lungs, minus bones, minus trachea - StartNewStep("Patient contours minus lungs, trachea [and bones]"); + this->StartNewStep("[Mediastinum] Patient contours minus lungs, trachea [and bones]"); typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); boolFilter->InPlaceOn(); @@ -237,7 +231,7 @@ GenerateOutputInformation() { output = boolFilter->GetOutput(); // Auto crop to gain support area - output = clitk::AutoCrop(output, GetBackgroundValue()); + output = clitk::AutoCrop(output, this->GetBackgroundValue()); this->template StopCurrentStep(output); //-------------------------------------------------------------------- @@ -246,69 +240,50 @@ GenerateOutputInformation() { // (because RelativePositionPropImageFilter only consider fg=1); // (label must be '1' because right is greater than left). (WE DO // NOT NEED TO SEPARATE ? ) - StartNewStep("Left/Right limits with lungs"); + this->StartNewStep("[Mediastinum] Left/Right limits with lungs"); // The following cannot be "inplace" because mask is the same than input ... MaskImagePointer right_lung = clitk::SetBackground(lung, lung, 2, 0, false); MaskImagePointer left_lung = clitk::SetBackground(lung, lung, 1, 0, false); - right_lung = clitk::ResizeImageLike(right_lung, output, GetBackgroundValue()); - left_lung = clitk::ResizeImageLike(left_lung, output, GetBackgroundValue()); - // writeImage(right_lung, "right.mhd"); - // writeImage(left_lung, "left.mhd"); - - typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; - typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); - relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); - relPosFilter->VerboseStepFlagOff(); - relPosFilter->WriteStepFlagOff(); - relPosFilter->SetInput(output); - relPosFilter->SetInputObject(left_lung); - relPosFilter->AddOrientationType(RelPosFilterType::AtRightTo); // warning left lung is at right ;) - relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("LR_lungs")); - relPosFilter->Update(); - output = relPosFilter->GetOutput(); - - relPosFilter = RelPosFilterType::New(); - relPosFilter->SetInput(output); - relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); - relPosFilter->VerboseStepFlagOff(); - relPosFilter->WriteStepFlagOff(); - relPosFilter->SetInput(output); - relPosFilter->SetInputObject(right_lung); - relPosFilter->AddOrientationType(RelPosFilterType::AtLeftTo); - relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("LR_lungs")); - relPosFilter->Update(); - output = relPosFilter->GetOutput(); - this->template StopCurrentStep(output); - - //-------------------------------------------------------------------- - // Step : superior limits - StartNewStep("Keep inferior to CricoidCartilag"); - // load Cricoid, get centroid, cut above (or below), lower bound - MaskImagePointer CricoidCartilag = GetAFDB()->template GetImage ("CricoidCartilag"); - MaskImagePointType p; - p[0] = p[1] = p[2] = 0.0; // to avoid warning - clitk::FindExtremaPointInAGivenDirection(CricoidCartilag, GetBackgroundValue(), 2, true, p); - output = clitk::CropImageRemoveGreaterThan(output, 2, p[2], true, GetBackgroundValue()); + left_lung = clitk::SetBackground(left_lung, left_lung, 2, 1, false); + right_lung = clitk::ResizeImageLike(right_lung, output, this->GetBackgroundValue()); + left_lung = clitk::ResizeImageLike(left_lung, output, this->GetBackgroundValue()); + this->GetAFDB()->template SetImage("RightLung", "seg/RightLung.mhd", + right_lung, true); + this->GetAFDB()->template SetImage("LeftLung", "seg/LeftLung.mhd", + left_lung, true); + this->GetAFDB()->Write(); this->template StopCurrentStep(output); //-------------------------------------------------------------------- // Step : AP limits from bones // Separate the bones in the ant-post middle - MaskImageConstPointer bones_ant; - MaskImageConstPointer bones_post; - MaskImagePointType middle_AntPost__position; + MaskImagePointer bones_ant; + MaskImagePointer bones_post; + MaskImagePointType middle_AntPost_position; if (GetUseBones()) { - StartNewStep("Ant/Post limits with bones"); - middle_AntPost__position[0] = middle_AntPost__position[2] = 0; - middle_AntPost__position[1] = bones->GetOrigin()[1]+(bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1])/2.0; + this->StartNewStep("[Mediastinum] Ant/Post limits with bones"); + + // To define ant and post part of the bones with a single horizontal line MaskImageIndexType index_bones_middle; - bones->TransformPhysicalPointToIndex(middle_AntPost__position, index_bones_middle); + + // Method1: cut in the middle (not optimal) + /* + middle_AntPost_position[0] = middle_AntPost_position[2] = 0; + middle_AntPost_position[1] = bones->GetOrigin()[1]+ + (bones->GetLargestPossibleRegion().GetSize()[1]*bones->GetSpacing()[1])/2.0; + DD(middle_AntPost_position); + bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle); + */ + // Method2: Use VertebralBody, take most ant point + MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage("VertebralBody"); + FindExtremaPointInAGivenDirection(VertebralBody, this->GetBackgroundValue(), + 1, true, middle_AntPost_position); + bones->TransformPhysicalPointToIndex(middle_AntPost_position, index_bones_middle); + // Split bone image first into two parts (ant and post), and crop // lateraly to get vertebral typedef itk::RegionOfInterestImageFilter ROIFilterType; @@ -319,8 +294,8 @@ GenerateOutputInformation() { MaskImageIndexType index = region.GetIndex(); // ANT part // crop LR to keep 1/4 center part - index[0] = size[0]/4+size[0]/8; - size[0] = size[0]/4; + // index[0] = size[0]/4+size[0]/8; + //size[0] = size[0]/4; // crop AP to keep first (ant) part size[1] = index_bones_middle[1]; //size[1]/2.0; region.SetSize(size); @@ -330,7 +305,7 @@ GenerateOutputInformation() { roiFilter->ReleaseDataFlagOff(); roiFilter->Update(); bones_ant = roiFilter->GetOutput(); - // writeImage(bones_ant, "b_ant.mhd"); + // writeImage(bones_ant, "b_ant.mhd"); // POST part roiFilter = ROIFilterType::New(); index[1] = bones->GetLargestPossibleRegion().GetIndex()[1] + size[1]-1; @@ -342,64 +317,60 @@ GenerateOutputInformation() { roiFilter->ReleaseDataFlagOff(); roiFilter->Update(); bones_post = roiFilter->GetOutput(); - // writeImage(bones_post, "b_post.mhd"); - - // Go ! - relPosFilter->SetCurrentStepNumber(0); - relPosFilter->ResetPipeline();// = RelPosFilterType::New(); - relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); - relPosFilter->VerboseStepFlagOff(); - relPosFilter->WriteStepFlagOff(); - relPosFilter->SetInput(output); - relPosFilter->SetInputObject(bones_post); - relPosFilter->AddOrientationType(RelPosFilterType::AntTo); - relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("bones")); - relPosFilter->Update(); - output = relPosFilter->GetOutput(); - // writeImage(output, "post.mhd"); - - relPosFilter->SetInput(relPosFilter->GetOutput()); - relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); - relPosFilter->VerboseStepFlagOff(); - relPosFilter->WriteStepFlagOff(); - relPosFilter->SetInput(output); - relPosFilter->SetInputObject(bones_ant); - relPosFilter->AddOrientationType(RelPosFilterType::PostTo); - relPosFilter->SetIntermediateSpacing(GetIntermediateSpacing()); - relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("bones")); - relPosFilter->Update(); - output = relPosFilter->GetOutput(); + // writeImage(bones_post, "b_post.mhd"); + + // Now, insert this image in the AFDB + this->GetAFDB()->template SetImage("Bones_Post", "seg/Bones_Post.mhd", + bones_post, true); + this->GetAFDB()->template SetImage("Bones_Ant", "seg/Bones_Ant.mhd", + bones_ant, true); + this->GetAFDB()->Write(); + this->template StopCurrentStep(output); } + //-------------------------------------------------------------------- + // Remove VertebralBody part + this->StartNewStep("[Mediastinum] Remove VertebralBody"); + MaskImagePointer VertebralBody = this->GetAFDB()->template GetImage("VertebralBody"); + clitk::AndNot(output, VertebralBody, this->GetBackgroundValue()); + this->template StopCurrentStep(output); + + //-------------------------------------------------------------------- + // Generic RelativePosition processes + output = this->ApplyRelativePositionList("Mediastinum", output); + + //-------------------------------------------------------------------- + // Step : SI limits It is better to do this limit *AFTER* the + // RelativePosition to avoid some issue due to superior boundaries. + this->StartNewStep("[Mediastinum] Keep inferior to CricoidCartilag"); + // load Cricoid, get centroid, cut above (or below), lower bound + MaskImagePointer CricoidCartilag = this->GetAFDB()->template GetImage ("CricoidCartilag"); + MaskImagePointType p; + p[0] = p[1] = p[2] = 0.0; // to avoid warning + clitk::FindExtremaPointInAGivenDirection(CricoidCartilag, + this->GetBackgroundValue(), 2, true, p); + output = clitk::CropImageRemoveGreaterThan(output, 2, p[2], true, this->GetBackgroundValue()); + this->template StopCurrentStep(output); + //-------------------------------------------------------------------- // Step: Get CCL - StartNewStep("Keep main connected component"); - output = clitk::Labelize(output, GetBackgroundValue(), false, 500); + this->StartNewStep("[Mediastinum] Keep main connected component"); + output = clitk::Labelize(output, this->GetBackgroundValue(), false, 500); // output = RemoveLabels(output, BG, param->GetLabelsToRemove()); - output = clitk::KeepLabels(output, GetBackgroundValue(), - GetForegroundValue(), 1, 1, 0); + output = clitk::KeepLabels(output, this->GetBackgroundValue(), + this->GetForegroundValue(), 1, 1, 0); this->template StopCurrentStep(output); //-------------------------------------------------------------------- // Step: Remove post part from VertebralBody - StartNewStep("Remove post part according to VertebralBody"); + this->StartNewStep("[Mediastinum] Remove post part according to VertebralBody"); RemovePostPartOfVertebralBody(); this->template StopCurrentStep(output); - //-------------------------------------------------------------------- - // Step: Remove ant part according to Sternum - StartNewStep("Remove ant part according to Sternum"); - MaskImagePointer Sternum = GetAFDB()->template GetImage ("Sternum"); - output = clitk::SliceBySliceRelativePosition(output, Sternum, 2, - GetFuzzyThreshold("ant_sternum"), - "PostTo", false, 3, true, false); - this->template StopCurrentStep(output); - //-------------------------------------------------------------------- // Step: Slice by Slice CCL - StartNewStep("Slice by Slice keep only one component"); + this->StartNewStep("[Mediastinum] Slice by Slice keep only one component"); typedef clitk::ExtractSliceFilter ExtractSliceFilterType; // typename ExtractSliceFilterType::Pointer ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New(); @@ -426,8 +397,8 @@ GenerateOutputInformation() { //-------------------------------------------------------------------- // Step 10 : AutoCrop - StartNewStep("AutoCrop"); - output = clitk::AutoCrop(output, GetBackgroundValue()); + this->StartNewStep("[Mediastinum] AutoCrop"); + output = clitk::AutoCrop(output, this->GetBackgroundValue()); this->template StopCurrentStep(output); // End, set the real size @@ -447,8 +418,8 @@ GenerateData() { this->GraftOutput(output); // Store image filenames into AFDB - GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename()); - WriteAFDB(); + this->GetAFDB()->SetImageFilename("Mediastinum", this->GetOutputMediastinumFilename()); + this->WriteAFDB(); } //-------------------------------------------------------------------- @@ -475,7 +446,7 @@ RemovePostPartOfVertebralBody() // Get VertebralBody mask image MaskImagePointer VertebralBody = - GetAFDB()->template GetImage ("VertebralBody"); + this->GetAFDB()->template GetImage ("VertebralBody"); // Consider vertebral body slice by slice std::vector vertebralSlices; @@ -486,7 +457,7 @@ RemovePostPartOfVertebralBody() for(uint i=0; i(vertebralSlices[i], - GetBackgroundValue(), + this->GetBackgroundValue(), 1, true, p); if (found) { vertebralAntPositionBySlice[i] = p; @@ -540,7 +511,7 @@ RemovePostPartOfVertebralBody() itk::ImageRegionIterator it(output, region); it.GoToBegin(); while (!it.IsAtEnd()) { - it.Set(GetBackgroundValue()); + it.Set(this->GetBackgroundValue()); ++it; } } diff --git a/segmentation/clitkExtractMediastinumGenericFilter.txx b/segmentation/clitkExtractMediastinumGenericFilter.txx index 0670ed7..679fb36 100644 --- a/segmentation/clitkExtractMediastinumGenericFilter.txx +++ b/segmentation/clitkExtractMediastinumGenericFilter.txx @@ -78,11 +78,10 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) f->SetDistanceMaxToAnteriorPartOfTheVertebralBody(mArgsInfo.maxAntSpine_arg); f->SetUseBones(mArgsInfo.useBones_flag); - f->SetIntermediateSpacing(mArgsInfo.spacing_arg); - f->SetFuzzyThreshold("LR_lungs", mArgsInfo.ft_LR_lungs_arg); - f->SetFuzzyThreshold("bones", mArgsInfo.ft_bones_arg); - f->SetFuzzyThreshold("inf_lungs", mArgsInfo.ft_inf_lungs_arg); - f->SetFuzzyThreshold("ant_sternum", mArgsInfo.ft_ant_sternum_arg); + // Set RelativePositionList filenames + for(uint i=0; iAddRelativePositionListFilename(mArgsInfo.relpos_arg[i]); + } } //-------------------------------------------------------------------- diff --git a/segmentation/clitkMorphoMath.ggo b/segmentation/clitkMorphoMath.ggo index c754f5e..7959bc9 100644 --- a/segmentation/clitkMorphoMath.ggo +++ b/segmentation/clitkMorphoMath.ggo @@ -22,7 +22,7 @@ option "fg" - "Foreground value" float no default="1" option "bg" - "Background value (0,1,3: filling value)" float no default="0" option "bound" b "0-1: Set borders to foreground/ 2:safe borders " flag off option "radius" r "Use binary ball element with given radius" int no multiple default="1" -option "radiusInMM" m "Use binary ball element with given radius in mm (rounded to nearest voxel value)" double no multiple default="1" +option "radiusInMM" m "Use binary ball element with given radius in mm (rounded to nearest voxel value), you can give one radius by dimension" double no multiple default="1" option "extend" - "Extend the image size according to radius" flag off diff --git a/segmentation/clitkRelativePositionList.h b/segmentation/clitkRelativePositionList.h new file mode 100644 index 0000000..20c6b77 --- /dev/null +++ b/segmentation/clitkRelativePositionList.h @@ -0,0 +1,109 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +#ifndef CLITKRELATIVEPOSITIONLIST_H +#define CLITKRELATIVEPOSITIONLIST_H + +// clitk +#include "clitkSegmentationUtils.h" +#include "clitkFilterWithAnatomicalFeatureDatabaseManagement.h" +#include "clitkRelativePosition_ggo.h" + +namespace clitk { + + /*-------------------------------------------------------------------- + Manage a list of RelativePosition operations, to be performed on + the same input image, with different objects and parameters. + ------------------------------------------------------------------*/ + + template + class ITK_EXPORT RelativePositionList: + public virtual clitk::FilterBase, + public clitk::FilterWithAnatomicalFeatureDatabaseManagement, + public itk::ImageToImageFilter + { + + public: + /** Standard class typedefs. */ + typedef itk::ImageToImageFilter Superclass; + typedef RelativePositionList Self; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(RelativePositionList, ImageToImageFilter); + + /** Some convenient typedefs. */ + typedef TImageType ImageType; + typedef typename ImageType::ConstPointer ImageConstPointer; + typedef typename ImageType::Pointer ImagePointer; + typedef typename ImageType::RegionType ImageRegionType; + typedef typename ImageType::PixelType ImagePixelType; + typedef typename ImageType::SizeType ImageSizeType; + typedef typename ImageType::IndexType ImageIndexType; + typedef typename ImageType::PointType ImagePointType; + typedef struct args_info_clitkRelativePosition ArgsInfoType; + typedef SliceBySliceRelativePositionFilter SliceRelPosFilterType; + typedef AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; + + /** ImageDimension constants */ + itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension); + FILTERBASE_INIT; + + itkGetConstMacro(BackgroundValue, ImagePixelType); + itkGetConstMacro(ForegroundValue, ImagePixelType); + itkSetMacro(BackgroundValue, ImagePixelType); + itkSetMacro(ForegroundValue, ImagePixelType); + + itkSetMacro(InputName, std::string); + itkGetConstMacro(InputName, std::string); + + void Read(std::string filename); + void SetFilterOptions(typename RelPosFilterType::Pointer filter, ArgsInfoType & options); + + protected: + RelativePositionList(); + virtual ~RelativePositionList() {} + + private: + RelativePositionList(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + void GenerateInputRequestedRegion(); + void GenerateOutputInformation(); + void GenerateData(); + + std::string m_InputName; + ImagePixelType m_BackgroundValue; + ImagePixelType m_ForegroundValue; + typename SliceRelPosFilterType::Pointer mFilter; + std::vector mArgsInfoList; + ImagePointer m_working_input; + + }; // end class + //-------------------------------------------------------------------- + +} // end namespace clitk +//-------------------------------------------------------------------- + +#include "clitkRelativePositionList.txx" + +#endif diff --git a/segmentation/clitkRelativePositionList.txx b/segmentation/clitkRelativePositionList.txx new file mode 100644 index 0000000..3ce5559 --- /dev/null +++ b/segmentation/clitkRelativePositionList.txx @@ -0,0 +1,233 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + + +//-------------------------------------------------------------------- +template +clitk::RelativePositionList:: +RelativePositionList() { +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::RelativePositionList:: +Read(std::string filename) { + /* + The goal here is to read a text file that contains options for the + RelativePosition filter. The text file contains options in the + same form of the config file of the clitkRelativePosition tool. In + this text file, each time a "object" option is found, a new set of + options is considered. + */ + + // Open the file + std::ifstream is; + openFileForReading(is, filename); + + // Read input -> the structure name that will be used as input + // (initial support) + skipComment(is); + std::string s; + is >> s; + if (s != "input") { + FATAL("while reading RelativePositionList file. This file must start with 'input = '"); + exit(0); + } + is >> s; + if (s != "=") { + FATAL("while reading RelativePositionList file. This file must start with 'input = '"); + exit(0); + } + std::string name; + is >> name; + skipComment(is); + + // Create a temporary filename + char n[] = "tmp_clitkRelativePosition_XXXXXX"; + mkstemp(n); // tmpnam(n); + std::string tmpfilename(n); + + // Loop on the file text ; Every time we see a "object" token, we + // split the file. + while (is) { + bool stop=false; + std::ostringstream ss; + // first part of ss is the last 'object = ' found, stored in s + ss << s << std::endl; + while (!stop) { + skipComment(is); + if (!is) stop = true; + else { + std::getline(is, s); + // DD(s); + if (s.find("object") != std::string::npos) stop=true; + else ss << s << std::endl; + if (!is) stop = true; + } + } + std::string text = ss.str(); + if (text.size() > 10) { // if too small, it is not a list of option + std::ofstream os; + openFileForWriting(os, tmpfilename); + os << text; + os.close(); + + // Create a struct to store options. I use two step to allow to + // fill the args values with de default and then check + // automatically the options. + ArgsInfoType args_info; + std::vector writable(tmpfilename.size() + 1); + std::copy(tmpfilename.begin(), tmpfilename.end(), writable.begin()); + char ** argv; + cmdline_parser_clitkRelativePosition2(0, argv, &args_info, 1, 1, 0); + args_info.input_given = 1; + args_info.input_arg = new char[1]; + args_info.output_given = 1; + args_info.output_arg = new char[1]; + cmdline_parser_clitkRelativePosition_configfile(&writable[0], &args_info, 0, 0, 1); + + // Store the args + mArgsInfoList.push_back(args_info); + + // Delete the temporary file + std::remove(tmpfilename.c_str()); + } + } + + // Set the main input name + SetInputName(name); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::RelativePositionList:: +GenerateInputRequestedRegion() +{ + // Call default + itk::ImageToImageFilter::GenerateInputRequestedRegion(); + // Get input pointers and set requested region to common region + ImagePointer input1 = dynamic_cast(itk::ProcessObject::GetInput(0)); + input1->SetRequestedRegion(input1->GetLargestPossibleRegion()); +} +//-------------------------------------------------------------------- + + + +//-------------------------------------------------------------------- +template +void +clitk::RelativePositionList:: +GenerateOutputInformation() { + + // Get input + m_working_input = dynamic_cast(itk::ProcessObject::GetInput(0)); + + // Loop on RelativePositionList of operations + std::string s = GetInputName(); + for(uint i=0; iSetDirection(2); + // Set SbS specific options + f->SetUniqueConnectedComponentBySliceFlag(mArgsInfoList[i].uniqueCCL_flag); + f->SetObjectCCLSelectionFlag(mArgsInfoList[i].uniqueObjectCCL_flag); + f->IgnoreEmptySliceObjectFlagOn(); + //f->SetObjectCCLSelectionDimension(0); + //f->SetObjectCCLSelectionDirection(-1); + //f->SetAutoCropFlag(false); + // Print if needed + if (mArgsInfoList[i].verboseOptions_flag) f->PrintOptions(); + } + else { + relPosFilter = clitk::AddRelativePositionConstraintToLabelImageFilter::New(); + SetFilterOptions(relPosFilter, mArgsInfoList[i]); + // Print if needed + if (mArgsInfoList[i].verboseOptions_flag) relPosFilter->PrintOptions(); + } + + // Set input + relPosFilter->SetInput(m_working_input); + + // Run the filter + relPosFilter->Update(); + m_working_input = relPosFilter->GetOutput(); + StopCurrentStep(m_working_input); + } +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::RelativePositionList:: +GenerateData() +{ + // Final Step -> set output + this->GraftOutput(m_working_input); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::RelativePositionList:: +SetFilterOptions(typename RelPosFilterType::Pointer filter, ArgsInfoType & options) { + + if (options.orientation_given != 1) { + DD("ERRROR DEBUG TODO no more than 1 orientation yet"); + exit(0); + } + + ImagePointer object = GetAFDB()->template GetImage(options.object_arg); + filter->SetInputObject(object); + filter->WriteStepFlagOff(); + filter->SetVerboseImageSizeFlag(GetVerboseImageSizeFlag()); + filter->SetFuzzyThreshold(options.threshold_arg); + filter->SetInverseOrientationFlag(options.inverse_flag); // MUST BE BEFORE AddOrientationTypeString + for(uint i=0; iAddOrientationTypeString(options.orientation_arg[i]); + filter->SetIntermediateSpacing(options.spacing_arg); + if (options.spacing_arg == -1) filter->IntermediateSpacingFlagOff(); + else filter->IntermediateSpacingFlagOn(); + filter->SetVerboseStepFlag(options.verboseStep_flag); + filter->SetAutoCropFlag(!options.noAutoCrop_flag); +} +//-------------------------------------------------------------------- diff --git a/segmentation/clitkStructuresExtractionFilter.h b/segmentation/clitkStructuresExtractionFilter.h new file mode 100644 index 0000000..3e7d303 --- /dev/null +++ b/segmentation/clitkStructuresExtractionFilter.h @@ -0,0 +1,120 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ===========================================================================**/ + +#ifndef CLITKSTRUCTURESEXTRACTIONFILTER_H +#define CLITKSTRUCTURESEXTRACTIONFILTER_H + +// clitk +#include "clitkFilterWithAnatomicalFeatureDatabaseManagement.h" +#include "clitkFilterBase.h" +#include "clitkRelativePositionList.h" + +namespace clitk { + + //-------------------------------------------------------------------- + /* + - Convenient class to add some capabilities to a filter : + FilterBase, FilterWithAnatomicalFeatureDatabaseManagement and + RelativePositionList + */ + //-------------------------------------------------------------------- + template + class StructuresExtractionFilter: + public virtual FilterBase, + public clitk::FilterWithAnatomicalFeatureDatabaseManagement, + public itk::ImageToImageFilter > + { + public: + // Standard class typedefs + typedef StructuresExtractionFilter Self; + typedef itk::ImageToImageFilter > Superclass; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + // Run-time type information (and related methods) + itkTypeMacro(StructuresExtractionFilter, ImageToImageFilter); + + /** Some convenient typedefs. */ + typedef TImageType ImageType; + typedef typename ImageType::ConstPointer ImageConstPointer; + typedef typename ImageType::Pointer ImagePointer; + typedef typename ImageType::RegionType ImageRegionType; + typedef typename ImageType::PixelType ImagePixelType; + typedef typename ImageType::SizeType ImageSizeType; + typedef typename ImageType::IndexType ImageIndexType; + typedef typename ImageType::PointType ImagePointType; + + typedef uchar MaskImagePixelType; + typedef itk::Image MaskImageType; + typedef typename MaskImageType::Pointer MaskImagePointer; + typedef typename MaskImageType::RegionType MaskImageRegionType; + typedef typename MaskImageType::SizeType MaskImageSizeType; + typedef typename MaskImageType::IndexType MaskImageIndexType; + typedef typename MaskImageType::PointType MaskImagePointType; + + typedef itk::Image MaskSliceType; + typedef typename MaskSliceType::Pointer MaskSlicePointer; + typedef typename MaskSliceType::PointType MaskSlicePointType; + typedef typename MaskSliceType::RegionType MaskSliceRegionType; + typedef typename MaskSliceType::SizeType MaskSliceSizeType; + typedef typename MaskSliceType::IndexType MaskSliceIndexType; + + /** ImageDimension constants */ + itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension); + FILTERBASE_INIT; + + itkGetConstMacro(BackgroundValue, MaskImagePixelType); + itkGetConstMacro(ForegroundValue, MaskImagePixelType); + itkSetMacro(BackgroundValue, MaskImagePixelType); + itkSetMacro(ForegroundValue, MaskImagePixelType); + + // RelativePositionList management + void AddRelativePositionListFilename(std::string s); + MaskImagePointer ApplyRelativePositionList(std::string name, MaskImageType * input); + + protected: + StructuresExtractionFilter(); + virtual ~StructuresExtractionFilter() {} + + virtual void GenerateData() {} + + MaskImagePixelType m_BackgroundValue; + MaskImagePixelType m_ForegroundValue; + + // RelativePositionList + std::vector mListOfRelativePositionListFilename; + typedef clitk::RelativePositionList RelPosListType; + typedef typename RelPosListType::Pointer RelPosListPointer; + std::map mMapOfRelativePositionList; + + private: + StructuresExtractionFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + }; // end class + //-------------------------------------------------------------------- + +} // end namespace clitk +//-------------------------------------------------------------------- + +#include "clitkStructuresExtractionFilter.txx" + +#endif // CLITKSTRUCTURESEXTRACTIONFILTER_H diff --git a/segmentation/clitkStructuresExtractionFilter.txx b/segmentation/clitkStructuresExtractionFilter.txx new file mode 100644 index 0000000..1562bea --- /dev/null +++ b/segmentation/clitkStructuresExtractionFilter.txx @@ -0,0 +1,81 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ===========================================================================**/ + +// clitk +#include "clitkStructuresExtractionFilter.h" + +//-------------------------------------------------------------------- +template +clitk::StructuresExtractionFilter:: +StructuresExtractionFilter(): + clitk::FilterBase(), + clitk::FilterWithAnatomicalFeatureDatabaseManagement(), + itk::ImageToImageFilter() +{ + SetBackgroundValue(0); + SetForegroundValue(1); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::StructuresExtractionFilter:: +AddRelativePositionListFilename(std::string s) { + mListOfRelativePositionListFilename.push_back(s); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +typename clitk::StructuresExtractionFilter::MaskImagePointer +clitk::StructuresExtractionFilter:: +ApplyRelativePositionList(std::string name, MaskImageType * input) +{ + // Create all RelativePositionList + for(int i=0; iSetAFDB(GetAFDB()); + rpl->Read(mListOfRelativePositionListFilename[i]); + std::string s = rpl->GetInputName(); + mMapOfRelativePositionList[s] = rpl; + } + + RelPosListPointer relpos; + if (mMapOfRelativePositionList.find(name) == mMapOfRelativePositionList.end()) { + std::cerr << "Warning: I do not find '" << name << "' in the RelativePositionList." << std::endl; + //DD("Not find !"); // do nothing + } + else { + relpos = mMapOfRelativePositionList[name]; + relpos->SetVerboseStepFlag(GetVerboseStepFlag()); + relpos->SetCurrentStepBaseId(GetCurrentStepBaseId()); + relpos->SetCurrentStepId(GetCurrentStepId()); + relpos->SetCurrentStepNumber(GetCurrentStepNumber()); + relpos->SetWriteStepFlag(GetWriteStepFlag()); + relpos->SetInput(input); + relpos->Update(); + input = relpos->GetOutput(); + SetCurrentStepNumber(relpos->GetCurrentStepNumber()); + } + return input; +} +//-------------------------------------------------------------------- + diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index f5a62b4..d785a49 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -157,10 +157,15 @@ IF (CLITK_BUILD_TOOLS) TARGET_LINK_LIBRARIES(clitkAutoCrop clitkCommon ${ITK_LIBRARIES} ) SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkAutoCrop) - WRAP_GGO(clitkDicomRTStruct2BinaryImage_GGO_C clitkDicomRTStruct2BinaryImage.ggo) - ADD_EXECUTABLE(clitkDicomRTStruct2BinaryImage clitkDicomRTStruct2BinaryImage.cxx ${clitkDicomRTStruct2BinaryImage_GGO_C}) - TARGET_LINK_LIBRARIES(clitkDicomRTStruct2BinaryImage clitkDicomRTStruct clitkCommon ${ITK_LIBRARIES} ) - SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDicomRTStruct2BinaryImage) + WRAP_GGO(clitkDicomRTStruct2Image_GGO_C clitkDicomRTStruct2Image.ggo) + ADD_EXECUTABLE(clitkDicomRTStruct2Image clitkDicomRTStruct2Image.cxx ${clitkDicomRTStruct2Image_GGO_C}) + TARGET_LINK_LIBRARIES(clitkDicomRTStruct2Image clitkDicomRTStruct clitkCommon ${ITK_LIBRARIES} ) + SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDicomRTStruct2Image) + + # WRAP_GGO(clitkImage2DicomRTStruct_GGO_C clitkImage2DicomRTStruct.ggo) + # ADD_EXECUTABLE(clitkImage2DicomRTStruct clitkImage2DicomRTStruct.cxx ${clitkImage2DicomRTStruct_GGO_C}) + # TARGET_LINK_LIBRARIES(clitkImage2DicomRTStruct clitkDicomRTStruct clitkCommon ${ITK_LIBRARIES} ) + # SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkImage2DicomRTStruct) WRAP_GGO(clitkImageLog_GGO_C clitkImageLog.ggo) ADD_EXECUTABLE(clitkImageLog clitkImageLog.cxx ${clitkImageLog_GGO_C}) diff --git a/tools/clitkDicomRTStruct2BinaryImage.cxx b/tools/clitkDicomRTStruct2Image.cxx similarity index 65% rename from tools/clitkDicomRTStruct2BinaryImage.cxx rename to tools/clitkDicomRTStruct2Image.cxx index 97eff97..cc23126 100644 --- a/tools/clitkDicomRTStruct2BinaryImage.cxx +++ b/tools/clitkDicomRTStruct2Image.cxx @@ -17,15 +17,15 @@ =========================================================================*/ -#include "clitkDicomRT_ROI_ConvertToImageFilter.h" +#include "clitkDicomRTStruct2ImageFilter.h" #include "clitkDicomRT_StructureSet.h" -#include "clitkDicomRTStruct2BinaryImage_ggo.h" +#include "clitkDicomRTStruct2Image_ggo.h" //-------------------------------------------------------------------- int main(int argc, char * argv[]) { // Init command line - GGO(clitkDicomRTStruct2BinaryImage, args_info); + GGO(clitkDicomRTStruct2Image, args_info); // Read and display information clitk::DicomRT_StructureSet::Pointer s = clitk::DicomRT_StructureSet::New(); @@ -35,17 +35,41 @@ int main(int argc, char * argv[]) { } // New filter to convert to binary image - clitk::DicomRT_ROI_ConvertToImageFilter filter; + clitk::DicomRTStruct2ImageFilter filter; filter.SetCropMaskEnabled(args_info.crop_flag); filter.SetImageFilename(args_info.image_arg); // Used to get spacing + origin if (args_info.roi_arg != -1) { - filter.SetROI(s->GetROI(args_info.roi_arg)); + filter.SetROI(s->GetROIFromROINumber(args_info.roi_arg)); filter.SetOutputImageFilename(args_info.output_arg); filter.Update(); } else { + clitk::DicomRT_StructureSet::ROIConstIteratorType iter; + for(iter = s->GetROIs().begin(); iter != s->GetROIs().end(); iter++) { + clitk::DicomRT_ROI::Pointer roi = iter->second; + // Create the filter + clitk::DicomRTStruct2ImageFilter filter; + filter.SetCropMaskEnabled(args_info.crop_flag); + filter.SetImageFilename(args_info.image_arg); // Used to get spacing + origin + std::string name = roi->GetName(); + int num = roi->GetROINumber(); + filter.SetROI(roi); + name.erase(remove_if(name.begin(), name.end(), isspace), name.end()); + std::string n = std::string(args_info.output_arg).append + (clitk::toString(num)).append + ("_").append + (name).append + (".mhd"); + if (args_info.verbose_flag) { + std::cout << num << " " << roi->GetName() << " num=" << num << " : " << n << std::endl; + } + filter.SetOutputImageFilename(n); + filter.Update(); + } + + /* for(unsigned int i=0; iGetListOfROI().size(); i++) { - clitk::DicomRT_ROI_ConvertToImageFilter filter; + clitk::DicomRTStruct2ImageFilter filter; filter.SetCropMaskEnabled(args_info.crop_flag); filter.SetImageFilename(args_info.image_arg); // Used to get spacing + origin std::string name = s->GetListOfROI()[i]->GetName(); @@ -62,7 +86,7 @@ int main(int argc, char * argv[]) { } filter.SetOutputImageFilename(n); filter.Update(); - } + }*/ } // This is the end my friend diff --git a/tools/clitkDicomRTStruct2BinaryImage.ggo b/tools/clitkDicomRTStruct2Image.ggo similarity index 100% rename from tools/clitkDicomRTStruct2BinaryImage.ggo rename to tools/clitkDicomRTStruct2Image.ggo diff --git a/tools/clitkImage2DicomRTStruct.cxx b/tools/clitkImage2DicomRTStruct.cxx new file mode 100644 index 0000000..d71491b --- /dev/null +++ b/tools/clitkImage2DicomRTStruct.cxx @@ -0,0 +1,47 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + Main authors : XX XX XX + + Authors belongs to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + - BSD http://www.opensource.org/licenses/bsd-license.php + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + + =========================================================================*/ + +#include "clitkImage2DicomRTStructFilter.h" +#include "clitkDicomRT_StructureSet.h" +#include "clitkImage2DicomRTStruct_ggo.h" + +//-------------------------------------------------------------------- +int main(int argc, char * argv[]) { + + // Init command line + GGO(clitkImage2DicomRTStruct, args_info); + + // Read initial 3D image + typedef float PixelType; + typedef itk::Image ImageType; + ImageType::Pointer input = clitk::readImage(args_info.input_arg, true); + + // Create a filter to convert image into dicomRTStruct + clitk::Image2DicomRTStructFilter filter; + filter.SetInput(input); + filter.Update(); + + // Write result + clitk::DicomRT_StructureSet::Pointer s = filter.GetDicomRTStruct(); + // s->Write(args_info.output_arg); + + // This is the end my friend + return 0; +} +//-------------------------------------------------------------------- diff --git a/tools/clitkImage2DicomRTStruct.ggo b/tools/clitkImage2DicomRTStruct.ggo new file mode 100644 index 0000000..68334b1 --- /dev/null +++ b/tools/clitkImage2DicomRTStruct.ggo @@ -0,0 +1,20 @@ +# file clitkImage2DicomRTStruct.ggo +package "clitk" +version "Convert (binary) image to DICOM RT Structure Set (contours)" + +option "config" - "Config file" string no +option "verbose" v "Verbose" flag off + +option "input" i "Input image file (binary image" string yes +option "output" o "Output DicomRT filename" string yes + + +# option "image" j "Used to read image info (spacing, origin)" string yes +# option "roi" r "ROI to binarize (if -1 = all roi)" int no default="-1" + +# option "crop" c "Crop binary mask" flag off + +#option "roi" r "ROI to print (ID)" int no +#option "contour" c "contour to print (ID)" int no +#option "offset" o "to display points as image offsets" flag off + diff --git a/tools/clitkRelativePosition.ggo b/tools/clitkRelativePosition.ggo index 2192012..bc3f511 100644 --- a/tools/clitkRelativePosition.ggo +++ b/tools/clitkRelativePosition.ggo @@ -20,7 +20,7 @@ section "Main options" option "orientation" r "L R A P I S (LeftRightAntPostInfSup)" string no multiple default="L" option "angle1" a "Angle 1 (deg)" double no default="0" option "angle2" b "Angle 2 (deg)" double no default="0" -option "spacing" s "Resample before (faster)" double no +option "spacing" s "Resample before (faster) (-1 if not resampling" double default = "-1" no option "threshold" t "Fuzzy threshold" double no default="0.6" option "inverse" n "Not flag : inverse of the orientation" flag off option "doNotRemoveObject" - "if flag is on, do not remove the object" flag off @@ -31,5 +31,6 @@ section "Slice by slice processing" option "sliceBySlice" - "Slice by slice relative position" flag off option "direction" d "If SbS, indicate the slice direction" int no default="2" option "uniqueCCL" u "Keep only one CC in each slice" flag off +option "uniqueObjectCCL" - "Keep only one CC in each object slice" flag off diff --git a/tools/clitkRelativePositionGenericFilter.txx b/tools/clitkRelativePositionGenericFilter.txx index 6e4814f..66c4be4 100644 --- a/tools/clitkRelativePositionGenericFilter.txx +++ b/tools/clitkRelativePositionGenericFilter.txx @@ -77,12 +77,16 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) f->IntermediateSpacingFlagOn(); f->SetIntermediateSpacing(mArgsInfo.spacing_arg); } + else { + f->IntermediateSpacingFlagOff(); + } f->SetFuzzyThreshold(mArgsInfo.threshold_arg); f->SetRemoveObjectFlag(!mArgsInfo.doNotRemoveObject_flag); f->SetAutoCropFlag(!mArgsInfo.noAutoCrop_flag); f->SetCombineWithOrFlag(mArgsInfo.combineWithOr_flag); f->SetInverseOrientationFlag(mArgsInfo.inverse_flag); + } //-------------------------------------------------------------------- @@ -112,7 +116,13 @@ UpdateWithInputImageType() // Set options only for SliceBySliceRelativePositionFilter filter->SetDirection(mArgsInfo.direction_arg); - filter->SetUniqueConnectedComponentBySlice(mArgsInfo.uniqueCCL_flag); + filter->SetUniqueConnectedComponentBySliceFlag(mArgsInfo.uniqueCCL_flag); + if (mArgsInfo.uniqueObjectCCL_flag) { + filter->UseTheLargestObjectCCLFlagOn(); + } + else { + filter->UseTheLargestObjectCCLFlagOff(); + } // Go ! filter->Update(); diff --git a/tools/clitkSetBackground.cxx b/tools/clitkSetBackground.cxx index 04b814b..a44b4eb 100644 --- a/tools/clitkSetBackground.cxx +++ b/tools/clitkSetBackground.cxx @@ -16,16 +16,6 @@ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html ===========================================================================**/ -/* ================================================= - * @file clitkSetBackground.cxx - * @author - * @date - * - * @brief - * - ===================================================*/ - - // clitk #include "clitkSetBackground_ggo.h" #include "clitkIO.h" diff --git a/tools/clitkSetBackgroundGenericFilter.txx b/tools/clitkSetBackgroundGenericFilter.txx index c6bbe18..8a825ad 100644 --- a/tools/clitkSetBackgroundGenericFilter.txx +++ b/tools/clitkSetBackgroundGenericFilter.txx @@ -18,16 +18,6 @@ #ifndef clitkSetBackgroundGenericFilter_txx #define clitkSetBackgroundGenericFilter_txx -/* ================================================= - * @file clitkSetBackgroundGenericFilter.txx - * @author - * @date - * - * @brief - * - ===================================================*/ - - namespace clitk { diff --git a/vv/vvStructureSetActor.cxx b/vv/vvStructureSetActor.cxx index 82ef95b..7fcfb01 100644 --- a/vv/vvStructureSetActor.cxx +++ b/vv/vvStructureSetActor.cxx @@ -83,7 +83,7 @@ vvROIActor * vvStructureSetActor::GetROIActor(int n) void vvStructureSetActor::CreateNewROIActor(int n, bool modeBG) { // Check - clitk::DicomRT_ROI * roi = mStructureSet->GetROI(n); + clitk::DicomRT_ROI * roi = mStructureSet->GetROIFromROINumber(n); if (roi == NULL) { std::cerr << "Error. No ROI number " << n << std::endl; exit(0); diff --git a/vv/vvToolStructureSetManager.cxx b/vv/vvToolStructureSetManager.cxx index b424359..f625e4c 100644 --- a/vv/vvToolStructureSetManager.cxx +++ b/vv/vvToolStructureSetManager.cxx @@ -202,11 +202,19 @@ void vvToolStructureSetManager::AddRoiInTreeWidget(clitk::DicomRT_ROI * roi, QTr //------------------------------------------------------------------------------ void vvToolStructureSetManager::UpdateStructureSetInTreeWidget(int index, clitk::DicomRT_StructureSet * s) { // Insert ROI + /* const std::vector & rois = s->GetListOfROI(); for(unsigned int i=0; iGetROIs().begin(); iter != s->GetROIs().end(); iter++) { + clitk::DicomRT_ROI::Pointer roi = iter->second; + if (mMapROIToTreeWidget.find(roi) == mMapROIToTreeWidget.end()) + AddRoiInTreeWidget(roi, mTree); // replace mTree with ss if several SS + } } //------------------------------------------------------------------------------ @@ -334,14 +342,14 @@ void vvToolStructureSetManager::AddImage(vvImage * binaryImage, std::string file int n = mCurrentStructureSet->AddBinaryImageAsNewROI(binaryImage, filename); mLoadedROIIndex.push_back(n); if (m_modeBG) - mCurrentStructureSet->GetROI(n)->SetBackgroundValueLabelImage(BG); + mCurrentStructureSet->GetROIFromROINumber(n)->SetBackgroundValueLabelImage(BG); else - mCurrentStructureSet->GetROI(n)->SetForegroundValueLabelImage(BG); + mCurrentStructureSet->GetROIFromROINumber(n)->SetForegroundValueLabelImage(BG); // Change color if (nGetNumberOfTableValues ()) { double * color = mDefaultLUTColor->GetTableValue(n % mDefaultLUTColor->GetNumberOfTableValues ()); - mCurrentStructureSet->GetROI(n)->SetDisplayColor(color[0], color[1], color[2]); + mCurrentStructureSet->GetROIFromROINumber(n)->SetDisplayColor(color[0], color[1], color[2]); } // Add a new roi actor