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