segmentation/clitkExtractLymphStation_7.txx

   1
   2 //--------------------------------------------------------------------
   3 template <class TImageType>
   4 void
   5 clitk::ExtractLymphStationsFilter<TImageType>::
   6 ExtractStation_7_SI_Limits()
   7 {
   8   // Get Inputs
   9   MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
  10   double m_CarinaZ = GetAFDB()->GetPoint3D("Carina", 2);
  11   DD(m_CarinaZ);
  12   double m_OriginOfRightMiddleLobeBronchusZ = GetAFDB()->GetPoint3D("OriginOfRightMiddleLobeBronchus", 2);
  13   DD(m_OriginOfRightMiddleLobeBronchusZ);
  14
  15   /* Crop support :
  16      Superior limit = carina
  17      Inferior limit = origin right middle lobe bronchus */
  18   StartNewStep("[Station7] Inf/Sup mediastinum limits with carina/bronchus");
  19   m_Working_Support =
  20     clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
  21                                                 m_OriginOfRightMiddleLobeBronchusZ,
  22                                                 m_CarinaZ, true,
  23                                                 GetBackgroundValue());
  24   /* Crop trachea
  25      Superior limit = carina
  26      Inferior limit = origin right middle lobe bronchus*/
  27   m_working_trachea =
  28     clitk::CropImageAlongOneAxis<MaskImageType>(Trachea, 2,
  29                                                 m_OriginOfRightMiddleLobeBronchusZ,
  30                                                 m_CarinaZ, true,
  31                                                 GetBackgroundValue());
  32
  33   StopCurrentStep<MaskImageType>(m_Working_Support);
  34   m_ListOfStations["7"] = m_Working_Support;
  35 }
  36 //--------------------------------------------------------------------
  37
  38
  39 //--------------------------------------------------------------------
  40 template <class TImageType>
  41 void
  42 clitk::ExtractLymphStationsFilter<TImageType>::
  43 ExtractStation_7_RL_Limits()
  44 {
  45   // ----------------------------------------------------------------
  46   // Separate trachea in two CCL
  47   StartNewStep("[Station7] Separate trachea under carina");
  48
  49   // Labelize and consider two main labels
  50   m_working_trachea = Labelize<MaskImageType>(m_working_trachea, 0, true, 1);
  51
  52   // Carina position must at the first slice that separate the two
  53   // main bronchus (not superiorly) Check that upper slice is composed
  54   // of at least two labels
  55   typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
  56   SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
  57   iter.SetFirstDirection(0);
  58   iter.SetSecondDirection(1);
  59   iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
  60   int maxLabel=0;
  61   while (!iter.IsAtReverseEndOfSlice()) {
  62     while (!iter.IsAtReverseEndOfLine()) {
  63       if (iter.Get() > maxLabel) maxLabel = iter.Get();
  64       --iter;
  65     }
  66     iter.PreviousLine();
  67   }
  68   if (maxLabel < 2) {
  69     clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
  70   }
  71
  72   // Compute 3D centroids of both parts to identify the left from the
  73   // right bronchus
  74   std::vector<ImagePointType> c;
  75   clitk::ComputeCentroids<MaskImageType>(m_working_trachea, GetBackgroundValue(), c);
  76   ImagePointType C1 = c[1];
  77   ImagePointType C2 = c[2];
  78
  79   ImagePixelType leftLabel;
  80   ImagePixelType rightLabel;
  81   if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
  82   else { leftLabel = 2; rightLabel = 1; }
  83
  84   StopCurrentStep<MaskImageType>(m_working_trachea);
  85
  86   //-----------------------------------------------------
  87   // Select LeftLabel (set one label to Backgroundvalue)
  88   m_LeftBronchus =
  89     SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
  90                                                 rightLabel, GetBackgroundValue(), false);
  91   m_RightBronchus  =
  92     SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
  93                                                 leftLabel, GetBackgroundValue(), false);
  94
  95   StartNewStep("[Station7] Limits with bronchus (slice by slice) : RightTo left bronchus");
  96   m_Working_Support =
  97     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
  98                                                        m_LeftBronchus, 2,
  99                                                        GetFuzzyThreshold(), "RightTo",
 100                                                        true, 4);
 101
 102   StartNewStep("[Station7] Limits with bronchus (slice by slice) : LeftTo right bronchus");
 103   m_Working_Support =
 104     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 105                                                        m_RightBronchus,
 106                                                        2, GetFuzzyThreshold(), "LeftTo",
 107                                                        true, 4);
 108
 109   StartNewStep("[Station7] Limits with bronchus (slice by slice) : not AntTo left bronchus");
 110   m_Working_Support =
 111     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 112                                                        m_LeftBronchus,
 113                                                        2, GetFuzzyThreshold(), "AntTo",
 114                                                        true, 4, true); // NOT
 115
 116   StartNewStep("[Station7] Limits with bronchus (slice by slice) : not AntTo right bronchus");
 117   m_Working_Support =
 118     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 119                                                        m_RightBronchus,
 120                                                        2, GetFuzzyThreshold(), "AntTo",
 121                                                        true, 4, true);
 122
 123   StartNewStep("[Station7] Limits with bronchus (slice by slice) : not PostTo left bronchus");
 124   m_Working_Support =
 125     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 126                                                        m_LeftBronchus,
 127                                                        2, GetFuzzyThreshold(), "PostTo",
 128                                                        true, 4, true);
 129
 130   StartNewStep("[Station7] Limits with bronchus (slice by slice) : not PostTo right bronchus");
 131   m_Working_Support =
 132     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 133                                                        m_RightBronchus,
 134                                                        2, GetFuzzyThreshold(), "PostTo",
 135                                                        true, 4, true);
 136   m_Station7 = m_Working_Support;
 137   StopCurrentStep<MaskImageType>(m_Station7);
 138 }
 139 //--------------------------------------------------------------------
 140
 141
 142 //--------------------------------------------------------------------
 143 template <class TImageType>
 144 void
 145 clitk::ExtractLymphStationsFilter<TImageType>::
 146 ExtractStation_7_Posterior_Limits()
 147 {
 148   StartNewStep("[Station7] Posterior limits -> must be AntTo post wall of the bronchi");
 149
 150   // Search for points that are the most left/post/ant and most
 151   // right/post/ant of the left and right bronchus
 152
 153   // extract, loop slices, label/keep, find extrema x 3
 154   FindExtremaPointsInBronchus(m_LeftBronchus, 0, 15,
 155                               m_RightMostInLeftBronchus,
 156                               m_AntMostInLeftBronchus,
 157                               m_PostMostInLeftBronchus);
 158   FindExtremaPointsInBronchus(m_RightBronchus, 1, 15,
 159                               m_LeftMostInRightBronchus,
 160                               m_AntMostInRightBronchus,
 161                               m_PostMostInRightBronchus);
 162   // DEBUG
 163   std::ofstream osrl; openFileForWriting(osrl, "osrl.txt"); osrl << "LANDMARKS1" << std::endl;
 164   std::ofstream osal; openFileForWriting(osal, "osal.txt"); osal << "LANDMARKS1" << std::endl;
 165   std::ofstream ospl; openFileForWriting(ospl, "ospl.txt"); ospl << "LANDMARKS1" << std::endl;
 166   std::ofstream osrr; openFileForWriting(osrr, "osrr.txt"); osrr << "LANDMARKS1" << std::endl;
 167   std::ofstream osar; openFileForWriting(osar, "osar.txt"); osar << "LANDMARKS1" << std::endl;
 168   std::ofstream ospr; openFileForWriting(ospr, "ospr.txt"); ospr << "LANDMARKS1" << std::endl;
 169
 170   for(uint i=0; i<m_RightMostInLeftBronchus.size(); i++) {
 171     osrl << i << " " << m_RightMostInLeftBronchus[i][0] << " " << m_RightMostInLeftBronchus[i][1]
 172          << " " << m_RightMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
 173     osal << i << " " << m_AntMostInLeftBronchus[i][0] << " " << m_AntMostInLeftBronchus[i][1]
 174          << " " << m_AntMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
 175     ospl << i << " " << m_PostMostInLeftBronchus[i][0] << " " << m_PostMostInLeftBronchus[i][1]
 176          << " " << m_PostMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
 177
 178     osrr << i << " " << m_LeftMostInRightBronchus[i][0] << " " << m_LeftMostInRightBronchus[i][1]
 179          << " " << m_LeftMostInRightBronchus[i][2] << " 0 0 " << std::endl;
 180     osar << i << " " << m_AntMostInRightBronchus[i][0] << " " << m_AntMostInRightBronchus[i][1]
 181          << " " << m_AntMostInRightBronchus[i][2] << " 0 0 " << std::endl;
 182     ospr << i << " " << m_PostMostInRightBronchus[i][0] << " " << m_PostMostInRightBronchus[i][1]
 183          << " " << m_PostMostInRightBronchus[i][2] << " 0 0 " << std::endl;
 184   }
 185   osrl.close();
 186   osal.close();
 187   ospl.close();
 188
 189   // Now uses these points to limit, slice by slice
 190   // http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
 191   /*
 192     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
 193     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
 194     This will equal zero if the point C is on the line formed by
 195     points A and B, and will have a different sign depending on the
 196     side. Which side this is depends on the orientation of your (x,y)
 197     coordinates, but you can plug test values for A,B and C into this
 198     formula to determine whether negative values are to the left or to
 199     the right.
 200     => to accelerate, start with formula, when change sign -> stop and fill
 201   */
 202   typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
 203   SliceIteratorType iter = SliceIteratorType(m_Working_Support,
 204                                              m_Working_Support->GetLargestPossibleRegion());
 205   iter.SetFirstDirection(0);
 206   iter.SetSecondDirection(1);
 207   iter.GoToBegin();
 208   int i=0;
 209   MaskImageType::PointType A;
 210   MaskImageType::PointType B;
 211   MaskImageType::PointType C;
 212   while (!iter.IsAtEnd()) {
 213     A = m_PostMostInLeftBronchus[i];
 214     B = m_PostMostInRightBronchus[i];
 215     C = A;
 216     C[1] -= 10; // I know I must keep this point
 217     double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 218     bool isPositive = s<0;
 219     while (!iter.IsAtEndOfSlice()) {
 220       while (!iter.IsAtEndOfLine()) {
 221         // Very slow, I know ... but image should be very small
 222         m_Working_Support->TransformIndexToPhysicalPoint(iter.GetIndex(), C);
 223         double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 224         if (s == 0) iter.Set(2);
 225         if (isPositive) {
 226           if (s > 0) iter.Set(GetBackgroundValue());
 227         }
 228         else {
 229           if (s < 0) iter.Set(GetBackgroundValue());
 230         }
 231         ++iter;
 232       }
 233       iter.NextLine();
 234     }
 235     iter.NextSlice();
 236     ++i;
 237   }
 238
 239   //-----------------------------------------------------
 240   // StartNewStep("[Station7] Anterior limits");
 241
 242
 243   // MISSING FROM NOW
 244
 245   // Station 4R, Station 4L, the right pulmonary artery, and/or the
 246   // left superior pulmonary vein
 247
 248
 249   //-----------------------------------------------------
 250   //-----------------------------------------------------
 251   // ALSO SUBSTRACT ARTERY/VEIN (initially in the support)
 252
 253
 254   // Set output
 255   m_Station7 = m_Working_Support;
 256 }
 257 //--------------------------------------------------------------------
 258
 259