]> Creatis software - clitk.git/blobdiff - segmentation/clitkExtractLymphStation_7.txx
correct options scheme
[clitk.git] / segmentation / clitkExtractLymphStation_7.txx
index e1cbe1f7750f47ae1319c5f8d8ef5c28c1580386..278efc8e076e3eb9f6e4410a92fb281863286e37 100644 (file)
@@ -1,46 +1,41 @@
+
 //--------------------------------------------------------------------
 template <class TImageType>
 void 
 clitk::ExtractLymphStationsFilter<TImageType>::
 ExtractStation_7_SI_Limits() 
 {
-  // Local variables
-  double m_CarinaZPositionInMM;
-  double m_MiddleLobeBronchusZPositionInMM;
-    
   // Get Inputs
-  m_Trachea = GetAFDB()->template GetImage <MaskImageType>("trachea");  
-  m_CarinaZPositionInMM = GetAFDB()->GetPoint3D("carina", 2);
-  DD(m_CarinaZPositionInMM);
-  m_MiddleLobeBronchusZPositionInMM = GetAFDB()->GetPoint3D("rightMiddleLobeBronchus", 2);
-  DD(m_MiddleLobeBronchusZPositionInMM);
+  MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");  
+  double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2);
+  DD(m_CarinaZ);
+  double m_OriginOfRightMiddleLobeBronchusZ = GetAFDB()->GetPoint3D("OriginOfRightMiddleLobeBronchus", 2);
+  DD(m_OriginOfRightMiddleLobeBronchusZ);
 
   /* Crop support :
-       Superior limit = carina
-       Inferior limit = origin right middle lobe bronchus */
-  StartNewStep("Inf/Sup mediastinum limits with carina/bronchus");
+     Superior limit = carina
+     Inferior limit = origin right middle lobe bronchus */
+  StartNewStep("[Station7] Inf/Sup mediastinum limits with carina/bronchus");
   m_Working_Support = 
-    clitk::CropImageAlongOneAxis<MaskImageType>(m_Support, 2, 
-                                                m_MiddleLobeBronchusZPositionInMM
-                                                m_CarinaZPositionInMM, true,
+    clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2, 
+                                                m_OriginOfRightMiddleLobeBronchusZ
+                                                m_CarinaZ, true,
                                                 GetBackgroundValue());
-  StopCurrentStep<MaskImageType>(m_Working_Support);
-
   /* Crop trachea
-       Superior limit = carina
-       Inferior limit = origin right middle lobe bronchus*/
-  StartNewStep("Inf/Sup trachea limits with carina/bronchus");
+     Superior limit = carina
+     Inferior limit = origin right middle lobe bronchus*/
   m_working_trachea = 
-    clitk::CropImageAlongOneAxis<MaskImageType>(m_Trachea, 2, 
-                                                m_MiddleLobeBronchusZPositionInMM
-                                                m_CarinaZPositionInMM, true,
+    clitk::CropImageAlongOneAxis<MaskImageType>(Trachea, 2, 
+                                                m_OriginOfRightMiddleLobeBronchusZ
+                                                m_CarinaZ, true,
                                                 GetBackgroundValue());
-  StopCurrentStep<MaskImageType>(m_working_trachea);
 
-  m_Station7 = m_Working_Support;
+  StopCurrentStep<MaskImageType>(m_Working_Support);
+  m_ListOfStations["7"] = m_Working_Support;
 }
 //--------------------------------------------------------------------
 
+
 //--------------------------------------------------------------------
 template <class TImageType>
 void 
