segmentation/clitkExtractLymphStation_8.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================*/
  18
  19 #include <itkBinaryDilateImageFilter.h>
  20 #include <itkMirrorPadImageFilter.h>
  21
  22 //--------------------------------------------------------------------
  23 template <class ImageType>
  24 void
  25 clitk::ExtractLymphStationsFilter<ImageType>::
  26 ExtractStation_8_SetDefaultValues()
  27 {
  28   SetDistanceMaxToAnteriorPartOfTheSpine(10);
  29   MaskImagePointType p;
  30   p[0] = 15; p[1] = 2; p[2] = 1;
  31   SetEsophagusDiltationForAnt(p);
  32   p[0] = 5; p[1] = 10; p[2] = 1;
  33   SetEsophagusDiltationForRight(p);
  34   SetFuzzyThresholdForS8(0.5);
  35   SetInjectedThresholdForS8(150);
  36 }
  37 //--------------------------------------------------------------------
  38
  39 //--------------------------------------------------------------------
  40 template <class ImageType>
  41 void
  42 clitk::ExtractLymphStationsFilter<ImageType>::
  43 ExtractStation_8_SI_Limits()
  44 {
  45   /*
  46     Station 8: paraeosphageal nodes
  47
  48     "Superiorly, Station 8 begins at the level of the carina and,
  49     therefore, abuts Station 3P (Fig. 3C). Inferiorly, the lower limit
  50     of Station 8 was not speciﬁed in the Mountain and Dresler
  51     classiﬁcation (1). We delineate this volume until it reaches the
  52     gastroesphogeal junction. "
  53
  54     => new classification IASCL 2009:
  55
  56     "Lower border: the diaphragm"
  57
  58     "Upper border: the upper border of the lower lobe bronchus on the
  59     left; the lower border of the bronchus intermedius on the right"
  60
  61   */
  62   StartNewStep("[Station8] Inf/Sup mediastinum limits with carina/diaphragm junction");
  63
  64   // Get Carina Z position
  65   MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
  66
  67   std::vector<MaskImagePointType> centroids;
  68   clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
  69   m_CarinaZ = centroids[1][2];
  70   // DD(m_CarinaZ);
  71   // add one slice to include carina ?
  72   m_CarinaZ += m_Mediastinum->GetSpacing()[2];
  73   // We dont need Carina structure from now
  74   Carina->Delete();
  75   GetAFDB()->SetDouble("CarinaZ", m_CarinaZ);
  76
  77   // Found most inferior part of the lung
  78   MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
  79   // It should be already croped, so I took the origin and add 10mm above
  80   m_DiaphragmInferiorLimit = Lungs->GetOrigin()[2]+10;
  81   //  Lungs->Delete(); // we don't need it, release memory -> it we want to release, also free in AFDB
  82   clitk::PrintMemory(GetVerboseMemoryFlag(), "after reading lungs");
  83   GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
  84   clitk::PrintMemory(GetVerboseMemoryFlag(), "after release lungs");
  85
  86   /* Crop support :
  87      Superior limit = carina
  88      Inferior limit = DiaphragmInferiorLimit (old=Gastroesphogeal junction) */
  89   m_Working_Support =
  90     clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
  91                                                 m_DiaphragmInferiorLimit, //GOjunctionZ,
  92                                                 m_CarinaZ, true,
  93                                                 GetBackgroundValue());
  94
  95   // Remove some structures that we know are excluded
  96   MaskImagePointer VertebralBody =
  97     GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
  98   MaskImagePointer Aorta =
  99     GetAFDB()->template GetImage <MaskImageType>("Aorta");
 100
 101   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
 102   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 103   boolFilter->InPlaceOn();
 104   boolFilter->SetInput1(m_Working_Support);
 105   boolFilter->SetInput2(VertebralBody);
 106   boolFilter->SetOperationType(BoolFilterType::AndNot);
 107   boolFilter->Update();
 108   /*
 109
 110     DO NOT REMOVE AORTA YET (LATER)
 111
 112     boolFilter->SetInput1(boolFilter->GetOutput());
 113     boolFilter->SetInput2(Aorta);
 114     boolFilter->SetOperationType(BoolFilterType::AndNot);
 115     boolFilter->Update();
 116   */
 117   m_Working_Support = boolFilter->GetOutput();
 118
 119   // Done.
 120   StopCurrentStep<MaskImageType>(m_Working_Support);
 121   m_ListOfStations["8"] = m_Working_Support;
 122 }
 123 //--------------------------------------------------------------------
 124
 125
 126 //--------------------------------------------------------------------
 127 template <class ImageType>
 128 void
 129 clitk::ExtractLymphStationsFilter<ImageType>::
 130 ExtractStation_8_Post_Limits()
 131 {
 132   /*
 133     Station 8: paraeosphageal nodes
 134
 135     Anteriorly, it is in contact with Station 7 and the
 136     left main stem bronchus in its superior aspect (Fig. 3C–H) and
 137     with the heart more inferiorly. Posteriorly, Station 8 abuts the
 138     descending aorta and the anterior aspect of the vertebral body
 139     until an imaginary horizontal line running 1 cm posterior to the
 140     anterior border of the vertebral body (Fig. 3C).
 141
 142     New classification IASCL 2009 :
 143
 144     "Nodes lying adjacent to the wall of the esophagus and to the
 145     right or left of the midline, excluding subcarinal nodes (S7)
 146     Upper"
 147
 148   */
 149
 150   StartNewStep("[Station8] Post limits with vertebral body");
 151   MaskImagePointer VertebralBody =
 152     GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
 153
 154   // Consider vertebral body slice by slice
 155   typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
 156   typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 157   typedef typename ExtractSliceFilterType::SliceType SliceType;
 158   std::vector<typename SliceType::Pointer> vertebralSlices;
 159   extractSliceFilter->SetInput(VertebralBody);
 160   extractSliceFilter->SetDirection(2);
 161   extractSliceFilter->Update();
 162   extractSliceFilter->GetOutputSlices(vertebralSlices);
 163
 164   // For each slice, compute the most anterior point
 165   std::map<int, typename SliceType::PointType> vertebralAntPositionBySlice;
 166   for(uint i=0; i<vertebralSlices.size(); i++) {
 167     typename SliceType::PointType p;
 168     bool found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vertebralSlices[i],
 169                                                                      GetBackgroundValue(),
 170                                                                      1, true, p);
 171     if (found) {
 172       vertebralAntPositionBySlice[i] = p;
 173     }
 174     else {
 175       // no Foreground in this slice (it appends generally at the
 176       // beginning and the end of the volume). Do nothing in this
 177       // case.
 178     }
 179   }
 180
 181   // Convert 2D points in slice into 3D points
 182   std::vector<MaskImagePointType> vertebralAntPositions;
 183   clitk::PointsUtils<MaskImageType>::Convert2DTo3DList(vertebralAntPositionBySlice,
 184                                                        VertebralBody,
 185                                                        vertebralAntPositions);
 186
 187   // DEBUG : write list of points
 188   clitk::WriteListOfLandmarks<MaskImageType>(vertebralAntPositions,
 189                                              "S8-vertebralMostAntPositions-points.txt");
 190
 191   // Cut support posteriorly 1cm the most anterior point of the
 192   // VertebralBody. Do nothing below and above the VertebralBody. To
 193   // do that compute several region, slice by slice and fill.
 194   MaskImageRegionType region;
 195   MaskImageSizeType size;
 196   MaskImageIndexType start;
 197   size[2] = 1; // a single slice
 198   start[0] = m_Working_Support->GetLargestPossibleRegion().GetIndex()[0];
 199   size[0] = m_Working_Support->GetLargestPossibleRegion().GetSize()[0];
 200   for(uint i=0; i<vertebralAntPositions.size(); i++) {
 201     typedef typename itk::ImageRegionIterator<MaskImageType> IteratorType;
 202     IteratorType iter =
 203       IteratorType(m_Working_Support, m_Working_Support->GetLargestPossibleRegion());
 204     MaskImageIndexType index;
 205     // Consider some cm posterior to most anterior positions (usually
 206     // 1 cm).
 207     vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheSpine();
 208     // Get index of this point
 209     m_Working_Support->TransformPhysicalPointToIndex(vertebralAntPositions[i], index);
 210     // Compute region (a single slice)
 211     start[2] = index[2];
 212     start[1] = m_Working_Support->GetLargestPossibleRegion().GetIndex()[1]+index[1];
 213     size[1] = m_Working_Support->GetLargestPossibleRegion().GetSize()[1]-start[1];
 214     region.SetSize(size);
 215     region.SetIndex(start);
 216     // Fill region
 217     if (m_Working_Support->GetLargestPossibleRegion().IsInside(start))  {
 218       itk::ImageRegionIterator<MaskImageType> it(m_Working_Support, region);
 219       it.GoToBegin();
 220       while (!it.IsAtEnd()) {
 221         it.Set(GetBackgroundValue());
 222         ++it;
 223       }
 224     }
 225   }
 226
 227   StopCurrentStep<MaskImageType>(m_Working_Support);
 228   m_ListOfStations["8"] = m_Working_Support;
 229 }
 230 //--------------------------------------------------------------------
 231
 232
 233 //--------------------------------------------------------------------
 234 template <class ImageType>
 235 void
 236 clitk::ExtractLymphStationsFilter<ImageType>::
 237 ExtractStation_8_Ant_Sup_Limits()
 238 {
 239   //--------------------------------------------------------------------
 240   StartNewStep("[Station8] Ant limits with S7 above Carina");
 241   /*
 242     Anteriorly, it is in contact with Station 7 and the
 243     left main stem bronchus in its superior aspect (Fig. 3C–H) and
 244     with the heart more inferiorly.
 245
 246     1) line post wall bronchi (S7), till originRMLB
 247     - LUL bronchus : to detect in trachea. But not needed here ??
 248     2) heart ! -> not delineated.
 249     ==> below S7, crop CT not to far away (ant), then try with threshold
 250
 251     1) ==> how to share with S7 ? Prepare both support at the same time !
 252     NEED vector of initial support for each station ? No -> map if it exist before, take it !!
 253
 254   */
 255
 256   // Ant limit from carina (start) to end of S7 = originRMLB
 257   // Get Trachea
 258   MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
 259
 260   MaskImagePointer m_Working_Trachea =
 261     clitk::CropImageAbove<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
 262                                          GetBackgroundValue());
 263
 264   // Seprate into two main bronchi
 265   MaskImagePointer RightBronchus;
 266   MaskImagePointer LeftBronchus;
 267
 268   // Labelize and consider the two first (main) labels
 269   m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
 270
 271   // Carina position must at the first slice that separate the two
 272   // main bronchus (not superiorly). We thus first check that the
 273   // upper slice is composed of at least two labels
 274   typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
 275   SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
 276   iter.SetFirstDirection(0);
 277   iter.SetSecondDirection(1);
 278   iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
 279   int maxLabel=0;
 280   while (!iter.IsAtReverseEndOfSlice()) {
 281     while (!iter.IsAtReverseEndOfLine()) {
 282       if (iter.Get() > maxLabel) maxLabel = iter.Get();
 283       --iter;
 284     }
 285     iter.PreviousLine();
 286   }
 287   if (maxLabel < 2) {
 288     clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
 289   }
 290
 291   // Compute 3D centroids of both parts to identify the left from the
 292   // right bronchus
 293   std::vector<ImagePointType> c;
 294   clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
 295   ImagePointType C1 = c[1];
 296   ImagePointType C2 = c[2];
 297
 298   ImagePixelType rightLabel;
 299   ImagePixelType leftLabel;
 300   if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
 301   else { rightLabel = 2; leftLabel = 1; }
 302
 303   // Select LeftLabel (set one label to Backgroundvalue)
 304   RightBronchus =
 305     clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
 306                                    GetBackgroundValue(), GetForegroundValue());
 307   /*
 308     SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
 309     leftLabel, GetBackgroundValue(), false);
 310   */
 311   LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
 312                                                 GetBackgroundValue(), GetForegroundValue());
 313   /*
 314     SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
 315     rightLabel, GetBackgroundValue(), false);
 316   */
 317
 318   // Crop images
 319   RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
 320   LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
 321
 322   // Insert int AFDB if need after
 323   GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
 324                                                RightBronchus, true);
 325   GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
 326                                                LeftBronchus, true);
 327
 328   // Now crop below OriginOfRightMiddleLobeBronchusZ
 329   // It is not done before to keep entire bronchi.
 330
 331   MaskImagePointer OriginOfRightMiddleLobeBronchus =
 332     GetAFDB()->template GetImage<MaskImageType>("OriginOfRightMiddleLobeBronchus");
 333   std::vector<MaskImagePointType> centroids;
 334   clitk::ComputeCentroids<MaskImageType>(OriginOfRightMiddleLobeBronchus, GetBackgroundValue(), centroids);
 335   m_OriginOfRightMiddleLobeBronchusZ = centroids[1][2];
 336   // add one slice to include carina ?
 337   m_OriginOfRightMiddleLobeBronchusZ += RightBronchus->GetSpacing()[2];
 338   // We dont need Carina structure from now
 339   OriginOfRightMiddleLobeBronchus->Delete();
 340
 341   RightBronchus =
 342     clitk::CropImageBelow<MaskImageType>(RightBronchus, 2,
 343                                          m_OriginOfRightMiddleLobeBronchusZ,
 344                                          true, // AutoCrop
 345                                          GetBackgroundValue());
 346   LeftBronchus =
 347     clitk::CropImageBelow<MaskImageType>(LeftBronchus, 2,
 348                                          m_OriginOfRightMiddleLobeBronchusZ,
 349                                          true, // AutoCrop
 350                                          GetBackgroundValue());
 351
 352   // Search for points that are the most left/post/ant and most
 353   // right/post/ant of the left and right bronchus
 354   // 15  = not 15 mm more distance than the middle point.
 355   FindExtremaPointsInBronchus(RightBronchus, 0, 10, m_LeftMostInRightBronchus,
 356                               m_AntMostInRightBronchus, m_PostMostInRightBronchus);
 357
 358   FindExtremaPointsInBronchus(LeftBronchus, 1, 10, m_RightMostInLeftBronchus,
 359                               m_AntMostInLeftBronchus, m_PostMostInLeftBronchus);
 360
 361   // DEBUG : write the list of points
 362   ListOfPointsType v;
 363   v.reserve(m_LeftMostInRightBronchus.size()+m_AntMostInRightBronchus.size()+
 364             m_PostMostInRightBronchus.size());
 365   v.insert(v.end(), m_LeftMostInRightBronchus.begin(), m_LeftMostInRightBronchus.end() );
 366   v.insert(v.end(), m_AntMostInRightBronchus.begin(), m_AntMostInRightBronchus.end() );
 367   v.insert(v.end(), m_PostMostInRightBronchus.begin(), m_PostMostInRightBronchus.end() );
 368   clitk::WriteListOfLandmarks<MaskImageType>(v, "S8-RightBronchus-points.txt");
 369
 370   v.clear();
 371   v.reserve(m_RightMostInLeftBronchus.size()+m_AntMostInLeftBronchus.size()+
 372             m_PostMostInLeftBronchus.size());
 373   v.insert(v.end(), m_RightMostInLeftBronchus.begin(), m_RightMostInLeftBronchus.end() );
 374   v.insert(v.end(), m_AntMostInLeftBronchus.begin(), m_AntMostInLeftBronchus.end() );
 375   v.insert(v.end(), m_PostMostInLeftBronchus.begin(), m_PostMostInLeftBronchus.end() );
 376   clitk::WriteListOfLandmarks<MaskImageType>(v, "S8-LeftBronchus-points.txt");
 377
 378   v.clear();
 379   v.reserve(m_PostMostInLeftBronchus.size()+m_PostMostInRightBronchus.size());
 380   v.insert(v.end(), m_PostMostInLeftBronchus.begin(), m_PostMostInLeftBronchus.end() );
 381   v.insert(v.end(), m_PostMostInRightBronchus.begin(), m_PostMostInRightBronchus.end() );
 382   clitk::WriteListOfLandmarks<MaskImageType>(v, "S8-RightLeftBronchus-points.txt");
 383
 384
 385   // Now uses these points to limit, slice by slice
 386   // line is mainly horizontal, so mainDirection=1
 387   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 388                                                                     m_PostMostInRightBronchus,
 389                                                                     m_PostMostInLeftBronchus,
 390                                                                     GetBackgroundValue(), 1, 10);
 391
 392   // Keep main 3D CCL :
 393   m_Working_Support = Labelize<MaskImageType>(m_Working_Support, 0, false, 10);
 394   m_Working_Support = KeepLabels<MaskImageType>(m_Working_Support,
 395                                                 GetBackgroundValue(),
 396                                                 GetForegroundValue(), 1, 1, true);
 397
 398   // Autocrop
 399   m_Working_Support = clitk::AutoCrop<MaskImageType>(m_Working_Support, GetBackgroundValue());
 400
 401   // End of step
 402   StopCurrentStep<MaskImageType>(m_Working_Support);
 403   //  m_ListOfStations["8"] = m_Working_Support;
 404
 405 }
 406
 407 //--------------------------------------------------------------------
 408 template <class ImageType>
 409 void
 410 clitk::ExtractLymphStationsFilter<ImageType>::
 411 ExtractStation_8_Ant_Inf_Limits()
 412 {
 413
 414   //--------------------------------------------------------------------
 415   StartNewStep("[Station8] Ant part: not post to Esophagus");
 416   /*
 417     Consider Esophagus, dilate it and remove ant part. It remains part
 418     on L & R, than can be partly removed by cutting what remains at
 419     right of vertebral body.
 420   */
 421
 422   // Get Esophagus
 423   m_Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
 424   clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Esophagus");
 425
 426   // In images from the original article, Atlas – UM, the oesophagus
 427   //was included in nodal stations 3p and 8.  Having said that, in the
 428   //description for station 8, it indicates that “the delineation of
 429   //station 8 is limited to the soft tissues surrounding the
 430   //oesophagus”.  In the recent article, The IASLC Lung Cancer Staging
 431   //Project (J Thorac Oncol 4:5, 568-77), the images drawn by
 432   //Dr. Aletta Frasier exclude this structure.  From an oncological
 433   //prospective, the oesophagus should be excluded from these nodal
 434   //stations.
 435
 436   /* NOT YET !! DO IT LATER
 437
 438      typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
 439      typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 440      boolFilter->InPlaceOn();
 441      boolFilter->SetInput1(m_Working_Support);
 442      boolFilter->SetInput2(m_Esophagus);
 443      boolFilter->SetOperationType(BoolFilterType::AndNot);
 444      boolFilter->Update();
 445      m_Working_Support = boolFilter->GetOutput();
 446
 447   */
 448
 449   // Crop Esophagus : keep only below the OriginOfRightMiddleLobeBronchusZ
 450   m_Esophagus =
 451     clitk::CropImageAbove<MaskImageType>(m_Esophagus, 2,
 452                                          m_OriginOfRightMiddleLobeBronchusZ,
 453                                          true, // AutoCrop
 454                                          GetBackgroundValue());
 455
 456   // Dilate to keep only not Anterior positions
 457   MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt();
 458
 459   //  m_Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(m_Esophagus);
 460
 461   m_Esophagus = clitk::Dilate<MaskImageType>(m_Esophagus,
 462                                              radiusInMM,
 463                                              GetBackgroundValue(),
 464                                              GetForegroundValue(), true);
 465
 466   // Remove Anterior part according to this dilatated esophagus. Note:
 467   // because we crop Esophagus with ORML, the support will also be
 468   // croped in the same way. Here it is a desired feature. If we dont
 469   // want, use SetIgnoreEmptySliceObject(true)
 470
 471   // In the new IASCL definition, it is not clear if sup limits is
 472   //  around carina or On the right, it is “the lower border of the
 473   //  bronchus intermedius”, indicated on the image set as a point
 474   //  (“lower border of the bronchus intermedius”)
 475
 476   typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
 477   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 478   relPosFilter->VerboseStepFlagOff();
 479   relPosFilter->WriteStepFlagOff();
 480   relPosFilter->SetInput(m_Working_Support);
 481   relPosFilter->SetInputObject(m_Esophagus);
 482   relPosFilter->AddOrientationTypeString("PostTo");
 483   //  relPosFilter->InverseOrientationFlagOff();
 484   relPosFilter->SetDirection(2); // Z axis
 485   relPosFilter->UniqueConnectedComponentBySliceOff();
 486   relPosFilter->SetIntermediateSpacing(3);
 487   relPosFilter->IntermediateSpacingFlagOn();
 488   relPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS8());
 489   relPosFilter->RemoveObjectFlagOff(); // Do not exclude here because it is dilated
 490   relPosFilter->CombineWithOrFlagOff(); // NO !
 491   relPosFilter->IgnoreEmptySliceObjectFlagOn();
 492   relPosFilter->Update();
 493   m_Working_Support = relPosFilter->GetOutput();
 494
 495   // AutoCrop (OR SbS ?)
 496   m_Working_Support = clitk::AutoCrop<MaskImageType>(m_Working_Support, GetBackgroundValue());
 497
 498   StopCurrentStep<MaskImageType>(m_Working_Support);
 499   m_ListOfStations["8"] = m_Working_Support;
 500 }
 501 //--------------------------------------------------------------------
 502
 503
 504 //--------------------------------------------------------------------
 505 template <class ImageType>
 506 void
 507 clitk::ExtractLymphStationsFilter<ImageType>::
 508 ExtractStation_8_Ant_Injected_Limits()
 509 {
 510
 511   //--------------------------------------------------------------------
 512   StartNewStep("[Station8] Ant part (remove high density, injected part)");
 513
 514   // Consider initial image, crop to current support
 515   ImagePointer working_input =
 516     clitk::ResizeImageLike<ImageType>(m_Input,
 517                                       m_Working_Support,
 518                                       (short)GetBackgroundValue());
 519
 520   // Threshold
 521   typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> BinarizeFilterType;
 522   typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
 523   binarizeFilter->SetInput(working_input);
 524   binarizeFilter->SetLowerThreshold(GetInjectedThresholdForS8());
 525   binarizeFilter->SetInsideValue(GetForegroundValue());
 526   binarizeFilter->SetOutsideValue(GetBackgroundValue());
 527   binarizeFilter->Update();
 528   MaskImagePointer injected = binarizeFilter->GetOutput();
 529
 530   // Combine with current support
 531   clitk::AndNot<MaskImageType>(m_Working_Support, injected, GetBackgroundValue());
 532
 533   StopCurrentStep<MaskImageType>(m_Working_Support);
 534   m_ListOfStations["8"] = m_Working_Support;
 535 }
 536 //--------------------------------------------------------------------
 537
 538
 539
 540 //--------------------------------------------------------------------
 541 template <class ImageType>
 542 void
 543 clitk::ExtractLymphStationsFilter<ImageType>::
 544 ExtractStation_8_LR_1_Limits()
 545 {
 546   //--------------------------------------------------------------------
 547   StartNewStep("[Station8] Left and Right (from Carina to PulmonaryTrunk): Right to LeftPulmonaryArtery");
 548
 549   /*
 550     We remove LeftPulmonaryArtery structure and what is at Left to
 551     this structure.
 552   */
 553   MaskImagePointer LeftPulmonaryArtery = GetAFDB()->template GetImage<MaskImageType>("LeftPulmonaryArtery");
 554
 555   // Relative Position : not at Left
 556   typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
 557   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 558   relPosFilter->VerboseStepFlagOff();
 559   relPosFilter->WriteStepFlagOff();
 560   relPosFilter->SetBackgroundValue(GetBackgroundValue());
 561   relPosFilter->SetInput(m_Working_Support);
 562   relPosFilter->SetInputObject(LeftPulmonaryArtery);
 563   relPosFilter->RemoveObjectFlagOn(); // remove the object too
 564   relPosFilter->AddOrientationTypeString("L");
 565   relPosFilter->InverseOrientationFlagOn(); // Not at Left
 566   relPosFilter->SetIntermediateSpacing(3);
 567   relPosFilter->IntermediateSpacingFlagOn();
 568   relPosFilter->SetFuzzyThreshold(0.7);
 569   relPosFilter->AutoCropFlagOn();
 570   relPosFilter->Update();
 571   m_Working_Support = relPosFilter->GetOutput();
 572
 573   // Release LeftPulmonaryArtery
 574   GetAFDB()->template ReleaseImage<MaskImageType>("LeftPulmonaryArtery");
 575
 576   StopCurrentStep<MaskImageType>(m_Working_Support);
 577   m_ListOfStations["8"] = m_Working_Support;
 578 }
 579 //--------------------------------------------------------------------
 580
 581
 582 //--------------------------------------------------------------------
 583 template <class ImageType>
 584 void
 585 clitk::ExtractLymphStationsFilter<ImageType>::
 586 ExtractStation_8_LR_2_Limits()
 587 {
 588   //--------------------------------------------------------------------
 589   StartNewStep("[Station8] Left and Right (from PulmTrunk to OriginMiddleLobeBronchus) Right to line from Aorta to PulmonaryTrunk");
 590
 591   /*
 592     We consider a line from Left part of Aorta to left part of
 593     PulmonaryTrunk and remove what is at Left.
 594   */
 595
 596   // Load the structures
 597   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 598   MaskImagePointer PulmonaryTrunk = GetAFDB()->template GetImage<MaskImageType>("PulmonaryTrunk");
 599
 600   // Resize like the PT and define the slices
 601   MaskImagePointType min, max;
 602   clitk::GetMinMaxPointPosition<MaskImageType>(PulmonaryTrunk, min, max);
 603   Aorta = clitk::CropImageAlongOneAxis<MaskImageType>(Aorta, 2, min[2], max[2], false, GetBackgroundValue());
 604   std::vector<MaskSlicePointer> slices_aorta;
 605   clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_aorta);
 606   std::vector<MaskSlicePointer> slices_PT;
 607   clitk::ExtractSlices<MaskImageType>(PulmonaryTrunk, 2, slices_PT);
 608
 609   // Find the points at left
 610   std::vector<MaskImagePointType> p_A;
 611   std::vector<MaskImagePointType> p_B;
 612   MaskImagePointType pA;
 613   MaskImagePointType pB;
 614   for(uint i=0; i<slices_PT.size(); i++) {
 615     typename MaskSliceType::PointType ps;
 616     // In Aorta (assume a single CCL)
 617     bool found =
 618       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices_aorta[i], GetBackgroundValue(), 0, false, ps);
 619     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(ps, Aorta, i, pA);
 620
 621     if (found) {
 622       // In PT : generally 2 CCL, we keep the most at left
 623       found =
 624         clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices_PT[i], GetBackgroundValue(), 0, false, ps);
 625       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(ps, PulmonaryTrunk, i, pB);
 626     }
 627
 628     if (found) {
 629       p_A.push_back(pA);
 630       p_B.push_back(pB);
 631     }
 632   }
 633   clitk::WriteListOfLandmarks<MaskImageType>(p_A, "S8-Aorta-Left-points.txt");
 634   clitk::WriteListOfLandmarks<MaskImageType>(p_B, "S8-PT-Left-points.txt");
 635
 636   // Remove part at Left
 637   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 638                                                                     p_A, p_B,
 639                                                                     GetBackgroundValue(),
 640                                                                     0, -10);
 641
 642   StopCurrentStep<MaskImageType>(m_Working_Support);
 643   m_ListOfStations["8"] = m_Working_Support;
 644 }
 645 //--------------------------------------------------------------------
 646
 647
 648 //--------------------------------------------------------------------
 649 template <class ImageType>
 650 void
 651 clitk::ExtractLymphStationsFilter<ImageType>::
 652 ExtractStation_8_Single_CCL_Limits()
 653 {
 654   //--------------------------------------------------------------------
 655   StartNewStep("[Station8 or 3P] Slice by slice, keep a single CCL (the closest to VertebralBody)");
 656
 657   // Consider slices
 658   std::vector<typename MaskSliceType::Pointer> slices;
 659   std::vector<typename MaskSliceType::Pointer> slices_vb;
 660   clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices);
 661   MaskImagePointer VertebralBody =
 662     GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
 663   VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
 664   clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, slices_vb);
 665
 666   for(uint i=0; i<slices.size(); i++) {
 667     // Decompose in labels
 668     slices[i] = Labelize<MaskSliceType>(slices[i], 0, true, 100);
 669     // Compute centroids coordinate
 670     std::vector<typename MaskSliceType::PointType> centroids;
 671     std::vector<typename MaskSliceType::PointType> c;
 672     clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), centroids);
 673     clitk::ComputeCentroids<MaskSliceType>(slices_vb[i], GetBackgroundValue(), c);
 674     if ((c.size() > 1) && (centroids.size() > 1)) {
 675       // keep the one which is the closer to vertebralbody centroid
 676       double distance=1000000;
 677       int index=0;
 678       for(uint j=1; j<centroids.size(); j++) {
 679         double d = centroids[j].EuclideanDistanceTo(c[1]);
 680         if (d<distance) {
 681           distance = d;
 682           index = j;
 683         }
 684       }
 685       // remove all others label
 686       slices[i] = KeepLabels<MaskSliceType>(slices[i], GetBackgroundValue(),
 687                                             GetForegroundValue(), index, index, true);
 688     }
 689   }
 690   m_Working_Support = clitk::JoinSlices<MaskImageType>(slices, m_Working_Support, 2);
 691
 692   StopCurrentStep<MaskImageType>(m_Working_Support);
 693   m_ListOfStations["8"] = m_Working_Support;
 694 }
 695 //--------------------------------------------------------------------
 696
 697
 698 //--------------------------------------------------------------------
 699 template <class ImageType>
 700 void
 701 clitk::ExtractLymphStationsFilter<ImageType>::
 702 ExtractStation_8_LR_Limits_old2()
 703 {
 704
 705   //--------------------------------------------------------------------
 706   StartNewStep("[Station8] Left and Right limits arround esophagus (below Carina)");
 707
 708   // Estract slices for current support for slice by slice processing
 709   std::vector<typename MaskSliceType::Pointer> slices;
 710   clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices);
 711
 712   // Dilate the Esophagus to consider a margins around
 713   MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt();
 714   m_Esophagus = clitk::Dilate<MaskImageType>(m_Esophagus,
 715                                              radiusInMM,
 716                                              GetBackgroundValue(),
 717                                              GetForegroundValue(), true);
 718
 719   // Remove what is outside the mediastinum in this enlarged Esophagus -> it allows to select
 720   // 'better' extrema points (not too post).
 721   MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
 722   clitk::AndNot<MaskImageType>(m_Esophagus, Lungs, GetBackgroundValue());
 723   GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
 724
 725   // Estract slices of Esophagus (resize like support before to have the same set of slices)
 726   MaskImagePointer EsophagusForSlice = clitk::ResizeImageLike<MaskImageType>(m_Esophagus, m_Working_Support, GetBackgroundValue());
 727
 728   std::vector<typename MaskSliceType::Pointer> eso_slices;
 729   clitk::ExtractSlices<MaskImageType>(EsophagusForSlice, 2, eso_slices);
 730
 731   // Estract slices of Vertebral (resize like support before to have the same set of slices)
 732   MaskImagePointer VertebralBody = GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
 733   VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
 734   std::vector<typename MaskSliceType::Pointer> vert_slices;
 735   clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vert_slices);
 736
 737   // Estract slices of Aorta (resize like support before to have the same set of slices)
 738   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 739   Aorta = clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
 740   std::vector<typename MaskSliceType::Pointer> aorta_slices;
 741   clitk::ExtractSlices<MaskImageType>(Aorta, 2, aorta_slices);
 742
 743   // Extract slices of Mediastinum (resize like support before to have the same set of slices)
 744   m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
 745   m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Working_Support, GetBackgroundValue());
 746   std::vector<typename MaskSliceType::Pointer> mediast_slices;
 747   clitk::ExtractSlices<MaskImageType>(m_Mediastinum, 2, mediast_slices);
 748
 749   // List of points
 750   std::vector<MaskImagePointType> p_RightMostAnt;
 751   std::vector<MaskImagePointType> p_RightMostPost;
 752   std::vector<MaskImagePointType> p_LeftMostAnt;
 753   std::vector<MaskImagePointType> p_LeftMostPost;
 754   std::vector<MaskImagePointType> p_AllPoints;
 755   std::vector<MaskImagePointType> p_LeftAorta;
 756   std::vector<MaskImagePointType> p_LeftEso;
 757
 758   /*
 759     In the following, we search for the LeftRight limits.  We search
 760     for the most Right points in Esophagus and in VertebralBody and
 761     consider a line between those to most right points. All points in
 762     the support which are most right to this line are discarded. Same
 763     for the left part. The underlying assumption is that the support
 764     is concave between Eso/VertebralBody. Esophagus is a bit
 765     dilatated. On VertebralBody we go right (or left) until we reach
 766     the lung (but no more 20 mm).
 767   */
 768
 769   // Loop slices
 770   MaskImagePointType p;
 771   MaskImagePointType pp;
 772   for(uint i=0; i<slices.size() ; i++) {
 773     // Declare all needed points (sp = slice point)
 774     typename MaskSliceType::PointType sp_maxRight_Eso;
 775     typename MaskSliceType::PointType sp_maxRight_Aorta;
 776     typename MaskSliceType::PointType sp_maxRight_Vertebra;
 777     typename MaskSliceType::PointType sp_maxLeft_Eso;
 778     typename MaskSliceType::PointType sp_maxLeft_Aorta;
 779     typename MaskSliceType::PointType sp_maxLeft_Vertebra;
 780
 781     // Right is at left on screen, coordinate decrease
 782     // Left is at right on screen, coordinate increase
 783
 784     // Find limit of Vertebral -> only at most Post part of current
 785     // slice support.  First found most ant point in VertebralBody
 786     typedef MaskSliceType SliceType;
 787     typename SliceType::PointType p_slice_ant;
 788     bool found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vert_slices[i], GetBackgroundValue(), 1, true, p_slice_ant);
 789     if (!found) {
 790       // It should not happen ! But sometimes, a contour is missing or
 791       // the VertebralBody is not delineated enough inferiorly ... in
 792       // those cases, we consider the first found slice.
 793       std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
 794       int j=i++;
 795       bool found = false;
 796       while (!found) {
 797         found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vert_slices[j], GetBackgroundValue(), 1, true, p_slice_ant);
 798         //clitkExceptionMacro("No foreground pixels in this VertebralBody slices ??");
 799         j++;
 800       }
 801     }
 802     p_slice_ant[1] += GetDistanceMaxToAnteriorPartOfTheSpine(); // Consider offset
 803
 804     // The, find most Right and Left points on that AP position
 805     typename SliceType::IndexType indexR;
 806     typename SliceType::IndexType indexL;
 807     vert_slices[i]->TransformPhysicalPointToIndex(p_slice_ant, indexR);
 808     indexL = indexR;
 809     // Check that is inside the mask
 810     indexR[1] = std::min(indexR[1], (long)vert_slices[i]->GetLargestPossibleRegion().GetSize()[1]-1);
 811     indexL[1] = indexR[1];
 812     while (vert_slices[i]->GetPixel(indexR) != GetBackgroundValue()) {
 813       indexR[0] --; // Go to the right
 814     }
 815     while (vert_slices[i]->GetPixel(indexL) != GetBackgroundValue()) {
 816       indexL[0] ++; // Go to the left
 817     }
 818     vert_slices[i]->TransformIndexToPhysicalPoint(indexR, sp_maxRight_Vertebra);
 819     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, VertebralBody, i, p);
 820     p_AllPoints.push_back(p);
 821     vert_slices[i]->TransformIndexToPhysicalPoint(indexL, sp_maxLeft_Vertebra);
 822     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, VertebralBody, i, p);
 823     p_AllPoints.push_back(p);
 824
 825     // Find last point out of the mediastinum on this line, Right :
 826     mediast_slices[i]->TransformPhysicalPointToIndex(sp_maxRight_Vertebra, indexR);
 827     double distance = 0.0;
 828     while (mediast_slices[i]->GetPixel(indexR) != GetBackgroundValue()) {
 829       indexR[0] --;
 830       distance += mediast_slices[i]->GetSpacing()[0];
 831     }
 832     if (distance < 30) { // Ok in this case, we found limit with lung
 833       mediast_slices[i]->TransformIndexToPhysicalPoint(indexR, sp_maxRight_Vertebra);
 834       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, m_Mediastinum, i, p);
 835     }
 836     else { // in that case, we are probably below the diaphragm, so we
 837            // add aribtrarly few mm
 838       sp_maxRight_Vertebra[0] -= 2; // Leave 2 mm around the VertebralBody
 839       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, m_Mediastinum, i, p);
 840     }
 841     p_RightMostPost.push_back(p);
 842     p_AllPoints.push_back(p);
 843
 844     // Find last point out of the mediastinum on this line, Left :
 845     mediast_slices[i]->TransformPhysicalPointToIndex(sp_maxLeft_Vertebra, indexL);
 846     distance = 0.0;
 847     while (mediast_slices[i]->GetPixel(indexL) != GetBackgroundValue()) {
 848       indexL[0] ++;
 849       distance += mediast_slices[i]->GetSpacing()[0];
 850     }
 851     if (distance < 30) { // Ok in this case, we found limit with lung
 852       mediast_slices[i]->TransformIndexToPhysicalPoint(indexL, sp_maxLeft_Vertebra);
 853       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, m_Mediastinum, i, p);
 854     }
 855     else { // in that case, we are probably below the diaphragm, so we
 856            // add aribtrarly few mm
 857       sp_maxLeft_Vertebra[0] += 2; // Leave 2 mm around the VertebralBody
 858       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, m_Mediastinum, i, p);
 859     }
 860     p_LeftMostPost.push_back(p);
 861     p_AllPoints.push_back(p);
 862
 863     // Find Eso slice centroid and do not consider what is post to
 864     // this centroid.
 865     std::vector<typename MaskSliceType::PointType> c;
 866     clitk::ComputeCentroids<MaskSliceType>(eso_slices[i], GetBackgroundValue(), c);
 867     if (c.size() >1) {
 868       eso_slices[i] =
 869         clitk::CropImageAbove<MaskSliceType>(eso_slices[i], 1, c[1][1], false, GetBackgroundValue());
 870       eso_slices[i] =
 871         clitk::ResizeImageLike<MaskSliceType>(eso_slices[i], aorta_slices[i], GetBackgroundValue());
 872       // writeImage<MaskSliceType>(eso_slices[i], "eso-slice-"+toString(i)+".mhd");
 873     }
 874
 875     // Find right limit of Esophagus and Aorta
 876     bool f =
 877       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(),
 878                                                               0, true, sp_maxRight_Eso);
 879     f = f &&
 880       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(),
 881                                                               0, true, sp_maxRight_Aorta);
 882     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Eso, EsophagusForSlice, i, p);
 883     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Aorta, Aorta, i, pp);
 884     pp[0] -= 2; // Add a margin of 2 mm to include the Aorta 'wall'
 885     p_AllPoints.push_back(p);
 886     if (f) {
 887       p_AllPoints.push_back(pp);
 888       MaskImagePointType A = p_RightMostPost.back();
 889       MaskImagePointType B = p;
 890       MaskImagePointType C = pp;
 891       double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 892       if (s>0) p_RightMostAnt.push_back(p);
 893       else p_RightMostAnt.push_back(pp);
 894     }
 895     else { // No more Esophagus in this slice : do nothing
 896       //      p_RightMostAnt.push_back(p);
 897       p_RightMostPost.pop_back();
 898     }
 899
 900     // --------------------------------------------------------------------------
 901     // Find the limit on the Left: most left point between Eso and
 902     // Eso. (Left is left on screen, coordinate increase)
 903
 904     // Find left limit of Esophagus
 905     clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 0, false, sp_maxLeft_Eso);
 906     f = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(), 0, false, sp_maxLeft_Aorta);
 907     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Eso, EsophagusForSlice, i, p);
 908     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Aorta, Aorta, i, pp);
 909     p_AllPoints.push_back(p);
 910     pp[0] += 2; // Add a margin of 2 mm to include the 'wall'
 911     if (f) { // not below Aorta
 912       p_AllPoints.push_back(pp);
 913       MaskImagePointType A = p_LeftMostPost.back();
 914       MaskImagePointType B = p;
 915       MaskImagePointType C = pp;
 916       double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 917       if (s<0) {
 918         p_LeftMostAnt.push_back(p); // Insert point most at Left, Eso
 919       }
 920       else {
 921         // in this case -> two lines !
 922         p_LeftMostAnt.push_back(pp);  // Insert point most at Left, Aorta (Vert to Aorta)
 923         // but also consider Aorta to Eso
 924         p_LeftAorta.push_back(pp);
 925         p_LeftEso.push_back(p);
 926       }
 927     }
 928     else { // No more Esophagus in this slice : do nothing
 929       p_LeftMostPost.pop_back();
 930       //p_LeftMostAnt.push_back(p);
 931     }
 932   } // End of slice loop
 933
 934   clitk::WriteListOfLandmarks<MaskImageType>(p_AllPoints, "S8-LR-Eso-Vert.txt");
 935   clitk::WriteListOfLandmarks<MaskImageType>(p_LeftEso, "S8-Left-Eso.txt");
 936   clitk::WriteListOfLandmarks<MaskImageType>(p_LeftAorta, "S8-Left-Aorta.txt");
 937
 938   // Now uses these points to limit, slice by slice
 939   // Line is mainly vertical, so mainDirection=0
 940   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 941                                                                     p_RightMostAnt, p_RightMostPost,
 942                                                                     GetBackgroundValue(), 0, 10);
 943   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 944                                                                     p_LeftMostAnt, p_LeftMostPost,
 945                                                                     GetBackgroundValue(), 0, -10);
 946   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 947                                                                     p_LeftEso,p_LeftAorta,
 948                                                                     GetBackgroundValue(), 0, -10);
 949   // END
 950   StopCurrentStep<MaskImageType>(m_Working_Support);
 951   m_ListOfStations["8"] = m_Working_Support;
 952   return;
 953 }
 954 //--------------------------------------------------------------------
 955
 956
 957 //--------------------------------------------------------------------
 958 template <class ImageType>
 959 void
 960 clitk::ExtractLymphStationsFilter<ImageType>::
 961 ExtractStation_8_LR_Limits()
 962 {
 963
 964   //--------------------------------------------------------------------
 965   StartNewStep("[Station8] Left and Right limits arround esophagus with lines from VertebralBody-Aorta-Esophagus");
 966
 967   // Estract slices for current support for slice by slice processing
 968   std::vector<typename MaskSliceType::Pointer> slices;
 969   clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices);
 970
 971   // Dilate the Esophagus to consider a margins around
 972   MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt();
 973   m_Esophagus = clitk::Dilate<MaskImageType>(m_Esophagus, radiusInMM,
 974                                              GetBackgroundValue(), GetForegroundValue(), true);
 975   //  m_Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(m_Esophagus);
 976
 977   // Remove what is outside the mediastinum in this enlarged Esophagus -> it allows to select
 978   // 'better' extrema points (not too post).
 979   MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
 980   clitk::AndNot<MaskImageType>(m_Esophagus, Lungs, GetBackgroundValue());
 981   GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
 982
 983   // Estract slices of Esophagus (resize like support before to have the same set of slices)
 984   MaskImagePointer EsophagusForSlice = clitk::ResizeImageLike<MaskImageType>(m_Esophagus, m_Working_Support, GetBackgroundValue());
 985   std::vector<typename MaskSliceType::Pointer> eso_slices;
 986   clitk::ExtractSlices<MaskImageType>(EsophagusForSlice, 2, eso_slices);
 987
 988   // Estract slices of Vertebral (resize like support before to have the same set of slices)
 989   MaskImagePointer VertebralBody = GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
 990   VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
 991   std::vector<typename MaskSliceType::Pointer> vert_slices;
 992   clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vert_slices);
 993
 994   // Estract slices of Aorta (resize like support before to have the same set of slices)
 995   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 996   Aorta = clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
 997   std::vector<typename MaskSliceType::Pointer> aorta_slices;
 998   clitk::ExtractSlices<MaskImageType>(Aorta, 2, aorta_slices);
 999
