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