clitkDicomRT_Contour.cxx
clitkDicomRT_ROI.cxx
clitkDicomRT_StructureSet.cxx
- clitkDicomRT_ROI_ConvertToImageFilter.cxx
+ clitkDicomRTStruct2ImageFilter.cxx
)
TARGET_LINK_LIBRARIES(clitkDicomRTStruct vtkHybrid)
#include <algorithm>
// clitk
-#include "clitkDicomRT_ROI_ConvertToImageFilter.h"
+#include "clitkDicomRTStruct2ImageFilter.h"
#include "clitkImageCommon.h"
// vtk
//--------------------------------------------------------------------
-clitk::DicomRT_ROI_ConvertToImageFilter::DicomRT_ROI_ConvertToImageFilter()
+clitk::DicomRTStruct2ImageFilter::DicomRTStruct2ImageFilter()
{
mROI = NULL;
mWriteOutput = false;
//--------------------------------------------------------------------
-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;
}
//--------------------------------------------------------------------
-void clitk::DicomRT_ROI_ConvertToImageFilter::SetCropMaskEnabled(bool b)
+void clitk::DicomRTStruct2ImageFilter::SetCropMaskEnabled(bool b)
{
mCropMask = b;
}
//--------------------------------------------------------------------
-void clitk::DicomRT_ROI_ConvertToImageFilter::SetOutputImageFilename(std::string s)
+void clitk::DicomRTStruct2ImageFilter::SetOutputImageFilename(std::string s)
{
mOutputFilename = s;
mWriteOutput = true;
//--------------------------------------------------------------------
-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) {
}
//--------------------------------------------------------------------
-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;
//--------------------------------------------------------------------
-vtkImageData * clitk::DicomRT_ROI_ConvertToImageFilter::GetOutput()
+vtkImageData * clitk::DicomRTStruct2ImageFilter::GetOutput()
{
assert(mBinaryImage);
return mBinaryImage;
=========================================================================*/
-#ifndef CLITKDICOMRT_ROI_CONVERTTOIMAGEFILTER_H
-#define CLITKDICOMRT_ROI_CONVERTTOIMAGEFILTER_H
+#ifndef CLITKDICOMRTSTRUCT2IMAGEFILTER_H
+#define CLITKDICOMRT_TRUCT2IMAGEFILTER_H
#include "clitkDicomRT_ROI.h"
#include "clitkImageCommon.h"
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
//--------------------------------------------------------------------
template <int Dimension>
-typename itk::Image<unsigned char,Dimension>::ConstPointer clitk::DicomRT_ROI_ConvertToImageFilter::GetITKOutput()
+typename itk::Image<unsigned char,Dimension>::ConstPointer clitk::DicomRTStruct2ImageFilter::GetITKOutput()
{
assert(mBinaryImage);
typedef itk::Image<unsigned char,Dimension> ConnectorImageType;
return connector->GetOutput();
}
//--------------------------------------------------------------------
-#endif // CLITKDICOMRT_ROI_CONVERTTOIMAGEFILTER_H
+#endif // CLITKDICOMRT_TRUCT2IMAGEFILTER_H
//--------------------------------------------------------------------
+//--------------------------------------------------------------------
+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<double> points;
+ points.resize(mNbOfPoints*3);
+ for(unsigned int i=0; i<nb; i++) {
+ double * p = mData->GetPoint(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;
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<unsigned int>(mNbOfPoints)*3);
// Organize values
vtkPolyData * clitk::DicomRT_Contour::GetMesh()
{
if (!mMeshIsUpToDate) {
- ComputeMesh();
+ ComputeMeshFromDataPoints();
}
return mMesh;
}
//--------------------------------------------------------------------
-void clitk::DicomRT_Contour::ComputeMesh()
+void clitk::DicomRT_Contour::SetMesh(vtkPolyData * mesh)
+{
+ mMesh = mesh;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+void clitk::DicomRT_Contour::ComputeMeshFromDataPoints()
{
// DD("ComputeMesh Contour");
mMesh = vtkSmartPointer<vtkPolyData>::New();
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<vtkPolyData>::New();
+ mMesh->Allocate(); //for cell structures
+ mPoints = vtkSmartPointer<vtkPoints>::New();
+ mMesh->SetPoints(mPoints);
+ vtkIdType ids[2];
+ for (unsigned int idx=0 ; idx<mNbOfPoints ; idx++) {
+ mMesh->GetPoints()->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;
+*/
+}
+//--------------------------------------------------------------------
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<vtkPoints> mData;
bool mMeshIsUpToDate;
///Z location of the contour
double mZ;
+
+ gdcm::Item * mItem;
private:
DicomRT_Contour();
#include "clitkDicomRT_ROI.h"
#include <vtkSmartPointer.h>
#include <vtkAppendPolyData.h>
+#include <vtkImageClip.h>
+#include <vtkMarchingSquares.h>
#if GDCM_MAJOR_VERSION == 2
#include "gdcmAttribute.h"
{
mName = "NoName";
mNumber = -1;
+ mImage = NULL;
mColor.resize(3);
mColor[0] = mColor[1] = mColor[2] = 0;
mMeshIsUpToDate = false;
mBackgroundValue = 0;
mForegroundValue = 1;
+ SetDicomUptodateFlag(false);
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
#if GDCM_MAJOR_VERSION == 2
-void clitk::DicomRT_ROI::Read(std::map<int, std::string> & 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 = {};
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<gdcm::SequenceOfItems> sqi2 = csq.GetValueAsSQ();
+ mContoursSequenceOfItems = csq.GetValueAsSQ();
+ gdcm::SmartPointer<gdcm::SequenceOfItems> & sqi2 = mContoursSequenceOfItems;
if( !sqi2 || !sqi2->GetNumberOfItems() )
{
}
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<int, std::string> & rois, gdcm::SQItem * item)
{
-
- // Change number if needed
-
- // TODO
-
// ROI number [Referenced ROI Number]
mNumber = atoi(item->GetEntryValue(0x3006,0x0084).c_str());
else {
std::cerr << "Warning. Could not read contour for structure <" << mName << ">, number" << mNumber << " ? I ignore it" << std::endl;
}
+ SetDicomUptodateFlag(true);
}
#endif
//--------------------------------------------------------------------
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<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
for(unsigned int i=0; i<mListOfContours.size(); i++) {
//--------------------------------------------------------------------
+#if GDCM_MAJOR_VERSION == 2
+
+//--------------------------------------------------------------------
+void clitk::DicomRT_ROI::UpdateDicomItem()
+{
+ if (GetDicomUptoDateFlag()) return;
+ DD("ROI::UpdateDicomItem");
+ DD(GetName());
+
+ // From now, only some item can be modified
+
+ // Set ROI Name 0x3006,0x26>
+ 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; i<mListOfContours.size(); i++) {
+ DD(i);
+ DicomRT_Contour::Pointer contour = mListOfContours[i];
+ contour->UpdateDicomItem();//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,
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<vtkImageClip> clipper = vtkSmartPointer<vtkImageClip>::New();
+ clipper->SetInput(image);
+ int* extent = image->GetExtent();
+ DDV(extent, 6);
+ // std::vector<int> extend;
+
+ // Prepare the marching squares
+ vtkSmartPointer<vtkMarchingSquares> squares = vtkSmartPointer<vtkMarchingSquares>::New();
+ squares->SetInput(clipper->GetOutput());
+
+ // Loop on slice
+ uint n = image->GetDimensions()[2];
+ DD(n);
+ for(uint i=0; i<n; i++) {
+ // Clip to the current slice
+ extent[4] = extent[5] = image->GetOrigin()[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());
+ }
+}
+//--------------------------------------------------------------------
itkNewMacro(Self);
void Print(std::ostream & os = std::cout) const;
-#if GDCM_MAJOR_VERSION == 2
- void Read(std::map<int, std::string> & rois, gdcm::Item const & item);
-#else
- void Read(std::map<int, std::string> & rois, gdcm::SQItem * item);
-#endif
void SetFromBinaryImage(vvImage * image, int n,
std::string name,
std::vector<double> color,
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<int, std::string> & rois, gdcm::SQItem * item);
+#endif
+
protected:
- void ComputeMesh();
std::string mName;
std::string mFilename;
int mNumber;
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<gdcm::SequenceOfItems> mContoursSequenceOfItems;
+#endif
private:
DicomRT_ROI();
#include "clitkDicomRT_StructureSet.h"
#include <vtksys/SystemTools.hxx>
#include "gdcmFile.h"
-#if GDCM_MAJOR_VERSION == 2
-#include "gdcmReader.h"
-#include "gdcmAttribute.h"
-#endif
//--------------------------------------------------------------------
clitk::DicomRT_StructureSet::DicomRT_StructureSet()
mName = "NoName";
mDate = "NoDate";
mTime = "NoTime";
+ mFile = NULL;
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-const std::vector<clitk::DicomRT_ROI::Pointer> & clitk::DicomRT_StructureSet::GetListOfROI() const
-{
- return mListOfROI;
-}
+// const std::vector<clitk::DicomRT_ROI::Pointer> & 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];
}
//--------------------------------------------------------------------
<< "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; i<mListOfROI.size(); i++) {
- mListOfROI[i]->Print(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<gdcm::SequenceOfItems> 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
mName = atname.GetValue();
mTime = time.GetValue();
+ // Temporary store the list of items
+ std::map<int, gdcm::Item*> mMapOfROIInfo;
+ std::map<int, gdcm::Item*> 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<gdcm::SequenceOfItems> roi_seq = ssroisq.GetValueAsSQ();
+ mROIInfoSequenceOfItems = ssroisq.GetValueAsSQ();
+ gdcm::SmartPointer<gdcm::SequenceOfItems> & roi_seq = mROIInfoSequenceOfItems;
assert(roi_seq); // TODO error message
for(unsigned int ridx = 0; ridx < roi_seq->GetNumberOfItems(); ++ridx)
{
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()) {
// 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<gdcm::SequenceOfItems> 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<int, gdcm::Item*>::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;
}
// 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 "
// 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++;
}
//--------------------------------------------------------------------
int clitk::DicomRT_StructureSet::AddBinaryImageAsNewROI(vvImage * im, std::string n)
{
- //DD("AddBinaryImageAsNewROI");
// Search max ROI number
int max = -1;
- for(unsigned int i=0; i<mListOfROI.size(); i++) {
- if (mListOfROI[i]->GetROINumber() > max)
- max = mListOfROI[i]->GetROINumber();
+ for(ROIConstIteratorType iter = mROIs.begin(); iter != mROIs.end(); iter++) {
+ // for(unsigned int i=0; i<mListOfROI.size(); i++) {
+ clitk::DicomRT_ROI::Pointer roi = iter->second;
+ 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
// 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;
}
//--------------------------------------------------------------------
#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 {
//--------------------------------------------------------------------
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro(Self);
+ typedef typename std::map<int, clitk::DicomRT_ROI::Pointer>::iterator ROIIteratorType;
+ typedef typename std::map<int, clitk::DicomRT_ROI::Pointer>::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<DicomRT_ROI::Pointer> & GetListOfROI() const;
- clitk::DicomRT_ROI * GetROI(int n);
+ clitk::DicomRT_ROI * GetROIFromROINumber(int n);
+ std::map<int, clitk::DicomRT_ROI::Pointer> & GetROIs() { return mROIs; }
const std::string & GetStudyID() const;
const std::string & GetStudyTime() const;
const std::string & GetStudyDate() const;
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;
std::string mName;
std::string mDate;
std::string mTime;
+
+ std::map<int, clitk::DicomRT_ROI::Pointer> mROIs;
std::map<int, std::string> mMapOfROIName;
- std::map<int, int> mMapOfROIIndex;
- std::vector<clitk::DicomRT_ROI::Pointer> mListOfROI;
+#if GDCM_MAJOR_VERSION == 2
+ gdcm::Reader * mReader;
+ gdcm::SmartPointer<gdcm::SequenceOfItems> mROIInfoSequenceOfItems;
+ gdcm::SmartPointer<gdcm::SequenceOfItems> mROIContoursSequenceOfItems;
+#endif
+ gdcm::File * mFile;
private:
DicomRT_StructureSet();
--- /dev/null
+/*=========================================================================
+ 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 PixelType>
+ class Image2DicomRTStructFilter: public clitk::FilterBase {
+
+ public:
+ Image2DicomRTStructFilter();
+ ~Image2DicomRTStructFilter();
+
+ typedef itk::Image<PixelType, 3> 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
+
--- /dev/null
+/*=========================================================================
+ 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 <iterator>
+#include <algorithm>
+
+// clitk
+#include "clitkImage2DicomRTStructFilter.h"
+#include "clitkImageCommon.h"
+#include "vvFromITK.h"
+
+// vtk
+#include <vtkPolyDataToImageStencil.h>
+#include <vtkSmartPointer.h>
+#include <vtkImageStencil.h>
+#include <vtkLinearExtrusionFilter.h>
+#include <vtkMetaImageWriter.h>
+#include <vtkImageData.h>
+
+// itk
+#include <itkImage.h>
+#include <itkVTKImageToImageFilter.h>
+
+//--------------------------------------------------------------------
+template<class PixelType>
+clitk::Image2DicomRTStructFilter<PixelType>::Image2DicomRTStructFilter()
+{
+ DD("Image2DicomRTStructFilter Const");
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class PixelType>
+clitk::Image2DicomRTStructFilter<PixelType>::~Image2DicomRTStructFilter()
+{
+ DD("Image2DicomRTStructFilter Destructor");
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class PixelType>
+void clitk::Image2DicomRTStructFilter<PixelType>::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<vtkPolyData> mesh = roi->GetMesh();
+ DD("done");
+
+ // Change the mesh (shift by 10);
+ // const vtkSmartPointer<vtkPoints> & points = mesh->GetPoints();
+ // for(uint i=0; i<mesh->GetNumberOfVerts (); 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");
+}
+//--------------------------------------------------------------------
+
+
+
+
itkSetMacro(CombineWithOrFlag, bool);
itkBooleanMacro(CombineWithOrFlag);
+ // I dont want to verify inputs information
+ virtual void VerifyInputInformation() { }
+
protected:
AddRelativePositionConstraintToLabelImageFilter();
virtual ~AddRelativePositionConstraintToLabelImageFilter() {}
/** ImageDimension constants */
itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension);
+ // I dont want to verify inputs information
+ virtual void VerifyInputInformation() { }
+
protected:
CropLikeImageFilter();
virtual ~CropLikeImageFilter() {}
// clitk
#include "clitkCommon.h"
+#include "clitkPasteImageFilter.h"
// itk
-#include "itkPasteImageFilter.h"
+//#include "itkPasteImageFilter.h"
//--------------------------------------------------------------------
template <class ImageType>
output->FillBuffer(GetBackgroundValue());
// Paste image inside
- typedef itk::PasteImageFilter<ImageType,ImageType> PasteFilterType;
+ typedef clitk::PasteImageFilter<ImageType,ImageType> PasteFilterType;
typename PasteFilterType::Pointer pasteFilter = PasteFilterType::New();
//pasteFilter->ReleaseDataFlagOn(); // change nothing ?
// pasteFilter->InPlaceOn(); // makt it seg fault
--- /dev/null
+/*=========================================================================
+ *
+ * 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
--- /dev/null
+/*=========================================================================
+ *
+ * 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
- BSD See included LICENSE.txt file
- CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
===========================================================================**/
-#ifndef __itkImageToVTKImageFilter_h\r
-#define __itkImageToVTKImageFilter_h\r
-#include "itkVTKImageExport.h"\r
-#include "vtkImageImport.h"\r
-#include "vtkImageData.h"\r
-\r
-namespace itk\r
-{\r
-\r
-/** \class ImageToVTKImageFilter\r
- * \brief Converts an ITK image into a VTK image and plugs a\r
- * itk data pipeline to a VTK datapipeline.\r
- *\r
- * This class puts together an itkVTKImageExporter and a vtkImageImporter.\r
- * It takes care of the details related to the connection of ITK and VTK\r
- * pipelines. The User will perceive this filter as an adaptor to which\r
- * an itk::Image can be plugged as input and a vtkImage is produced as\r
- * output.\r
- *\r
- * \ingroup ImageFilters\r
- */\r
-template <class TInputImage >\r
-class ITK_EXPORT ImageToVTKImageFilter : public ProcessObject\r
-{\r
-public:\r
- /** Standard class typedefs. */\r
- typedef ImageToVTKImageFilter Self;\r
- typedef ProcessObject Superclass;\r
- typedef SmartPointer<Self> Pointer;\r
- typedef SmartPointer<const Self> ConstPointer;\r
-\r
- /** Method for creation through the object factory. */\r
- itkNewMacro(Self);\r
-\r
- /** Run-time type information (and related methods). */\r
- itkTypeMacro(ImageToVTKImageFilter, ProcessObject);\r
-\r
- /** Some typedefs. */\r
- typedef TInputImage InputImageType;\r
- typedef typename InputImageType::ConstPointer InputImagePointer;\r
- typedef VTKImageExport< InputImageType> ExporterFilterType;\r
- typedef typename ExporterFilterType::Pointer ExporterFilterPointer;\r
-\r
- /** Get the output in the form of a vtkImage.\r
- This call is delegated to the internal vtkImageImporter filter */\r
- vtkImageData * GetOutput() const;\r
-\r
- /** Set the input in the form of an itk::Image */\r
- void SetInput( const InputImageType * );\r
-\r
- /** Return the internal VTK image importer filter.\r
- This is intended to facilitate users the access\r
- to methods in the importer */\r
- vtkImageImport * GetImporter() const;\r
-\r
- /** Return the internal ITK image exporter filter.\r
- This is intended to facilitate users the access\r
- to methods in the exporter */\r
- ExporterFilterType * GetExporter() const;\r
-\r
- /** This call delegate the update to the importer */\r
- void Update();\r
-\r
-protected:\r
- ImageToVTKImageFilter();\r
- virtual ~ImageToVTKImageFilter();\r
-\r
-private:\r
- ImageToVTKImageFilter(const Self&); //purposely not implemented\r
- void operator=(const Self&); //purposely not implemented\r
-\r
- ExporterFilterPointer m_Exporter;\r
- vtkImageImport * m_Importer;\r
-\r
-};\r
-\r
-} // end namespace itk\r
-\r
-#ifndef ITK_MANUAL_INSTANTIATION\r
-#include "itkImageToVTKImageFilter.txx"\r
-#endif\r
-\r
-#endif\r
-\r
-\r
-\r
+#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 TInputImage >
+class ITK_EXPORT ImageToVTKImageFilter : public ProcessObject
+{
+public:
+ /** Standard class typedefs. */
+ typedef ImageToVTKImageFilter Self;
+ typedef ProcessObject 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(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
+
// itk
#include <itkImageDuplicator.h>
-//--------------------------------------------------------------------
-template<class PointType>
-class comparePointsX {
-public:
- bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
-};
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template<class PairType>
-class comparePointsWithAngle {
-public:
- bool operator() (PairType i, PairType j) { return (i.second < j.second); }
-};
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template<int Dim>
-void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
- std::vector<itk::Point<double, Dim-1> > previous;
- HypercubeCorners<Dim-1>(previous);
- out.resize(previous.size()*2);
- for(unsigned int i=0; i<out.size(); i++) {
- itk::Point<double, Dim> p;
- if (i<previous.size()) p[0] = 0;
- else p[0] = 1;
- for(int j=0; j<Dim-1; j++)
- {
- p[j+1] = previous[i%previous.size()][j];
- }
- out[i] = p;
- }
-}
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template<>
-void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
- out.resize(2);
- out[0][0] = 0;
- out[1][0] = 1;
-}
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template<class ImageType>
-void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
- std::vector<typename ImageType::PointType> & 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; i<dim; i++)
- max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
- image->TransformIndexToPhysicalPoint(min_i, min_c);
- image->TransformIndexToPhysicalPoint(max_i, max_c);
-
- // Get corners coordinates
- HypercubeCorners<ImageType::ImageDimension>(bounds);
- for(unsigned int i=0; i<bounds.size(); i++) {
- for(unsigned int j=0; j<dim; j++) {
- if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
- if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
- }
- }
-}
-//--------------------------------------------------------------------
-
-
//--------------------------------------------------------------------
template <class ImageType>
void
//--------------------------------------------------------------------
+//--------------------------------------------------------------------
+template <class TImageType>
+void
+clitk::ExtractLymphStationsFilter<TImageType>::
+ExtractStation_2RL()
+{
+ if ((!CheckForStation("2R")) && (!CheckForStation("2L"))) return;
+
+ 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<MaskImageType>(m_ListOfStations["2R"], "seg/Station2R.mhd");
+ writeImage<MaskImageType>(m_ListOfStations["2L"], "seg/Station2L.mhd");
+ GetAFDB()->SetImageFilename("Station2R", "seg/Station2R.mhd");
+ GetAFDB()->SetImageFilename("Station2L", "seg/Station2L.mhd");
+ WriteAFDB();
+
+}
+//--------------------------------------------------------------------
+
+
//--------------------------------------------------------------------
template <class ImageType>
void
template <class ImageType>
void
clitk::ExtractLymphStationsFilter<ImageType>::
-ExtractStation_2RL_Ant_Limits2()
+ExtractStation_2RL_Ant_Limits2()
{
- // -----------------------------------------------------
- /* 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<MaskImageType>("BrachioCephalicArtery");
- MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
- MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
- MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
- MaskImagePointer Thyroid = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
- MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
- MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
-
- // Resize all structures like support
- BrachioCephalicArtery =
- clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, m_Working_Support, GetBackgroundValue());
- CommonCarotidArtery =
- clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, m_Working_Support, GetBackgroundValue());
- SubclavianArtery =
- clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, m_Working_Support, GetBackgroundValue());
- Thyroid =
- clitk::ResizeImageLike<MaskImageType>(Thyroid, m_Working_Support, GetBackgroundValue());
- Aorta =
- clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
- BrachioCephalicVein =
- clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, m_Working_Support, GetBackgroundValue());
- Trachea =
- clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
-
- // Extract slices
- std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
- clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
- std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
- clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
- std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
- clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
- std::vector<MaskSlicePointer> slices_SubclavianArtery;
- clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
- std::vector<MaskSlicePointer> slices_Thyroid;
- clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
- std::vector<MaskSlicePointer> slices_Aorta;
- clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
- std::vector<MaskSlicePointer> slices_Trachea;
- clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
- unsigned int n= slices_BrachioCephalicArtery.size();
- // Get the boundaries of one slice
- std::vector<MaskSlicePointType> bounds;
- ComputeImageBoundariesCoordinates<MaskSliceType>(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<MaskImagePointType> p3D;
-
- vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
- for(unsigned int i=0; i<n; i++) {
- // Labelize the slices
- slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
- GetBackgroundValue(), true, 1);
- slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
- GetBackgroundValue(), true, 1);
- slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
- GetBackgroundValue(), true, 1);
- slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
- GetBackgroundValue(), true, 1);
- slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
- GetBackgroundValue(), true, 1);
- slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
- GetBackgroundValue(), true, 1);
-
- // Search centroids
- std::vector<MaskSlicePointType> points2D;
- std::vector<MaskSlicePointType> centroids1;
- std::vector<MaskSlicePointType> centroids2;
- std::vector<MaskSlicePointType> centroids3;
- std::vector<MaskSlicePointType> centroids4;
- std::vector<MaskSlicePointType> centroids5;
- std::vector<MaskSlicePointType> centroids6;
- ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
- ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
- ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
- ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
- ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
- ComputeCentroids<MaskSliceType>(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<centroids1.size(); j++) points2D.push_back(centroids1[j]);
- for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
- for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
- for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
- for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
- for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
-
- // Sort by angle according to trachea centroid and vertical line,
- // in polar coordinates :
- // http://en.wikipedia.org/wiki/Polar_coordinate_system
- std::vector<MaskSlicePointType> centroids_trachea;
- ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
- typedef std::pair<MaskSlicePointType, double> PointAngleType;
- std::vector<PointAngleType> angles;
- for(unsigned int j=0; j<points2D.size(); j++) {
- //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
- double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
- double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
- double angle = 0;
- if (x>0) 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<PointAngleType>());
- for(unsigned int j=0; j<angles.size(); j++) {
- points2D[j] = angles[j].first;
- }
- // DDV(points2D, points2D.size());
-
- /* When vessels are far away, we try to replace the line segment
- with a curved line that join the two vessels but stay
- approximately at the same distance from the trachea centroids
- than the vessels.
-
- For that:
- - let a and c be two successive vessels centroids
- - id distance(a,c) < threshold, next point
-
- TODO HERE
-
- - compute mid position between two successive points -
- compute dist to trachea centroid for the 3 pts - if middle too
- low, add one point
- */
- std::vector<MaskSlicePointType> toadd;
- unsigned int index = 0;
- double dmax = 5;
- while (index<points2D.size()-1) {
- MaskSlicePointType a = points2D[index];
- MaskSlicePointType c = points2D[index+1];
-
- double dac = a.EuclideanDistanceTo(c);
- if (dac>dmax) {
-
- 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<ImagePointType> points3D;
- clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, m_Working_Support, points3D);
- for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
-
- // Build the mesh from the contour's points
- vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
- append->AddInput(mesh);
- }
-
- // DEBUG: write points3D
- clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
-
- // Build the final 3D mesh form the list 2D mesh
- append->Update();
- vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
-
- // Debug, write the mesh
- /*
- vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
- w->SetInput(mesh);
- w->SetFileName("bidon.vtk");
- w->Write();
- */
+ // 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<MaskImageType>(binarizedContour,
+ m_Working_Support,
+ GetForegroundValue());
- // Compute a single binary 3D image from the list of contours
- clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
- clitk::MeshToBinaryImageFilter<MaskImageType>::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());
-
// remove from support
typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
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();
//--------------------------------------------------------------------
+//--------------------------------------------------------------------
+template <class TImageType>
+void
+clitk::ExtractLymphStationsFilter<TImageType>::
+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<MaskImageType>(m_Working_Support);
+
+ ExtractStation_3A_AntPost_S5();
+ ExtractStation_3A_AntPost_S6();
+ ExtractStation_3A_AntPost_Superiorly();
+ ExtractStation_3A_Remove_Structures();
+
+ 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");
+
+ ExtractStation_3A_PostToBones();
+
+ // Keep a single CCL
+ m_ListOfStations["3A"] =
+ clitk::SliceBySliceKeepMainCCL<MaskImageType>(m_ListOfStations["3A"],
+ GetBackgroundValue(),
+ GetForegroundValue());
+
+ // Store image filenames into AFDB
+ writeImage<MaskImageType>(m_ListOfStations["3A"], "seg/Station3A.mhd");
+ GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd");
+ WriteAFDB();
+ StopSubStep();
+}
+//--------------------------------------------------------------------
+
+
+
//--------------------------------------------------------------------
template <class ImageType>
void
clitk::ExtractLymphStationsFilter<ImageType>::
-ExtractStation_3A_SI_Limits()
+ExtractStation_3A_AntPost_S5()
{
- // Apex of the chest or Sternum & Carina.
- StartNewStep("[Station 3A] Inf/Sup limits with Sternum and Carina");
+ StartNewStep("[Station 3A] Post limits around S5");
+
+ // First remove post to SVC
+ MaskImagePointer SVC = GetAFDB()->template GetImage <MaskImageType>("SVC");
+
+ // Trial in 3D -> difficulties superiorly. Stay slice by slice.
+ // Slice by slice not post to SVC. Use initial spacing
+ m_Working_Support =
+ clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, SVC, 2,
+ GetFuzzyThreshold("3A", "SVC"),
+ "NotPostTo", true,
+ SVC->GetSpacing()[0], false, false);
- // 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 <MaskImageType>("Aorta");
+ Aorta = clitk::ResizeImageLike<MaskImageType>(Aorta, S5, GetBackgroundValue());
- // Get Sternum and search for the upper position
- MaskImagePointer Sternum = GetAFDB()->template GetImage<MaskImageType>("Sternum");
+ // 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<MaskSlicePointer> slices;
+ clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
+ std::vector<MaskImagePointType> points;
+ for(uint i=0; i<slices.size(); i++) {
+ // Select most ant CCL
+ slices[i] = clitk::Labelize<MaskSliceType>(slices[i], GetBackgroundValue(), false, 1);
+ std::vector<MaskSlicePointType> c;
+ clitk::ComputeCentroids<MaskSliceType>(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<MaskSliceType, MaskSliceType>(slices[i], slices[i], l,
+ GetBackgroundValue(), true);
+
+ // Detect the most ant point
+ MaskSlicePointType p;
+ MaskImagePointType pA;
+ clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(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<MaskImageType>::Convert2DTo3D(p, Aorta, i, pA);
+ points.push_back(pA);
+ }
+
+ // DEBUG
+ // MaskImagePointer o = clitk::JoinSlices<MaskImageType>(slices, Aorta, 2);
+ // writeImage<MaskImageType>(o, "o.mhd");
+ // clitk::WriteListOfLandmarks<MaskImageType>(points, "Ant-Aorta.txt");
- // Search most sup point
- MaskImagePointType ps = Sternum->GetOrigin(); // initialise to avoid warning
- clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Sternum, GetBackgroundValue(), 2, false, ps);
- double m_SternumZ = ps[2]+Sternum->GetSpacing()[2]; // One more slice, because it is below this point
+ // Remove Post/Left to this point
+ m_Working_Support =
+ clitk::SliceBySliceSetBackgroundFromPoints<MaskImageType>(m_Working_Support,
+ GetBackgroundValue(), 2,
+ points,
+ true, // Set BG if X greater than point[x], and
+ true); // if Y greater than point[y]
+
+ StopCurrentStep<MaskImageType>(m_Working_Support);
+ m_ListOfStations["3A"] = m_Working_Support;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+ExtractStation_3A_AntPost_S6()
+{
+ StartNewStep("[Station 3A] Post limits around S6");
- //* Crop support :
+ // Consider Aorta
+ MaskImagePointer Aorta = GetAFDB()->template GetImage <MaskImageType>("Aorta");
+
+ // Limits the support to S6
+ MaskImagePointer S6 = m_ListOfSupports["S6"];
+ Aorta = clitk::ResizeImageLike<MaskImageType>(Aorta, S6, GetBackgroundValue());
+
+ // Extend 1cm anteriorly
+ MaskImagePointType radius; // in mm
+ radius[0] = 10;
+ radius[1] = 10;
+ radius[2] = 0; // required
+ Aorta = clitk::Dilate<MaskImageType>(Aorta, radius, GetBackgroundValue(), GetForegroundValue(), false);
+
+ // Not Post to
m_Working_Support =
- clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
- m_CarinaZ, m_SternumZ, true,
- GetBackgroundValue());
+ clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, Aorta, 2,
+ GetFuzzyThreshold("3A", "Aorta"),
+ "NotPostTo", true,
+ Aorta->GetSpacing()[0], false, false);
+
+ StopCurrentStep<MaskImageType>(m_Working_Support);
+ m_ListOfStations["3A"] = m_Working_Support;
+}
+//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+ExtractStation_3A_AntPost_Superiorly()
+{
+ StartNewStep("[Station 3A] Post limits superiorly");
+
+ /*
+ MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage <MaskImageType>("BrachioCephalicVein");
+ MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage <MaskImageType>("BrachioCephalicArtery");
+ MaskImagePointer CommonCarotidArteryLeft = GetAFDB()->template GetImage <MaskImageType>("CommonCarotidArteryLeft");
+ MaskImagePointer CommonCarotidArteryRight = GetAFDB()->template GetImage <MaskImageType>("CommonCarotidArteryRight");
+ MaskImagePointer SubclavianArteryLeft = GetAFDB()->template GetImage <MaskImageType>("SubclavianArteryLeft");
+ MaskImagePointer SubclavianArteryRight = GetAFDB()->template GetImage <MaskImageType>("SubclavianArteryRight");
+
+ // Not Post to
+#define RP(STRUCTURE) \
+ m_Working_Support = \
+ clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, STRUCTURE, 2, \
+ 0.5, \
+ "NotPostTo", true, \
+ STRUCTURE->GetSpacing()[0], false, false);
+
+ // RP(BrachioCephalicVein);
+ RP(BrachioCephalicArtery);
+ RP(CommonCarotidArteryRight);
+ RP(CommonCarotidArteryLeft);
+ RP(SubclavianArteryRight);
+ RP(SubclavianArteryLeft);
+ */
+
+ // Get or compute the binary mask that separate Ant/Post part
+ // according to vessels
+ MaskImagePointer binarizedContour = FindAntPostVessels2();
+ binarizedContour = clitk::ResizeImageLike<MaskImageType>(binarizedContour,
+ m_Working_Support,
+ GetBackgroundValue());
+
+ // remove from support
+ typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> 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<MaskImageType>(m_Working_Support);
m_ListOfStations["3A"] = m_Working_Support;
}
template <class ImageType>
void
clitk::ExtractLymphStationsFilter<ImageType>::
-ExtractStation_3A_Ant_Limits()
+ExtractStation_3A_Remove_Structures()
{
- StartNewStep("[Station 3A] Ant limits with Sternum");
+ 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", "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 <MaskImageType>("BrachioCephalicVein");
+ BrachioCephalicVein = clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein,
+ m_Working_Support,
+ GetBackgroundValue());
+ std::vector<MaskSlicePointer> slices;
+ std::vector<MaskSlicePointer> slices_BCV;
+ clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices);
+ clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BCV);
+ for(uint i=0; i<slices.size(); i++) {
+ // Labelize slices_BCV
+ slices_BCV[i] = Labelize<MaskSliceType>(slices_BCV[i], 0, true, 1);
+
+ // Compute centroids
+ std::vector<typename MaskSliceType::PointType> centroids;
+ ComputeCentroids<MaskSliceType>(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<MaskSliceType, MaskSliceType>(slices_BCV[i],
+ slices_BCV[i],
+ label,
+ GetBackgroundValue(),
+ true);
+ }
+
+ // Remove from the support
+ clitk::AndNot<MaskSliceType>(slices[i], slices_BCV[i], GetBackgroundValue());
+ }
+
+ // Joint
+ m_Working_Support = clitk::JoinSlices<MaskImageType>(slices, m_Working_Support, 2);
- // Get Sternum, keep posterior part.
- MaskImagePointer Sternum = GetAFDB()->template GetImage<MaskImageType>("Sternum");
- m_Working_Support =
- clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, Sternum, 2,
- GetFuzzyThreshold("3A", "Sternum"), "PostTo",
- false, 3, true, false);
StopCurrentStep<MaskImageType>(m_Working_Support);
m_ListOfStations["3A"] = m_Working_Support;
}
template <class ImageType>
void
clitk::ExtractLymphStationsFilter<ImageType>::
-ExtractStation_3A_Post_Limits()
+ExtractStation_3A_PostToBones()
{
- StartNewStep("[Station 3A] Post limits with SubclavianArtery");
+ StartNewStep("[Station 3A] Post limits with bones");
- // Get Sternum, keep posterior part.
- MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
+ // limits with bones
+ MaskImagePointer Bones = GetAFDB()->template GetImage<MaskImageType>("Bones");
m_Working_Support =
- clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, SubclavianArtery, 2,
- GetFuzzyThreshold("3A", "SubclavianArtery"), "AntTo",
+ clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, Bones, 2,
+ GetFuzzyThreshold("3A", "Bones"), "NotAntTo",
false, 3, true, false);
+
StopCurrentStep<MaskImageType>(m_Working_Support);
m_ListOfStations["3A"] = m_Working_Support;
}
//--------------------------------------------------------------------
-
-
+
//--------------------------------------------------------------------
-template <class ImageType>
+template <class TImageType>
void
-clitk::ExtractLymphStationsFilter<ImageType>::
-ExtractStation_3P_SI_Limits()
+clitk::ExtractLymphStationsFilter<TImageType>::
+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<MaskImageType>("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<MaskImageType>(m_Working_Support, 2,
- m_CarinaZ,
- m_ApexOfTheChest, true,
- GetBackgroundValue());
+ StartNewStep("Station 3P");
+ StartSubStep();
- StopCurrentStep<MaskImageType>(m_Working_Support);
+ // 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<MaskImageType>(m_Working_Support);
+
+ // LR limits inferiorly
+ ExtractStation_3P_LR_inf_Limits();
+
+ // 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
+
+ ExtractStation_8_Single_CCL_Limits(); // YES 8 !
+ ExtractStation_3P_Remove_Structures(); // after CCL
+
+ // Store image filenames into AFDB
+ writeImage<MaskImageType>(m_ListOfStations["3P"], "seg/Station3P.mhd");
+ GetAFDB()->SetImageFilename("Station3P", "seg/Station3P.mhd");
+ WriteAFDB();
+ StopSubStep();
}
//--------------------------------------------------------------------
clitk::ExtractLymphStationsFilter<ImageType>::
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", "SubclavianArteryLeft");
+ Remove_Structures("3P", "SubclavianArteryRight");
Remove_Structures("3P", "Esophagus");
- Remove_Structures("3P", "Azygousvein");
+ Remove_Structures("3P", "AzygousVein");
Remove_Structures("3P", "Thyroid");
Remove_Structures("3P", "VertebralArtery");
//--------------------------------------------------------------------
-//--------------------------------------------------------------------
-template <class ImageType>
-void
-clitk::ExtractLymphStationsFilter<ImageType>::
-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<MaskImageType>("Trachea");
-
- // Crop like current support (need by SliceBySliceSetBackgroundFromLineSeparation after)
- Trachea =
- clitk::ResizeImageLike<MaskImageType>(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<MaskImagePointType> p_A;
- std::vector<MaskImagePointType> p_B;
- std::vector<typename MaskSliceType::Pointer> slices;
- clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices);
- MaskImagePointType p;
- typename MaskSliceType::PointType sp;
- for(uint i=0; i<slices.size(); i++) {
- // First only consider the main CCL (trachea, not bronchus)
- slices[i] = Labelize<MaskSliceType>(slices[i], 0, true, 100);
- slices[i] = KeepLabels<MaskSliceType>(slices[i], GetBackgroundValue(),
- GetForegroundValue(), 1, 1, true);
- // Find most posterior point
- clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i], GetBackgroundValue(),
- 1, false, sp);
- clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp, Trachea, i, p);
- p_A.push_back(p);
- p[0] -= 20;
- p_B.push_back(p);
- }
- clitk::WriteListOfLandmarks<MaskImageType>(p_A, "trachea-post.txt");
-
- // Remove Ant part above those lines
- clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
- p_A, p_B,
- GetBackgroundValue(),
- 1, 10);
- // End, release memory
- GetAFDB()->template ReleaseImage<MaskImageType>("Trachea");
- StopCurrentStep<MaskImageType>(m_Working_Support);
- m_ListOfStations["3P"] = m_Working_Support;
-}
-//--------------------------------------------------------------------
-
-
-//--------------------------------------------------------------------
-template <class ImageType>
-void
-clitk::ExtractLymphStationsFilter<ImageType>::
-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<MaskImageType>("VertebralBody");
-
- // Crop like current support (need by SliceBySliceSetBackgroundFromLineSeparation after)
- VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
-
- writeImage<MaskImageType>(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<MaskImagePointType> p_A;
- std::vector<MaskImagePointType> p_B;
- std::vector<typename MaskSliceType::Pointer> slices;
- clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, slices);
- MaskImagePointType p;
- typename MaskSliceType::PointType sp;
- for(uint i=0; i<slices.size(); i++) {
- // Find most anterior point
- bool found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(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<MaskImageType>::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<MaskImageType>(p_A, "vb-ant.txt");
-
-
- // Remove Ant part above those lines
- clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
- p_A, p_B,
- GetBackgroundValue(),
- 1, -10);
- writeImage<MaskImageType>(m_Working_Support, "a.mhd");
-
- StopCurrentStep<MaskImageType>(m_Working_Support);
- m_ListOfStations["3P"] = m_Working_Support;
-}
-//--------------------------------------------------------------------
-
-
//--------------------------------------------------------------------
template <class ImageType>
void
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<MaskImageType>("VertebralArtery");
-
- clitk::AndNot<MaskImageType>(m_Working_Support, VertebralArtery);
+ "On the right side, the limit is defined by the air–soft-tissue
+ interface. On the left side, it is defined by the air–tissue
+ interface superiorly (Fig. 2A–C) and the aorta inferiorly
+ (Figs. 2D–I and 3A–C)."
- StopCurrentStep<MaskImageType>(m_Working_Support);
- m_ListOfStations["3P"] = m_Working_Support;
+ 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 (superior part) ");
+
+
}
//--------------------------------------------------------------------
MaskImagePointer AzygousVein = GetAFDB()->template GetImage<MaskImageType>("AzygousVein");
MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
+ DD("ici");
+ writeImage<MaskImageType>(m_Working_Support, "ws.mhd");
+
typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
relPosFilter->VerboseStepFlagOff();
relPosFilter->SetInputObject(AzygousVein);
relPosFilter->AddOrientationTypeString("R");
relPosFilter->SetInverseOrientationFlag(true);
- // relPosFilter->SetIntermediateSpacing(3);
+ relPosFilter->SetIntermediateSpacing(3);
relPosFilter->SetIntermediateSpacingFlag(false);
relPosFilter->SetFuzzyThreshold(0.7);
relPosFilter->AutoCropFlagOn();
relPosFilter->Update();
m_Working_Support = relPosFilter->GetOutput();
+ DD("la");
writeImage<MaskImageType>(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->SetIntermediateSpacing(3);
relPosFilter->SetIntermediateSpacingFlag(false);
relPosFilter->SetFuzzyThreshold(0.7);
relPosFilter->AutoCropFlagOn();
A, B, 2, 0, false);
// Get the CarinaZ position
- m_CarinaZ = FindCarinaSlicePosition();
+ double m_CarinaZ = FindCarina();
// Crop support
m_Working_Support =
clitk::ExtractLymphStationsFilter<TImageType>::
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<MaskImageType>(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<MaskImageType>(m_ListOfStations["8"], "seg/Station8.mhd");
+ GetAFDB()->SetImageFilename("Station8", "seg/Station8.mhd");
+ WriteAFDB();
+ StopSubStep();
}
//--------------------------------------------------------------------
// Get initial Mediastinum
m_Working_Support = m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("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<MaskImageType>(m_Working_Support, 2,
m_CarinaZ, true, GetBackgroundValue());
MaskImagePointer m_Support_Inferior_to_Carina =
clitk::CropImageRemoveGreaterThan<MaskImageType>(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<MaskImageType>(m_Support_Inferior_to_Carina, "seg/Support_Inf_Carina.mhd");
+ GetAFDB()->SetImageFilename("Support_Inf_Carina", "seg/Support_Inf_Carina.mhd");
+ writeImage<MaskImageType>(m_Support_Superior_to_Carina, "seg/Support_Sup_Carina.mhd");
+ 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<MaskImageType>(m_Support_Inferior_to_Carina);
+ m_ListOfSupports["S8"] = clitk::Clone<MaskImageType>(m_Support_Inferior_to_Carina);
+ m_ListOfSupports["S9"] = clitk::Clone<MaskImageType>(m_Support_Inferior_to_Carina);
+ m_ListOfSupports["S10"] = clitk::Clone<MaskImageType>(m_Support_Inferior_to_Carina);
+ m_ListOfSupports["S11"] = clitk::Clone<MaskImageType>(m_Support_Inferior_to_Carina);
+
+ // Store image filenames into AFDB
+ writeImage<MaskImageType>(m_ListOfSupports["S1R"], "seg/Support_S1R.mhd");
+ GetAFDB()->SetImageFilename("Support_S1R", "seg/Support_S1R.mhd");
+ writeImage<MaskImageType>(m_ListOfSupports["S1L"], "seg/Support_S1L.mhd");
+ GetAFDB()->SetImageFilename("Support_S1L", "seg/Support_S1L.mhd");
+
+ writeImage<MaskImageType>(m_ListOfSupports["S2L"], "seg/Support_S2L.mhd");
+ GetAFDB()->SetImageFilename("Support_S2L", "seg/Support_S2L.mhd");
+ writeImage<MaskImageType>(m_ListOfSupports["S2R"], "seg/Support_S2R.mhd");
+ GetAFDB()->SetImageFilename("Support_S2R", "seg/Support_S2R.mhd");
+
+ writeImage<MaskImageType>(m_ListOfSupports["S3P"], "seg/Support_S3P.mhd");
+ GetAFDB()->SetImageFilename("Support_S3P", "seg/Support_S3P.mhd");
+ writeImage<MaskImageType>(m_ListOfSupports["S3A"], "seg/Support_S3A.mhd");
+ GetAFDB()->SetImageFilename("Support_S3A", "seg/Support_S3A.mhd");
+
+ writeImage<MaskImageType>(m_ListOfSupports["S4L"], "seg/Support_S4L.mhd");
+ GetAFDB()->SetImageFilename("Support_S4L", "seg/Support_S4L.mhd");
+ writeImage<MaskImageType>(m_ListOfSupports["S4R"], "seg/Support_S4R.mhd");
+ GetAFDB()->SetImageFilename("Support_S4R", "seg/Support_S4R.mhd");
+
+ writeImage<MaskImageType>(m_ListOfSupports["S5"], "seg/Support_S5.mhd");
+ GetAFDB()->SetImageFilename("Support_S5", "seg/Support_S5.mhd");
+ writeImage<MaskImageType>(m_ListOfSupports["S6"], "seg/Support_S6.mhd");
+ GetAFDB()->SetImageFilename("Support_S6", "seg/Support_S6.mhd");
+
+ writeImage<MaskImageType>(m_ListOfSupports["S7"], "seg/Support_S7.mhd");
+ GetAFDB()->SetImageFilename("Support_S7", "seg/Support_S7.mhd");
+
+ writeImage<MaskImageType>(m_ListOfSupports["S8"], "seg/Support_S8.mhd");
+ GetAFDB()->SetImageFilename("Support_S8", "seg/Support_S8.mhd");
+
+ writeImage<MaskImageType>(m_ListOfSupports["S9"], "seg/Support_S9.mhd");
+ GetAFDB()->SetImageFilename("Support_S9", "seg/Support_S9.mhd");
+
+ writeImage<MaskImageType>(m_ListOfSupports["S10"], "seg/Support_S10.mhd");
+ GetAFDB()->SetImageFilename("Support_S10", "seg/Support_S10.mhd");
+
+ writeImage<MaskImageType>(m_ListOfSupports["S11"], "seg/Support_S11.mhd");
+ GetAFDB()->SetImageFilename("Support_S11", "seg/Support_S11.mhd");
+ WriteAFDB();
+}
+//--------------------------------------------------------------------
- // Consider only Superior to Carina
- m_Working_Support = m_Support_Superior_to_Carina;
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+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 = GetAFDB()->template GetImage <MaskImageType>("Sternum");
+ MaskImagePointType p;
+ p[0] = p[1] = p[2] = 0.0; // to avoid warning
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Sternum, GetBackgroundValue(), 2, false, p);
+ MaskImagePointer S1RL =
+ clitk::CropImageRemoveLowerThan<MaskImageType>(m_Working_Support, 2,
+ p[2], true, GetBackgroundValue());
+ m_Working_Support =
+ clitk::CropImageRemoveGreaterThan<MaskImageType>(m_Working_Support, 2,
+ p[2], true, GetBackgroundValue());
+ m_ListOfSupports["S1RL"] = S1RL;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+Support_LeftRight_S1R_S1L()
+{
+ // Step S1RL : Left-Right
+ StartNewStep("[Support] Left-Right S1R S1L");
+ std::vector<ImagePointType> A;
+ std::vector<ImagePointType> B;
+ // Search for centroid positions of trachea
+ MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
+ MaskImagePointer S1RL = m_ListOfSupports["S1RL"];
+ Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, S1RL, GetBackgroundValue());
+ std::vector<MaskSlicePointer> slices;
+ clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices);
+ for(uint i=0; i<slices.size(); i++) {
+ slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
+ std::vector<typename MaskSliceType::PointType> c;
+ clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
+ ImagePointType a,b;
+ clitk::PointsUtils<MaskImageType>::Convert2DTo3D(c[1], Trachea, i, a);
+ A.push_back(a);
+ b = a;
+ b[1] += 50;
+ B.push_back(b);
+ }
+ clitk::WriteListOfLandmarks<MaskImageType>(A, "S1LR_A.txt");
+ clitk::WriteListOfLandmarks<MaskImageType>(B, "S1LR_B.txt");
+ // Clone support
+ MaskImagePointer S1R = clitk::Clone<MaskImageType>(S1RL);
+ MaskImagePointer S1L = clitk::Clone<MaskImageType>(S1RL);
+ // Right part
+ clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(S1R, A, B,
+ GetBackgroundValue(), 0, 10);
+ S1R = clitk::AutoCrop<MaskImageType>(S1R, GetBackgroundValue());
+ m_ListOfSupports["S1R"] = S1R;
+ // Left part
+ clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(S1L, A, B,
+ GetBackgroundValue(), 0, -10);
+ S1L = clitk::AutoCrop<MaskImageType>(S1L, GetBackgroundValue());
+ m_ListOfSupports["S1L"] = S1L;
+}
+//--------------------------------------------------------------------
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+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 = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
+ MaskImagePointType p;
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
+ // I add slightly more than a slice
+ double CaudalMarginOfLeftBrachiocephalicVeinZ=p[2]+ 1.1*m_Working_Support->GetSpacing()[2];
+ GetAFDB()->SetDouble("CaudalMarginOfLeftBrachiocephalicVeinZ", CaudalMarginOfLeftBrachiocephalicVeinZ);
+ MaskImagePointer S2R =
+ clitk::CropImageRemoveLowerThan<MaskImageType>(m_Working_Support, 2,
+ CaudalMarginOfLeftBrachiocephalicVeinZ, true,
+ GetBackgroundValue());
+ m_ListOfSupports["S2R"] = S2R;
+
+ // S2L : Top Of Aortic Arch
+ MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
+ // I add slightly more than a slice
+ double TopOfAorticArchZ=p[2]+ 1.1*m_Working_Support->GetSpacing()[2];
+ GetAFDB()->SetDouble("TopOfAorticArchZ", TopOfAorticArchZ);
+ MaskImagePointer S2L =
+ clitk::CropImageRemoveLowerThan<MaskImageType>(m_Working_Support, 2,
+ TopOfAorticArchZ, true,
+ GetBackgroundValue());
+ m_ListOfSupports["S2L"] = S2L;
+}
+//--------------------------------------------------------------------
- // // LeftRight cut along trachea
- // MaskImagePointer Trachea = GetAFDB()->GetImage("Trachea");
- // // build a ant-post line for each slice
- // MaskImagePointer m_Support_SupRight =
- // clitk::CropImageRemoveLowerThan<MaskImageType>(m_Working_Support, 2,
- // m_CarinaZ, true, GetBackgroundValue());
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+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;
+ GetAFDB()->template ReleaseImage<MaskImageType>("Trachea");
+}
+//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+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<MaskImageType>(m_Working_Support);
+ MaskImagePointer S2R = m_ListOfSupports["S2R"];
+ MaskImagePointer S2L = m_ListOfSupports["S2L"];
+ clitk::AndNot<MaskImageType>(S4RL, S2R, GetBackgroundValue());
+ clitk::AndNot<MaskImageType>(S4RL, S2L, GetBackgroundValue());
+ S4RL = clitk::AutoCrop<MaskImageType>(S4RL, GetBackgroundValue());
+ // Copy, stop 4R at AzygousVein and 4L at LeftPulmonaryArtery
+ MaskImagePointer S4R = clitk::Clone<MaskImageType>(S4RL);
+ MaskImagePointer S4L = clitk::Clone<MaskImageType>(S4RL);
+
+ // Get AzygousVein and limit according to LowerBorderAzygousVein
+ MaskImagePointer LowerBorderAzygousVein
+ = GetAFDB()->template GetImage<MaskImageType>("LowerBorderAzygousVein");
+ std::vector<MaskImagePointType> c;
+ clitk::ComputeCentroids<MaskImageType>(LowerBorderAzygousVein, GetBackgroundValue(), c);
+ S4R = clitk::CropImageRemoveLowerThan<MaskImageType>(S4RL, 2,
+ c[1][2], true, GetBackgroundValue());
+ S4R = clitk::AutoCrop<MaskImageType>(S4R, GetBackgroundValue());
+ m_ListOfSupports["S4R"] = S4R;
- // Store image filenames into AFDB
- m_ListOfSupports["S1R"] = m_Working_Support;
- writeImage<MaskImageType>(m_ListOfSupports["S1R"], "seg/Support_S1R.mhd");
- GetAFDB()->SetImageFilename("Support_S1R", "seg/Support_S1R.mhd");
- WriteAFDB();
+ // Limit according to LeftPulmonaryArtery
+ MaskImagePointer LeftPulmonaryArtery
+ = GetAFDB()->template GetImage<MaskImageType>("LeftPulmonaryArtery");
+ MaskImagePointType p;
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(LeftPulmonaryArtery, GetBackgroundValue(),
+ 2, false, p);
+ S4L = clitk::CropImageRemoveLowerThan<MaskImageType>(S4RL, 2,
+ p[2], true, GetBackgroundValue());
+ S4L = clitk::AutoCrop<MaskImageType>(S4L, GetBackgroundValue());
+ m_ListOfSupports["S4L"] = S4L;
}
//--------------------------------------------------------------------
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+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 <class ImageType>
+typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
+clitk::ExtractLymphStationsFilter<ImageType>::
+LimitsWithTrachea(MaskImageType * input, int extremaDirection, int lineDirection,
+ double offset)
+{
+ MaskImagePointType min, max;
+ GetMinMaxBoundary<MaskImageType>(input, min, max);
+ return LimitsWithTrachea(input, extremaDirection, lineDirection, offset, max[2]);
+}
+template <class ImageType>
+typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
+clitk::ExtractLymphStationsFilter<ImageType>::
+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 = GetAFDB()->template GetImage<MaskImageType>("Trachea");
+
+ // Find extrema post positions
+ std::vector<MaskImagePointType> tracheaLeftPositionsA;
+ std::vector<MaskImagePointType> tracheaLeftPositionsB;
+
+ // Crop Trachea only on the Sup-Inf axes, without autocrop
+ // Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, input, GetBackgroundValue());
+ MaskImagePointType min, max;
+ GetMinMaxBoundary<MaskImageType>(input, min, max);
+ Trachea = clitk::CropImageAlongOneAxis<MaskImageType>(Trachea, 2, min[2], max[2],
+ false, GetBackgroundValue());
+
+ // Select the main CCL (because of bronchus)
+ Trachea = SliceBySliceKeepMainCCL<MaskImageType>(Trachea, GetBackgroundValue(), GetForegroundValue());
+
+ // Slice by slice, build the separation line
+ clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(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<tracheaLeftPositionsA.size(); i++) {
+ if (tracheaLeftPositionsA[i][2] > 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<MaskImageType>(input,
+ tracheaLeftPositionsA,
+ tracheaLeftPositionsB,
+ GetBackgroundValue(),
+ extremaDirection, offset);
+ MaskImagePointer output = clitk::AutoCrop<MaskImageType>(input, GetBackgroundValue());
+ return output;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+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 <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+Support_S3P()
+{
+ StartNewStep("[Support] Ant limits of S3P and Post limits of S1RL, S2RL, S4RL");
+
+ // Initial S3P support
+ MaskImagePointer S3P = clitk::Clone<MaskImageType>(m_ListOfSupports["Support_Superior_to_Carina"]);
+
+ // Stop at Lung Apex
+ double m_ApexOfTheChest = FindApexOfTheChest();
+ S3P =
+ clitk::CropImageRemoveGreaterThan<MaskImageType>(S3P, 2,
+ m_ApexOfTheChest, true,
+ GetBackgroundValue());
+ // Ant limits with Trachea
+ S3P = LimitsWithTrachea(S3P, 1, 0, 10);
+ m_ListOfSupports["S3P"] = S3P;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+Support_S3A()
+{
+ StartNewStep("[Support] Sup-Inf and Post limits for S3A");
+
+ // Initial S3A support
+ MaskImagePointer S3A = clitk::Clone<MaskImageType>(m_ListOfSupports["Support_Superior_to_Carina"]);
+
+ // Stop at Lung Apex or like S2/S1 (upper border Sternum - manubrium) ?
+
+ //double m_ApexOfTheChest = FindApexOfTheChest();
+
+ MaskImagePointer Sternum = GetAFDB()->template GetImage <MaskImageType>("Sternum");
+ MaskImagePointType p;
+ p[0] = p[1] = p[2] = 0.0; // to avoid warning
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Sternum, GetBackgroundValue(), 2, false, p);
+
+ S3A =
+ clitk::CropImageRemoveGreaterThan<MaskImageType>(S3A, 2,
+ //m_ApexOfTheChest
+ p[2], true,
+ GetBackgroundValue());
+ // Ant limits with Trachea
+ S3A = LimitsWithTrachea(S3A, 1, 0, -10);
+ m_ListOfSupports["S3A"] = S3A;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+Support_S5()
+{
+ StartNewStep("[Support] Sup-Inf limits S5 with aorta");
+
+ // Initial S5 support
+ MaskImagePointer S5 = clitk::Clone<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("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 = GetAFDB()->template GetImage<MaskImageType>("PulmonaryTrunk");
+ MaskImagePointType p;
+ p[0] = p[1] = p[2] = 0.0; // to avoid warning
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(PulmonaryTrunk, GetBackgroundValue(), 2, false, p);
+
+ // Cut Sup/Inf
+ S5 = clitk::CropImageAlongOneAxis<MaskImageType>(S5, 2, p[2], sup, true, GetBackgroundValue());
+
+ m_ListOfSupports["S5"] = S5;
+}
+//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+Support_S6()
+{
+ StartNewStep("[Support] Sup-Inf limits S6 with aorta");
+
+ // Initial S6 support like S3A
+ MaskImagePointer S6 = clitk::Clone<MaskImageType>(m_ListOfSupports["S3A"]);
+
+ // Inf Sup limits with Aorta
+ double sup = FindSuperiorBorderOfAorticArch();
+ double inf = FindInferiorBorderOfAorticArch();
+
+ // Cut Sup/Inf
+ S6 = clitk::CropImageAlongOneAxis<MaskImageType>(S6, 2, inf, sup, true, GetBackgroundValue());
+
+ m_ListOfSupports["S6"] = S6;
+}
+//--------------------------------------------------------------------
+
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 "nosupport" - "Do not recompute the station supports if already available" flag off
+
+
+section "Options for Station 3A"
+option "S3A_ft_SVC" - "Fuzzy Threshold for NotPost to SVC" double default="0.5" no
+option "S3A_ft_Bones" - "Fuzzy Threshold for Post to Bones" double default="0.5" no
+option "S3A_t_Bones" - "Binary HU Threshold for bones" double default="250" no
+option "S3A_ft_Aorta" - "Fuzzy Threshold for NotPost to Aorta" double default="0.5" no
+option "S3A_ft_SubclavianArtery" - "Fuzzy Threshold for Ant to SubclavianArtery" double default="0.2" no
+
section "Options for Station 8"
option "S8_esophagusDilatationForAnt" - "Dilatation of esophagus, in mm, for 'anterior' limits (default=15,2,1)" double no multiple
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
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();
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 FindAntPostVessels();
+ MaskImagePointer FindAntPostVessels2();
// Global parameters
typedef std::map<std::string, double> FuzzyThresholdByStructureType;
std::map<std::string, FuzzyThresholdByStructureType> m_FuzzyThreshold;
+ typedef std::map<std::string, double> ThresholdByStructureType;
+ std::map<std::string, ThresholdByStructureType> 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;
// 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_AntPost_S5();
+ void ExtractStation_3A_AntPost_S6();
+ void ExtractStation_3A_AntPost_Superiorly();
+ void ExtractStation_3A_Remove_Structures();
+ void ExtractStation_3A_PostToBones();
// Station 2RL
void ExtractStation_2RL();
void ExtractStation_2RL_SeparateRL();
vtkSmartPointer<vtkPolyData> Build3DMeshFrom2DContour(const std::vector<ImagePointType> & 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 7
void ExtractStation_7();
void ExtractStation_7_SetDefaultValues();
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;
#include "clitkAutoCropFilter.h"
#include "clitkSegmentationUtils.h"
#include "clitkSliceBySliceRelativePositionFilter.h"
+#include "clitkMeshToBinaryImageFilter.h"
// itk
#include <itkStatisticsLabelMapFilter.h>
// itk ENST
#include "RelativePositionPropImageFilter.h"
+// vtk
+#include <vtkAppendPolyData.h>
+#include <vtkPolyDataWriter.h>
+#include <vtkCellArray.h>
+
//--------------------------------------------------------------------
template <class TImageType>
clitk::ExtractLymphStationsFilter<TImageType>::
this->SetNumberOfRequiredInputs(1);
SetBackgroundValue(0);
SetForegroundValue(1);
+ ComputeStationsSupportsFlagOn();
// Default values
ExtractStation_8_SetDefaultValues();
m_Input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
m_Mediastinum = GetAFDB()->template GetImage <MaskImageType>("Mediastinum");
+ // Clean some computer landmarks to force the recomputation
+ GetAFDB()->RemoveTag("AntPostVesselsSeparation");
+
// Global supports for stations
- StartNewStep("Supports for stations");
- StartSubStep();
- ExtractStationSupports();
- StopSubStep();
+ if (GetComputeStationsSupportsFlag()) {
+ StartNewStep("Supports for stations");
+ StartSubStep();
+ GetAFDB()->RemoveTag("CarinaZ");
+ GetAFDB()->RemoveTag("ApexOfTheChestZ");
+ GetAFDB()->RemoveTag("ApexOfTheChest");
+ GetAFDB()->RemoveTag("RightBronchus");
+ GetAFDB()->RemoveTag("LeftBronchus");
+ GetAFDB()->RemoveTag("SuperiorBorderOfAorticArchZ");
+ GetAFDB()->RemoveTag("SuperiorBorderOfAorticArch");
+ GetAFDB()->RemoveTag("InferiorBorderOfAorticArchZ");
+ GetAFDB()->RemoveTag("InferiorBorderOfAorticArch");
+ ExtractStationSupports();
+ StopSubStep();
+ }
+ else {
+ m_ListOfSupports["S1R"] = GetAFDB()->template GetImage<MaskImageType>("Support_S1R");
+ m_ListOfSupports["S1L"] = GetAFDB()->template GetImage<MaskImageType>("Support_S1L");
+ m_ListOfSupports["S2R"] = GetAFDB()->template GetImage<MaskImageType>("Support_S2R");
+ m_ListOfSupports["S2L"] = GetAFDB()->template GetImage<MaskImageType>("Support_S2L");
+ m_ListOfSupports["S4R"] = GetAFDB()->template GetImage<MaskImageType>("Support_S4R");
+ m_ListOfSupports["S4L"] = GetAFDB()->template GetImage<MaskImageType>("Support_S4L");
+
+ m_ListOfSupports["S3A"] = GetAFDB()->template GetImage<MaskImageType>("Support_S3A");
+ m_ListOfSupports["S3P"] = GetAFDB()->template GetImage<MaskImageType>("Support_S3P");
+ m_ListOfSupports["S5"] = GetAFDB()->template GetImage<MaskImageType>("Support_S5");
+ m_ListOfSupports["S6"] = GetAFDB()->template GetImage<MaskImageType>("Support_S6");
+ m_ListOfSupports["S7"] = GetAFDB()->template GetImage<MaskImageType>("Support_S7");
+ m_ListOfSupports["S8"] = GetAFDB()->template GetImage<MaskImageType>("Support_S8");
+ m_ListOfSupports["S9"] = GetAFDB()->template GetImage<MaskImageType>("Support_S9");
+ m_ListOfSupports["S10"] = GetAFDB()->template GetImage<MaskImageType>("Support_S10");
+ m_ListOfSupports["S11"] = GetAFDB()->template GetImage<MaskImageType>("Support_S11");
+ }
// Extract Station8
- StartNewStep("Station 8");
- StartSubStep();
ExtractStation_8();
- 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();
+
+ // HERE
+
+ // Extract Station2RL
+ StartNewStep("Station 2RL");
+ StartSubStep();
+ ExtractStation_2RL();
StopSubStep();
// Extract Station7
//--------------------------------------------------------------------
+
//--------------------------------------------------------------------
-template <class TImageType>
-void
-clitk::ExtractLymphStationsFilter<TImageType>::
-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<MaskImageType>(m_ListOfStations["3P"], "seg/Station3P.mhd");
- GetAFDB()->SetImageFilename("Station3P", "seg/Station3P.mhd");
- WriteAFDB();
+template<class PointType>
+class comparePointsX {
+public:
+ bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
+};
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class PairType>
+class comparePointsWithAngle {
+public:
+ bool operator() (PairType i, PairType j) { return (i.second < j.second); }
+};
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<int Dim>
+void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
+ std::vector<itk::Point<double, Dim-1> > previous;
+ HypercubeCorners<Dim-1>(previous);
+ out.resize(previous.size()*2);
+ for(unsigned int i=0; i<out.size(); i++) {
+ itk::Point<double, Dim> p;
+ if (i<previous.size()) p[0] = 0;
+ else p[0] = 1;
+ for(int j=0; j<Dim-1; j++)
+ {
+ p[j+1] = previous[i%previous.size()][j];
+ }
+ out[i] = p;
}
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TImageType>
-void
-clitk::ExtractLymphStationsFilter<TImageType>::
-ExtractStation_3A()
-{
- if (CheckForStation("3A")) {
- ExtractStation_3A_SI_Limits();
- ExtractStation_3A_Ant_Limits();
- ExtractStation_3A_Post_Limits();
- // Store image filenames into AFDB
- writeImage<MaskImageType>(m_ListOfStations["3A"], "seg/Station3A.mhd");
- GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd");
- WriteAFDB();
- }
+template<>
+void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
+ out.resize(2);
+ out[0][0] = 0;
+ out[1][0] = 1;
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
-template <class TImageType>
-void
-clitk::ExtractLymphStationsFilter<TImageType>::
-ExtractStation_2RL()
+template<class ImageType>
+void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
+ std::vector<typename ImageType::PointType> & bounds)
{
- 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<MaskImageType>(m_ListOfStations["2R"], "seg/Station2R.mhd");
- writeImage<MaskImageType>(m_ListOfStations["2L"], "seg/Station2L.mhd");
- GetAFDB()->SetImageFilename("Station2R", "seg/Station2R.mhd");
- GetAFDB()->SetImageFilename("Station2L", "seg/Station2L.mhd");
- WriteAFDB();
+ // 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; i<dim; i++)
+ max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
+ image->TransformIndexToPhysicalPoint(min_i, min_c);
+ image->TransformIndexToPhysicalPoint(max_i, max_c);
+
+ // Get corners coordinates
+ HypercubeCorners<ImageType::ImageDimension>(bounds);
+ for(unsigned int i=0; i<bounds.size(); i++) {
+ for(unsigned int j=0; j<dim; j++) {
+ if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
+ if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
+ }
}
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
+//--------------------------------------------------------------------
+template <class TImageType>
+void
+clitk::ExtractLymphStationsFilter<TImageType>::
+SetThreshold(std::string station, std::string tag, double value)
+{
+ m_Threshold[station][tag] = value;
+}
+//--------------------------------------------------------------------
+
+
//--------------------------------------------------------------------
template <class TImageType>
double
//--------------------------------------------------------------------
+//--------------------------------------------------------------------
+template <class TImageType>
+double
+clitk::ExtractLymphStationsFilter<TImageType>::
+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 <class TImageType>
void
template <class TImageType>
double
clitk::ExtractLymphStationsFilter<TImageType>::
-FindCarinaSlicePosition()
+FindCarina()
{
double z;
try {
clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
// We dont need Carina structure from now
- Carina->Delete();
+ GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
// Put inside the AFDB
GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
//--------------------------------------------------------------------
+//--------------------------------------------------------------------
+template <class TImageType>
+double
+clitk::ExtractLymphStationsFilter<TImageType>::
+FindApexOfTheChest()
+{
+ double z;
+ try {
+ z = GetAFDB()->GetDouble("ApexOfTheChestZ");
+ }
+ catch(clitk::ExceptionObject e) {
+ DD("FindApexOfTheChestPosition");
+ MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
+ MaskImagePointType p;
+ p[0] = p[1] = p[2] = 0.0; // to avoid warning
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
+
+ // We dont need Lungs structure from now
+ GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
+
+ // Put inside the AFDB
+ GetAFDB()->SetPoint3D("ApexOfTheChest", p);
+ p[2] -= 5; // We consider 5 mm lower
+ GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
+ WriteAFDB();
+ z = p[2];
+ }
+ return z;
+}
+//--------------------------------------------------------------------
+
+
//--------------------------------------------------------------------
template <class TImageType>
void
MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
// Get the Carina position
- m_CarinaZ = FindCarinaSlicePosition();
+ double m_CarinaZ = FindCarina();
// Consider only inferiorly to the Carina
MaskImagePointer m_Working_Trachea =
}
//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+template <class TImageType>
+double
+clitk::ExtractLymphStationsFilter<TImageType>::
+FindSuperiorBorderOfAorticArch()
+{
+ double z;
+ try {
+ z = GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
+ }
+ catch(clitk::ExceptionObject e) {
+ DD("FindSuperiorBorderOfAorticArch");
+ MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
+ MaskImagePointType p;
+ p[0] = p[1] = p[2] = 0.0; // to avoid warning
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
+ p[2] += Aorta->GetSpacing()[2]; // the slice above
+
+ // We dont need Lungs structure from now
+ GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
+
+ // Put inside the AFDB
+ GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
+ GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
+ WriteAFDB();
+ z = p[2];
+ }
+ return z;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class TImageType>
+double
+clitk::ExtractLymphStationsFilter<TImageType>::
+FindInferiorBorderOfAorticArch()
+{
+ double z;
+ try {
+ z = GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
+ }
+ catch(clitk::ExceptionObject e) {
+ DD("FindInferiorBorderOfAorticArch");
+ MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
+ std::vector<MaskSlicePointer> slices;
+ clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
+ bool found=false;
+ uint i = slices.size()-1;
+ while (!found) {
+ slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
+ std::vector<typename MaskSliceType::PointType> c;
+ clitk::ComputeCentroids<MaskSliceType>(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
+ GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
+
+ // Put inside the AFDB
+ GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
+ GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
+ WriteAFDB();
+ z = lower[2];
+ }
+ return z;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
+clitk::ExtractLymphStationsFilter<ImageType>::
+FindAntPostVessels()
+{
+ // -----------------------------------------------------
+ /* 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 {
+ DD("FindAntPostVessels try to get");
+ binarizedContour = GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
+ }
+ catch(clitk::ExceptionObject e) {
+ DD("not found");
+ 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 = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
+ MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
+ MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
+ MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
+ MaskImagePointer Thyroid = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
+ MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
+ MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
+
+ // Create a temporay support
+ // From first slice of BrachioCephalicVein to end of 3A
+ MaskImagePointer support = GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
+ MaskImagePointType p;
+ p[0] = p[1] = p[2] = 0.0; // to avoid warning
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
+ double inf = p [2];
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
+ GetBackgroundValue(), 2, false, p);
+ double sup = p [2];
+ support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
+ false, GetBackgroundValue());
+
+ // Resize all structures like support
+ BrachioCephalicArtery =
+ clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
+ CommonCarotidArtery =
+ clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
+ SubclavianArtery =
+ clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
+ Thyroid =
+ clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
+ Aorta =
+ clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
+ BrachioCephalicVein =
+ clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
+ Trachea =
+ clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
+
+ // Extract slices
+ std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
+ clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
+ std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
+ clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
+ std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
+ clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
+ std::vector<MaskSlicePointer> slices_SubclavianArtery;
+ clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
+ std::vector<MaskSlicePointer> slices_Thyroid;
+ clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
+ std::vector<MaskSlicePointer> slices_Aorta;
+ clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
+ std::vector<MaskSlicePointer> slices_Trachea;
+ clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
+ unsigned int n= slices_BrachioCephalicArtery.size();
+
+ // Get the boundaries of one slice
+ std::vector<MaskSlicePointType> bounds;
+ ComputeImageBoundariesCoordinates<MaskSliceType>(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<MaskImagePointType> p3D;
+
+ vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
+ for(unsigned int i=0; i<n; i++) {
+ // Labelize the slices
+ slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
+ GetBackgroundValue(), true, 1);
+ slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
+ GetBackgroundValue(), true, 1);
+ slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
+ GetBackgroundValue(), true, 1);
+ slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
+ GetBackgroundValue(), true, 1);
+ slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
+ GetBackgroundValue(), true, 1);
+ slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
+ GetBackgroundValue(), true, 1);
+
+ // Search centroids
+ std::vector<MaskSlicePointType> points2D;
+ std::vector<MaskSlicePointType> centroids1;
+ std::vector<MaskSlicePointType> centroids2;
+ std::vector<MaskSlicePointType> centroids3;
+ std::vector<MaskSlicePointType> centroids4;
+ std::vector<MaskSlicePointType> centroids5;
+ std::vector<MaskSlicePointType> centroids6;
+ ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
+ ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
+ ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
+ ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
+ ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
+ ComputeCentroids<MaskSliceType>(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<centroids1.size(); j++) points2D.push_back(centroids1[j]);
+ for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
+ for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
+ for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
+ for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
+ for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
+
+ // Sort by angle according to trachea centroid and vertical line,
+ // in polar coordinates :
+ // http://en.wikipedia.org/wiki/Polar_coordinate_system
+ std::vector<MaskSlicePointType> centroids_trachea;
+ ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
+ typedef std::pair<MaskSlicePointType, double> PointAngleType;
+ std::vector<PointAngleType> angles;
+ for(unsigned int j=0; j<points2D.size(); j++) {
+ //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
+ double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
+ double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
+ double angle = 0;
+ if (x>0) 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<PointAngleType>());
+ for(unsigned int j=0; j<angles.size(); j++) {
+ points2D[j] = angles[j].first;
+ }
+ // DDV(points2D, points2D.size());
+
+ /* When vessels are far away, we try to replace the line segment
+ with a curved line that join the two vessels but stay
+ approximately at the same distance from the trachea centroids
+ than the vessels.
+
+ For that:
+ - let a and c be two successive vessels centroids
+ - id distance(a,c) < threshold, next point
+
+ TODO HERE
+
+ - compute mid position between two successive points -
+ compute dist to trachea centroid for the 3 pts - if middle too
+ low, add one point
+ */
+ std::vector<MaskSlicePointType> toadd;
+ unsigned int index = 0;
+ double dmax = 5;
+ while (index<points2D.size()-1) {
+ MaskSlicePointType a = points2D[index];
+ MaskSlicePointType c = points2D[index+1];
+
+ double dac = a.EuclideanDistanceTo(c);
+ if (dac>dmax) {
+
+ 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<ImagePointType> points3D;
+ clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
+ for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
+
+ // Build the mesh from the contour's points
+ vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
+ append->AddInput(mesh);
+ }
+
+ // DEBUG: write points3D
+ clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
+
+ // Build the final 3D mesh form the list 2D mesh
+ append->Update();
+ vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
+
+ // Debug, write the mesh
+ /*
+ vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
+ w->SetInput(mesh);
+ w->SetFileName("bidon.vtk");
+ w->Write();
+ */
+
+ // Compute a single binary 3D image from the list of contours
+ clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
+ clitk::MeshToBinaryImageFilter<MaskImageType>::New();
+ filter->SetMesh(mesh);
+ filter->SetLikeImage(support);
+ filter->Update();
+ binarizedContour = filter->GetOutput();
+
+ // Crop
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
+ inf = p[2];
+ DD(p);
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
+ sup = p[2];
+ DD(p);
+ binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
+ false, GetBackgroundValue());
+ // Update the AFDB
+ writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
+ GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
+ 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<MaskImageType> 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<MaskImageType>(m_Working_Support);
+ //m_ListOfStations["2R"] = m_Working_Support;
+ //m_ListOfStations["2L"] = m_Working_Support;
+}
+//--------------------------------------------------------------------
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
+clitk::ExtractLymphStationsFilter<ImageType>::
+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 {
+ DD("FindAntPostVessels try to get");
+ binarizedContour = GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
+ }
+ catch(clitk::ExceptionObject e) {
+ DD("not found");
+ 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<std::string, MaskImagePointer> MapOfStructures;
+ typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
+ MapOfStructures["BrachioCephalicArtery"] = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
+ MapOfStructures["BrachioCephalicVein"] = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
+ MapOfStructures["CommonCarotidArteryLeft"] = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryLeft");
+ MapOfStructures["CommonCarotidArteryRight"] = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryRight");
+ MapOfStructures["SubclavianArteryLeft"] = GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryLeft");
+ MapOfStructures["SubclavianArteryRight"] = GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryRight");
+ MapOfStructures["Thyroid"] = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
+ MapOfStructures["Aorta"] = GetAFDB()->template GetImage<MaskImageType>("Aorta");
+ MapOfStructures["Trachea"] = GetAFDB()->template GetImage<MaskImageType>("Trachea");
+
+ std::vector<std::string> ListOfStructuresNames;
+
+ // Create a temporay support
+ // From first slice of BrachioCephalicVein to end of 3A
+ MaskImagePointer support = GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
+ MaskImagePointType p;
+ p[0] = p[1] = p[2] = 0.0; // to avoid warning
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
+ GetBackgroundValue(), 2, true, p);
+ double inf = p[2];
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
+ GetBackgroundValue(), 2, false, p);
+ double sup = p[2]+support->GetSpacing()[2];//one slice more ?
+ support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
+ false, GetBackgroundValue());
+
+ // Resize all structures like support
+ for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
+ it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
+ }
+
+ // Extract slices
+ typedef std::vector<MaskSlicePointer> SlicesType;
+ std::map<std::string, SlicesType> MapOfSlices;
+ for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
+ SlicesType s;
+ clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
+ MapOfSlices[it->first] = s;
+ }
+
+ unsigned int n= MapOfSlices["Trachea"].size();
+
+ // Get the boundaries of one slice
+ std::vector<MaskSlicePointType> bounds;
+ ComputeImageBoundariesCoordinates<MaskSliceType>(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<MaskImagePointType> p3D;
+ vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
+ for(unsigned int i=0; i<n; i++) {
+
+ // Labelize the slices
+ for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
+ MaskSlicePointer & s = MapOfSlices[it->first][i];
+ s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
+ }
+
+ // Search centroids
+ std::vector<MaskSlicePointType> points2D;
+ typedef std::vector<MaskSlicePointType> CentroidsType;
+ std::map<std::string, CentroidsType> MapOfCentroids;
+ for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
+ std::string structure = it->first;
+ MaskSlicePointer & s = MapOfSlices[structure][i];
+ CentroidsType c;
+ clitk::ComputeCentroids<MaskSliceType>(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<c.size(); j++) points2D.push_back(c[j]);
+ }
+ }
+
+ // Sort by angle according to trachea centroid and vertical line,
+ // in polar coordinates :
+ // http://en.wikipedia.org/wiki/Polar_coordinate_system
+ // std::vector<MaskSlicePointType> centroids_trachea;
+ //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
+ CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
+ typedef std::pair<MaskSlicePointType, double> PointAngleType;
+ std::vector<PointAngleType> angles;
+ for(unsigned int j=0; j<points2D.size(); j++) {
+ //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
+ double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
+ double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
+ double angle = 0;
+ if (x>0) 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<PointAngleType>());
+ for(unsigned int j=0; j<angles.size(); j++) {
+ points2D[j] = angles[j].first;
+ }
+ // DDV(points2D, points2D.size());
+
+ /* When vessels are far away, we try to replace the line segment
+ with a curved line that join the two vessels but stay
+ approximately at the same distance from the trachea centroids
+ than the vessels.
+
+ For that:
+ - let a and c be two successive vessels centroids
+ - id distance(a,c) < threshold, next point
+
+ TODO HERE
+
+ - compute mid position between two successive points -
+ compute dist to trachea centroid for the 3 pts - if middle too
+ low, add one point
+ */
+ std::vector<MaskSlicePointType> toadd;
+ unsigned int index = 0;
+ double dmax = 5;
+ while (index<points2D.size()-1) {
+ MaskSlicePointType a = points2D[index];
+ MaskSlicePointType c = points2D[index+1];
+
+ double dac = a.EuclideanDistanceTo(c);
+ if (dac>dmax) {
+
+ 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<ImagePointType> points3D;
+ clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
+ for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
+
+ // Build the mesh from the contour's points
+ vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
+ append->AddInput(mesh);
+ // if (i ==n-1) { // last slice
+ // clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
+ // vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
+ // append->AddInput(mesh);
+ // }
+ }
+
+ // DEBUG: write points3D
+ clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
+
+ // Build the final 3D mesh form the list 2D mesh
+ append->Update();
+ vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
+
+ // Debug, write the mesh
+ /*
+ vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
+ w->SetInput(mesh);
+ w->SetFileName("bidon.vtk");
+ w->Write();
+ */
+
+ // Compute a single binary 3D image from the list of contours
+ clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
+ clitk::MeshToBinaryImageFilter<MaskImageType>::New();
+ filter->SetMesh(mesh);
+ filter->SetLikeImage(support);
+ filter->Update();
+ binarizedContour = filter->GetOutput();
+
+ // Crop
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
+ GetForegroundValue(), 2, true, p);
+ inf = p[2];
+ clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
+ GetForegroundValue(), 2, false, p);
+ sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
+ binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
+ false, GetBackgroundValue());
+
+ // Update the AFDB
+ writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
+ GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
+ 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<MaskImageType> 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<MaskImageType>(m_Working_Support);
+ //m_ListOfStations["2R"] = m_Working_Support;
+ //m_ListOfStations["2L"] = m_Working_Support;
+}
+//--------------------------------------------------------------------
+
+
+
#endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX
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);
for(uint i=0; i<mArgsInfo.station_given; i++)
f->AddComputeStation(mArgsInfo.station_arg[i]);
+ // Station 3A
+ f->SetFuzzyThreshold("3A", "SVC", mArgsInfo.S3A_ft_SVC_arg);
+ f->SetFuzzyThreshold("3A", "Bones", mArgsInfo.S3A_ft_Bones_arg);
+ f->SetThreshold("3A", "Bones", mArgsInfo.S3A_t_Bones_arg);
+ f->SetFuzzyThreshold("3A", "Aorta", mArgsInfo.S3A_ft_Aorta_arg);
+ f->SetFuzzyThreshold("3A", "SubclavianArtery", mArgsInfo.S3A_ft_SubclavianArtery_arg);
+
// Station 7
f->SetFuzzyThreshold("7", "Bronchi", mArgsInfo.S7_ft_Bronchi_arg);
f->SetFuzzyThreshold("7", "LeftSuperiorPulmonaryVein", mArgsInfo.S7_ft_LeftSuperiorPulmonaryVein_arg);
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);
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
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})
=========================================================================*/
-#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();
}
// 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; i<s->GetListOfROI().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();
}
filter.SetOutputImageFilename(n);
filter.Update();
- }
+ }*/
}
// This is the end my friend
--- /dev/null
+/*=========================================================================
+ 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<PixelType, 3> ImageType;
+ ImageType::Pointer input = clitk::readImage<ImageType>(args_info.input_arg, true);
+
+ // Create a filter to convert image into dicomRTStruct
+ clitk::Image2DicomRTStructFilter<PixelType> 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;
+}
+//--------------------------------------------------------------------
--- /dev/null
+# 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
+
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);
//------------------------------------------------------------------------------
void vvToolStructureSetManager::UpdateStructureSetInTreeWidget(int index, clitk::DicomRT_StructureSet * s) {
// Insert ROI
+ /*
const std::vector<clitk::DicomRT_ROI::Pointer> & rois = s->GetListOfROI();
for(unsigned int i=0; i<rois.size(); i++) {
if (mMapROIToTreeWidget.find(rois[i]) == mMapROIToTreeWidget.end())
AddRoiInTreeWidget(rois[i], mTree); // replace mTree with ss if several SS
}
+ */
+ clitk::DicomRT_StructureSet::ROIConstIteratorType iter;
+ for(iter = s->GetROIs().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
+ }
}
//------------------------------------------------------------------------------
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 (n<mDefaultLUTColor->GetNumberOfTableValues ()) {
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
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff