X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=segmentation%2FclitkExtractLymphStation_4RL.txx;h=462ae9d34e142a7d8a40f1ec7168166730dabfce;hb=9a79f8cf7020e0a7ed5b1b97293bfa13a3a27b58;hp=a4478a6c4621a1b69804e9e8a9ef4e732b520ef3;hpb=5ff46b355b778ae675ca60085430c062ea6a3945;p=clitk.git

diff --git a/segmentation/clitkExtractLymphStation_4RL.txx b/segmentation/clitkExtractLymphStation_4RL.txx
index a4478a6..462ae9d 100644
--- a/segmentation/clitkExtractLymphStation_4RL.txx
+++ b/segmentation/clitkExtractLymphStation_4RL.txx
@@ -4,9 +4,33 @@ void
 clitk::ExtractLymphStationsFilter<TImageType>::
 ExtractStation_4RL_SI_Limits() 
 {
+
+  /*
+    Station 4R: right lower paratracheal nodes From superior to
+    inferior, the delineation of Station 4R starts at the top of the
+    aortic arch (Fig. 2D) and ends at the upper lobe bronchus or where the
+    right pulmonary artery crosses the midline of the mediastinum
+    (Fig. 3E,F). On the left side, Station 4R is defined by the midline of
+    the trachea (Fig. 2D). On the right side, it is contained within the
+    pleural envelope in the upper part, medial to the superior vena cava
+    and the arch of the azygos vein in the intermediate section (Fig. 2I
+    and 3A,B) and the right upper lobe pulmonary vein in its very caudal
+    part. Anteriorly, it is limited most supe- riorly by the right
+    brachiocephalic vein (Fig. 2D–H), fol- lowed by the superior vena cava
+    and the arch or ascending section of the aorta (Figs. 2I and 3A–E). In
+    between the superior vena cava and the aorta, we recommend delineating
+    Station 4R so that it extends halfway between the two vessels where it
+    will contact Station 3A or 6 (Figs. 2H,I and 3A–D). Posteriorly,
+    Station 4R is defined at its superior extent by an imaginary horizontal
+    line running along the posterior wall of the trachea
+    (Fig. 2E). Inferiorly, it remains anterior to the right main stem
+    bronchus, filling the soft- tissue space between the vessels.
+  */
+
+  StartNewStep("[Station 4R]Inf/Sup mediastinum limits with aortic arch/upperLBronchus");
   /* SupInf limits : 
-  - top of aortic arch
-  - ends at the upper lobe bronchus or where the right pulmonary artery crosses the midline of the mediastinum
+     - top of aortic arch
+     - ends at the upper lobe bronchus or where the right pulmonary artery crosses the midline of the mediastinum
   */
 
   // Local variables
@@ -15,24 +39,21 @@ ExtractStation_4RL_SI_Limits()
   double m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM;
 
   // Get Inputs
-  m_TopOfAorticArchInMM = GetAFDB()->GetPoint3D("topOfAorticArch", 2);
+  m_TopOfAorticArchInMM = GetAFDB()->GetPoint3D("TopOfAorticArch", 2);
   DD(m_TopOfAorticArchInMM);
-  m_UpperLobeBronchusZPositionInMM = GetAFDB()->GetPoint3D("rightUpperLobeBronchus", 2);
+  m_UpperLobeBronchusZPositionInMM = GetAFDB()->GetPoint3D("RightUpperLobeBronchus", 2);
   DD(m_UpperLobeBronchusZPositionInMM);
-  m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM = GetAFDB()->GetPoint3D("rightPulmoArteryCrossesMidMediastinum", 2);
+  m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM = GetAFDB()->GetPoint3D("RightPulmoArteryCrossesMidMediastinum", 2);
   DD(m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM);
 
   /* Crop support */
-  StartNewStep("Inf/Sup mediastinum limits with aortic arch/upperLBronchus");
   double inf = std::max(m_UpperLobeBronchusZPositionInMM, m_RightPulmoArteyrCrossesMidMediastinumZPositionInMM);
   m_Working_Support = 
-    clitk::CropImageAlongOneAxis<MaskImageType>(m_Support, 2, 
+    clitk::CropImageAlongOneAxis<MaskImageType>(m_Mediastinum, 2, 
                                                 inf,
                                                 m_TopOfAorticArchInMM, true,
                                                 GetBackgroundValue());
   StopCurrentStep<MaskImageType>(m_Working_Support);
-
-  m_Station4RL = m_Working_Support;
 }
 //--------------------------------------------------------------------
 
@@ -51,50 +72,249 @@ ExtractStation_4RL_LR_Limits()
   // - very caudal part : right upper lobe pulmonary vein"
   // AAV ??
   
-  // Constraint at right from the SVC
-  MaskImagePointer svc = GetAFDB()->template GetImage<MaskImageType>("SVC");
-  typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
-  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
-  relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
-  relPosFilter->VerboseStepOff();
-  relPosFilter->WriteStepOff();
-  relPosFilter->SetInput(m_Working_Support); 
-  relPosFilter->SetInputObject(svc); 
-  relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
-  relPosFilter->SetIntermediateSpacing(2);
-  relPosFilter->SetFuzzyThreshold(0.3);
-  relPosFilter->Update();
-  m_Working_Support = relPosFilter->GetOutput();
-  m_Station4RL = m_Working_Support;
-  
   // Left -> midline of the trachea
   // slice by slice, find X coord of 2D centroid (?)
   // check with previous line in order to not move too fast
+
+
+  // ----------------------------------------------------------
+  StartNewStep("[Station 4R] Left limits with midline of trachea ");
   
-  // skeleton ? -> need path description ? follow from slice to slice
-  // OR CENTROID at each slice ?
-  
-  // Crop trachea
-  // Extract list of slice from trachea
-  // Loop slice -> Get centroid crop along line (BB limit) -> two supports
+  /*
+    Two possible approaches
+    1) use trachea skeleton (clitkExtractAirwaysTreeInfo)
+    -> but need to analyse the tree to remove "false" bifurcation
+    -> need to track from top to bottom, with each bronchus
+    -> how to stay on the main path ?
+       
+    2) analyse slice by slice the trachea, labelize and take the first
+    or two first connected components -> find centroids. 
+    -> not really "smooth" when bronchus slicing lead to "splated" region
+       
+    ==> we choose 2
+  */
 
   // Crop the trachea like the current support
+  MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");    
   MaskImagePointer crop_trachea = 
-    clitk::ResizeImageLike<MaskImageType>(m_Trachea, m_Working_Support, GetBackgroundValue());
+    clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
   writeImage<MaskImageType>(crop_trachea, "croptrachea.mhd");
 
   // Extract all the slices
-  typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
-  typedef typename ExtractSliceFilterType::SliceType SliceType;
-  typename ExtractSliceFilterType::Pointer 
-    extractSliceFilter = ExtractSliceFilterType::New();
-  extractSliceFilter->SetInput(crop_trachea);
-  extractSliceFilter->SetDirection(2);
-  extractSliceFilter->Update();
-  std::vector<typename SliceType::Pointer> trachea_slices;
-  extractSliceFilter->GetOutputSlices(trachea_slices);
+  std::vector<typename MaskSliceType::Pointer> bronchi_slices;
+  clitk::ExtractSlices<MaskImageType>(crop_trachea, 2, bronchi_slices);
+
+  // List of midpoints
+  std::vector<MaskImagePointType> midpoints;
+
+  // Add mid points below carina, from foot to head
+  for(uint i=0; i<m_LeftMostInRightBronchus.size(); i++) {
+    MaskImagePointType p;
+    p[0] = (m_LeftMostInRightBronchus[i][0]+m_RightMostInLeftBronchus[i][0])/2.0;
+    p[1] = (m_LeftMostInRightBronchus[i][1]+m_RightMostInLeftBronchus[i][1])/2.0;
+    p[2] = m_LeftMostInRightBronchus[i][2];
+    midpoints.push_back(p);
+  }
+
+  // Crop the trachea above the carina. Find largest Z position
+  MaskImageIndexType p = Trachea->GetLargestPossibleRegion().GetIndex();
+  for(uint i=0; i<3; i++) {
+    p[i] += Trachea->GetLargestPossibleRegion().GetSize()[i];
+  }
+  MaskImagePointType q;
+  Trachea->TransformIndexToPhysicalPoint(p, q);
+  double maxZ = q[2];
+  double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2);
+  MaskImagePointer m_above_carina = 
+    clitk::CropImageAlongOneAxis<MaskImageType>(Trachea, 2, 
+                                                m_CarinaZ,
+                                                maxZ, true,
+                                                GetBackgroundValue());
+  writeImage<MaskImageType>(m_above_carina, "above.mhd");
+
+  // Extract all the slices
+  std::vector<typename MaskSliceType::Pointer> trachea_slices;
+  clitk::ExtractSlices<MaskImageType>(m_above_carina, 2, trachea_slices);
+
+  // Find centroid of the trachea in each slice
+  std::vector<MaskImagePointType> points;
+  typedef typename MaskSliceType::PointType SlicePointType;
+  SlicePointType previous;
+  // Start from patient top (head)
+  for(uint i=0; i<trachea_slices.size(); i++) {
+    // Labelize 
+    trachea_slices[i] = Labelize<MaskSliceType>(trachea_slices[i], GetBackgroundValue(), true, 10);
+    // Get centroid
+    std::vector<SlicePointType> c;
+    clitk::ComputeCentroids<MaskSliceType>(trachea_slices[i], GetBackgroundValue(), c);
+    // Keep first one (first connected component label) and convert to 3D
+    MaskImagePointType p;
+    p[2] = m_above_carina->GetOrigin()[2] + i*m_above_carina->GetSpacing()[2];
+    p[0] = c[1][0];
+    p[1] = c[1][1];
+    midpoints.push_back(p);
+  }
+
+  // DEBUG POINTS
+  std::ofstream osp;
+  openFileForWriting(osp, "mp.txt");
+  osp << "LANDMARKS1" << std::endl;
+  for(uint i=0; i<midpoints.size(); i++) {
+    osp << i << " " << midpoints[i][0] << " " 
+	<< midpoints[i][1] << " " << midpoints[i][2] << " 0 0 " << std::endl;
+  }
+  osp.close();
+
+  // Create and allocate left and right part of the support (we create
+  // images because we will then crop them)
+  m_LeftSupport = clitk::NewImageLike<MaskImageType>(m_Working_Support);
+  m_RightSupport = clitk::NewImageLike<MaskImageType>(m_Working_Support);
+
+  // Loop on current support, slice by slice
+  typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
+  SliceIteratorType iter = SliceIteratorType(m_Working_Support, 
+                                             m_Working_Support->GetLargestPossibleRegion());
+  SliceIteratorType iterL = SliceIteratorType(m_LeftSupport, 
+					      m_LeftSupport->GetLargestPossibleRegion());
+  SliceIteratorType iterR = SliceIteratorType(m_RightSupport, 
+					      m_RightSupport->GetLargestPossibleRegion());
+  iter.SetFirstDirection(0);
+  iter.SetSecondDirection(1);
+  iter.GoToBegin();
+  iterL.SetFirstDirection(0);
+  iterL.SetSecondDirection(1);
+  iterL.GoToBegin();
+  iterR.SetFirstDirection(0);
+  iterR.SetSecondDirection(1);
+  iterR.GoToBegin();
+
+  int slice=0;
+  // Assert starting of image has the same Z than first midpoints
+  while (midpoints[slice][2] != m_Working_Support->GetOrigin()[2]) {
+    slice++;
+    if ((uint)slice >= midpoints.size()) {
+      clitkExceptionMacro("Bug while searching for first midpoint to use");
+    }
+  }
+
+  // Start loop
+  while (!iter.IsAtEnd()) {
+    ImageIndexType index;
+    m_Working_Support->TransformPhysicalPointToIndex(midpoints[slice], index);
+    while (!iter.IsAtEndOfSlice()) {
+      // convert into index
+      while (!iter.IsAtEndOfLine()) {
+	// if indexcourant <index this is Right part. Left otherwise
+	if (iter.GetIndex()[0]<index[0]) {
+	  iterR.Set(iter.Get());
+	  iterL.Set(GetBackgroundValue());
+	}
+	else {
+	  iterL.Set(iter.Get());
+	  iterR.Set(GetBackgroundValue());
+	}
+        ++iter;
+	++iterL;
+	++iterR;
+      }
+      iter.NextLine();
+      iterL.NextLine();
+      iterR.NextLine();
+
+    }
+    iter.NextSlice();
+    iterL.NextSlice();
+    iterR.NextSlice();
+    ++slice;
+  }
+
+  // Auto crop !
+  m_LeftSupport = clitk::AutoCrop<MaskImageType>(m_LeftSupport, GetBackgroundValue());               
+  m_RightSupport = clitk::AutoCrop<MaskImageType>(m_RightSupport, GetBackgroundValue());               
+  writeImage<MaskImageType>(m_LeftSupport, "lsac.mhd");
+  writeImage<MaskImageType>(m_RightSupport, "rsac.mhd");
+
+  StopCurrentStep<MaskImageType>(m_RightSupport);
+
+  // -------------------------------------------------------------------------
+  // Constraint at right from the SVC. Must be after MidTrachea limits
+  // (if not, bronchus can be split, midposition is false)
+  StartNewStep("[Station 4R] R limits with SVC ");
+  MaskImagePointer svc = GetAFDB()->template GetImage<MaskImageType>("SVC");
+  typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
+  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
+  /*  relPosFilter->SetVerboseStep(false);
+  relPosFilter->SetInput(m_RightSupport);  // only right here ...
+  relPosFilter->SetInputObject(svc); 
+  relPosFilter->SetOrientationType(RelPosFilterType::RightTo);
+  relPosFilter->SetIntermediateSpacing(2);
+  relPosFilter->SetFuzzyThreshold(0.3);
+  relPosFilter->Update();
+  m_RightSupport = relPosFilter->GetOutput();
+  */
 