1000   // Extract slices of Mediastinum (resize like support before to have the same set of slices)
1001   m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
1002   m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Working_Support, GetBackgroundValue());
1003   std::vector<typename MaskSliceType::Pointer> mediast_slices;
1004   clitk::ExtractSlices<MaskImageType>(m_Mediastinum, 2, mediast_slices);
1005
1006   // List of points
1007   std::vector<MaskImagePointType> p_MostLeftVertebralBody;
1008   std::vector<MaskImagePointType> p_MostRightVertebralBody;
1009   std::vector<MaskImagePointType> p_MostLeftAorta;
1010   std::vector<MaskImagePointType> p_MostLeftEsophagus;
1011
1012   /*
1013     In the following, we search for the LeftRight limits.  We search
1014     for the most Right points in Esophagus and in VertebralBody and
1015     consider a line between those to most right points. All points in
1016     the support which are most right to this line are discarded. Same
1017     for the left part. The underlying assumption is that the support
1018     is concave between Eso/VertebralBody. Esophagus is a bit
1019     dilatated. On VertebralBody we go right (or left) until we reach
1020     the lung (but no more 20 mm).
1021   */
1022
1023   // Temporary 3D point
1024   MaskImagePointType p;
1025
1026   // Loop slices
1027   for(uint i=0; i<slices.size() ; i++) {
1028     // Declare all needed 2D points (sp = slice point)
1029     typename MaskSliceType::PointType sp_MostRightVertebralBody;
1030     typename MaskSliceType::PointType sp_MostLeftVertebralBody;
1031     typename MaskSliceType::PointType sp_MostLeftAorta;
1032     typename MaskSliceType::PointType sp_temp;
1033     typename MaskSliceType::PointType sp_MostLeftEsophagus;
1034
1035     // Right is at left on screen, coordinate decrease
1036     // Left is at right on screen, coordinate increase
1037
1038     /* -------------------------------------------------------------------------------------
1039        Find most Left point in VertebralBody. Consider only the
1040        horizontal line which is 'DistanceMaxToAnteriorPartOfTheSpine'
1041        away from the most ant point.
1042      */
1043     typename MaskSliceType::PointType sp_MostAntVertebralBody;
1044     bool found = false;
1045     int j=i;
1046     while (!found) {
1047       found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vert_slices[j], GetBackgroundValue(), 1, true, sp_MostAntVertebralBody);
1048       if (!found) {
1049         // It should not happen ! But sometimes, a contour is missing or
1050         // the VertebralBody is not delineated enough inferiorly ... in
1051         // those cases, we consider the first found slice.
1052         std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
1053         j++;
1054       }
1055     }
1056     sp_MostAntVertebralBody[1] += GetDistanceMaxToAnteriorPartOfTheSpine(); // Consider offset
1057
1058     // Crop the vertebralbody below this most post line
1059     vert_slices[j] =
1060       clitk::CropImageAbove<MaskSliceType>(vert_slices[j], 1, sp_MostAntVertebralBody[1], false, GetBackgroundValue());
1061     vert_slices[j] =
1062       clitk::ResizeImageLike<MaskSliceType>(vert_slices[j], aorta_slices[i], GetBackgroundValue());
1063     //    writeImage<MaskSliceType>(vert_slices[i], "vert-slice-"+toString(i)+".mhd");
1064
1065     // Find first point not in mediastinum
1066     clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vert_slices[j], GetBackgroundValue(), 0, false, sp_MostLeftVertebralBody);
1067     sp_MostLeftVertebralBody = clitk::FindExtremaPointInAGivenLine<MaskSliceType>(mediast_slices[i], 0, false, sp_MostLeftVertebralBody, GetBackgroundValue(), 30);
1068     sp_MostLeftVertebralBody[0] += 2; // 2mm margin
1069     clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vert_slices[j], GetBackgroundValue(), 0, true, sp_MostRightVertebralBody);
1070     sp_MostRightVertebralBody = clitk::FindExtremaPointInAGivenLine<MaskSliceType>(mediast_slices[i], 0, true, sp_MostRightVertebralBody, GetBackgroundValue(),30);
1071     sp_MostRightVertebralBody[0] -= 2; // 2 mm margin
1072
1073     // Convert 2D points into 3D
1074     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostRightVertebralBody, VertebralBody, i, p);
1075     p_MostRightVertebralBody.push_back(p);
1076     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftVertebralBody, VertebralBody, i, p);
1077     p_MostLeftVertebralBody.push_back(p);
1078
1079     /* -------------------------------------------------------------------------------------
1080        Find most Left point in Esophagus
1081      */
1082     found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 0, false, sp_MostLeftEsophagus);
1083     if (!found) { // No more Esophagus, I remove the previous point
1084
1085       // if (p_MostLeftEsophagus.size() < 1)  {
1086       p_MostLeftVertebralBody.pop_back();
1087       // }
1088       // else {
1089       //   // Consider the previous point
1090       //   p = p_MostLeftEsophagus.back();
1091       //   p_MostLeftEsophagus.push_back(p);
1092       //   sp_MostLeftEsophagus = sp_temp; // Retrieve previous 2D position
1093       //   found = true;
1094       // }
1095     }
1096     else {
1097       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftEsophagus, EsophagusForSlice, i, p);
1098       p_MostLeftEsophagus.push_back(p);
1099       // sp_temp = sp_MostLeftEsophagus; // Store previous 2D position
1100     }
1101
1102     /* -------------------------------------------------------------------------------------
1103        Find most Left point in Aorta
1104     */
1105     if (found) {
1106       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(), 0, false, sp_MostLeftAorta);
1107       sp_MostLeftAorta = clitk::FindExtremaPointInAGivenLine<MaskSliceType>(mediast_slices[i], 0, false, sp_MostLeftAorta, GetBackgroundValue(), 10);
1108       typename MaskSliceType::PointType temp=sp_MostLeftEsophagus;
1109       temp[0] -= 10;
1110       if (clitk::IsOnTheSameLineSide(sp_MostLeftAorta, sp_MostLeftVertebralBody, sp_MostLeftEsophagus, temp)) {
1111         // sp_MostLeftAorta is on the same side than temp (at Right) -> ignore Aorta
1112         sp_MostLeftAorta = sp_MostLeftEsophagus;
1113       }
1114       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftAorta, Aorta, i, p);
1115       p_MostLeftAorta.push_back(p);
1116     }
1117
1118   } // End of slice loop
1119
1120   clitk::WriteListOfLandmarks<MaskImageType>(p_MostLeftVertebralBody, "S8-MostLeft-VB-points.txt");
1121   clitk::WriteListOfLandmarks<MaskImageType>(p_MostRightVertebralBody, "S8-MostRight-VB-points.txt");
1122   clitk::WriteListOfLandmarks<MaskImageType>(p_MostLeftAorta, "S8-MostLeft-Aorta-points.txt");
1123   clitk::WriteListOfLandmarks<MaskImageType>(p_MostLeftEsophagus, "S8-MostLeft-eso-points.txt");
1124
1125   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
1126                                                                     p_MostLeftVertebralBody, p_MostLeftAorta,
1127                                                                     GetBackgroundValue(), 0, -10);
1128
1129   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
1130                                                                     p_MostLeftAorta, p_MostLeftEsophagus,
1131                                                                     GetBackgroundValue(), 0, -10);
1132
1133   // END
1134   StopCurrentStep<MaskImageType>(m_Working_Support);
1135   m_ListOfStations["8"] = m_Working_Support;
1136   return;
1137 }
1138 //--------------------------------------------------------------------
1139
1140
1141 //--------------------------------------------------------------------
1142 template <class ImageType>
1143 void
1144 clitk::ExtractLymphStationsFilter<ImageType>::
1145 ExtractStation_8_Remove_Structures()
1146 {
1147
1148   //--------------------------------------------------------------------
1149   StartNewStep("[Station8] remove some structures");
1150
1151   Remove_Structures("Aorta");
1152   Remove_Structures("Esophagus");
1153
1154   // END
1155   StopCurrentStep<MaskImageType>(m_Working_Support);
1156   m_ListOfStations["8"] = m_Working_Support;
1157   return;
1158 }
1159 //--------------------------------------------------------------------
1160
1161
1162 //--------------------------------------------------------------------
1163 template <class ImageType>
1164 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1165 clitk::ExtractLymphStationsFilter<ImageType>::
1166 EnlargeEsophagusDilatationRadiusInferiorly(MaskImagePointer & Esophagus)
1167 {
1168   // Check if Esophagus is delineated at least until Diaphragm. Now,
1169   // because we use AutoCrop, Origin[2] gives this max inferior
1170   // position.
1171
1172   DD("BUGGY DONT USE");
1173   exit(0);
1174
1175   if (Esophagus->GetOrigin()[2] > m_DiaphragmInferiorLimit) {
1176     // crop first slice without mask
1177     MaskImagePointType pt;
1178     clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Esophagus, GetBackgroundValue(), 2, true, pt);
1179     DD(pt);
1180     Esophagus =
1181       clitk::CropImageBelow<MaskImageType>(Esophagus, 2,
1182                                            pt[2],
1183                                            false, // AutoCrop
1184                                            GetBackgroundValue());
1185     writeImage<MaskImageType>(Esophagus, "crop-eso.mhd");
1186
1187     std::cout << "Warning Esophagus is not delineated until Diaphragm. I mirror-pad it."
1188               << std::endl;
1189     double extraSize = Esophagus->GetOrigin()[2]-m_DiaphragmInferiorLimit;
1190
1191     // Pad with few more slices
1192     typedef itk::MirrorPadImageFilter<MaskImageType, MaskImageType> PadFilterType;
1193     typename PadFilterType::Pointer padFilter = PadFilterType::New();
1194     padFilter->SetInput(Esophagus);
1195     MaskImageSizeType b;
1196     b[0] = 0; b[1] = 0; b[2] = (uint)ceil(extraSize/Esophagus->GetSpacing()[2])+1;
1197     padFilter->SetPadLowerBound(b);
1198     padFilter->Update();
1199     Esophagus = padFilter->GetOutput();
1200   }
1201   return Esophagus;
1202 }
1203 //--------------------------------------------------------------------
1204
1205
1206 //--------------------------------------------------------------------
1207 template <class ImageType>
1208 void
1209 clitk::ExtractLymphStationsFilter<ImageType>::
1210 ExtractStation_8_LR_Limits_old()
1211 {
1212   /*
1213     Station 8: paraeosphageal nodes
1214
1215     Laterally, it is within the pleural envelope and again abuts the
1216     descending aorta on the left. Reasonably, the delineation of
1217     Station 8 is limited to the soft tissue surrounding the esophagus
1218     (Fig.  3C–H).
1219   */
1220
1221   StartNewStep("[Station8] Right limits (around esophagus)");
1222   // Get Esophagus
1223   MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
1224
1225   // Autocrop to get first slice with starting Esophagus
1226   Esophagus = clitk::AutoCrop<MaskImageType>(Esophagus, GetBackgroundValue());
1227
1228   // Dilate
1229   // LR dilatation -> large to keep point inside
1230   // AP dilatation -> few mm
1231   // SI dilatation -> enough to cover Diaphragm (old=GOjunctionZ)
1232   MaskImagePointType radiusInMM = GetEsophagusDiltationForRight();
1233   Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
1234   Esophagus = clitk::Dilate<MaskImageType>(Esophagus,
1235                                            radiusInMM,
1236                                            GetBackgroundValue(),
1237                                            GetForegroundValue(), true);
1238   writeImage<MaskImageType>(Esophagus, "dilateEso2.mhd");
1239
1240   writeImage<MaskImageType>(m_Working_Support, "before-relpos2.mhd");
1241
1242   // Remove Right (Left on the screen) part according to this
1243   // dilatated esophagus
1244   typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
1245   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
1246   relPosFilter->VerboseStepFlagOff();
1247   relPosFilter->WriteStepFlagOff();
1248   relPosFilter->SetInput(m_Working_Support);
1249   relPosFilter->SetInputObject(Esophagus);
1250   relPosFilter->AddOrientationTypeString("NotLeftTo");
1251   //  relPosFilter->InverseOrientationFlagOn(); // Not Left to
1252   relPosFilter->SetDirection(2); // Z axis
1253   relPosFilter->UniqueConnectedComponentBySliceOff(); // important
1254   relPosFilter->SetIntermediateSpacing(4);
1255   relPosFilter->IntermediateSpacingFlagOn();
1256   relPosFilter->SetFuzzyThreshold(0.9); // remove few part only
1257   relPosFilter->RemoveObjectFlagOff();
1258   relPosFilter->Update();
1259   m_Working_Support = relPosFilter->GetOutput();
1260
1261   // Get a single 3D CCL
1262   m_Working_Support = Labelize<MaskImageType>(m_Working_Support, 0, false, 10);
1263   m_Working_Support = KeepLabels<MaskImageType>(m_Working_Support,
1264                                                 GetBackgroundValue(),
1265                                                 GetForegroundValue(), 1, 1, true);
1266
1267
1268   /*
1269   // Re-Add post to Esophagus -> sometimes previous relpos remove some
1270   // correct part below esophagus.
1271   MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
1272   EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
1273   writeImage<MaskImageType>(Esophagus, "e-again.mhd");
1274   typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
1275   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
1276   relPosFilter->VerboseStepOff();
1277   relPosFilter->WriteStepOff();
1278   relPosFilter->SetInput(m_Working_Support);
1279   relPosFilter->SetInputObject(Esophagus);
1280   relPosFilter->SetOrientationTypeString("P");
1281   relPosFilter->InverseOrientationFlagOff();
1282   relPosFilter->SetDirection(2); // Z axis
1283   relPosFilter->UniqueConnectedComponentBySliceOff(); // important
1284   relPosFilter->SetIntermediateSpacing(4);
1285   relPosFilter->IntermediateSpacingFlagOn();
1286   relPosFilter->CombineWithOrFlagOn();
1287   relPosFilter->SetFuzzyThreshold(0.9); // remove few part only
1288   relPosFilter->RemoveObjectFlagOff();
1289   relPosFilter->Update();
1290   m_Working_Support = relPosFilter->GetOutput();
1291   */
1292
1293   StopCurrentStep<MaskImageType>(m_Working_Support);
1294   m_ListOfStations["8"] = m_Working_Support;
1295 }
1296 //--------------------------------------------------------------------
1297
1298
1299 //--------------------------------------------------------------------
1300 template <class TImageType>
1301 void
1302 clitk::ExtractLymphStationsFilter<TImageType>::
1303 FindExtremaPointsInBronchus(MaskImagePointer input,
1304                             int direction,
1305                             double distance_max_from_center_point,
1306                             ListOfPointsType & LR,
1307                             ListOfPointsType & Ant,
1308                             ListOfPointsType & Post)
1309 {
1310
1311   // Other solution ==> with auto bounding box ! (but pb to prevent to
1312   // be too distant from the center point
1313
1314   // Extract slices
1315   std::vector<typename MaskSliceType::Pointer> slices;
1316   clitk::ExtractSlices<MaskImageType>(input, 2, slices);
1317
1318   // Loop on slices
1319   bool found;
1320   for(uint i=0; i<slices.size(); i++) {
1321     /*
1322     // Keep main CCL
1323     slices[i] = Labelize<MaskSliceType>(slices[i], 0, true, 10);
1324     slices[i] = KeepLabels<MaskSliceType>(slices[i],
1325     GetBackgroundValue(),
1326     GetForegroundValue(), 1, 1, true);
1327     */
1328
1329     // ------- Find rightmost or leftmost point  -------
1330     MaskSliceType::PointType LRMost;
1331     found =
1332       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i],
1333                                                               GetBackgroundValue(),
1334                                                               0, // axis XY
1335                                                               (direction==0?false:true),  // right or left according to direction
1336                                                               LRMost);
1337     // ------- Find postmost point  -------
1338     MaskSliceType::PointType postMost;
1339     found =
1340       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i],
1341                                                               GetBackgroundValue(),
1342                                                               1, false, LRMost,
1343                                                               distance_max_from_center_point,
1344                                                               postMost);
1345     // ------- Find antmost point  -------
1346     MaskSliceType::PointType antMost;
1347     found =
1348       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i],
1349                                                               GetBackgroundValue(),
1350                                                               1, true, LRMost,
1351                                                               distance_max_from_center_point,
1352                                                               antMost);
1353     // Only add point if found
1354     if (found)  {
1355       // ------- Convert 2D to 3D points --------
1356       MaskImageType::PointType p;
1357       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(LRMost, input, i, p);
1358       LR.push_back(p);
1359       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(antMost, input, i, p);
1360       Ant.push_back(p);
1361       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(postMost, input, i, p);
1362       Post.push_back(p);
1363     }
1364   }
1365 }
1366 //--------------------------------------------------------------------