segmentation/clitkExtractLymphStation_7.txx

   1 //--------------------------------------------------------------------
   2 template <class TImageType>
   3 void
   4 clitk::ExtractLymphStationsFilter<TImageType>::
   5 ExtractStation_7_SI_Limits()
   6 {
   7   // Local variables
   8   double m_CarinaZPositionInMM;
   9   double m_MiddleLobeBronchusZPositionInMM;
  10
  11   // Get Inputs
  12   m_Trachea = GetAFDB()->template GetImage <MaskImageType>("trachea");
  13   m_CarinaZPositionInMM = GetAFDB()->GetPoint3D("carina", 2);
  14   DD(m_CarinaZPositionInMM);
  15   m_MiddleLobeBronchusZPositionInMM = GetAFDB()->GetPoint3D("rightMiddleLobeBronchus", 2);
  16   DD(m_MiddleLobeBronchusZPositionInMM);
  17
  18   /* Crop support :
  19        Superior limit = carina
  20        Inferior limit = origin right middle lobe bronchus */
  21   StartNewStep("Inf/Sup mediastinum limits with carina/bronchus");
  22   m_Working_Support =
  23     clitk::CropImageAlongOneAxis<MaskImageType>(m_Support, 2,
  24                                                 m_MiddleLobeBronchusZPositionInMM,
  25                                                 m_CarinaZPositionInMM, true,
  26                                                 GetBackgroundValue());
  27   StopCurrentStep<MaskImageType>(m_Working_Support);
  28
  29   /* Crop trachea
  30        Superior limit = carina
  31        Inferior limit = origin right middle lobe bronchus*/
  32   StartNewStep("Inf/Sup trachea limits with carina/bronchus");
  33   m_working_trachea =
  34     clitk::CropImageAlongOneAxis<MaskImageType>(m_Trachea, 2,
  35                                                 m_MiddleLobeBronchusZPositionInMM,
  36                                                 m_CarinaZPositionInMM, true,
  37                                                 GetBackgroundValue());
  38   StopCurrentStep<MaskImageType>(m_working_trachea);
  39
  40   m_Station7 = m_Working_Support;
  41 }
  42 //--------------------------------------------------------------------
  43
  44 //--------------------------------------------------------------------
  45 template <class TImageType>
  46 void
  47 clitk::ExtractLymphStationsFilter<TImageType>::
  48 ExtractStation_7_RL_Limits()
  49 {
  50   // ----------------------------------------------------------------
  51   // Separate trachea in two CCL
  52   StartNewStep("Separate trachea under carina");
  53
  54   // Labelize and consider two main labels
  55   m_working_trachea = Labelize<MaskImageType>(m_working_trachea, 0, true, 1);
  56
  57   // Carina position must at the first slice that separate the two main bronchus (not superiorly)
  58   // Check that upper slice is composed of at least two labels
  59   typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
  60   SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion());
  61   iter.SetFirstDirection(0);
  62   iter.SetSecondDirection(1);
  63   iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
  64   int maxLabel=0;
  65   while (!iter.IsAtReverseEndOfSlice()) {
  66     while (!iter.IsAtReverseEndOfLine()) {
  67       //  DD(iter.GetIndex());
  68       if (iter.Get() > maxLabel) maxLabel = iter.Get();
  69       --iter;
  70     }
  71     iter.PreviousLine();
  72   }
  73   if (maxLabel < 2) {
  74     clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort");
  75   }
  76
  77   // Compute centroid of both parts to identify the left from the right bronchus
  78   typedef long LabelType;
  79   static const unsigned int Dim = ImageType::ImageDimension;
  80   typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
  81   typedef itk::LabelMap< LabelObjectType > LabelMapType;
  82   typedef itk::LabelImageToLabelMapFilter<MaskImageType, LabelMapType> ImageToMapFilterType;
  83   typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
  84   typedef itk::ShapeLabelMapFilter<LabelMapType, MaskImageType> ShapeFilterType;
  85   typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
  86   imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
  87   imageToLabelFilter->SetInput(m_working_trachea);
  88   statFilter->SetInput(imageToLabelFilter->GetOutput());
  89   statFilter->Update();
  90   typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
  91
  92   ImagePointType C1 = labelMap->GetLabelObject(1)->GetCentroid();
  93   ImagePointType C2 = labelMap->GetLabelObject(2)->GetCentroid();
  94
  95   ImagePixelType leftLabel;
  96   ImagePixelType rightLabel;
  97   if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; }
  98   else { leftLabel = 2; rightLabel = 1; }
  99   DD(leftLabel);
 100   DD(rightLabel);
 101
 102   StopCurrentStep<MaskImageType>(m_working_trachea);
 103
 104   //-----------------------------------------------------
 105   StartNewStep("Left limits with bronchus (slice by slice)");
 106   // Select LeftLabel (set one label to Backgroundvalue)
 107   m_LeftBronchus =
 108     SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
 109                                                 rightLabel, GetBackgroundValue());
 110   m_RightBronchus  =
 111     SetBackground<MaskImageType, MaskImageType>(m_working_trachea, m_working_trachea,
 112                                                 leftLabel, GetBackgroundValue());
 113   writeImage<MaskImageType>(m_LeftBronchus, "left.mhd");
 114   writeImage<MaskImageType>(m_RightBronchus, "right.mhd");
 115
 116   m_Working_Support =
 117     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 118                                                        m_LeftBronchus, 2,
 119                                                        GetFuzzyThreshold(), "RightTo",
 120                                                        true, 4);
 121   m_Working_Support =
 122     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 123                                                        m_RightBronchus,
 124                                                        2, GetFuzzyThreshold(), "LeftTo",
 125                                                        true, 4);
 126   m_Working_Support =
 127     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 128                                                        m_LeftBronchus,
 129                                                        2, GetFuzzyThreshold(), "AntTo",
 130                                                        true, 4, true); // NOT
 131   m_Working_Support =
 132     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 133                                                        m_RightBronchus,
 134                                                        2, GetFuzzyThreshold(), "AntTo",
 135                                                        true, 4, true);
 136   m_Working_Support =
 137     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 138                                                        m_LeftBronchus,
 139                                                        2, GetFuzzyThreshold(), "PostTo",
 140                                                        true, 4, true);
 141   m_Working_Support =
 142     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support,
 143                                                        m_RightBronchus,
 144                                                        2, GetFuzzyThreshold(), "PostTo",
 145                                                        true, 4, true);
 146   m_Station7 = m_Working_Support;
 147   StopCurrentStep<MaskImageType>(m_Station7);
 148 }
 149 //--------------------------------------------------------------------
 150
 151
 152 //--------------------------------------------------------------------
 153 template <class TImageType>
 154 void
 155 clitk::ExtractLymphStationsFilter<TImageType>::
 156 ExtractStation_7_Posterior_Limits()
 157 {
 158   StartNewStep("Posterior limits");
 159
 160   // Left Bronchus slices
 161   typedef clitk::ExtractSliceFilter<MaskImageType> ExtractSliceFilterType;
 162   typedef typename ExtractSliceFilterType::SliceType SliceType;
 163   typename ExtractSliceFilterType::Pointer
 164     extractSliceFilter = ExtractSliceFilterType::New();
 165   extractSliceFilter->SetInput(m_LeftBronchus);
 166   extractSliceFilter->SetDirection(2);
 167   extractSliceFilter->Update();
 168   std::vector<typename SliceType::Pointer> leftBronchusSlices;
 169   extractSliceFilter->GetOutputSlices(leftBronchusSlices);
 170
 171   // Right Bronchus slices
 172   extractSliceFilter = ExtractSliceFilterType::New();
 173   extractSliceFilter->SetInput(m_RightBronchus);
 174   extractSliceFilter->SetDirection(2);
 175   extractSliceFilter->Update();
 176   std::vector<typename SliceType::Pointer> rightBronchusSlices ;
 177   extractSliceFilter->GetOutputSlices(rightBronchusSlices);
 178
 179   assert(leftBronchusSlices.size() == rightBronchusSlices.size());
 180
 181   std::vector<MaskImageType::PointType> leftPoints;
 182   std::vector<MaskImageType::PointType> rightPoints;
 183   for(uint i=0; i<leftBronchusSlices.size(); i++) {
 184     // Keep main CCL
 185     leftBronchusSlices[i] = Labelize<SliceType>(leftBronchusSlices[i], 0, true, 10);
 186     leftBronchusSlices[i] = KeepLabels<SliceType>(leftBronchusSlices[i],
 187                                                   GetBackgroundValue(),
 188                                                   GetForegroundValue(), 1, 1, true);
 189     rightBronchusSlices[i] = Labelize<SliceType>(rightBronchusSlices[i], 0, true, 10);
 190     rightBronchusSlices[i] = KeepLabels<SliceType>(rightBronchusSlices[i],
 191                                                    GetBackgroundValue(),
 192                                                    GetForegroundValue(), 1, 1, true);
 193     double distance_max_from_center_point = 15;
 194
 195     // ------- Find point in left Bronchus -------
 196     // find right most point in left  = rightMost
 197     SliceType::PointType a;
 198     SliceType::PointType rightMost =
 199       clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i],
 200                                                           GetBackgroundValue(),
 201                                                           0, false, a, 0);
 202     // find post most point in left, not far away from rightMost
 203     SliceType::PointType p =
 204       clitk::FindExtremaPointInAGivenDirection<SliceType>(leftBronchusSlices[i],
 205                                                           GetBackgroundValue(),
 206                                                           1, false, rightMost,
 207                                                           distance_max_from_center_point);
 208     MaskImageType::PointType pp;
 209     pp[0] = p[0]; pp[1] = p[1];
 210     pp[2] = i*m_LeftBronchus->GetSpacing()[2] + m_LeftBronchus->GetOrigin()[2];
 211     leftPoints.push_back(pp);
 212
 213     // ------- Find point in right Bronchus -------
 214     // find left most point in right  = leftMost
 215     SliceType::PointType leftMost =
 216       clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i],
 217                                                           GetBackgroundValue(),
 218                                                           0, true, a, 0);
 219     // find post most point in left, not far away from leftMost
 220     p = clitk::FindExtremaPointInAGivenDirection<SliceType>(rightBronchusSlices[i],
 221                                                             GetBackgroundValue(),
 222                                                             1, false, leftMost,
 223                                                             distance_max_from_center_point);
 224     pp[0] = p[0]; pp[1] = p[1];
 225     pp[2] = i*m_RightBronchus->GetSpacing()[2] + m_RightBronchus->GetOrigin()[2];
 226     rightPoints.push_back(pp);
 227   }
 228
 229   // DEBUG
 230   std::ofstream osl;
 231   openFileForWriting(osl, "left.txt");
 232   osl << "LANDMARKS1" << std::endl;
 233   std::ofstream osr;
 234   openFileForWriting(osr, "right.txt");
 235   osr << "LANDMARKS1" << std::endl;
 236   for(uint i=0; i<leftBronchusSlices.size(); i++) {
 237     osl << i << " " << leftPoints[i][0] << " " << leftPoints[i][1]
 238         << " " << leftPoints[i][2] << " 0 0 " << std::endl;
 239     osr << i << " " << rightPoints[i][0] << " " << rightPoints[i][1]
 240         << " " << rightPoints[i][2] << " 0 0 " << std::endl;
 241   }
 242   osl.close();
 243   osr.close();
 244
 245   // Now uses these points to limit, slice by slice
 246   // http://www.gamedev.net/community/forums/topic.asp?topic_id=542870
 247   /*
 248     Assuming the points are (Ax,Ay) (Bx,By) and (Cx,Cy), you need to compute:
 249     (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax)
 250     This will equal zero if the point C is on the line formed by
 251     points A and B, and will have a different sign depending on the
 252     side. Which side this is depends on the orientation of your (x,y)
 253     coordinates, but you can plug test values for A,B and C into this
 254     formula to determine whether negative values are to the left or to
 255     the right.
 256     => to accelerate, start with formula, when change sign -> stop and fill
 257   */
 258   typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
 259   SliceIteratorType iter = SliceIteratorType(m_Working_Support,
 260                                              m_Working_Support->GetLargestPossibleRegion());
 261   iter.SetFirstDirection(0);
 262   iter.SetSecondDirection(1);
 263   iter.GoToBegin();
 264   int i=0;
 265   MaskImageType::PointType A;
 266   MaskImageType::PointType B;
 267   MaskImageType::PointType C;
 268   while (!iter.IsAtEnd()) {
 269     A = leftPoints[i];
 270     B = rightPoints[i];
 271     C = A;
 272     C[1] -= 10; // I know I must keep this point
 273     double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 274     bool isPositive = s<0;
 275     while (!iter.IsAtEndOfSlice()) {
 276       while (!iter.IsAtEndOfLine()) {
 277         // Very slow, I know ... but image should be very small
 278         m_Working_Support->TransformIndexToPhysicalPoint(iter.GetIndex(), C);
 279         double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]);
 280         if (s == 0) iter.Set(2);
 281         if (isPositive) {
 282           if (s > 0) iter.Set(GetBackgroundValue());
 283         }
 284         else {
 285           if (s < 0) iter.Set(GetBackgroundValue());
 286         }
 287         ++iter;
 288       }
 289       iter.NextLine();
 290     }
 291     iter.NextSlice();
 292     ++i;
 293   }
 294
 295   //-----------------------------------------------------
 296   // StartNewStep("Anterior limits");
 297
 298
 299   // MISSING FROM NOW
 300
 301   // Station 4R, Station 4L, the right pulmonary artery, and/or the
 302   // left superior pulmonary vein
 303
 304
 305   //-----------------------------------------------------
 306   //-----------------------------------------------------
 307   // ALSO SUBSTRACT ARTERY/VEIN (initially in the support)
 308
 309
 310   // Set output
 311   m_Station7 = m_Working_Support;
 312 }
 313 //--------------------------------------------------------------------