segmentation/clitkExtractLymphStation_4RL.txx

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