X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLymphStation_2RL.txx;h=260e50a42ef4ed97d153d3c8f61d05bfe6f260ac;hb=3c86758765bc9bcba20d439424bcf97091b5af6f;hp=7e8c834b4aaa6d1d52510f7676881922cf20251f;hpb=ea2fc85081a73d762517dc9663bcc6d5bf15241d;p=clitk.git

diff --git a/segmentation/clitkExtractLymphStation_2RL.txx b/segmentation/clitkExtractLymphStation_2RL.txx
index 7e8c834..260e50a 100644
--- a/segmentation/clitkExtractLymphStation_2RL.txx
+++ b/segmentation/clitkExtractLymphStation_2RL.txx
@@ -11,82 +11,6 @@
 // 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(uint 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 uint 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(uint 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(uint i=0; i<bounds.size(); i++) {
-    for(uint 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 
@@ -103,6 +27,35 @@ ExtractStation_2RL_SetDefaultValues()
 //--------------------------------------------------------------------
 
 
+//--------------------------------------------------------------------
+template <class TImageType>
+void 
+clitk::ExtractLymphStationsFilter<TImageType>::
+ExtractStation_2RL()
+{
+  if ((!CheckForStation("2R")) && (!CheckForStation("2L"))) return;
+
+  ExtractStation_2RL_SI_Limits();
+  ExtractStation_2RL_Post_Limits();
+  
+  ExtractStation_2RL_Ant_Limits2();
+  //ExtractStation_2RL_Ant_Limits(); 
+  
+  ExtractStation_2RL_LR_Limits(); 
+  ExtractStation_2RL_Remove_Structures(); 
+  ExtractStation_2RL_SeparateRL(); 
+  
+  // Store image filenames into AFDB 
+  writeImage<MaskImageType>(m_ListOfStations["2R"], "seg/Station2R.mhd");
+  writeImage<MaskImageType>(m_ListOfStations["2L"], "seg/Station2L.mhd");
+  GetAFDB()->SetImageFilename("Station2R", "seg/Station2R.mhd"); 
+  GetAFDB()->SetImageFilename("Station2L", "seg/Station2L.mhd"); 
+  WriteAFDB(); 
+
+}
+//--------------------------------------------------------------------
+
+
 //--------------------------------------------------------------------
 template <class ImageType>
 void 
@@ -299,299 +252,23 @@ ExtractStation_2RL_Ant_Limits()
 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);
-  uint 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(uint 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(uint j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
-    for(uint j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
-    for(uint j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
-    for(uint j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
-    for(uint j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
-    for(uint 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(uint 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(uint 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;
-    uint 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(uint 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(); 
@@ -600,10 +277,7 @@ ExtractStation_2RL_Ant_Limits2()
   boolFilter->SetInput2(binarizedContour);
   boolFilter->SetBackgroundValue1(GetBackgroundValue());
   boolFilter->SetBackgroundValue2(GetBackgroundValue());
-  if (isInside)
-    boolFilter->SetOperationType(BoolFilterType::And);
-  else
-    boolFilter->SetOperationType(BoolFilterType::AndNot);
+  boolFilter->SetOperationType(BoolFilterType::And);
   boolFilter->Update();
   m_Working_Support = boolFilter->GetOutput();