segmentation/clitkExtractLymphStation_8.txx

   1
   2 #include <itkBinaryDilateImageFilter.h>
   3 #include <itkMirrorPadImageFilter.h>
   4
   5 //--------------------------------------------------------------------
   6 template <class ImageType>
   7 void
   8 clitk::ExtractLymphStationsFilter<ImageType>::
   9 ExtractStation_8_SetDefaultValues()
  10 {
  11   SetDistanceMaxToAnteriorPartOfTheSpine(10);
  12   MaskImagePointType p;
  13   p[0] = 15; p[1] = 2; p[2] = 1;
  14   SetEsophagusDiltationForAnt(p);
  15   p[0] = 5; p[1] = 10; p[2] = 1;
  16   SetEsophagusDiltationForRight(p);
  17   SetFuzzyThresholdForS8(0.5);
  18 }
  19 //--------------------------------------------------------------------
  20
  21 //--------------------------------------------------------------------
  22 template <class ImageType>
  23 void
  24 clitk::ExtractLymphStationsFilter<ImageType>::
  25 ExtractStation_8_SI_Limits()
  26 {
  27   /*
  28     Station 8: paraeosphageal nodes
  29
  30     "Superiorly, Station 8 begins at the level of the carina and,
  31     therefore, abuts Station 3P (Fig. 3C). Inferiorly, the lower limit
  32     of Station 8 was not speciﬁed in the Mountain and Dresler
  33     classiﬁcation (1). We delineate this volume until it reaches the
  34     gastroesphogeal junction. "
  35
  36     => new classification IASCL 2009:
  37
  38     "Lower border: the diaphragm"
  39
  40     "Upper border: the upper border of the lower lobe bronchus on the
  41     left; the lower border of the bronchus intermedius on the right"
  42
  43   */
  44   StartNewStep("[Station8] Inf/Sup mediastinum limits with carina/diaphragm junction");
  45
  46   // Get Carina Z position
  47   MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
  48
  49   // CHANGE TO -> OriginOfRightMiddleLobeBronchus ???
  50
  51   std::vector<MaskImagePointType> centroids;
  52   clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
  53   m_CarinaZ = centroids[1][2];
  54   // add one slice to include carina ?
  55   m_CarinaZ += m_Mediastinum->GetSpacing()[2];
  56   // We dont need Carina structure from now
  57   Carina->Delete();
  58   GetAFDB()->SetDouble("CarinaZ", m_CarinaZ);
  59
  60   // Get left lower lobe bronchus (L), take the upper border
  61   // Get right bronchus intermedius (RML), take the lower border
  62
  63   /*
  64     double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2) +
  65     m_Mediastinum->GetSpacing()[2]; // add one slice to include carina
  66     // double GOjunctionZ = GetAFDB()->GetPoint3D("GOjunction", 2);
  67     */
  68
  69   // Found most inferior part of the lung
  70   MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
  71   // It should be already croped, so I took the origin and add 10mm above
  72   m_DiaphragmInferiorLimit = Lungs->GetOrigin()[2]+10;
  73   //  Lungs->Delete(); // we don't need it, release memory -> it we want to release, also free in AFDB
  74   clitk::PrintMemory(GetVerboseMemoryFlag(), "after reading lungs");
  75   GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
  76   clitk::PrintMemory(GetVerboseMemoryFlag(), "after release lungs");
  77
  78   /* Crop support :
  79      Superior limit = carina
  80      Inferior limit = DiaphragmInferiorLimit (old=Gastroesphogeal junction) */
  81   m_Working_Support =
  82     clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
  83                                                 m_DiaphragmInferiorLimit, //GOjunctionZ,
  84                                                 m_CarinaZ, true,
  85                                                 GetBackgroundValue());
  86
  87   // Remove some structures that we know are excluded
  88   MaskImagePointer VertebralBody =
  89     GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
  90   MaskImagePointer Aorta =
  91     GetAFDB()->template GetImage <MaskImageType>("Aorta");
  92
  93   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
  94   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
  95   boolFilter->InPlaceOn();
  96   boolFilter->SetInput1(m_Working_Support);
  97   boolFilter->SetInput2(VertebralBody);
  98   boolFilter->SetOperationType(BoolFilterType::AndNot);
  99   boolFilter->Update();
 100   boolFilter->SetInput1(boolFilter->GetOutput());
 101   boolFilter->SetInput2(Aorta);
 102   boolFilter->SetOperationType(BoolFilterType::AndNot);
 103   boolFilter->Update();
 104   m_Working_Support = boolFilter->GetOutput();
 105
 106   // Done.
 107   StopCurrentStep<MaskImageType>(m_Working_Support);
 108   m_ListOfStations["8"] = m_Working_Support;
 109 }
 110 //--------------------------------------------------------------------
 111
 112
 113 //--------------------------------------------------------------------
 114 template <class ImageType>
 115 void
 116 clitk::ExtractLymphStationsFilter<ImageType>::
 117 ExtractStation_8_Post_Limits()
 118 {
 119   /*
 120     Station 8: paraeosphageal nodes
 121
 122     Anteriorly, it is in contact with Station 7 and the
 123     left main stem bronchus in its superior aspect (Fig. 3C–H) and
 124     with the heart more inferiorly. Posteriorly, Station 8 abuts the
 125     descending aorta and the anterior aspect of the vertebral body
 126     until an imaginary horizontal line running 1 cm posterior to the
 127     anterior border of the vertebral body (Fig. 3C).
 128
 129     New classification IASCL 2009 :
 130
 131     "Nodes lying adjacent to the wall of the esophagus and to the
 132     right or left of the midline, excluding subcarinal nodes (S7)
 133     Upper"
 134
 135   */
 136
 137   StartNewStep("[Station8] Post limits with vertebral body");
 138   MaskImagePointer VertebralBody =
 139     GetAFDB()->template GetImage <MaskImageType>("VertebralBody");
 140
 141   // Consider vertebral body slice by slice
 142   typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
 143   typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 144   typedef typename ExtractSliceFilterType::SliceType SliceType;
 145   std::vector<typename SliceType::Pointer> vertebralSlices;
 146   extractSliceFilter->SetInput(VertebralBody);
 147   extractSliceFilter->SetDirection(2);
 148   extractSliceFilter->Update();
 149   extractSliceFilter->GetOutputSlices(vertebralSlices);
 150
 151   // For each slice, compute the most anterior point
 152   std::map<int, typename SliceType::PointType> vertebralAntPositionBySlice;
 153   for(uint i=0; i<vertebralSlices.size(); i++) {
 154     typename SliceType::PointType p;
 155     bool found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vertebralSlices[i],
 156                                                                      GetBackgroundValue(),
 157                                                                      1, true, p);
 158     if (found) {
 159       vertebralAntPositionBySlice[i] = p;
 160     }
 161     else {
 162       // no Foreground in this slice (it appends generally at the
 163       // beginning and the end of the volume). Do nothing in this
 164       // case.
 165     }
 166   }
 167
 168   // Convert 2D points in slice into 3D points
 169   std::vector<MaskImagePointType> vertebralAntPositions;
 170   clitk::PointsUtils<MaskImageType>::Convert2DTo3DList(vertebralAntPositionBySlice,
 171                                                        VertebralBody,
 172                                                        vertebralAntPositions);
 173
 174   // DEBUG : write list of points
 175   clitk::WriteListOfLandmarks<MaskImageType>(vertebralAntPositions,
 176                                              "vertebralMostAntPositions.txt");
 177
 178   // Cut support posteriorly 1cm the most anterior point of the
 179   // VertebralBody. Do nothing below and above the VertebralBody. To
 180   // do that compute several region, slice by slice and fill.
 181   MaskImageRegionType region;
 182   MaskImageSizeType size;
 183   MaskImageIndexType start;
 184   size[2] = 1; // a single slice
 185   start[0] = m_Working_Support->GetLargestPossibleRegion().GetIndex()[0];
 186   size[0] = m_Working_Support->GetLargestPossibleRegion().GetSize()[0];
 187   for(uint i=0; i<vertebralAntPositions.size(); i++) {
 188     typedef typename itk::ImageRegionIterator<MaskImageType> IteratorType;
 189     IteratorType iter =
 190       IteratorType(m_Working_Support, m_Working_Support->GetLargestPossibleRegion());
 191     MaskImageIndexType index;
 192     // Consider some cm posterior to most anterior positions (usually
 193     // 1 cm).
 194     vertebralAntPositions[i][1] += GetDistanceMaxToAnteriorPartOfTheSpine();
 195     // Get index of this point
 196     m_Working_Support->TransformPhysicalPointToIndex(vertebralAntPositions[i], index);
 197     // Compute region (a single slice)
 198     start[2] = index[2];
 199     start[1] = m_Working_Support->GetLargestPossibleRegion().GetIndex()[1]+index[1];
 200     size[1] = m_Working_Support->GetLargestPossibleRegion().GetSize()[1]-start[1];
 201     region.SetSize(size);
 202     region.SetIndex(start);
 203     // Fill region
 204     if (m_Working_Support->GetLargestPossibleRegion().IsInside(start))  {
 205       itk::ImageRegionIterator<MaskImageType> it(m_Working_Support, region);
 206       it.GoToBegin();
 207       while (!it.IsAtEnd()) {
 208         it.Set(GetBackgroundValue());
 209         ++it;
 210       }
 211     }
 212   }
 213
 214   StopCurrentStep<MaskImageType>(m_Working_Support);
 215   m_ListOfStations["8"] = m_Working_Support;
 216 }
 217 //--------------------------------------------------------------------
 218
 219
 220 //--------------------------------------------------------------------
 221 template <class ImageType>
 222 void
 223 clitk::ExtractLymphStationsFilter<ImageType>::
 224 ExtractStation_8_Ant_Limits()
 225 {
 226   //--------------------------------------------------------------------
 227   StartNewStep("[Station8] Ant limits with S7 above Carina");
 228   /*
 229     Anteriorly, it is in contact with Station 7 and the
 230     left main stem bronchus in its superior aspect (Fig. 3C–H) and
 231     with the heart more inferiorly.
 232
 233     1) line post wall bronchi (S7), till originRMLB
 234     - LUL bronchus : to detect in trachea. But not needed here ??
 235     2) heart ! -> not delineated.
 236     ==> below S7, crop CT not to far away (ant), then try with threshold
 237
 238     1) ==> how to share with S7 ? Prepare both support at the same time !
 239     NEED vector of initial support for each station ? No -> map if it exist before, take it !!
 240
 241   */
 242
 243   // Ant limit from carina (start) to end of S7 = originRMLB
 244   // Get Trachea
 245   MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
 246
 247   MaskImagePointer m_Working_Trachea =
 248     clitk::CropImageAbove<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
 249                                          GetBackgroundValue());
 250
 251   // Seprate into two main bronchi
 252   MaskImagePointer LeftBronchus;
 253   MaskImagePointer RightBronchus;
 254
 255   // Labelize and consider the two first (main) labels
 256   m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
 257
 258   // Carina position must at the first slice that separate the two
 259   // main bronchus (not superiorly). We thus first check that the
 260   // upper slice is composed of at least two labels
 261   typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
 262   SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
 263   iter.SetFirstDirection(0);
 264   iter.SetSecondDirection(1);
 265   iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
 266   int maxLabel=0;
 267   while (!iter.IsAtReverseEndOfSlice()) {
 268     while (!iter.IsAtReverseEndOfLine()) {
 269       if (iter.Get() > maxLabel) maxLabel = iter.Get();
 270       --iter;
 271     }
 272     iter.PreviousLine();
 273   }
 274   if (maxLabel < 2) {
 275     clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
 276   }
 277
 278   // Compute 3D centroids of both parts to identify the left from the
 279   // right bronchus
 280   std::vector<ImagePointType> c;
 281   clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
 282   ImagePointType C1 = c[1];
 283   ImagePointType C2 = c[2];
 284
 285   ImagePixelType leftLabel;
 286   ImagePixelType rightLabel;
 287   if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
 288   else { leftLabel = 2; rightLabel = 1; }
 289
 290   // Select LeftLabel (set one label to Backgroundvalue)
 291   LeftBronchus =
 292     SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
 293                                                 rightLabel, GetBackgroundValue(), false);
 294   RightBronchus  =
 295     SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
 296                                                 leftLabel, GetBackgroundValue(), false);
 297   // Crop images
 298   LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
 299   RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
 300
 301   // Insert int AFDB if need after
 302   GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
 303                                                LeftBronchus, true);
 304   GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
 305                                                RightBronchus, true);
 306
 307   // Now crop below OriginOfRightMiddleLobeBronchusZ
 308   // It is not done before to keep entire bronchi.
 309
 310   MaskImagePointer OriginOfRightMiddleLobeBronchus =
 311     GetAFDB()->template GetImage<MaskImageType>("OriginOfRightMiddleLobeBronchus");
 312   std::vector<MaskImagePointType> centroids;
 313   clitk::ComputeCentroids<MaskImageType>(OriginOfRightMiddleLobeBronchus, GetBackgroundValue(), centroids);
 314   m_OriginOfRightMiddleLobeBronchusZ = centroids[1][2];
 315   // add one slice to include carina ?
 316   m_OriginOfRightMiddleLobeBronchusZ += LeftBronchus->GetSpacing()[2];
 317   // We dont need Carina structure from now
 318   OriginOfRightMiddleLobeBronchus->Delete();
 319
 320   LeftBronchus =
 321     clitk::CropImageBelow<MaskImageType>(LeftBronchus, 2,
 322                                          m_OriginOfRightMiddleLobeBronchusZ,
 323                                          true, // AutoCrop
 324                                          GetBackgroundValue());
 325   RightBronchus =
 326     clitk::CropImageBelow<MaskImageType>(RightBronchus, 2,
 327                                          m_OriginOfRightMiddleLobeBronchusZ,
 328                                          true, // AutoCrop
 329                                          GetBackgroundValue());
 330
 331   // Search for points that are the most left/post/ant and most
 332   // right/post/ant of the left and right bronchus
 333   // 15  = not 15 mm more distance than the middle point.
 334   FindExtremaPointsInBronchus(LeftBronchus, 0, 10, m_RightMostInLeftBronchus,
 335                               m_AntMostInLeftBronchus, m_PostMostInLeftBronchus);
 336   FindExtremaPointsInBronchus(RightBronchus, 1, 10, m_LeftMostInRightBronchus,
 337                               m_AntMostInRightBronchus, m_PostMostInRightBronchus);
 338
 339   // DEBUG : write the list of points
 340   ListOfPointsType v;
 341   v.reserve(m_RightMostInLeftBronchus.size()+m_AntMostInLeftBronchus.size()+
 342             m_PostMostInLeftBronchus.size());
 343   v.insert(v.end(), m_RightMostInLeftBronchus.begin(), m_RightMostInLeftBronchus.end() );
 344   v.insert(v.end(), m_AntMostInLeftBronchus.begin(), m_AntMostInLeftBronchus.end() );
 345   v.insert(v.end(), m_PostMostInLeftBronchus.begin(), m_PostMostInLeftBronchus.end() );
 346   clitk::WriteListOfLandmarks<MaskImageType>(v, "LeftBronchusPoints.txt");
 347
 348   v.clear();
 349   v.reserve(m_LeftMostInRightBronchus.size()+m_AntMostInRightBronchus.size()+
 350             m_PostMostInRightBronchus.size());
 351   v.insert(v.end(), m_LeftMostInRightBronchus.begin(), m_LeftMostInRightBronchus.end() );
 352   v.insert(v.end(), m_AntMostInRightBronchus.begin(), m_AntMostInRightBronchus.end() );
 353   v.insert(v.end(), m_PostMostInRightBronchus.begin(), m_PostMostInRightBronchus.end() );
 354   clitk::WriteListOfLandmarks<MaskImageType>(v, "RightBronchusPoints.txt");
 355
 356
 357   // Now uses these points to limit, slice by slice
 358   // line is mainly horizontal, so mainDirection=1
 359   writeImage<MaskImageType>(m_Working_Support, "before.mhd");
 360   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 361                                                                     m_PostMostInLeftBronchus,
 362                                                                     m_PostMostInRightBronchus,
 363                                                                     GetBackgroundValue(), 1, 10);
 364   writeImage<MaskImageType>(m_Working_Support, "after.mhd");
 365
 366   // Keep main 3D CCL :
 367   m_Working_Support = Labelize<MaskImageType>(m_Working_Support, 0, false, 10);
 368   m_Working_Support = KeepLabels<MaskImageType>(m_Working_Support,
 369                                                 GetBackgroundValue(),
 370                                                 GetForegroundValue(), 1, 1, true);
 371
 372   // Autocrop
 373   m_Working_Support = clitk::AutoCrop<MaskImageType>(m_Working_Support, GetBackgroundValue());
 374
 375   // End of step
 376   StopCurrentStep<MaskImageType>(m_Working_Support);
 377   //  m_ListOfStations["8"] = m_Working_Support;
 378
 379 }
 380
 381 //--------------------------------------------------------------------
 382 template <class ImageType>
 383 void
 384 clitk::ExtractLymphStationsFilter<ImageType>::
 385 ExtractStation_8_LR_Limits()
 386 {
 387
 388   //--------------------------------------------------------------------
 389   StartNewStep("[Station8] Ant part (not post to Esophagus)");
 390   /*
 391     Consider Esophagus, dilate it and remove ant part. It remains part
 392     on L & R, than can be partly removed by cutting what remains at
 393     right of vertebral body.
 394   */
 395
 396   // Get Esophagus
 397   MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
 398   clitk::PrintMemory(GetVerboseMemoryFlag(), "after read Esophagus");
 399
 400   // In images from the original article, Atlas – UM, the oesophagus
 401   //was included in nodal stations 3p and 8.  Having said that, in the
 402   //description for station 8, it indicates that “the delineation of
 403   //station 8 is limited to the soft tissues surrounding the
 404   //oesophagus”.  In the recent article, The IASLC Lung Cancer Staging
 405   //Project (J Thorac Oncol 4:5, 568-77), the images drawn by
 406   //Dr. Aletta Frasier exclude this structure.  From an oncological
 407   //prospective, the oesophagus should be excluded from these nodal
 408   //stations.
 409   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
 410   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
 411   boolFilter->InPlaceOn();
 412   boolFilter->SetInput1(m_Working_Support);
 413   boolFilter->SetInput2(Esophagus);
 414   boolFilter->SetOperationType(BoolFilterType::AndNot);
 415   boolFilter->Update();
 416   m_Working_Support = boolFilter->GetOutput();
 417
 418   // Crop Esophagus : keep only below the OriginOfRightMiddleLobeBronchusZ
 419   Esophagus =
 420     clitk::CropImageAbove<MaskImageType>(Esophagus, 2,
 421                                          m_OriginOfRightMiddleLobeBronchusZ,
 422                                          true, // AutoCrop
 423                                          GetBackgroundValue());
 424
 425   // Dilate to keep only not Anterior positions
 426   MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt();
 427   Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
 428   Esophagus = clitk::Dilate<MaskImageType>(Esophagus,
 429                                            radiusInMM,
 430                                            GetBackgroundValue(),
 431                                            GetForegroundValue(), true);
 432   writeImage<MaskImageType>(Esophagus, "enlarged-eso.mhd");
 433
 434   // Remove Anterior part according to this dilatated esophagus. Note:
 435   // because we crop Esophagus with ORML, the support will also be
 436   // croped in the same way. Here it is a desired feature. If we dont
 437   // want, use SetIgnoreEmptySliceObject(true)
 438
 439   // In the new IASCL definition, it is not clear if sup limits is
 440   //  around carina or On the right, it is “the lower border of the
 441   //  bronchus intermedius”, indicated on the image set as a point
 442   //  (“lower border of the bronchus intermedius”)
 443
 444   typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
 445   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 446   relPosFilter->VerboseStepFlagOff();
 447   relPosFilter->SetInput(m_Working_Support);
 448   relPosFilter->SetInputObject(Esophagus);
 449   relPosFilter->AddOrientationTypeString("P");
 450   relPosFilter->InverseOrientationFlagOff();
 451   relPosFilter->SetDirection(2); // Z axis
 452   relPosFilter->UniqueConnectedComponentBySliceOff();
 453   relPosFilter->SetIntermediateSpacing(3);
 454   relPosFilter->ResampleBeforeRelativePositionFilterOn();
 455   relPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS8());
 456   relPosFilter->RemoveObjectFlagOff(); // Do not exclude here because it is dilated
 457   relPosFilter->Update();
 458   m_Working_Support = relPosFilter->GetOutput();
 459
 460   // AutoCrop (OR SbS ?)
 461   m_Working_Support = clitk::AutoCrop<MaskImageType>(m_Working_Support, GetBackgroundValue());
 462
 463   StopCurrentStep<MaskImageType>(m_Working_Support);
 464   //--------------------------------------------------------------------
 465   StartNewStep("[Station8] Left and Right limits arround esophagus (below Carina)");
 466
 467   // Estract slices for current support for slice by slice processing
 468   std::vector<typename MaskSliceType::Pointer> slices;
 469   clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices);
 470
 471   // Estract slices of Esophagus (resize like support before to have the same set of slices)
 472   MaskImagePointer EsophagusForSlice = clitk::ResizeImageLike<MaskImageType>(Esophagus, m_Working_Support, GetBackgroundValue());
 473   std::vector<typename MaskSliceType::Pointer> eso_slices;
 474   clitk::ExtractSlices<MaskImageType>(EsophagusForSlice, 2, eso_slices);
 475
 476   // Estract slices of Vertebral (resize like support before to have the same set of slices)
 477   MaskImagePointer VertebralBody = GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
 478   VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
 479   std::vector<typename MaskSliceType::Pointer> vert_slices;
 480   clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vert_slices);
 481
 482   // Estract slices of Aorta (resize like support before to have the same set of slices)
 483   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 484   Aorta = clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
 485   std::vector<typename MaskSliceType::Pointer> aorta_slices;
 486   clitk::ExtractSlices<MaskImageType>(Aorta, 2, aorta_slices);
 487
 488   // Extract slices of Mediastinum (resize like support before to have the same set of slices)
 489   m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
 490   m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Working_Support, GetBackgroundValue());
 491   std::vector<typename MaskSliceType::Pointer> mediast_slices;
 492   clitk::ExtractSlices<MaskImageType>(m_Mediastinum, 2, mediast_slices);
 493
 494   // List of points
 495   std::vector<MaskImagePointType> p_RightMostAnt;
 496   std::vector<MaskImagePointType> p_RightMostPost;
 497   std::vector<MaskImagePointType> p_LeftMostAnt;
 498   std::vector<MaskImagePointType> p_LeftMostPost;
 499   std::vector<MaskImagePointType> p_AllPoints;
 500
 501   /*
 502     In the following, we search for the LeftRight limits.  We search
 503     for the most Right points in Esophagus and in VertebralBody and
 504     consider a line between those to most right points. All points in
 505     the support which are most right to this line are discarded. Same
 506     for the left part. The underlying assumption is that the support
 507     is concave between Eso/VertebralBody. Esophagus is a bit
 508     dilatated. On VertebralBody we go right (or left) until we reach
 509     the lung (but no more 20 mm).
 510    */
 511
 512   // Loop slices
 513   MaskImagePointType p;
 514   MaskImagePointType pp;
 515   for(uint i=0; i<slices.size() ; i++) {
 516     // Declare all needed points (sp = slice point)
 517     typename MaskSliceType::PointType sp_maxRight_Eso;
 518     typename MaskSliceType::PointType sp_maxRight_Aorta;
 519     typename MaskSliceType::PointType sp_maxRight_Vertebra;
 520     typename MaskSliceType::PointType sp_maxLeft_Eso;
 521     typename MaskSliceType::PointType sp_maxLeft_Aorta;
 522     typename MaskSliceType::PointType sp_maxLeft_Vertebra;
 523
 524     // Right is at left on screen, coordinate decrease
 525     // Left is at right on screen, coordinate increase
 526
 527     // Find limit of Vertebral -> only at most Post part of current
 528     // slice support.  First found most ant point in VertebralBody
 529     typedef MaskSliceType SliceType;
 530     typename SliceType::PointType p_slice_ant;
 531     bool found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vert_slices[i], GetBackgroundValue(), 1, true, p_slice_ant);
 532     if (!found) {
 533       // It should not happen ! But sometimes, a contour is missing or
 534       // the VertebralBody is not delineated enough inferiorly ... in
 535       // those cases, we consider the first found slice.
 536       std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
 537       int j=i++;
 538       bool found = false;
 539       while (!found) {
 540         found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vert_slices[j], GetBackgroundValue(), 1, true, p_slice_ant);
 541         //clitkExceptionMacro("No foreground pixels in this VertebralBody slices ??");
 542         j++;
 543       }
 544     }
 545     p_slice_ant[1] += GetDistanceMaxToAnteriorPartOfTheSpine(); // Consider offset
 546
 547     // The, find most Right and Left points on that AP position
 548     typename SliceType::IndexType indexR;
 549     typename SliceType::IndexType indexL;
 550     vert_slices[i]->TransformPhysicalPointToIndex(p_slice_ant, indexR);
 551     indexL = indexR;
 552     // Check that is inside the mask
 553     indexR[1] = std::min(indexR[1], (long)vert_slices[i]->GetLargestPossibleRegion().GetSize()[1]-1);
 554     indexL[1] = indexR[1];
 555     while (vert_slices[i]->GetPixel(indexR) != GetBackgroundValue()) {
 556       indexR[0] --; // Go to the right
 557     }
 558     while (vert_slices[i]->GetPixel(indexL) != GetBackgroundValue()) {
 559       indexL[0] ++; // Go to the left
 560     }
 561     vert_slices[i]->TransformIndexToPhysicalPoint(indexR, sp_maxRight_Vertebra);
 562     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, VertebralBody, i, p);
 563     p_AllPoints.push_back(p);
 564     vert_slices[i]->TransformIndexToPhysicalPoint(indexL, sp_maxLeft_Vertebra);
 565     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, VertebralBody, i, p);
 566     p_AllPoints.push_back(p);
 567
 568     // Find last point out of the mediastinum on this line, Right :
 569     mediast_slices[i]->TransformPhysicalPointToIndex(sp_maxRight_Vertebra, indexR);
 570     double distance = 0.0;
 571     while (mediast_slices[i]->GetPixel(indexR) != GetBackgroundValue()) {
 572       indexR[0] --;
 573       distance += mediast_slices[i]->GetSpacing()[0];
 574     }
 575     if (distance < 30) { // Ok in this case, we found limit with lung
 576       mediast_slices[i]->TransformIndexToPhysicalPoint(indexR, sp_maxRight_Vertebra);
 577       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, m_Mediastinum, i, p);
 578     }
 579     else { // in that case, we are probably below the diaphragm, so we
 580            // add aribtrarly few mm
 581       sp_maxRight_Vertebra[0] -= 2; // Leave 2 mm around the VertebralBody
 582       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, m_Mediastinum, i, p);
 583     }
 584     p_RightMostPost.push_back(p);
 585     p_AllPoints.push_back(p);
 586
 587     // Find last point out of the mediastinum on this line, Left :
 588     mediast_slices[i]->TransformPhysicalPointToIndex(sp_maxLeft_Vertebra, indexL);
 589     distance = 0.0;
 590     while (mediast_slices[i]->GetPixel(indexL) != GetBackgroundValue()) {
 591       indexL[0] ++;
 592       distance += mediast_slices[i]->GetSpacing()[0];
 593     }
 594     if (distance < 30) { // Ok in this case, we found limit with lung
 595       mediast_slices[i]->TransformIndexToPhysicalPoint(indexL, sp_maxLeft_Vertebra);
 596       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, m_Mediastinum, i, p);
 597     }
 598     else { // in that case, we are probably below the diaphragm, so we
 599            // add aribtrarly few mm
 600       sp_maxLeft_Vertebra[0] += 2; // Leave 2 mm around the VertebralBody
 601       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, m_Mediastinum, i, p);
 602     }
 603     p_LeftMostPost.push_back(p);
 604     p_AllPoints.push_back(p);
 605
 606     // Find right limit of Esophagus and Aorta
 607     clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 0, true, sp_maxRight_Eso);
 608     bool f = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(), 0, true, sp_maxRight_Aorta);
 609     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Eso, EsophagusForSlice, i, p);
 610     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Aorta, Aorta, i, pp);
 611     pp[0] -= 2; // Add a margin of 2 mm to include the Aorta 'wall'
 612     p_AllPoints.push_back(p);
 613     if (f) {
 614       p_AllPoints.push_back(pp);
 615       MaskImagePointType A = p_RightMostPost.back();
 616       MaskImagePointType B = p;
 617       MaskImagePointType C = pp;
 618       double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 619       if (s>0) p_RightMostAnt.push_back(p);
 620       else p_RightMostAnt.push_back(pp);
 621     }
 622     else {
 623       p_RightMostAnt.push_back(p);
 624     }
 625
 626     // --------------------------------------------------------------------------
 627     // Find the limit on the Left: most left point between Eso and
 628     // Vertebra. (Left is left on screen, coordinate increase)
 629
 630     // Find left limit of Esophagus
 631     clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 0, false, sp_maxLeft_Eso);
 632     f = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(), 0, false, sp_maxLeft_Aorta);
 633     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Eso, EsophagusForSlice, i, p);
 634     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Aorta, Aorta, i, pp);
 635     p_AllPoints.push_back(p);
 636     pp[0] += 2; // Add a margin of 2 mm to include the 'wall'
 637     if (f) { // not below Aorta
 638       p_AllPoints.push_back(pp);
 639       MaskImagePointType A = p_LeftMostPost.back();
 640       MaskImagePointType B = p;
 641       MaskImagePointType C = pp;
 642       double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 643       if (s<0) p_LeftMostAnt.push_back(p); // Insert point most at Left
 644       else p_LeftMostAnt.push_back(pp);
 645     }
 646     else {
 647       p_LeftMostAnt.push_back(p);
 648     }
 649   } // End of slice loop
 650
 651   clitk::WriteListOfLandmarks<MaskImageType>(p_AllPoints, "LR-Eso-Vert.txt");
 652
 653   // Now uses these points to limit, slice by slice
 654   // Line is mainly vertical, so mainDirection=0
 655   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 656                                                                     p_RightMostAnt, p_RightMostPost,
 657                                                                     GetBackgroundValue(), 0, 10);
 658   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 659                                                                     p_LeftMostAnt, p_LeftMostPost,
 660                                                                     GetBackgroundValue(), 0, -10);
 661   // END
 662   StopCurrentStep<MaskImageType>(m_Working_Support);
 663   m_ListOfStations["8"] = m_Working_Support;
 664   return;
 665 }
 666 //--------------------------------------------------------------------
 667
 668
 669 //--------------------------------------------------------------------
 670 template <class ImageType>
 671 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
 672 clitk::ExtractLymphStationsFilter<ImageType>::
 673 EnlargeEsophagusDilatationRadiusInferiorly(MaskImagePointer & Esophagus)
 674 {
 675   // Check if Esophagus is delineated at least until GOjunction. Now,
 676   // because we use AutoCrop, Origin[2] gives this max inferior
 677   // position.
 678   //  double GOjunctionZ = GetAFDB()->GetPoint3D("GOjunction", 2);
 679
 680   if (Esophagus->GetOrigin()[2] > m_DiaphragmInferiorLimit) {
 681     std::cout << "Warning Esophagus is not delineated until Diaphragm. I mirror-pad it."
 682               << std::endl;
 683     double extraSize = Esophagus->GetOrigin()[2]-m_DiaphragmInferiorLimit;
 684
 685     // Pad with few more slices
 686     typedef itk::MirrorPadImageFilter<MaskImageType, MaskImageType> PadFilterType;
 687     typename PadFilterType::Pointer padFilter = PadFilterType::New();
 688     padFilter->SetInput(Esophagus);
 689     MaskImageSizeType b;
 690     b[0] = 0; b[1] = 0; b[2] = (uint)ceil(extraSize/Esophagus->GetSpacing()[2])+1;
 691     padFilter->SetPadLowerBound(b);
 692     padFilter->Update();
 693     Esophagus = padFilter->GetOutput();
 694   }
 695   return Esophagus;
 696 }
 697 //--------------------------------------------------------------------
 698
 699
 700 //--------------------------------------------------------------------
 701 template <class ImageType>
 702 void
 703 clitk::ExtractLymphStationsFilter<ImageType>::
 704 ExtractStation_8_LR_Limits_old()
 705 {
 706   /*
 707     Station 8: paraeosphageal nodes
 708
 709     Laterally, it is within the pleural envelope and again abuts the
 710     descending aorta on the left. Reasonably, the delineation of
 711     Station 8 is limited to the soft tissue surrounding the esophagus
 712     (Fig.  3C–H).
 713   */
 714
 715   StartNewStep("[Station8] Right limits (around esophagus)");
 716   // Get Esophagus
 717   MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
 718
 719   // Autocrop to get first slice with starting Esophagus
 720   Esophagus = clitk::AutoCrop<MaskImageType>(Esophagus, GetBackgroundValue());
 721
 722   // Dilate
 723   // LR dilatation -> large to keep point inside
 724   // AP dilatation -> few mm
 725   // SI dilatation -> enough to cover Diaphragm (old=GOjunctionZ)
 726   MaskImagePointType radiusInMM = GetEsophagusDiltationForRight();
 727   Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
 728   Esophagus = clitk::Dilate<MaskImageType>(Esophagus,
 729                                            radiusInMM,
 730                                            GetBackgroundValue(),
 731                                            GetForegroundValue(), true);
 732   writeImage<MaskImageType>(Esophagus, "dilateEso2.mhd");
 733
 734   writeImage<MaskImageType>(m_Working_Support, "before-relpos2.mhd");
 735
 736   // Remove Right (Left on the screen) part according to this
 737   // dilatated esophagus
 738   typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
 739   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 740   relPosFilter->VerboseStepFlagOff();
 741   relPosFilter->WriteStepFlagOff();
 742   relPosFilter->SetInput(m_Working_Support);
 743   relPosFilter->SetInputObject(Esophagus);
 744   relPosFilter->AddOrientationTypeString("L");
 745   relPosFilter->InverseOrientationFlagOn(); // Not Left to
 746   relPosFilter->SetDirection(2); // Z axis
 747   relPosFilter->UniqueConnectedComponentBySliceOff(); // important
 748   relPosFilter->SetIntermediateSpacing(4);
 749   relPosFilter->ResampleBeforeRelativePositionFilterOn();
 750   relPosFilter->SetFuzzyThreshold(0.9); // remove few part only
 751   relPosFilter->RemoveObjectFlagOff();
 752   relPosFilter->Update();
 753   m_Working_Support = relPosFilter->GetOutput();
 754
 755   // Get a single 3D CCL
 756   m_Working_Support = Labelize<MaskImageType>(m_Working_Support, 0, false, 10);
 757   m_Working_Support = KeepLabels<MaskImageType>(m_Working_Support,
 758                                                 GetBackgroundValue(),
 759                                                 GetForegroundValue(), 1, 1, true);
 760
 761
 762   /*
 763   // Re-Add post to Esophagus -> sometimes previous relpos remove some
 764   // correct part below esophagus.
 765   MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
 766   EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
 767   writeImage<MaskImageType>(Esophagus, "e-again.mhd");
 768   typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
 769   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 770   relPosFilter->VerboseStepOff();
 771   relPosFilter->WriteStepOff();
 772   relPosFilter->SetInput(m_Working_Support);
 773   relPosFilter->SetInputObject(Esophagus);
 774   relPosFilter->SetOrientationTypeString("P");
 775   relPosFilter->InverseOrientationFlagOff();
 776   relPosFilter->SetDirection(2); // Z axis
 777   relPosFilter->UniqueConnectedComponentBySliceOff(); // important
 778   relPosFilter->SetIntermediateSpacing(4);
 779   relPosFilter->ResampleBeforeRelativePositionFilterOn();
 780   relPosFilter->CombineWithOrFlagOn();
 781   relPosFilter->SetFuzzyThreshold(0.9); // remove few part only
 782   relPosFilter->RemoveObjectFlagOff();
 783   relPosFilter->Update();
 784   m_Working_Support = relPosFilter->GetOutput();
 785   */
 786
 787   StopCurrentStep<MaskImageType>(m_Working_Support);
 788   m_ListOfStations["8"] = m_Working_Support;
 789 }
 790 //--------------------------------------------------------------------