@@ -49,13 +44,14 @@ ExtractStation_7_RL_Limits()
 {
   // ----------------------------------------------------------------
   // Separate trachea in two CCL
-  StartNewStep("Separate trachea under carina");
+  StartNewStep("[Station7] Separate trachea under carina");
 
   // Labelize and consider two main labels
   m_working_trachea = Labelize<MaskImageType>(m_working_trachea, 0, true, 1);
 
-  // Carina position must at the first slice that separate the two main bronchus (not superiorly) 
-  // Check that upper slice is composed of at least two labels
+  // Carina position must at the first slice that separate the two
+  // main bronchus (not superiorly) Check that upper slice is composed
+  // of at least two labels
   typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
   SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
   iter.SetFirstDirection(0);
@@ -64,7 +60,6 @@ ExtractStation_7_RL_Limits()
   int maxLabel=0;
   while (!iter.IsAtReverseEndOfSlice()) {
     while (!iter.IsAtReverseEndOfLine()) {    
-      //  DD(iter.GetIndex());
       if (iter.Get() > maxLabel) maxLabel = iter.Get();
       --iter;
     }
@@ -74,70 +69,65 @@ ExtractStation_7_RL_Limits()
     clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
   }
 
-  // Compute centroid of both parts to identify the left from the right bronchus
-  typedef long LabelType;
-  static const unsigned int Dim = ImageType::ImageDimension;
-  typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
-  typedef itk::LabelMap< LabelObjectType > LabelMapType;
-  typedef itk::LabelImageToLabelMapFilter<MaskImageType, LabelMapType> ImageToMapFilterType;
-  typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); 
-  typedef itk::ShapeLabelMapFilter<LabelMapType, MaskImageType> ShapeFilterType; 
-  typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
-  imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
-  imageToLabelFilter->SetInput(m_working_trachea);
-  statFilter->SetInput(imageToLabelFilter->GetOutput());
-  statFilter->Update();
-  typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
-
-  ImagePointType C1 = labelMap->GetLabelObject(1)->GetCentroid();
-  ImagePointType C2 = labelMap->GetLabelObject(2)->GetCentroid();
+  // Compute 3D centroids of both parts to identify the left from the
+  // right bronchus
+  std::vector<ImagePointType> c;
+  clitk::ComputeCentroids<MaskImageType>(m_working_trachea, GetBackgroundValue(), c);
+  ImagePointType C1 = c[1];
+  ImagePointType C2 = c[2];
 
   ImagePixelType leftLabel;
   ImagePixelType rightLabel;  
   if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
   else { leftLabel = 2; rightLabel = 1; }
-  DD(leftLabel);
-  DD(rightLabel);
 
   StopCurrentStep<MaskImageType>(m_working_trachea);
 
   //-----------------------------------------------------
-  StartNewStep("Left limits with bronchus (slice by slice)");  
   // Select LeftLabel (set one label to Backgroundvalue)
   m_LeftBronchus = 
     SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea, 
-                                                rightLabel, GetBackgroundValue());
+                                                rightLabel, GetBackgroundValue(), false);
   m_RightBronchus  = 
     SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea, 
-                                                leftLabel, GetBackgroundValue());
-  writeImage<MaskImageType>(m_LeftBronchus, "left.mhd");
-  writeImage<MaskImageType>(m_RightBronchus, "right.mhd");
+                                                leftLabel, GetBackgroundValue(), false);
 
+  StartNewStep("[Station7] Limits with bronchus (slice by slice) : RightTo left bronchus");  
   m_Working_Support = 
     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, 
                                                       m_LeftBronchus, 2, 
                                                        GetFuzzyThreshold(), "RightTo", 
                                                        true, 4);
+
+  StartNewStep("[Station7] Limits with bronchus (slice by slice) : LeftTo right bronchus");  
   m_Working_Support = 
     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, 
                                                       m_RightBronchus, 
                                                       2, GetFuzzyThreshold(), "LeftTo", 
                                                        true, 4);
+
+  StartNewStep("[Station7] Limits with bronchus (slice by slice) : not AntTo left bronchus");  
   m_Working_Support = 
     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, 
                                                       m_LeftBronchus, 
                                                       2, GetFuzzyThreshold(), "AntTo", 
                                                        true, 4, true); // NOT
+
+  StartNewStep("[Station7] Limits with bronchus (slice by slice) : not AntTo right bronchus");  
   m_Working_Support = 
     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, 
                                                       m_RightBronchus, 
                                                       2, GetFuzzyThreshold(), "AntTo", 
                                                        true, 4, true);
+
+  StartNewStep("[Station7] Limits with bronchus (slice by slice) : not PostTo left bronchus");  
   m_Working_Support = 
     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, 
                                                       m_LeftBronchus, 
                                                       2, GetFuzzyThreshold(), "PostTo", 
                                                        true, 4, true);
