segmentation/clitkExtractLymphStation_4RL.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================*/
  18 //--------------------------------------------------------------------
  19 template <class TImageType>
  20 void
  21 clitk::ExtractLymphStationsFilter<TImageType>::
  22 ExtractStation_4RL_SI_Limits()
  23 {
  24
  25   /*
  26     Station 4R: right lower paratracheal nodes From superior to
  27     inferior, the delineation of Station 4R starts at the top of the
  28     aortic arch (Fig. 2D) and ends at the upper lobe bronchus or where the
  29     right pulmonary artery crosses the midline of the mediastinum
  30     (Fig. 3E,F). On the left side, Station 4R is deﬁned by the midline of
  31     the trachea (Fig. 2D). On the right side, it is contained within the
  32     pleural envelope in the upper part, medial to the superior vena cava
  33     and the arch of the azygos vein in the intermediate section (Fig. 2I
  34     and 3A,B) and the right upper lobe pulmonary vein in its very caudal
  35     part. Anteriorly, it is limited most supe- riorly by the right
  36     brachiocephalic vein (Fig. 2D–H), fol- lowed by the superior vena cava
  37     and the arch or ascending section of the aorta (Figs. 2I and 3A–E). In
  38     between the superior vena cava and the aorta, we recommend delineating
  39     Station 4R so that it extends halfway between the two vessels where it
  40     will contact Station 3A or 6 (Figs. 2H,I and 3A–D). Posteriorly,
  41     Station 4R is deﬁned at its superior extent by an imaginary horizontal
  42     line running along the posterior wall of the trachea
  43     (Fig. 2E). Inferiorly, it remains anterior to the right main stem
  44     bronchus, ﬁlling the soft- tissue space between the vessels.
  45   */
  46
  47   StartNewStep("[Station 4R]Inf/Sup mediastinum limits with aortic arch/upperLBronchus");
  48   /* SupInf limits :
  49      - top of aortic arch
  50      - ends at the upper lobe bronchus or where the right pulmonary artery crosses the midline of the mediastinum
  51   */
  52
  53   // Local variables
  54   double m_TopOfAorticArchInMM;
  55   double m_UpperLobeBronchusZPositionInMM;
  56   double m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM;
  57
  58   // Get Inputs
  59   m_TopOfAorticArchInMM = GetAFDB()->GetPoint3D("TopOfAorticArch", 2);
  60   DD(m_TopOfAorticArchInMM);
  61   m_UpperLobeBronchusZPositionInMM = GetAFDB()->GetPoint3D("RightUpperLobeBronchus", 2);
  62   DD(m_UpperLobeBronchusZPositionInMM);
  63   m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM = GetAFDB()->GetPoint3D("RightPulmoArteryCrossesMidMediastinum", 2);
  64   DD(m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM);
  65
  66   /* Crop support */
  67   double inf = std::max(m_UpperLobeBronchusZPositionInMM, m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM);
  68   m_Working_Support =
  69     clitk::CropImageAlongOneAxis<MaskImageType>(m_Mediastinum, 2,
  70                                                 inf,
  71                                                 m_TopOfAorticArchInMM, true,
  72                                                 GetBackgroundValue());
  73   StopCurrentStep<MaskImageType>(m_Working_Support);
  74 }
  75 //--------------------------------------------------------------------
  76
  77
  78 //--------------------------------------------------------------------
  79 template <class TImageType>
  80 void
  81 clitk::ExtractLymphStationsFilter<TImageType>::
  82 ExtractStation_4RL_LR_Limits()
  83 {
  84   // 4R first
  85
  86   // Left : midline of the trachea
  87   // Right : "- upper part : contained within the pleural envelope
  88   //- intermediate section : medial to the superior vena cava and the arch of the azygos vein
  89   // - very caudal part : right upper lobe pulmonary vein"
  90   // AAV ??
  91
  92   // Left -> midline of the trachea
  93   // slice by slice, find X coord of 2D centroid (?)
  94   // check with previous line in order to not move too fast
  95
  96
  97   // ----------------------------------------------------------
  98   StartNewStep("[Station 4R] Left limits with midline of trachea ");
  99
 100   /*
 101     Two possible approaches
 102     1) use trachea skeleton (clitkExtractAirwaysTreeInfo)
 103     -> but need to analyse the tree to remove "false" bifurcation
 104     -> need to track from top to bottom, with each bronchus
 105     -> how to stay on the main path ?
 106
 107     2) analyse slice by slice the trachea, labelize and take the first
 108     or two first connected components -> find centroids.
 109     -> not really "smooth" when bronchus slicing lead to "splated" region
 110
 111     ==> we choose 2
 112   */
 113
 114   // Crop the trachea like the current support
 115   MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
 116   MaskImagePointer crop_trachea =
 117     clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
 118   writeImage<MaskImageType>(crop_trachea, "croptrachea.mhd");
 119
 120   // Extract all the slices
 121   std::vector<typename MaskSliceType::Pointer> bronchi_slices;
 122   clitk::ExtractSlices<MaskImageType>(crop_trachea, 2, bronchi_slices);
 123
 124   // List of midpoints
 125   std::vector<MaskImagePointType> midpoints;
 126
 127   // Add mid points below carina, from foot to head
 128   for(uint i=0; i<m_LeftMostInRightBronchus.size(); i++) {
 129     MaskImagePointType p;
 130     p[0] = (m_LeftMostInRightBronchus[i][0]+m_RightMostInLeftBronchus[i][0])/2.0;
 131     p[1] = (m_LeftMostInRightBronchus[i][1]+m_RightMostInLeftBronchus[i][1])/2.0;
 132     p[2] = m_LeftMostInRightBronchus[i][2];
 133     midpoints.push_back(p);
 134   }
 135
 136   // Crop the trachea above the carina. Find largest Z position
 137   MaskImageIndexType p = Trachea->GetLargestPossibleRegion().GetIndex();
 138   for(uint i=0; i<3; i++) {
 139     p[i] += Trachea->GetLargestPossibleRegion().GetSize()[i];
 140   }
 141   MaskImagePointType q;
 142   Trachea->TransformIndexToPhysicalPoint(p, q);
 143   double maxZ = q[2];
 144   double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2);
 145   MaskImagePointer m_above_carina =
 146     clitk::CropImageAlongOneAxis<MaskImageType>(Trachea, 2,
 147                                                 m_CarinaZ,
 148                                                 maxZ, true,
 149                                                 GetBackgroundValue());
 150   writeImage<MaskImageType>(m_above_carina, "above.mhd");
 151
 152   // Extract all the slices
 153   std::vector<typename MaskSliceType::Pointer> trachea_slices;
 154   clitk::ExtractSlices<MaskImageType>(m_above_carina, 2, trachea_slices);
 155
 156   // Find centroid of the trachea in each slice
 157   std::vector<MaskImagePointType> points;
 158   typedef typename MaskSliceType::PointType SlicePointType;
 159   SlicePointType previous;
 160   // Start from patient top (head)
 161   for(uint i=0; i<trachea_slices.size(); i++) {
 162     // Labelize
 163     trachea_slices[i] = Labelize<MaskSliceType>(trachea_slices[i], GetBackgroundValue(), true, 10);
 164     // Get centroid
 165     std::vector<SlicePointType> c;
 166     clitk::ComputeCentroids<MaskSliceType>(trachea_slices[i], GetBackgroundValue(), c);
 167     // Keep first one (first connected component label) and convert to 3D
 168     MaskImagePointType p;
 169     p[2] = m_above_carina->GetOrigin()[2] + i*m_above_carina->GetSpacing()[2];
 170     p[0] = c[1][0];
 171     p[1] = c[1][1];
 172     midpoints.push_back(p);
 173   }
 174
 175   // DEBUG POINTS
 176   std::ofstream osp;
 177   openFileForWriting(osp, "mp.txt");
 178   osp << "LANDMARKS1" << std::endl;
 179   for(uint i=0; i<midpoints.size(); i++) {
 180     osp << i << " " << midpoints[i][0] << " "
 181         << midpoints[i][1] << " " << midpoints[i][2] << " 0 0 " << std::endl;
 182   }
 183   osp.close();
 184
 185   // Create and allocate left and right part of the support (we create
 186   // images because we will then crop them)
 187   m_LeftSupport = clitk::NewImageLike<MaskImageType>(m_Working_Support);
 188   m_RightSupport = clitk::NewImageLike<MaskImageType>(m_Working_Support);
 189
 190   // Loop on current support, slice by slice
 191   typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
 192   SliceIteratorType iter = SliceIteratorType(m_Working_Support,
 193                                              m_Working_Support->GetLargestPossibleRegion());
 194   SliceIteratorType iterL = SliceIteratorType(m_LeftSupport,
 195                                               m_LeftSupport->GetLargestPossibleRegion());
 196   SliceIteratorType iterR = SliceIteratorType(m_RightSupport,
 197                                               m_RightSupport->GetLargestPossibleRegion());
 198   iter.SetFirstDirection(0);
 199   iter.SetSecondDirection(1);
 200   iter.GoToBegin();
 201   iterL.SetFirstDirection(0);
 202   iterL.SetSecondDirection(1);
 203   iterL.GoToBegin();
 204   iterR.SetFirstDirection(0);
 205   iterR.SetSecondDirection(1);
 206   iterR.GoToBegin();
 207
 208   int slice=0;
 209   // Assert starting of image has the same Z than first midpoints
 210   while (midpoints[slice][2] != m_Working_Support->GetOrigin()[2]) {
 211     slice++;
 212     if ((uint)slice >= midpoints.size()) {
 213       clitkExceptionMacro("Bug while searching for first midpoint to use");
 214     }
 215   }
 216
 217   // Start loop
 218   while (!iter.IsAtEnd()) {
 219     ImageIndexType index;
 220     m_Working_Support->TransformPhysicalPointToIndex(midpoints[slice], index);
 221     while (!iter.IsAtEndOfSlice()) {
 222       // convert into index
 223       while (!iter.IsAtEndOfLine()) {
 224         // if indexcourant <index this is Right part. Left otherwise
 225         if (iter.GetIndex()[0]<index[0]) {
 226           iterR.Set(iter.Get());
 227           iterL.Set(GetBackgroundValue());
 228         }
 229         else {
 230           iterL.Set(iter.Get());
 231           iterR.Set(GetBackgroundValue());
 232         }
 233         ++iter;
 234         ++iterL;
 235         ++iterR;
 236       }
 237       iter.NextLine();
 238       iterL.NextLine();
 239       iterR.NextLine();
 240
 241     }
 242     iter.NextSlice();
 243     iterL.NextSlice();
 244     iterR.NextSlice();
 245     ++slice;
 246   }
 247
 248   // Auto crop !
 249   m_LeftSupport = clitk::AutoCrop<MaskImageType>(m_LeftSupport, GetBackgroundValue());
 250   m_RightSupport = clitk::AutoCrop<MaskImageType>(m_RightSupport, GetBackgroundValue());
 251   writeImage<MaskImageType>(m_LeftSupport, "lsac.mhd");
 252   writeImage<MaskImageType>(m_RightSupport, "rsac.mhd");
 253
 254   StopCurrentStep<MaskImageType>(m_RightSupport);
 255
 256   // -------------------------------------------------------------------------
 257   // Constraint at right from the SVC. Must be after MidTrachea limits
 258   // (if not, bronchus can be split, midposition is false)
 259   StartNewStep("[Station 4R] R limits with SVC ");
 260   MaskImagePointer svc = GetAFDB()->template GetImage<MaskImageType>("SVC");
 261   typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
 262   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 263   /*  relPosFilter->SetVerboseStep(false);
 264   relPosFilter->SetInput(m_RightSupport);  // only right here ...
 265   relPosFilter->SetInputObject(svc);
 266   relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
 267   relPosFilter->SetIntermediateSpacing(2);
 268   relPosFilter->SetFuzzyThreshold(0.3);
 269   relPosFilter->Update();
 270   m_RightSupport = relPosFilter->GetOutput();
 271   */
 272
 273   /*
 274     ==> TODO RIGHT TO MEDIAL TO SVC :
 275     get centroid, cut in X direction to get medial
 276
 277     ==> REDO OPERATOR to find extrma points ?
 278     centroid or skeleton ?
 279
 280   */
 281
 282   relPosFilter = RelPosFilterType::New();
 283   relPosFilter->VerboseStepFlagOff();
 284   relPosFilter->SetInput(m_RightSupport);  // only right here ...
 285   relPosFilter->SetInputObject(svc);
 286   relPosFilter->AddOrientationType(RelPosFilterType::AntTo);
 287   relPosFilter->InverseOrientationFlagOn();
 288   relPosFilter->SetIntermediateSpacing(2); // this is important to put it low
 289   relPosFilter->SetFuzzyThreshold(0.6);
 290   relPosFilter->Update();
 291   m_RightSupport = relPosFilter->GetOutput();
 292
 293   // AutoCrop
 294   m_RightSupport = clitk::AutoCrop<MaskImageType>(m_RightSupport, GetBackgroundValue());
 295
 296 }
 297 //--------------------------------------------------------------------
 298
 299
 300 //--------------------------------------------------------------------
 301 template <class TImageType>
 302 void
 303 clitk::ExtractLymphStationsFilter<TImageType>::
 304 ExtractStation_4RL_AP_Limits()
 305 {
 306
 307   /*
 308     post of S4R
 309     - sup part -> cut line post wall trachea
 310     - inf part (where ???) -> remains anterior to the right main stem bronchus
 311     ==> bo
 312   */
 313
 314   // Post (not Ant) to Aorta
 315   MaskImagePointer aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 316   // Crop according to current support
 317   aorta = clitk::ResizeImageLike<MaskImageType>(aorta, m_RightSupport, GetBackgroundValue());
 318   typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
 319   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 320   relPosFilter = RelPosFilterType::New();
 321   relPosFilter->VerboseStepFlagOff();
 322   relPosFilter->SetInput(m_RightSupport);  // only right here ...
 323   relPosFilter->SetInputObject(aorta);
 324   relPosFilter->AddOrientationType(RelPosFilterType::PostTo);
 325   //  relPosFilter->NotFlagOn();
 326   relPosFilter->SetIntermediateSpacing(2); // this is important to put it low
 327   relPosFilter->SetFuzzyThreshold(0.6);
 328   relPosFilter->Update();
 329   m_RightSupport = relPosFilter->GetOutput();
 330
 331   m_ListOfStations["4R"] = m_RightSupport;
 332   StopCurrentStep<MaskImageType>(m_RightSupport);
 333
 334   // POST -> horizontal lines bronchus --> points dectection in S7 to store.
 335 }
 336 //--------------------------------------------------------------------
 337