segmentation/clitkExtractLymphStation_8.txx

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