+
+  StartNewStep("[Station7] Limits with bronchus (slice by slice) : not PostTo right bronchus");  
   m_Working_Support = 
     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, 
                                                       m_RightBronchus, 
@@ -155,92 +145,46 @@ void
 clitk::ExtractLymphStationsFilter<TImageType>::
 ExtractStation_7_Posterior_Limits() 
 {
-  StartNewStep("Posterior limits");  
-
-  // Left Bronchus slices
-  typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
-  typedef typename ExtractSliceFilterType::SliceType SliceType;
-  typename ExtractSliceFilterType::Pointer 
-    extractSliceFilter = ExtractSliceFilterType::New();
-  extractSliceFilter->SetInput(m_LeftBronchus);
-  extractSliceFilter->SetDirection(2);
-  extractSliceFilter->Update();
-  std::vector<typename SliceType::Pointer> leftBronchusSlices;
-  extractSliceFilter->GetOutputSlices(leftBronchusSlices);
-  
-  // Right Bronchus slices
-  extractSliceFilter = ExtractSliceFilterType::New();
-  extractSliceFilter->SetInput(m_RightBronchus);
-  extractSliceFilter->SetDirection(2);
-  extractSliceFilter->Update();
-  std::vector<typename SliceType::Pointer> rightBronchusSlices ;
-  extractSliceFilter->GetOutputSlices(rightBronchusSlices);
-  
-  assert(leftBronchusSlices.size() == rightBronchusSlices.size());
-  
-  std::vector<MaskImageType::PointType> leftPoints;
-  std::vector<MaskImageType::PointType> rightPoints;
-  for(uint i=0; i<leftBronchusSlices.size(); i++) {
-    // Keep main CCL
-    leftBronchusSlices[i] = Labelize<SliceType>(leftBronchusSlices[i], 0, true, 10);
-    leftBronchusSlices[i] = KeepLabels<SliceType>(leftBronchusSlices[i], 
-                                                  GetBackgroundValue(), 
-                                                  GetForegroundValue(), 1, 1, true);
-    rightBronchusSlices[i] = Labelize<SliceType>(rightBronchusSlices[i], 0, true, 10);
-    rightBronchusSlices[i] = KeepLabels<SliceType>(rightBronchusSlices[i], 
-                                                   GetBackgroundValue(), 
-                                                   GetForegroundValue(), 1, 1, true);
-    double distance_max_from_center_point = 15;
-
-    // ------- Find point in left Bronchus ------- 
-    // find right most point in left  = rightMost
-    SliceType::PointType a;
-    SliceType::PointType rightMost = 
-      clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i], 
-                                                          GetBackgroundValue(), 
-                                                          0, false, a, 0);
-    // find post most point in left, not far away from rightMost
-    SliceType::PointType p = 
-      clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i], 
-                                                          GetBackgroundValue(), 
-                                                          1, false, rightMost, 
-                                                          distance_max_from_center_point);
-    MaskImageType::PointType pp;
-    pp[0] = p[0]; pp[1] = p[1];
-    pp[2] = i*m_LeftBronchus->GetSpacing()[2] + m_LeftBronchus->GetOrigin()[2];
-    leftPoints.push_back(pp);
-
-    // ------- Find point in right Bronchus ------- 
-    // find left most point in right  = leftMost
-    SliceType::PointType leftMost = 
-      clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i], 
-                                                          GetBackgroundValue(), 
-                                                          0, true, a, 0);
-    // find post most point in left, not far away from leftMost
-    p = clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i], 
-                                                            GetBackgroundValue(), 
-                                                            1, false, leftMost, 
-                                                            distance_max_from_center_point);
-    pp[0] = p[0]; pp[1] = p[1];
-    pp[2] = i*m_RightBronchus->GetSpacing()[2] + m_RightBronchus->GetOrigin()[2];
-    rightPoints.push_back(pp);
-  }
+  StartNewStep("[Station7] Posterior limits -> must be AntTo post wall of the bronchi");  
 
