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("PostTo");
 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 or 3P] Slice by slice, keep a single CCL (the closest to VertebralBody)");
 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     if ((c.size() > 1) && (centroids.size() > 1)) {
 658       // keep the one which is the closer to vertebralbody centroid
 659       double distance=1000000;
 660       int index=0;
 661       for(uint j=1; j<centroids.size(); j++) {
 662         double d = centroids[j].EuclideanDistanceTo(c[1]);
 663         if (d<distance) {
 664           distance = d;
 665           index = j;
 666         }
 667       }
 668       // remove all others label
 669       slices[i] = KeepLabels<MaskSliceType>(slices[i], GetBackgroundValue(),
 670                                             GetForegroundValue(), index, index, true);
 671     }
 672   }
 673   m_Working_Support = clitk::JoinSlices<MaskImageType>(slices, m_Working_Support, 2);
 674
 675   StopCurrentStep<MaskImageType>(m_Working_Support);
 676   m_ListOfStations["8"] = m_Working_Support;
 677 }
 678 //--------------------------------------------------------------------
 679
 680
 681 //--------------------------------------------------------------------
 682 template <class ImageType>
 683 void
 684 clitk::ExtractLymphStationsFilter<ImageType>::
 685 ExtractStation_8_LR_Limits_old2()
 686 {
 687
 688   //--------------------------------------------------------------------
 689   StartNewStep("[Station8] Left and Right limits arround esophagus (below Carina)");
 690
 691   // Estract slices for current support for slice by slice processing
 692   std::vector<typename MaskSliceType::Pointer> slices;
 693   clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices);
 694
 695   // Dilate the Esophagus to consider a margins around
 696   MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt();
 697   m_Esophagus = clitk::Dilate<MaskImageType>(m_Esophagus,
 698                                              radiusInMM,
 699                                              GetBackgroundValue(),
 700                                              GetForegroundValue(), true);
 701
 702   // Remove what is outside the mediastinum in this enlarged Esophagus -> it allows to select
 703   // 'better' extrema points (not too post).
 704   MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
 705   clitk::AndNot<MaskImageType>(m_Esophagus, Lungs, GetBackgroundValue());
 706   GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
 707
 708   // Estract slices of Esophagus (resize like support before to have the same set of slices)
 709   MaskImagePointer EsophagusForSlice = clitk::ResizeImageLike<MaskImageType>(m_Esophagus, m_Working_Support, GetBackgroundValue());
 710
 711   std::vector<typename MaskSliceType::Pointer> eso_slices;
 712   clitk::ExtractSlices<MaskImageType>(EsophagusForSlice, 2, eso_slices);
 713
 714   // Estract slices of Vertebral (resize like support before to have the same set of slices)
 715   MaskImagePointer VertebralBody = GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
 716   VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
 717   std::vector<typename MaskSliceType::Pointer> vert_slices;
 718   clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vert_slices);
 719
 720   // Estract slices of Aorta (resize like support before to have the same set of slices)
 721   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 722   Aorta = clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
 723   std::vector<typename MaskSliceType::Pointer> aorta_slices;
 724   clitk::ExtractSlices<MaskImageType>(Aorta, 2, aorta_slices);
 725
 726   // Extract slices of Mediastinum (resize like support before to have the same set of slices)
 727   m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
 728   m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Working_Support, GetBackgroundValue());
 729   std::vector<typename MaskSliceType::Pointer> mediast_slices;
 730   clitk::ExtractSlices<MaskImageType>(m_Mediastinum, 2, mediast_slices);
 731
 732   // List of points
 733   std::vector<MaskImagePointType> p_RightMostAnt;
 734   std::vector<MaskImagePointType> p_RightMostPost;
 735   std::vector<MaskImagePointType> p_LeftMostAnt;
 736   std::vector<MaskImagePointType> p_LeftMostPost;
 737   std::vector<MaskImagePointType> p_AllPoints;
 738   std::vector<MaskImagePointType> p_LeftAorta;
 739   std::vector<MaskImagePointType> p_LeftEso;
 740
 741   /*
 742     In the following, we search for the LeftRight limits.  We search
 743     for the most Right points in Esophagus and in VertebralBody and
 744     consider a line between those to most right points. All points in
 745     the support which are most right to this line are discarded. Same
 746     for the left part. The underlying assumption is that the support
 747     is concave between Eso/VertebralBody. Esophagus is a bit
 748     dilatated. On VertebralBody we go right (or left) until we reach
 749     the lung (but no more 20 mm).
 750   */
 751
 752   // Loop slices
 753   MaskImagePointType p;
 754   MaskImagePointType pp;
 755   for(uint i=0; i<slices.size() ; i++) {
 756     // Declare all needed points (sp = slice point)
 757     typename MaskSliceType::PointType sp_maxRight_Eso;
 758     typename MaskSliceType::PointType sp_maxRight_Aorta;
 759     typename MaskSliceType::PointType sp_maxRight_Vertebra;
 760     typename MaskSliceType::PointType sp_maxLeft_Eso;
 761     typename MaskSliceType::PointType sp_maxLeft_Aorta;
 762     typename MaskSliceType::PointType sp_maxLeft_Vertebra;
 763
 764     // Right is at left on screen, coordinate decrease
 765     // Left is at right on screen, coordinate increase
 766
 767     // Find limit of Vertebral -> only at most Post part of current
 768     // slice support.  First found most ant point in VertebralBody
 769     typedef MaskSliceType SliceType;
 770     typename SliceType::PointType p_slice_ant;
 771     bool found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vert_slices[i], GetBackgroundValue(), 1, true, p_slice_ant);
 772     if (!found) {
 773       // It should not happen ! But sometimes, a contour is missing or
 774       // the VertebralBody is not delineated enough inferiorly ... in
 775       // those cases, we consider the first found slice.
 776       std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
 777       int j=i++;
 778       bool found = false;
 779       while (!found) {
 780         found = clitk::FindExtremaPointInAGivenDirection<SliceType>(vert_slices[j], GetBackgroundValue(), 1, true, p_slice_ant);
 781         //clitkExceptionMacro("No foreground pixels in this VertebralBody slices ??");
 782         j++;
 783       }
 784     }
 785     p_slice_ant[1] += GetDistanceMaxToAnteriorPartOfTheSpine(); // Consider offset
 786
 787     // The, find most Right and Left points on that AP position
 788     typename SliceType::IndexType indexR;
 789     typename SliceType::IndexType indexL;
 790     vert_slices[i]->TransformPhysicalPointToIndex(p_slice_ant, indexR);
 791     indexL = indexR;
 792     // Check that is inside the mask
 793     indexR[1] = std::min(indexR[1], (long)vert_slices[i]->GetLargestPossibleRegion().GetSize()[1]-1);
 794     indexL[1] = indexR[1];
 795     while (vert_slices[i]->GetPixel(indexR) != GetBackgroundValue()) {
 796       indexR[0] --; // Go to the right
 797     }
 798     while (vert_slices[i]->GetPixel(indexL) != GetBackgroundValue()) {
 799       indexL[0] ++; // Go to the left
 800     }
 801     vert_slices[i]->TransformIndexToPhysicalPoint(indexR, sp_maxRight_Vertebra);
 802     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, VertebralBody, i, p);
 803     p_AllPoints.push_back(p);
 804     vert_slices[i]->TransformIndexToPhysicalPoint(indexL, sp_maxLeft_Vertebra);
 805     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, VertebralBody, i, p);
 806     p_AllPoints.push_back(p);
 807
 808     // Find last point out of the mediastinum on this line, Right :
 809     mediast_slices[i]->TransformPhysicalPointToIndex(sp_maxRight_Vertebra, indexR);
 810     double distance = 0.0;
 811     while (mediast_slices[i]->GetPixel(indexR) != GetBackgroundValue()) {
 812       indexR[0] --;
 813       distance += mediast_slices[i]->GetSpacing()[0];
 814     }
 815     if (distance < 30) { // Ok in this case, we found limit with lung
 816       mediast_slices[i]->TransformIndexToPhysicalPoint(indexR, sp_maxRight_Vertebra);
 817       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, m_Mediastinum, i, p);
 818     }
 819     else { // in that case, we are probably below the diaphragm, so we
 820            // add aribtrarly few mm
 821       sp_maxRight_Vertebra[0] -= 2; // Leave 2 mm around the VertebralBody
 822       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Vertebra, m_Mediastinum, i, p);
 823     }
 824     p_RightMostPost.push_back(p);
 825     p_AllPoints.push_back(p);
 826
 827     // Find last point out of the mediastinum on this line, Left :
 828     mediast_slices[i]->TransformPhysicalPointToIndex(sp_maxLeft_Vertebra, indexL);
 829     distance = 0.0;
 830     while (mediast_slices[i]->GetPixel(indexL) != GetBackgroundValue()) {
 831       indexL[0] ++;
 832       distance += mediast_slices[i]->GetSpacing()[0];
 833     }
 834     if (distance < 30) { // Ok in this case, we found limit with lung
 835       mediast_slices[i]->TransformIndexToPhysicalPoint(indexL, sp_maxLeft_Vertebra);
 836       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, m_Mediastinum, i, p);
 837     }
 838     else { // in that case, we are probably below the diaphragm, so we
 839            // add aribtrarly few mm
 840       sp_maxLeft_Vertebra[0] += 2; // Leave 2 mm around the VertebralBody
 841       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Vertebra, m_Mediastinum, i, p);
 842     }
 843     p_LeftMostPost.push_back(p);
 844     p_AllPoints.push_back(p);
 845
 846     // Find Eso slice centroid and do not consider what is post to
 847     // this centroid.
 848     std::vector<typename MaskSliceType::PointType> c;
 849     clitk::ComputeCentroids<MaskSliceType>(eso_slices[i], GetBackgroundValue(), c);
 850     if (c.size() >1) {
 851       eso_slices[i] =
 852         clitk::CropImageAbove<MaskSliceType>(eso_slices[i], 1, c[1][1], false, GetBackgroundValue());
 853       eso_slices[i] =
 854         clitk::ResizeImageLike<MaskSliceType>(eso_slices[i], aorta_slices[i], GetBackgroundValue());
 855       // writeImage<MaskSliceType>(eso_slices[i], "eso-slice-"+toString(i)+".mhd");
 856     }
 857
 858     // Find right limit of Esophagus and Aorta
 859     bool f =
 860       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(),
 861                                                               0, true, sp_maxRight_Eso);
 862     f = f &&
 863       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(),
 864                                                               0, true, sp_maxRight_Aorta);
 865     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Eso, EsophagusForSlice, i, p);
 866     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxRight_Aorta, Aorta, i, pp);
 867     pp[0] -= 2; // Add a margin of 2 mm to include the Aorta 'wall'
 868     p_AllPoints.push_back(p);
 869     if (f) {
 870       p_AllPoints.push_back(pp);
 871       MaskImagePointType A = p_RightMostPost.back();
 872       MaskImagePointType B = p;
 873       MaskImagePointType C = pp;
 874       double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 875       if (s>0) p_RightMostAnt.push_back(p);
 876       else p_RightMostAnt.push_back(pp);
 877     }
 878     else { // No more Esophagus in this slice : do nothing
 879       //      p_RightMostAnt.push_back(p);
 880       p_RightMostPost.pop_back();
 881     }
 882
 883     // --------------------------------------------------------------------------
 884     // Find the limit on the Left: most left point between Eso and
 885     // Eso. (Left is left on screen, coordinate increase)
 886
 887     // Find left limit of Esophagus
 888     clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 0, false, sp_maxLeft_Eso);
 889     f = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(), 0, false, sp_maxLeft_Aorta);
 890     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Eso, EsophagusForSlice, i, p);
 891     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_maxLeft_Aorta, Aorta, i, pp);
 892     p_AllPoints.push_back(p);
 893     pp[0] += 2; // Add a margin of 2 mm to include the 'wall'
 894     if (f) { // not below Aorta
 895       p_AllPoints.push_back(pp);
 896       MaskImagePointType A = p_LeftMostPost.back();
 897       MaskImagePointType B = p;
 898       MaskImagePointType C = pp;
 899       double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 900       if (s<0) {
 901         p_LeftMostAnt.push_back(p); // Insert point most at Left, Eso
 902       }
 903       else {
 904         // in this case -> two lines !
 905         p_LeftMostAnt.push_back(pp);  // Insert point most at Left, Aorta (Vert to Aorta)
 906         // but also consider Aorta to Eso
 907         p_LeftAorta.push_back(pp);
 908         p_LeftEso.push_back(p);
 909       }
 910     }
 911     else { // No more Esophagus in this slice : do nothing
 912       p_LeftMostPost.pop_back();
 913       //p_LeftMostAnt.push_back(p);
 914     }
 915   } // End of slice loop
 916
 917   clitk::WriteListOfLandmarks<MaskImageType>(p_AllPoints, "S8-LR-Eso-Vert.txt");
 918   clitk::WriteListOfLandmarks<MaskImageType>(p_LeftEso, "S8-Left-Eso.txt");
 919   clitk::WriteListOfLandmarks<MaskImageType>(p_LeftAorta, "S8-Left-Aorta.txt");
 920
 921   // Now uses these points to limit, slice by slice
 922   // Line is mainly vertical, so mainDirection=0
 923   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 924                                                                     p_RightMostAnt, p_RightMostPost,
 925                                                                     GetBackgroundValue(), 0, 10);
 926   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 927                                                                     p_LeftMostAnt, p_LeftMostPost,
 928                                                                     GetBackgroundValue(), 0, -10);
 929   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 930                                                                     p_LeftEso,p_LeftAorta,
 931                                                                     GetBackgroundValue(), 0, -10);
 932   // END
 933   StopCurrentStep<MaskImageType>(m_Working_Support);
 934   m_ListOfStations["8"] = m_Working_Support;
 935   return;
 936 }
 937 //--------------------------------------------------------------------
 938
 939
 940 //--------------------------------------------------------------------
 941 template <class ImageType>
 942 void
 943 clitk::ExtractLymphStationsFilter<ImageType>::
 944 ExtractStation_8_LR_Limits()
 945 {
 946
 947   //--------------------------------------------------------------------
 948   StartNewStep("[Station8] Left and Right limits arround esophagus with lines from VertebralBody-Aorta-Esophagus");
 949
 950   // Estract slices for current support for slice by slice processing
 951   std::vector<typename MaskSliceType::Pointer> slices;
 952   clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices);
 953
 954   // Dilate the Esophagus to consider a margins around
 955   MaskImagePointType radiusInMM = GetEsophagusDiltationForAnt();
 956   m_Esophagus = clitk::Dilate<MaskImageType>(m_Esophagus, radiusInMM,
 957                                              GetBackgroundValue(), GetForegroundValue(), true);
 958   //  m_Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(m_Esophagus);
 959
 960   // Remove what is outside the mediastinum in this enlarged Esophagus -> it allows to select
 961   // 'better' extrema points (not too post).
 962   MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
 963   clitk::AndNot<MaskImageType>(m_Esophagus, Lungs, GetBackgroundValue());
 964   GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
 965
 966   // Estract slices of Esophagus (resize like support before to have the same set of slices)
 967   MaskImagePointer EsophagusForSlice = clitk::ResizeImageLike<MaskImageType>(m_Esophagus, m_Working_Support, GetBackgroundValue());
 968   std::vector<typename MaskSliceType::Pointer> eso_slices;
 969   clitk::ExtractSlices<MaskImageType>(EsophagusForSlice, 2, eso_slices);
 970
 971   // Estract slices of Vertebral (resize like support before to have the same set of slices)
 972   MaskImagePointer VertebralBody = GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
 973   VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
 974   std::vector<typename MaskSliceType::Pointer> vert_slices;
 975   clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, vert_slices);
 976
 977   // Estract slices of Aorta (resize like support before to have the same set of slices)
 978   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 979   Aorta = clitk::ResizeImageLike<MaskImageType>(Aorta, m_Working_Support, GetBackgroundValue());
 980   std::vector<typename MaskSliceType::Pointer> aorta_slices;
 981   clitk::ExtractSlices<MaskImageType>(Aorta, 2, aorta_slices);
 982
 983   // Extract slices of Mediastinum (resize like support before to have the same set of slices)
 984   m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
 985   m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Working_Support, GetBackgroundValue());
 986   std::vector<typename MaskSliceType::Pointer> mediast_slices;
 987   clitk::ExtractSlices<MaskImageType>(m_Mediastinum, 2, mediast_slices);
 988
 989   // List of points
 990   std::vector<MaskImagePointType> p_MostLeftVertebralBody;
 991   std::vector<MaskImagePointType> p_MostRightVertebralBody;
 992   std::vector<MaskImagePointType> p_MostLeftAorta;
 993   std::vector<MaskImagePointType> p_MostLeftEsophagus;
 994
 995   /*
 996     In the following, we search for the LeftRight limits.  We search
 997     for the most Right points in Esophagus and in VertebralBody and
 998     consider a line between those to most right points. All points in
 999     the support which are most right to this line are discarded. Same
1000     for the left part. The underlying assumption is that the support
1001     is concave between Eso/VertebralBody. Esophagus is a bit
1002     dilatated. On VertebralBody we go right (or left) until we reach
1003     the lung (but no more 20 mm).
1004   */
1005
1006   // Temporary 3D point
1007   MaskImagePointType p;
1008
1009   // Loop slices
1010   for(uint i=0; i<slices.size() ; i++) {
1011     // Declare all needed 2D points (sp = slice point)
1012     typename MaskSliceType::PointType sp_MostRightVertebralBody;
1013     typename MaskSliceType::PointType sp_MostLeftVertebralBody;
1014     typename MaskSliceType::PointType sp_MostLeftAorta;
1015     typename MaskSliceType::PointType sp_temp;
1016     typename MaskSliceType::PointType sp_MostLeftEsophagus;
1017
1018     // Right is at left on screen, coordinate decrease
1019     // Left is at right on screen, coordinate increase
1020
1021     /* -------------------------------------------------------------------------------------
1022        Find most Left point in VertebralBody. Consider only the
1023        horizontal line which is 'DistanceMaxToAnteriorPartOfTheSpine'
1024        away from the most ant point.
1025      */
1026     typename MaskSliceType::PointType sp_MostAntVertebralBody;
1027     bool found = false;
1028     int j=i;
1029     while (!found) {
1030       found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vert_slices[j], GetBackgroundValue(), 1, true, sp_MostAntVertebralBody);
1031       if (!found) {
1032         // It should not happen ! But sometimes, a contour is missing or
1033         // the VertebralBody is not delineated enough inferiorly ... in
1034         // those cases, we consider the first found slice.
1035         std::cerr << "No foreground pixels in this VertebralBody slices !?? I try with the previous/next slice" << std::endl;
1036         j++;
1037       }
1038     }
1039     sp_MostAntVertebralBody[1] += GetDistanceMaxToAnteriorPartOfTheSpine(); // Consider offset
1040
1041     // Crop the vertebralbody below this most post line
1042     vert_slices[j] =
1043       clitk::CropImageAbove<MaskSliceType>(vert_slices[j], 1, sp_MostAntVertebralBody[1], false, GetBackgroundValue());
1044     vert_slices[j] =
1045       clitk::ResizeImageLike<MaskSliceType>(vert_slices[j], aorta_slices[i], GetBackgroundValue());
1046     //    writeImage<MaskSliceType>(vert_slices[i], "vert-slice-"+toString(i)+".mhd");
1047
1048     // Find first point not in mediastinum
1049     clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vert_slices[j], GetBackgroundValue(), 0, false, sp_MostLeftVertebralBody);
1050     sp_MostLeftVertebralBody = clitk::FindExtremaPointInAGivenLine<MaskSliceType>(mediast_slices[i], 0, false, sp_MostLeftVertebralBody, GetBackgroundValue(), 30);
1051     sp_MostLeftVertebralBody[0] += 2; // 2mm margin
1052     clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(vert_slices[j], GetBackgroundValue(), 0, true, sp_MostRightVertebralBody);
1053     sp_MostRightVertebralBody = clitk::FindExtremaPointInAGivenLine<MaskSliceType>(mediast_slices[i], 0, true, sp_MostRightVertebralBody, GetBackgroundValue(),30);
1054     sp_MostRightVertebralBody[0] -= 2; // 2 mm margin
1055
1056     // Convert 2D points into 3D
1057     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostRightVertebralBody, VertebralBody, i, p);
1058     p_MostRightVertebralBody.push_back(p);
1059     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftVertebralBody, VertebralBody, i, p);
1060     p_MostLeftVertebralBody.push_back(p);
1061
1062     /* -------------------------------------------------------------------------------------
1063        Find most Left point in Esophagus
1064      */
1065     found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(eso_slices[i], GetBackgroundValue(), 0, false, sp_MostLeftEsophagus);
1066     if (!found) { // No more Esophagus, I remove the previous point
1067
1068       // if (p_MostLeftEsophagus.size() < 1)  {
1069       p_MostLeftVertebralBody.pop_back();
1070       // }
1071       // else {
1072       //   // Consider the previous point
1073       //   p = p_MostLeftEsophagus.back();
1074       //   p_MostLeftEsophagus.push_back(p);
1075       //   sp_MostLeftEsophagus = sp_temp; // Retrieve previous 2D position
1076       //   found = true;
1077       // }
1078     }
1079     else {
1080       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftEsophagus, EsophagusForSlice, i, p);
1081       p_MostLeftEsophagus.push_back(p);
1082       // sp_temp = sp_MostLeftEsophagus; // Store previous 2D position
1083     }
1084
1085     /* -------------------------------------------------------------------------------------
1086        Find most Left point in Aorta
1087     */
1088     if (found) {
1089       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(aorta_slices[i], GetBackgroundValue(), 0, false, sp_MostLeftAorta);
1090       sp_MostLeftAorta = clitk::FindExtremaPointInAGivenLine<MaskSliceType>(mediast_slices[i], 0, false, sp_MostLeftAorta, GetBackgroundValue(), 10);
1091       typename MaskSliceType::PointType temp=sp_MostLeftEsophagus;
1092       temp[0] -= 10;
1093       if (clitk::IsOnTheSameLineSide(sp_MostLeftAorta, sp_MostLeftVertebralBody, sp_MostLeftEsophagus, temp)) {
1094         // sp_MostLeftAorta is on the same side than temp (at Right) -> ignore Aorta
1095         sp_MostLeftAorta = sp_MostLeftEsophagus;
1096       }
1097       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp_MostLeftAorta, Aorta, i, p);
1098       p_MostLeftAorta.push_back(p);
1099     }
1100
1101   } // End of slice loop
1102
1103   clitk::WriteListOfLandmarks<MaskImageType>(p_MostLeftVertebralBody, "S8-MostLeft-VB-points.txt");
1104   clitk::WriteListOfLandmarks<MaskImageType>(p_MostRightVertebralBody, "S8-MostRight-VB-points.txt");
1105   clitk::WriteListOfLandmarks<MaskImageType>(p_MostLeftAorta, "S8-MostLeft-Aorta-points.txt");
1106   clitk::WriteListOfLandmarks<MaskImageType>(p_MostLeftEsophagus, "S8-MostLeft-eso-points.txt");
1107
1108   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
1109                                                                     p_MostLeftVertebralBody, p_MostLeftAorta,
1110                                                                     GetBackgroundValue(), 0, -10);
1111
1112   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
1113                                                                     p_MostLeftAorta, p_MostLeftEsophagus,
1114                                                                     GetBackgroundValue(), 0, -10);
1115
1116   // END
1117   StopCurrentStep<MaskImageType>(m_Working_Support);
1118   m_ListOfStations["8"] = m_Working_Support;
1119   return;
1120 }
1121 //--------------------------------------------------------------------
1122
1123
1124 //--------------------------------------------------------------------
1125 template <class ImageType>
1126 void
1127 clitk::ExtractLymphStationsFilter<ImageType>::
1128 ExtractStation_8_Remove_Structures()
1129 {
1130
1131   //--------------------------------------------------------------------
1132   StartNewStep("[Station8] remove some structures");
1133
1134   Remove_Structures("Aorta");
1135   Remove_Structures("Esophagus");
1136
1137   // END
1138   StopCurrentStep<MaskImageType>(m_Working_Support);
1139   m_ListOfStations["8"] = m_Working_Support;
1140   return;
1141 }
1142 //--------------------------------------------------------------------
1143
1144
1145 //--------------------------------------------------------------------
1146 template <class ImageType>
1147 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1148 clitk::ExtractLymphStationsFilter<ImageType>::
1149 EnlargeEsophagusDilatationRadiusInferiorly(MaskImagePointer & Esophagus)
1150 {
1151   // Check if Esophagus is delineated at least until Diaphragm. Now,
1152   // because we use AutoCrop, Origin[2] gives this max inferior
1153   // position.
1154
1155   DD("BUGGY DONT USE");
1156   exit(0);
1157
1158   if (Esophagus->GetOrigin()[2] > m_DiaphragmInferiorLimit) {
1159     // crop first slice without mask
1160     MaskImagePointType pt;
1161     clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Esophagus, GetBackgroundValue(), 2, true, pt);
1162     DD(pt);
1163     Esophagus =
1164       clitk::CropImageBelow<MaskImageType>(Esophagus, 2,
1165                                            pt[2],
1166                                            false, // AutoCrop
1167                                            GetBackgroundValue());
1168     writeImage<MaskImageType>(Esophagus, "crop-eso.mhd");
1169
1170     std::cout << "Warning Esophagus is not delineated until Diaphragm. I mirror-pad it."
1171               << std::endl;
1172     double extraSize = Esophagus->GetOrigin()[2]-m_DiaphragmInferiorLimit;
1173
1174     // Pad with few more slices
1175     typedef itk::MirrorPadImageFilter<MaskImageType, MaskImageType> PadFilterType;
1176     typename PadFilterType::Pointer padFilter = PadFilterType::New();
1177     padFilter->SetInput(Esophagus);
1178     MaskImageSizeType b;
1179     b[0] = 0; b[1] = 0; b[2] = (uint)ceil(extraSize/Esophagus->GetSpacing()[2])+1;
1180     padFilter->SetPadLowerBound(b);
1181     padFilter->Update();
1182     Esophagus = padFilter->GetOutput();
1183   }
1184   return Esophagus;
1185 }
1186 //--------------------------------------------------------------------
1187
1188
1189 //--------------------------------------------------------------------
1190 template <class ImageType>
1191 void
1192 clitk::ExtractLymphStationsFilter<ImageType>::
1193 ExtractStation_8_LR_Limits_old()
1194 {
1195   /*
1196     Station 8: paraeosphageal nodes
1197
1198     Laterally, it is within the pleural envelope and again abuts the
1199     descending aorta on the left. Reasonably, the delineation of
1200     Station 8 is limited to the soft tissue surrounding the esophagus
1201     (Fig.  3C–H).
1202   */
1203
1204   StartNewStep("[Station8] Right limits (around esophagus)");
1205   // Get Esophagus
1206   MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
1207
1208   // Autocrop to get first slice with starting Esophagus
1209   Esophagus = clitk::AutoCrop<MaskImageType>(Esophagus, GetBackgroundValue());
1210
1211   // Dilate
1212   // LR dilatation -> large to keep point inside
1213   // AP dilatation -> few mm
1214   // SI dilatation -> enough to cover Diaphragm (old=GOjunctionZ)
1215   MaskImagePointType radiusInMM = GetEsophagusDiltationForRight();
1216   Esophagus = EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
1217   Esophagus = clitk::Dilate<MaskImageType>(Esophagus,
1218                                            radiusInMM,
1219                                            GetBackgroundValue(),
1220                                            GetForegroundValue(), true);
1221   writeImage<MaskImageType>(Esophagus, "dilateEso2.mhd");
1222
1223   writeImage<MaskImageType>(m_Working_Support, "before-relpos2.mhd");
1224
1225   // Remove Right (Left on the screen) part according to this
1226   // dilatated esophagus
1227   typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
1228   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
1229   relPosFilter->VerboseStepFlagOff();
1230   relPosFilter->WriteStepFlagOff();
1231   relPosFilter->SetInput(m_Working_Support);
1232   relPosFilter->SetInputObject(Esophagus);
1233   relPosFilter->AddOrientationTypeString("NotLeftTo");
1234   //  relPosFilter->InverseOrientationFlagOn(); // Not Left to
1235   relPosFilter->SetDirection(2); // Z axis
1236   relPosFilter->UniqueConnectedComponentBySliceOff(); // important
1237   relPosFilter->SetIntermediateSpacing(4);
1238   relPosFilter->IntermediateSpacingFlagOn();
1239   relPosFilter->SetFuzzyThreshold(0.9); // remove few part only
1240   relPosFilter->RemoveObjectFlagOff();
1241   relPosFilter->Update();
1242   m_Working_Support = relPosFilter->GetOutput();
1243
1244   // Get a single 3D CCL
1245   m_Working_Support = Labelize<MaskImageType>(m_Working_Support, 0, false, 10);
1246   m_Working_Support = KeepLabels<MaskImageType>(m_Working_Support,
1247                                                 GetBackgroundValue(),
1248                                                 GetForegroundValue(), 1, 1, true);
1249
1250
1251   /*
1252   // Re-Add post to Esophagus -> sometimes previous relpos remove some
1253   // correct part below esophagus.
1254   MaskImagePointer Esophagus = GetAFDB()->template GetImage<MaskImageType>("Esophagus");
1255   EnlargeEsophagusDilatationRadiusInferiorly(Esophagus);
1256   writeImage<MaskImageType>(Esophagus, "e-again.mhd");
1257   typedef clitk::SliceBySliceRelativePositionFilter<MaskImageType> RelPosFilterType;
1258   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
1259   relPosFilter->VerboseStepOff();
1260   relPosFilter->WriteStepOff();
1261   relPosFilter->SetInput(m_Working_Support);
1262   relPosFilter->SetInputObject(Esophagus);
1263   relPosFilter->SetOrientationTypeString("P");
1264   relPosFilter->InverseOrientationFlagOff();
1265   relPosFilter->SetDirection(2); // Z axis
1266   relPosFilter->UniqueConnectedComponentBySliceOff(); // important
1267   relPosFilter->SetIntermediateSpacing(4);
1268   relPosFilter->IntermediateSpacingFlagOn();
1269   relPosFilter->CombineWithOrFlagOn();
1270   relPosFilter->SetFuzzyThreshold(0.9); // remove few part only
1271   relPosFilter->RemoveObjectFlagOff();
1272   relPosFilter->Update();
1273   m_Working_Support = relPosFilter->GetOutput();
1274   */
1275
1276   StopCurrentStep<MaskImageType>(m_Working_Support);
1277   m_ListOfStations["8"] = m_Working_Support;
1278 }
1279 //--------------------------------------------------------------------
1280
1281
1282 //--------------------------------------------------------------------
1283 template <class TImageType>
1284 void
1285 clitk::ExtractLymphStationsFilter<TImageType>::
1286 FindExtremaPointsInBronchus(MaskImagePointer input,
1287                             int direction,
1288                             double distance_max_from_center_point,
1289                             ListOfPointsType & LR,
1290                             ListOfPointsType & Ant,
1291                             ListOfPointsType & Post)
1292 {
1293
1294   // Other solution ==> with auto bounding box ! (but pb to prevent to
1295   // be too distant from the center point
1296
1297   // Extract slices
1298   std::vector<typename MaskSliceType::Pointer> slices;
1299   clitk::ExtractSlices<MaskImageType>(input, 2, slices);
1300
1301   // Loop on slices
1302   bool found;
1303   for(uint i=0; i<slices.size(); i++) {
1304     /*
1305     // Keep main CCL
1306     slices[i] = Labelize<MaskSliceType>(slices[i], 0, true, 10);
1307     slices[i] = KeepLabels<MaskSliceType>(slices[i],
1308     GetBackgroundValue(),
1309     GetForegroundValue(), 1, 1, true);
1310     */
1311
1312     // ------- Find rightmost or leftmost point  -------
1313     MaskSliceType::PointType LRMost;
1314     found =
1315       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i],
1316                                                               GetBackgroundValue(),
1317                                                               0, // axis XY
1318                                                               (direction==0?false:true),  // right or left according to direction
1319                                                               LRMost);
1320     // ------- Find postmost point  -------
1321     MaskSliceType::PointType postMost;
1322     found =
1323       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i],
1324                                                               GetBackgroundValue(),
1325                                                               1, false, LRMost,
1326                                                               distance_max_from_center_point,
1327                                                               postMost);
1328     // ------- Find antmost point  -------
1329     MaskSliceType::PointType antMost;
1330     found =
1331       clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i],
1332                                                               GetBackgroundValue(),
1333                                                               1, true, LRMost,
1334                                                               distance_max_from_center_point,
1335                                                               antMost);
1336     // Only add point if found
1337     if (found)  {
1338       // ------- Convert 2D to 3D points --------
1339       MaskImageType::PointType p;
1340       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(LRMost, input, i, p);
1341       LR.push_back(p);
1342       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(antMost, input, i, p);
1343       Ant.push_back(p);
1344       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(postMost, input, i, p);
1345       Post.push_back(p);
1346     }
1347   }
1348 }
1349 //--------------------------------------------------------------------