segmentation/clitkExtractLymphStation_2RL.txx

   1
   2
   3 // vtk
   4 #include <vtkAppendPolyData.h>
   5 #include <vtkPolyDataWriter.h>
   6 #include <vtkCellArray.h>
   7
   8 // clitk
   9 #include "clitkMeshToBinaryImageFilter.h"
  10
  11 // itk
  12 #include <itkImageDuplicator.h>
  13
  14 //--------------------------------------------------------------------
  15 template<class PointType>
  16 class comparePointsX {
  17 public:
  18   bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
  19 };
  20 //--------------------------------------------------------------------
  21
  22
  23 //--------------------------------------------------------------------
  24 template<class PairType>
  25 class comparePointsWithAngle {
  26 public:
  27   bool operator() (PairType i, PairType j) { return (i.second < j.second); }
  28 };
  29 //--------------------------------------------------------------------
  30
  31
  32 //--------------------------------------------------------------------
  33 template<int Dim>
  34 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
  35   std::vector<itk::Point<double, Dim-1> > previous;
  36   HypercubeCorners<Dim-1>(previous);
  37   out.resize(previous.size()*2);
  38   for(unsigned int i=0; i<out.size(); i++) {
  39     itk::Point<double, Dim> p;
  40     if (i<previous.size()) p[0] = 0;
  41     else p[0] = 1;
  42     for(int j=0; j<Dim-1; j++)
  43       {
  44         p[j+1] = previous[i%previous.size()][j];
  45       }
  46     out[i] = p;
  47   }
  48 }
  49 //--------------------------------------------------------------------
  50
  51
  52 //--------------------------------------------------------------------
  53 template<>
  54 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
  55   out.resize(2);
  56   out[0][0] = 0;
  57   out[1][0] = 1;
  58 }
  59 //--------------------------------------------------------------------
  60
  61
  62 //--------------------------------------------------------------------
  63 template<class ImageType>
  64 void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
  65                                        std::vector<typename ImageType::PointType> & bounds)
  66 {
  67   // Get image max/min coordinates
  68   const unsigned int dim=ImageType::ImageDimension;
  69   typedef typename ImageType::PointType PointType;
  70   typedef typename ImageType::IndexType IndexType;
  71   PointType min_c, max_c;
  72   IndexType min_i, max_i;
  73   min_i = image->GetLargestPossibleRegion().GetIndex();
  74   for(unsigned int i=0; i<dim; i++)
  75     max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
  76   image->TransformIndexToPhysicalPoint(min_i, min_c);
  77   image->TransformIndexToPhysicalPoint(max_i, max_c);
  78
  79   // Get corners coordinates
  80   HypercubeCorners<ImageType::ImageDimension>(bounds);
  81   for(unsigned int i=0; i<bounds.size(); i++) {
  82     for(unsigned int j=0; j<dim; j++) {
  83       if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
  84       if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
  85     }
  86   }
  87 }
  88 //--------------------------------------------------------------------
  89
  90
  91 //--------------------------------------------------------------------
  92 template <class ImageType>
  93 void
  94 clitk::ExtractLymphStationsFilter<ImageType>::
  95 ExtractStation_2RL_SetDefaultValues()
  96 {
  97   SetFuzzyThreshold("2RL", "CommonCarotidArtery", 0.7);
  98   SetFuzzyThreshold("2RL", "BrachioCephalicArtery", 0.7);
  99   SetFuzzyThreshold("2RL", "BrachioCephalicVein", 0.3);
 100   SetFuzzyThreshold("2RL", "Aorta", 0.7);
 101   SetFuzzyThreshold("2RL", "SubclavianArteryRight", 0.5);
 102   SetFuzzyThreshold("2RL", "SubclavianArteryLeft", 0.8);
 103 }
 104 //--------------------------------------------------------------------
 105
 106
 107 //--------------------------------------------------------------------
 108 template <class ImageType>
 109 void
 110 clitk::ExtractLymphStationsFilter<ImageType>::
 111 ExtractStation_2RL_SI_Limits()
 112 {
 113   // Apex of the chest or Sternum & Carina.
 114   StartNewStep("[Station 2RL] Inf/Sup limits with Sternum and TopOfAorticArch/CaudalMarginOfLeftBrachiocephalicVein");
 115
 116   /* Rod says: "For the inferior border, unlike in Atlas – UM, there
 117    is now a difference between 2R and 2L.  2R stops at the
 118    intersection of the caudal margin of the innominate vein with the
 119    trachea.  2L extends less inferiorly to the superior border of the
 120    aortic arch." */
 121
 122   /* Get TopOfAorticArch and CaudalMarginOfLeftBrachiocephalicVein
 123      - TopOfAorticArch -> can be obtain from Aorta -> most sup part.
 124      - CaudalMarginOfLeftBrachiocephalicVein -> must inf part of BrachicephalicVein
 125    */
 126   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 127   MaskImagePointType p = Aorta->GetOrigin(); // initialise to avoid warning
 128   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
 129   double TopOfAorticArchZ=p[2];
 130   GetAFDB()->SetDouble("TopOfAorticArchZ", TopOfAorticArchZ);
 131
 132   MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
 133   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
 134   double CaudalMarginOfLeftBrachiocephalicVeinZ=p[2];
 135   GetAFDB()->SetDouble("CaudalMarginOfLeftBrachiocephalicVeinZ", CaudalMarginOfLeftBrachiocephalicVeinZ);
 136
 137   // First, cut on the most inferior part. Add one slice because this
 138   // inf slice should not be included.
 139   double inf = std::min(CaudalMarginOfLeftBrachiocephalicVeinZ, TopOfAorticArchZ) + m_Working_Support->GetSpacing()[2];
 140
 141   // Get Sternum and search for the upper position
 142   MaskImagePointer Sternum = GetAFDB()->template GetImage<MaskImageType>("Sternum");
 143
 144   // Search most sup point
 145   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Sternum, GetBackgroundValue(), 2, false, p);
 146   double m_SternumZ = p[2];
 147   // +Sternum->GetSpacing()[2]; // One more slice, because it is below this point // NOT HERE ?? WHY DIFFERENT FROM 3A ??
 148
 149   //* Crop support :
 150   m_Working_Support =
 151     clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
 152                                                 inf, m_SternumZ, true,
 153                                                 GetBackgroundValue());
 154
 155   StopCurrentStep<MaskImageType>(m_Working_Support);
 156   m_ListOfStations["2R"] = m_Working_Support;
 157   m_ListOfStations["2L"] = m_Working_Support;
 158 }
 159 //--------------------------------------------------------------------
 160
 161
 162 //--------------------------------------------------------------------
 163 template <class ImageType>
 164 void
 165 clitk::ExtractLymphStationsFilter<ImageType>::
 166 ExtractStation_2RL_Ant_Limits()
 167 {
 168   // -----------------------------------------------------
 169   /* Rod says: "The anterior border, as with the Atlas – UM, is
 170     posterior to the vessels (right subclavian vein, left
 171     brachiocephalic vein, right brachiocephalic vein, left subclavian
 172     artery, left common carotid artery and brachiocephalic trunk).
 173     These vessels are not included in the nodal station.  The anterior
 174     border is drawn to the midpoint of the vessel and an imaginary
 175     line joins the middle of these vessels.  Between the vessels,
 176     station 2 is in contact with station 3a." */
 177
 178   // -----------------------------------------------------
 179   StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery");
 180
 181   // Get CommonCarotidArtery
 182   MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
 183
 184   // Remove Ant to CommonCarotidArtery
 185   typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
 186   typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 187   sliceRelPosFilter->VerboseStepFlagOff();
 188   sliceRelPosFilter->WriteStepFlagOff();
 189   sliceRelPosFilter->SetInput(m_Working_Support);
 190   sliceRelPosFilter->SetInputObject(CommonCarotidArtery);
 191   sliceRelPosFilter->SetDirection(2);
 192   sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "CommonCarotidArtery"));
 193   sliceRelPosFilter->AddOrientationTypeString("NotAntTo");
 194   sliceRelPosFilter->IntermediateSpacingFlagOn();
 195   sliceRelPosFilter->SetIntermediateSpacing(2);
 196   sliceRelPosFilter->UniqueConnectedComponentBySliceOff();
 197   sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff();
 198   sliceRelPosFilter->AutoCropFlagOn();
 199   sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 200   sliceRelPosFilter->RemoveObjectFlagOff();
 201   sliceRelPosFilter->Update();
 202   m_Working_Support = sliceRelPosFilter->GetOutput();
 203
 204   // End
 205   StopCurrentStep<MaskImageType>(m_Working_Support);
 206   m_ListOfStations["2R"] = m_Working_Support;
 207   m_ListOfStations["2L"] = m_Working_Support;
 208
 209   // -----------------------------------------------------
 210   // Remove Ant to H line from the Ant most part of the
 211   // CommonCarotidArtery until we reach the first slice of
 212   // BrachioCephalicArtery
 213   StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery, H line");
 214
 215   // First, find the first slice of BrachioCephalicArtery
 216   MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
 217   MaskImagePointType p = BrachioCephalicArtery->GetOrigin(); // initialise to avoid warning
 218   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicArtery, GetBackgroundValue(), 2, false, p);
 219   double TopOfBrachioCephalicArteryZ=p[2] + BrachioCephalicArtery->GetSpacing()[2]; // Add one slice
 220
 221   // Remove CommonCarotidArtery below this point
 222   CommonCarotidArtery = clitk::CropImageRemoveLowerThan<MaskImageType>(CommonCarotidArtery, 2, TopOfBrachioCephalicArteryZ, true, GetBackgroundValue());
 223
 224   // Find most Ant points
 225   std::vector<MaskImagePointType> ccaAntPositionA;
 226   std::vector<MaskImagePointType> ccaAntPositionB;
 227   clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(CommonCarotidArtery,
 228                                                                                GetBackgroundValue(), 2,
 229                                                                                1, true, // Ant
 230                                                                                0, // Horizontal line
 231                                                                                -3, // margin
 232                                                                                ccaAntPositionA,
 233                                                                                ccaAntPositionB);
 234   // Cut ant to this line
 235   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 236                                                                     ccaAntPositionA,
 237                                                                     ccaAntPositionB,
 238                                                                     GetBackgroundValue(), 1, 10);
 239
 240   // End
 241   StopCurrentStep<MaskImageType>(m_Working_Support);
 242   m_ListOfStations["2R"] = m_Working_Support;
 243   m_ListOfStations["2L"] = m_Working_Support;
 244
 245   // -----------------------------------------------------
 246   // Ant limit with the BrachioCephalicArtery
 247   StartNewStep("[Station 2RL] Ant limits with BrachioCephalicArtery line");
 248
 249   // Remove Ant to BrachioCephalicArtery
 250   m_Working_Support =
 251     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, BrachioCephalicArtery, 2,
 252                                                        GetFuzzyThreshold("2RL", "BrachioCephalicArtery"), "NotAntTo", false, 2, true, false);
 253   // End
 254   StopCurrentStep<MaskImageType>(m_Working_Support);
 255   m_ListOfStations["2R"] = m_Working_Support;
 256   m_ListOfStations["2L"] = m_Working_Support;
 257
 258   // -----------------------------------------------------
 259   // Ant limit with the BrachioCephalicArtery H line
 260   StartNewStep("[Station 2RL] Ant limits with BrachioCephalicArtery, Horizontal line");
 261
 262   // Find most Ant points
 263   std::vector<MaskImagePointType> bctAntPositionA;
 264   std::vector<MaskImagePointType> bctAntPositionB;
 265   clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(BrachioCephalicArtery,
 266                                                                                GetBackgroundValue(), 2,
 267                                                                                1, true, // Ant
 268                                                                                0, // Horizontal line
 269                                                                                -1, // margin
 270                                                                                bctAntPositionA,
 271                                                                                bctAntPositionB);
 272   // Cut ant to this line
 273   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 274                                                                     bctAntPositionA,
 275                                                                     bctAntPositionB,
 276                                                                     GetBackgroundValue(), 1, 10);
 277  // End
 278   StopCurrentStep<MaskImageType>(m_Working_Support);
 279   m_ListOfStations["2R"] = m_Working_Support;
 280   m_ListOfStations["2L"] = m_Working_Support;
 281
 282   // -----------------------------------------------------
 283   // Ant limit with the BrachioCephalicVein
 284   StartNewStep("[Station 2RL] Ant limits with BrachioCephalicVein");
 285
 286   // Remove Ant to BrachioCephalicVein
 287   MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
 288   m_Working_Support =
 289     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, BrachioCephalicVein, 2,
 290                                                        GetFuzzyThreshold("2RL", "BrachioCephalicVein"), "NotAntTo", false, 2, true, false);
 291   // End
 292   StopCurrentStep<MaskImageType>(m_Working_Support);
 293   m_ListOfStations["2R"] = m_Working_Support;
 294   m_ListOfStations["2L"] = m_Working_Support;
 295 }
 296 //--------------------------------------------------------------------
 297
 298
 299 //--------------------------------------------------------------------
 300 template <class ImageType>
 301 void
 302 clitk::ExtractLymphStationsFilter<ImageType>::
 303 ExtractStation_2RL_Ant_Limits2()
 304 {
 305   // -----------------------------------------------------
 306   /* Rod says: "The anterior border, as with the Atlas – UM, is
 307     posterior to the vessels (right subclavian vein, left
 308     brachiocephalic vein, right brachiocephalic vein, left subclavian
 309     artery, left common carotid artery and brachiocephalic trunk).
 310     These vessels are not included in the nodal station.  The anterior
 311     border is drawn to the midpoint of the vessel and an imaginary
 312     line joins the middle of these vessels.  Between the vessels,
 313     station 2 is in contact with station 3a." */
 314
 315   // -----------------------------------------------------
 316   StartNewStep("[Station 2RL] Ant limits with vessels centroids");
 317
 318   /* Here, we consider the vessels form a kind of anterior barrier. We
 319      link all vessels centroids and remove what is post to it.
 320     - select the list of structure
 321             vessel1 = BrachioCephalicArtery
 322             vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
 323             vessel3 = CommonCarotidArtery
 324             vessel4 = SubclavianArtery
 325             other   = Thyroid
 326             other   = Aorta
 327      - crop images as needed
 328      - slice by slice, choose the CCL and extract centroids
 329      - slice by slice, sort according to polar angle wrt Trachea centroid.
 330      - slice by slice, link centroids and close contour
 331      - remove outside this contour
 332      - merge with support
 333   */
 334
 335   // Read structures
 336   MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
 337   MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
 338   MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
 339   MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
 340   MaskImagePointer Thyroid = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
 341   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 342   MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
 343
 344   // Resize all structures like support
 345   BrachioCephalicArtery =
 346     clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, m_Working_Support, GetBackgroundValue());
 347   CommonCarotidArtery =
 348     clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, m_Working_Support, GetBackgroundValue());
 349   SubclavianArtery =
 350     clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, m_Working_Support, GetBackgroundValue());
 351   Thyroid =
 352     clitk::ResizeImageLike<MaskImageType>(Thyroid, m_Working_Support, GetBackgroundValue());
 353   Aorta =
 354     clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
 355   BrachioCephalicVein =
 356     clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, m_Working_Support, GetBackgroundValue());
 357   Trachea =
 358     clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
 359
 360   // Extract slices
 361   std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
 362   clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
 363   std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
 364   clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
 365   std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
 366   clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
 367   std::vector<MaskSlicePointer> slices_SubclavianArtery;
 368   clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
 369   std::vector<MaskSlicePointer> slices_Thyroid;
 370   clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
 371   std::vector<MaskSlicePointer> slices_Aorta;
 372   clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
 373   std::vector<MaskSlicePointer> slices_Trachea;
 374   clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
 375   unsigned int n= slices_BrachioCephalicArtery.size();
 376
 377   // Get the boundaries of one slice
 378   std::vector<MaskSlicePointType> bounds;
 379   ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
 380
 381   // For all slices, for all structures, find the centroid and build the contour
 382   // List of 3D points (for debug)
 383   std::vector<MaskImagePointType> p3D;
 384
 385   vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
 386   for(unsigned int i=0; i<n; i++) {
 387     // Labelize the slices
 388     slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
 389                                                             GetBackgroundValue(), true, 1);
 390     slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
 391                                                          GetBackgroundValue(), true, 1);
 392     slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
 393                                                              GetBackgroundValue(), true, 1);
 394     slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
 395                                                             GetBackgroundValue(), true, 1);
 396     slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
 397                                                 GetBackgroundValue(), true, 1);
 398     slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
 399                                               GetBackgroundValue(), true, 1);
 400
 401     // Search centroids
 402     std::vector<MaskSlicePointType> points2D;
 403     std::vector<MaskSlicePointType> centroids1;
 404     std::vector<MaskSlicePointType> centroids2;
 405     std::vector<MaskSlicePointType> centroids3;
 406     std::vector<MaskSlicePointType> centroids4;
 407     std::vector<MaskSlicePointType> centroids5;
 408     std::vector<MaskSlicePointType> centroids6;
 409     ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
 410     ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
 411     ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
 412     ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
 413     ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
 414     ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
 415
 416     // BrachioCephalicVein -> when it is separated into two CCL, we
 417     // only consider the most at Right one
 418     if (centroids6.size() > 2) {
 419       if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
 420       else centroids6.erase(centroids6.begin()+1);
 421     }
 422
 423     // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
 424     // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
 425     if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
 426       centroids6.clear();
 427     }
 428
 429     for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
 430     for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
 431     for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
 432     for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
 433     for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
 434     for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
 435
 436     // Sort by angle according to trachea centroid and vertical line,
 437     // in polar coordinates :
 438     // http://en.wikipedia.org/wiki/Polar_coordinate_system
 439     std::vector<MaskSlicePointType> centroids_trachea;
 440     ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
 441     typedef std::pair<MaskSlicePointType, double> PointAngleType;
 442     std::vector<PointAngleType> angles;
 443     for(unsigned int j=0; j<points2D.size(); j++) {
 444       //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
 445       double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
 446       double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
 447       double angle = 0;
 448       if (x>0) angle = atan(y/x);
 449       if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
 450       if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
 451       if (x==0) {
 452         if (y>0) angle = M_PI/2.0;
 453         if (y<0) angle = -M_PI/2.0;
 454         if (y==0) angle = 0;
 455       }
 456       angle = clitk::rad2deg(angle);
 457       // Angle is [-180;180] wrt the X axis. We change the X axis to
 458       // be the vertical line, because we want to sort from Right to
 459       // Left from Post to Ant.
 460       if (angle>0)
 461         angle = (270-angle);
 462       if (angle<0) {
 463         angle = -angle-90;
 464         if (angle<0) angle = 360-angle;
 465       }
 466       PointAngleType a;
 467       a.first = points2D[j];
 468       a.second = angle;
 469       angles.push_back(a);
 470     }
 471
 472     // Do nothing if less than 2 points --> n
 473     if (points2D.size() < 3) { //continue;
 474       continue;
 475     }
 476
 477     // Sort points2D according to polar angles
 478     std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
 479     for(unsigned int j=0; j<angles.size(); j++) {
 480       points2D[j] = angles[j].first;
 481     }
 482     //    DDV(points2D, points2D.size());
 483
 484     /* When vessels are far away, we try to replace the line segment
 485        with a curved line that join the two vessels but stay
 486        approximately at the same distance from the trachea centroids
 487        than the vessels.
 488
 489        For that:
 490        - let a and c be two successive vessels centroids
 491        - id distance(a,c) < threshold, next point
 492
 493        TODO HERE
 494
 495        - compute mid position between two successive points -
 496        compute dist to trachea centroid for the 3 pts - if middle too
 497        low, add one point
 498     */
 499     std::vector<MaskSlicePointType> toadd;
 500     unsigned int index = 0;
 501     double dmax = 5;
 502     while (index<points2D.size()-1) {
 503       MaskSlicePointType a = points2D[index];
 504       MaskSlicePointType c = points2D[index+1];
 505
 506       double dac = a.EuclideanDistanceTo(c);
 507       if (dac>dmax) {
 508
 509         MaskSlicePointType b;
 510         b[0] = a[0]+(c[0]-a[0])/2.0;
 511         b[1] = a[1]+(c[1]-a[1])/2.0;
 512
 513         // Compute distance to trachea centroid
 514         MaskSlicePointType m = centroids_trachea[1];
 515         double da = m.EuclideanDistanceTo(a);
 516         double db = m.EuclideanDistanceTo(b);
 517         //double dc = m.EuclideanDistanceTo(c);
 518
 519         // Mean distance, find point on the line from trachea centroid
 520         // to b
 521         double alpha = (da+db)/2.0;
 522         MaskSlicePointType v;
 523         double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
 524         v[0] = (b[0]-m[0])/n;
 525         v[1] = (b[1]-m[1])/n;
 526         MaskSlicePointType r;
 527         r[0] = m[0]+alpha*v[0];
 528         r[1] = m[1]+alpha*v[1];
 529         points2D.insert(points2D.begin()+index+1, r);
 530       }
 531       else {
 532         index++;
 533       }
 534     }
 535     //    DDV(points2D, points2D.size());
 536
 537     // Add some points to close the contour
 538     // - H line towards Right
 539     MaskSlicePointType p = points2D[0];
 540     p[0] = bounds[0][0];
 541     points2D.insert(points2D.begin(), p);
 542     // - corner Right/Post
 543     p = bounds[0];
 544     points2D.insert(points2D.begin(), p);
 545     // - H line towards Left
 546     p = points2D.back();
 547     p[0] = bounds[2][0];
 548     points2D.push_back(p);
 549     // - corner Right/Post
 550     p = bounds[2];
 551     points2D.push_back(p);
 552     // Close contour with the first point
 553     points2D.push_back(points2D[0]);
 554     //    DDV(points2D, points2D.size());
 555
 556     // Build 3D points from the 2D points
 557     std::vector<ImagePointType> points3D;
 558     clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, m_Working_Support, points3D);
 559     for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
 560
 561     // Build the mesh from the contour's points
 562     vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
 563     append->AddInput(mesh);
 564   }
 565
 566   // DEBUG: write points3D
 567   clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
 568
 569   // Build the final 3D mesh form the list 2D mesh
 570   append->Update();
 571   vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
 572
 573   // Debug, write the mesh
 574   /*
 575     vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
 576     w->SetInput(mesh);
 577     w->SetFileName("bidon.vtk");
 578     w->Write();
 579   */
 580
 581   // Compute a single binary 3D image from the list of contours
 582   clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
 583     clitk::MeshToBinaryImageFilter<MaskImageType>::New();
 584   filter->SetMesh(mesh);
 585   filter->SetLikeImage(m_Working_Support);
 586   filter->Update();
 587   MaskImagePointer binarizedContour = filter->GetOutput();
 588
 589   // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
 590   ImagePointType p = p3D[2]; // This is the first centroid of the first slice
 591   p[1] += 50; // 50 mm Post from this point must be kept
 592   ImageIndexType index;
 593   binarizedContour->TransformPhysicalPointToIndex(p, index);
 594   bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
 595
 596   // remove from support
 597   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
 598   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 599   boolFilter->InPlaceOn();
 600   boolFilter->SetInput1(m_Working_Support);
 601   boolFilter->SetInput2(binarizedContour);
 602   boolFilter->SetBackgroundValue1(GetBackgroundValue());
 603   boolFilter->SetBackgroundValue2(GetBackgroundValue());
 604   if (isInside)
 605     boolFilter->SetOperationType(BoolFilterType::And);
 606   else
 607     boolFilter->SetOperationType(BoolFilterType::AndNot);
 608   boolFilter->Update();
 609   m_Working_Support = boolFilter->GetOutput();
 610
 611   // End
 612   StopCurrentStep<MaskImageType>(m_Working_Support);
 613   m_ListOfStations["2R"] = m_Working_Support;
 614   m_ListOfStations["2L"] = m_Working_Support;
 615 }
 616 //--------------------------------------------------------------------
 617
 618
 619 //--------------------------------------------------------------------
 620 template <class ImageType>
 621 void
 622 clitk::ExtractLymphStationsFilter<ImageType>::
 623 ExtractStation_2RL_Post_Limits()
 624 {
 625   StartNewStep("[Station 2RL] Post limits with post wall of Trachea");
 626
 627   // Get Trachea
 628   MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
 629
 630   // Resize like the current support (to have the same number of slices)
 631   Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
 632
 633   // Find extrema post positions
 634   std::vector<MaskImagePointType> tracheaPostPositionsA;
 635   std::vector<MaskImagePointType> tracheaPostPositionsB;
 636   clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(Trachea,
 637                                                                                GetBackgroundValue(), 2,
 638                                                                                1, false, // Post
 639                                                                                0, // Horizontal line
 640                                                                                1,
 641                                                                                tracheaPostPositionsA,
 642                                                                                tracheaPostPositionsB);
 643   // Cut post to this line
 644   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 645                                                                     tracheaPostPositionsA,
 646                                                                     tracheaPostPositionsB,
 647                                                                     GetBackgroundValue(), 1, -10);
 648   // END
 649   StopCurrentStep<MaskImageType>(m_Working_Support);
 650   m_ListOfStations["2R"] = m_Working_Support;
 651   m_ListOfStations["2L"] = m_Working_Support;
 652 }
 653 //--------------------------------------------------------------------
 654
 655
 656 //--------------------------------------------------------------------
 657 // Build a vtk mesh from a list of slice number/closed-contours
 658 template <class ImageType>
 659 vtkSmartPointer<vtkPolyData>
 660 clitk::ExtractLymphStationsFilter<ImageType>::
 661 Build3DMeshFrom2DContour(const std::vector<ImagePointType> & points)
 662 {
 663   // create a contour, polydata.
 664   vtkSmartPointer<vtkPolyData> mesh = vtkSmartPointer<vtkPolyData>::New();
 665   mesh->Allocate(); //for cell structures
 666   mesh->SetPoints(vtkPoints::New());
 667   vtkIdType ids[2];
 668   int point_number = points.size();
 669   for (unsigned int i=0; i<points.size(); i++) {
 670     mesh->GetPoints()->InsertNextPoint(points[i][0],points[i][1],points[i][2]);
 671     ids[0]=i;
 672     ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0
 673     mesh->GetLines()->InsertNextCell(2,ids);
 674   }
 675   // Return
 676   return mesh;
 677 }
 678 //--------------------------------------------------------------------
 679
 680
 681 //--------------------------------------------------------------------
 682 template <class ImageType>
 683 void
 684 clitk::ExtractLymphStationsFilter<ImageType>::
 685 ExtractStation_2RL_LR_Limits()
 686 {
 687   // ---------------------------------------------------------------------------
 688   StartNewStep("[Station 2RL] Left/Right limits with Aorta");
 689   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 690   //  DD(GetFuzzyThreshold("2RL", "BrachioCephalicVein"));
 691   m_Working_Support =
 692     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, Aorta, 2,
 693                                                        GetFuzzyThreshold("2RL", "Aorta"),
 694                                                        "RightTo", false, 2, true, false);
 695   // END
 696   StopCurrentStep<MaskImageType>(m_Working_Support);
 697   m_ListOfStations["2R"] = m_Working_Support;
 698   m_ListOfStations["2L"] = m_Working_Support;
 699
 700   // ---------------------------------------------------------------------------
 701   StartNewStep("[Station 2RL] Left/Right limits with SubclavianArtery (Right)");
 702
 703   // SliceBySliceRelativePosition + select CCL most at Right
 704   MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
 705   typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
 706   typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 707   sliceRelPosFilter->VerboseStepFlagOff();
 708   sliceRelPosFilter->WriteStepFlagOff();
 709   sliceRelPosFilter->SetInput(m_Working_Support);
 710   sliceRelPosFilter->SetInputObject(SubclavianArtery);
 711   sliceRelPosFilter->SetDirection(2);
 712   sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "SubclavianArteryRight"));
 713   sliceRelPosFilter->AddOrientationTypeString("NotRightTo");
 714   sliceRelPosFilter->IntermediateSpacingFlagOn();
 715   sliceRelPosFilter->SetIntermediateSpacing(2);
 716   sliceRelPosFilter->UniqueConnectedComponentBySliceOff();
 717   sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff();
 718
 719   sliceRelPosFilter->CCLSelectionFlagOn(); // select one CCL by slice
 720   sliceRelPosFilter->SetCCLSelectionDimension(0); // select according to X (0) axis
 721   sliceRelPosFilter->SetCCLSelectionDirection(-1); // select most at Right
 722   sliceRelPosFilter->CCLSelectionIgnoreSingleCCLFlagOn(); // ignore if only one CCL
 723
 724   sliceRelPosFilter->AutoCropFlagOn();
 725   sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 726   sliceRelPosFilter->RemoveObjectFlagOff();
 727   sliceRelPosFilter->Update();
 728   m_Working_Support = sliceRelPosFilter->GetOutput();
 729
 730   // END
 731   StopCurrentStep<MaskImageType>(m_Working_Support);
 732   m_ListOfStations["2R"] = m_Working_Support;
 733   m_ListOfStations["2L"] = m_Working_Support;
 734
 735
 736   // ---------------------------------------------------------------------------
 737   StartNewStep("[Station 2RL] Left/Right limits with SubclavianArtery (Left)");
 738
 739   // SliceBySliceRelativePosition + select CCL most at Right
 740    sliceRelPosFilter = SliceRelPosFilterType::New();
 741   sliceRelPosFilter->VerboseStepFlagOff();
 742   sliceRelPosFilter->WriteStepFlagOff();
 743   sliceRelPosFilter->SetInput(m_Working_Support);
 744   sliceRelPosFilter->SetInputObject(SubclavianArtery);
 745   sliceRelPosFilter->SetDirection(2);
 746   sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "SubclavianArteryLeft"));
 747   sliceRelPosFilter->AddOrientationTypeString("NotLeftTo");
 748   sliceRelPosFilter->IntermediateSpacingFlagOn();
 749   sliceRelPosFilter->SetIntermediateSpacing(2);
 750   sliceRelPosFilter->UniqueConnectedComponentBySliceOff();
 751   sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff();
 752
 753   sliceRelPosFilter->CCLSelectionFlagOn(); // select one CCL by slice
 754   sliceRelPosFilter->SetCCLSelectionDimension(0); // select according to X (0) axis
 755   sliceRelPosFilter->SetCCLSelectionDirection(+1); // select most at Left
 756   sliceRelPosFilter->CCLSelectionIgnoreSingleCCLFlagOff(); // do not ignore if only one CCL
 757
 758   sliceRelPosFilter->AutoCropFlagOn();
 759   sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 760   sliceRelPosFilter->RemoveObjectFlagOff();
 761   sliceRelPosFilter->Update();
 762   m_Working_Support = sliceRelPosFilter->GetOutput();
 763
 764   // END
 765   StopCurrentStep<MaskImageType>(m_Working_Support);
 766   m_ListOfStations["2R"] = m_Working_Support;
 767   m_ListOfStations["2L"] = m_Working_Support;
 768 }
 769 //--------------------------------------------------------------------
 770
 771 //--------------------------------------------------------------------
 772 template <class ImageType>
 773 void
 774 clitk::ExtractLymphStationsFilter<ImageType>::
 775 ExtractStation_2RL_Remove_Structures()
 776 {
 777   Remove_Structures("2RL", "BrachioCephalicVein");
 778   Remove_Structures("2RL", "CommonCarotidArtery");
 779   Remove_Structures("2RL", "SubclavianArtery");
 780   Remove_Structures("2RL", "Thyroid");
 781   Remove_Structures("2RL", "Aorta");
 782
 783   // END
 784   StopCurrentStep<MaskImageType>(m_Working_Support);
 785   m_ListOfStations["2R"] = m_Working_Support;
 786   m_ListOfStations["2L"] = m_Working_Support;
 787 }
 788 //--------------------------------------------------------------------
 789
 790
 791 //--------------------------------------------------------------------
 792 template <class ImageType>
 793 void
 794 clitk::ExtractLymphStationsFilter<ImageType>::
 795 ExtractStation_2RL_SeparateRL()
 796 {
 797   // ---------------------------------------------------------------------------
 798   StartNewStep("[Station 2RL] Separate 2R/2L according to Trachea");
 799
 800   /*Rod says:
 801
 802    "For station 2 there is a shift, dividing 2R from 2L, from midline
 803    to the left paratracheal border."
 804
 805    Algo:
 806     - Consider Trachea SliceBySlice
 807     - find extrema at Left
 808     - add margins towards Right
 809     - remove what is at Left of this line
 810    */
 811
 812   // Get Trachea
 813   MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
 814
 815   // Resize like the current support (to have the same number of slices)
 816   Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
 817
 818   // Find extrema post positions
 819   std::vector<MaskImagePointType> tracheaLeftPositionsA;
 820   std::vector<MaskImagePointType> tracheaLeftPositionsB;
 821   clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(Trachea,
 822                                                                                GetBackgroundValue(), 2,
 823                                                                                0, false, // Left
 824                                                                                1, // Vertical line
 825                                                                                1, // margins
 826                                                                                tracheaLeftPositionsA,
 827                                                                                tracheaLeftPositionsB);
 828   // Copy support for R and L
 829   typedef itk::ImageDuplicator<MaskImageType> DuplicatorType;
 830   DuplicatorType::Pointer duplicator = DuplicatorType::New();
 831   duplicator->SetInputImage(m_Working_Support);
 832   duplicator->Update();
 833   MaskImageType::Pointer m_Working_Support2 = duplicator->GetOutput();
 834
 835   // Cut post to this line for Right part
 836   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 837                                                                     tracheaLeftPositionsA,
 838                                                                     tracheaLeftPositionsB,
 839                                                                     GetBackgroundValue(), 0, -10);
 840   writeImage<MaskImageType>(m_Working_Support, "R.mhd");
 841
 842   // Cut post to this line for Left part
 843   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support2,
 844                                                                     tracheaLeftPositionsA,
 845                                                                     tracheaLeftPositionsB,
 846                                                                     GetBackgroundValue(), 0, +10);
 847   writeImage<MaskImageType>(m_Working_Support2, "L.mhd");
 848
 849   // END
 850   StopCurrentStep<MaskImageType>(m_Working_Support);
 851   m_ListOfStations["2R"] = m_Working_Support;
 852   m_ListOfStations["2L"] = m_Working_Support2;
 853 }
 854 //--------------------------------------------------------------------