+  /*
+    ==> TODO RIGHT TO MEDIAL TO SVC : 
+    get centroid, cut in X direction to get medial
 
+    ==> REDO OPERATOR to find extrma points ?
+    centroid or skeleton ? 
+
+  */
+
+  relPosFilter = RelPosFilterType::New();
+  relPosFilter->VerboseStepFlagOff();
+  relPosFilter->SetInput(m_RightSupport);  // only right here ...
+  relPosFilter->SetInputObject(svc); 
+  relPosFilter->AddOrientationType(RelPosFilterType::AntTo);
+  relPosFilter->InverseOrientationFlagOn();
+  relPosFilter->SetIntermediateSpacing(2); // this is important to put it low
+  relPosFilter->SetFuzzyThreshold(0.6);
+  relPosFilter->Update();
+  m_RightSupport = relPosFilter->GetOutput();
+
+  // AutoCrop
+  m_RightSupport = clitk::AutoCrop<MaskImageType>(m_RightSupport, GetBackgroundValue());               
 
 }
 //--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class TImageType>
+void 
+clitk::ExtractLymphStationsFilter<TImageType>::
+ExtractStation_4RL_AP_Limits() 
+{
+
+  /*
+    post of S4R
+    - sup part -> cut line post wall trachea
+    - inf part (where ???) -> remains anterior to the right main stem bronchus
+    ==> bo
+  */
+
+  // Post (not Ant) to Aorta
+  MaskImagePointer aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
+  // Crop according to current support
+  aorta = clitk::ResizeImageLike<MaskImageType>(aorta, m_RightSupport, GetBackgroundValue());
+  typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
+  typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
+  relPosFilter = RelPosFilterType::New();
+  relPosFilter->VerboseStepFlagOff();
+  relPosFilter->SetInput(m_RightSupport);  // only right here ...
+  relPosFilter->SetInputObject(aorta); 
+  relPosFilter->AddOrientationType(RelPosFilterType::PostTo);
+  //  relPosFilter->NotFlagOn();
+  relPosFilter->SetIntermediateSpacing(2); // this is important to put it low
+  relPosFilter->SetFuzzyThreshold(0.6);
+  relPosFilter->Update();
+  m_RightSupport = relPosFilter->GetOutput();
+
+  m_ListOfStations["4R"] = m_RightSupport;
+  StopCurrentStep<MaskImageType>(m_RightSupport);
+
+  // POST -> horizontal lines bronchus --> points dectection in S7 to store.
+}
+//--------------------------------------------------------------------
+