// 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("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();
+ }
+}
+//--------------------------------------------------------------------
+
+
//--------------------------------------------------------------------
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 ImageType>
+template <class TImageType>
void
-clitk::ExtractLymphStationsFilter<ImageType>::
-ExtractStation_3A_SI_Limits()
+clitk::ExtractLymphStationsFilter<TImageType>::
+ExtractStation_3A()
{
- // Apex of the chest or Sternum & Carina.
- StartNewStep("[Station 3A] Inf/Sup limits with Sternum and Carina");
-
- // Get Carina position (has been determined in Station8)
- m_CarinaZ = GetAFDB()->GetDouble("CarinaZ");
+ if (!CheckForStation("3A")) return;
+
+ // 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);
- // Get Sternum and search for the upper position
- MaskImagePointer Sternum = GetAFDB()->template GetImage<MaskImageType>("Sternum");
+ ExtractStation_3A_AntPost_S5();
+ ExtractStation_3A_AntPost_S6();
+ ExtractStation_3A_AntPost_Superiorly();
- // 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
+ ExtractStation_3A_Remove_Structures();
- //* Crop support :
- m_Working_Support =
- clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
- m_CarinaZ, m_SternumZ, true,
- GetBackgroundValue());
+ 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"); ?
- StopCurrentStep<MaskImageType>(m_Working_Support);
- m_ListOfStations["3A"] = m_Working_Support;
+
+ // ExtractStation_3A_Ant_Limits(); --> No, already in support; to remove
+ // ExtractStation_3A_Post_Limits(); --> No, more complex, see Vessels etc
+
+ // Store image filenames into AFDB
+ writeImage<MaskImageType>(m_ListOfStations["3A"], "seg/Station3A.mhd");
+ GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd");
+ WriteAFDB();
}
//--------------------------------------------------------------------
StartNewStep("[Station 3A] Post limits with SubclavianArtery");
// Get Sternum, keep posterior part.
- MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
+ MaskImagePointer SubclavianArteryLeft =
+ GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryLeft");
+ MaskImagePointer SubclavianArteryRight =
+ GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryRight");
+
m_Working_Support =
- clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, SubclavianArtery, 2,
+ clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, SubclavianArteryLeft, 2,
GetFuzzyThreshold("3A", "SubclavianArtery"), "AntTo",
false, 3, true, false);
+ m_Working_Support =
+ clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, SubclavianArteryRight, 2,
+ GetFuzzyThreshold("3A", "SubclavianArtery"), "AntTo",
+ false, 3, true, false);
+ StopCurrentStep<MaskImageType>(m_Working_Support);
+ m_ListOfStations["3A"] = m_Working_Support;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLymphStationsFilter<ImageType>::
+ExtractStation_3A_AntPost_S5()
+{
+ 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.
+ /*
+ typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
+ typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
+ relPosFilter->VerboseStepFlagOff();
+ relPosFilter->WriteStepFlagOff();
+ relPosFilter->SetBackgroundValue(GetBackgroundValue());
+ relPosFilter->SetInput(m_Working_Support);
+ relPosFilter->SetInputObject(SVC);
+ relPosFilter->AddOrientationTypeString("PostTo");
+ relPosFilter->SetInverseOrientationFlag(true);
+ relPosFilter->SetIntermediateSpacing(4);
+ relPosFilter->SetIntermediateSpacingFlag(false);
+ relPosFilter->SetFuzzyThreshold(0.5);
+ // relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
+ relPosFilter->Update();
+ m_Working_Support = relPosFilter->GetOutput();
+ */
+
+ // 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);
+
+ // 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());
+
+ // 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");
+
+ // 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");
+
+ // 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::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_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");
+ // 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
+ 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 (centroids.size() > 1) {
+ // Only keep the one most post
+ typename MaskSliceType::PixelType label;
+ if (centroids[1][1] > centroids[2][1]) {
+ label = 1;
+ }
+ else {
+ label = 2;
+ }
+
+ HERE
+
+ slices_BCV[i] = clitk::SetBackground<MaskSliceType>(slices_BCV[i], slices_BCV[i],
+ label,
+ GetBackgroundValue(), true);
+ }
+ // Remove from the support
+ slices[i] = clitk::AndNot(slices[i], slices_BCV[i], GetBackgroundValue());
+ }
+
+ // Joint
+ m_Working_Support = clitk::JoinSlices<MaskImageType>(slices, m_Working_Support, 2);
+
+ 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());
-
- 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();
}
//--------------------------------------------------------------------
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 =
// 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_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);
+ 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;
// 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_SI_Limits();
+ void ExtractStation_3A_Ant_Limits();
+ void ExtractStation_3A_Post_Limits();
// 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");
ExtractStation_3P();
StopSubStep();
- // Extract Station2RL
- /*
- StartNewStep("Station 2RL");
- StartSubStep();
- ExtractStation_2RL();
- StopSubStep();
- */
-
// Extract Station3A
StartNewStep("Station 3A");
StartSubStep();
ExtractStation_3A();
StopSubStep();
+
+ // HERE
+
+ // Extract Station2RL
+ StartNewStep("Station 2RL");
+ StartSubStep();
+ ExtractStation_2RL();
+ StopSubStep();
+
// Extract Station7
StartNewStep("Station 7");
StartSubStep();
//--------------------------------------------------------------------
+
//--------------------------------------------------------------------
-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>
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", "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);
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff