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