+  // Search for points that are the most left/post/ant and most
+  // right/post/ant of the left and right bronchus
+
+  // extract, loop slices, label/keep, find extrema x 3
+  FindExtremaPointsInBronchus(m_LeftBronchus, 0, 15,
+                             m_RightMostInLeftBronchus, 
+                             m_AntMostInLeftBronchus, 
+                             m_PostMostInLeftBronchus);
+  FindExtremaPointsInBronchus(m_RightBronchus, 1, 15,
+                             m_LeftMostInRightBronchus, 
+                             m_AntMostInRightBronchus, 
+                             m_PostMostInRightBronchus);
   // DEBUG
-  std::ofstream osl;
-  openFileForWriting(osl, "left.txt");
-  osl << "LANDMARKS1" << std::endl;
-  std::ofstream osr;
-  openFileForWriting(osr, "right.txt");
-  osr << "LANDMARKS1" << std::endl;
-  for(uint i=0; i<leftBronchusSlices.size(); i++) {
-    osl << i << " " << leftPoints[i][0] << " " << leftPoints[i][1] 
-        << " " << leftPoints[i][2] << " 0 0 " << std::endl;
-    osr << i << " " << rightPoints[i][0] << " " << rightPoints[i][1] 
-        << " " << rightPoints[i][2] << " 0 0 " << std::endl;
+  std::ofstream osrl; openFileForWriting(osrl, "osrl.txt"); osrl << "LANDMARKS1" << std::endl;
+  std::ofstream osal; openFileForWriting(osal, "osal.txt"); osal << "LANDMARKS1" << std::endl;
+  std::ofstream ospl; openFileForWriting(ospl, "ospl.txt"); ospl << "LANDMARKS1" << std::endl;
+  std::ofstream osrr; openFileForWriting(osrr, "osrr.txt"); osrr << "LANDMARKS1" << std::endl;
+  std::ofstream osar; openFileForWriting(osar, "osar.txt"); osar << "LANDMARKS1" << std::endl;
+  std::ofstream ospr; openFileForWriting(ospr, "ospr.txt"); ospr << "LANDMARKS1" << std::endl;
+
+  for(uint i=0; i<m_RightMostInLeftBronchus.size(); i++) {
+    osrl << i << " " << m_RightMostInLeftBronchus[i][0] << " " << m_RightMostInLeftBronchus[i][1] 
+        << " " << m_RightMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
+    osal << i << " " << m_AntMostInLeftBronchus[i][0] << " " << m_AntMostInLeftBronchus[i][1] 
+        << " " << m_AntMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
+    ospl << i << " " << m_PostMostInLeftBronchus[i][0] << " " << m_PostMostInLeftBronchus[i][1] 
+        << " " << m_PostMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
+
+    osrr << i << " " << m_LeftMostInRightBronchus[i][0] << " " << m_LeftMostInRightBronchus[i][1] 
+        << " " << m_LeftMostInRightBronchus[i][2] << " 0 0 " << std::endl;
+    osar << i << " " << m_AntMostInRightBronchus[i][0] << " " << m_AntMostInRightBronchus[i][1] 
+        << " " << m_AntMostInRightBronchus[i][2] << " 0 0 " << std::endl;
+    ospr << i << " " << m_PostMostInRightBronchus[i][0] << " " << m_PostMostInRightBronchus[i][1] 
+        << " " << m_PostMostInRightBronchus[i][2] << " 0 0 " << std::endl;
   }
-  osl.close();
-  osr.close();
+  osrl.close();
+  osal.close();
+  ospl.close();
 
   // Now uses these points to limit, slice by slice 
   // http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
@@ -266,8 +210,8 @@ ExtractStation_7_Posterior_Limits()
   MaskImageType::PointType B;
   MaskImageType::PointType C;
   while (!iter.IsAtEnd()) {
-    A = leftPoints[i];
-    B = rightPoints[i];
+    A = m_PostMostInLeftBronchus[i];
+    B = m_PostMostInRightBronchus[i];
     C = A;
     C[1] -= 10; // I know I must keep this point
     double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
@@ -293,7 +237,7 @@ ExtractStation_7_Posterior_Limits()
   }
 
   //-----------------------------------------------------
-  // StartNewStep("Anterior limits");  
+  // StartNewStep("[Station7] Anterior limits");  
  
 
   // MISSING FROM NOW 
@@ -311,3 +255,5 @@ ExtractStation_7_Posterior_Limits()
   m_Station7 = m_Working_Support;
 }
 //--------------------------------------------------------------------
+
+