segmentation/clitkExtractLymphStation_Supports.txx

   1
   2 #include <itkBinaryDilateImageFilter.h>
   3 #include <itkMirrorPadImageFilter.h>
   4
   5 //--------------------------------------------------------------------
   6 template <class ImageType>
   7 void
   8 clitk::ExtractLymphStationsFilter<ImageType>::
   9 ExtractStationSupports()
  10 {
  11   // Get initial Mediastinum
  12   m_Working_Support = m_Mediastinum = this->GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
  13
  14   // Remove some computed structures
  15   this->GetAFDB()->RemoveTag("CarinaZ");
  16   this->GetAFDB()->RemoveTag("ApexOfTheChestZ");
  17
  18   // Superior and inferior limits.
  19   Support_SI_Limit("inferior", "Sup_to_Carina", "inferior", "Carina", 0);
  20   Support_SI_Limit("superior", "Inf_to_Carina", "inferior", "Carina", m_Working_Support->GetSpacing()[2]);
  21
  22   // Initialise all others supports
  23   // m_ListOfSupports["S1R"] = m_ListOfSupports["Sup_to_Carina"];
  24   // m_ListOfSupports["S1L"] = m_ListOfSupports["Sup_to_Carina"];
  25   // m_ListOfSupports["S2R"] = m_ListOfSupports["Sup_to_Carina"];
  26   // m_ListOfSupports["S2L"] = m_ListOfSupports["Sup_to_Carina"];
  27   // m_ListOfSupports["S3A"] = m_ListOfSupports["Sup_to_Carina"];
  28   // m_ListOfSupports["S3P"] = m_ListOfSupports["Sup_to_Carina"];
  29   // m_ListOfSupports["S4R"] = m_ListOfSupports["Sup_to_Carina"];
  30   // m_ListOfSupports["S4L"] = m_ListOfSupports["Sup_to_Carina"];
  31   // m_ListOfSupports["S5"] = m_Mediastinum; // Not above Carina
  32   // m_ListOfSupports["S6"] = m_Mediastinum; // Not above Carina
  33
  34   // Read all support limits in a file and apply them
  35   ReadSupportLimits(GetSupportLimitsFilename());
  36   for(unsigned int i=0; i<m_ListOfSupportLimits.size(); i++) {
  37     SupportLimitsType s = m_ListOfSupportLimits[i];
  38     Support_SI_Limit(s.station_limit, s.station, s.structure_limit,
  39                      s.structure, s.offset*m_Working_Support->GetSpacing()[2]);
  40   }
  41
  42   // S1RL
  43   Support_LeftRight_S1R_S1L();
  44
  45   // S2RL
  46   Support_LeftRight_S2R_S2L();
  47
  48   // S4RL
  49   Support_LeftRight_S4R_S4L();
  50
  51   // Post limits of S1,S2,S4
  52   Support_Post_S1S2S4();
  53
  54   // S3P
  55   StartNewStep("[Support] Ant limits of S3P with trachea");
  56   m_ListOfSupports["S3P"] = LimitsWithTrachea(m_ListOfSupports["S3P"], 1, 0, 10);
  57
  58   // S3A
  59   StartNewStep("[Support] Ant limits of S3A with trachea");
  60   m_ListOfSupports["S3A"] = LimitsWithTrachea(m_ListOfSupports["S3A"], 1, 0, -10);
  61
  62   // I will do it later
  63   // Below Carina S7,8,9,10
  64   m_ListOfSupports["S7"] = clitk::Clone<MaskImageType>(m_ListOfSupports["Inf_to_Carina"]);
  65   m_ListOfSupports["S8"] = clitk::Clone<MaskImageType>(m_ListOfSupports["Inf_to_Carina"]);
  66   m_ListOfSupports["S9"] = clitk::Clone<MaskImageType>(m_ListOfSupports["Inf_to_Carina"]);
  67   m_ListOfSupports["S10"] = clitk::Clone<MaskImageType>(m_ListOfSupports["Inf_to_Carina"]);
  68   m_ListOfSupports["S11"] = clitk::Clone<MaskImageType>(m_ListOfSupports["Inf_to_Carina"]);
  69
  70   // Store image filenames into AFDB
  71   WriteImageSupport("S1R"); WriteImageSupport("S1L");
  72   WriteImageSupport("S2R"); WriteImageSupport("S2L");
  73   WriteImageSupport("S3A"); WriteImageSupport("S3P");
  74   WriteImageSupport("S4R"); WriteImageSupport("S4L");
  75   WriteImageSupport("S5");
  76   WriteImageSupport("S6");
  77   WriteImageSupport("S7");
  78   WriteImageSupport("S8");
  79   WriteImageSupport("S9");
  80   WriteImageSupport("S10");
  81   WriteImageSupport("S11");
  82   WriteAFDB();
  83 }
  84 //--------------------------------------------------------------------
  85
  86
  87 //--------------------------------------------------------------------
  88 template <class ImageType>
  89 void
  90 clitk::ExtractLymphStationsFilter<ImageType>::
  91 Support_SI_Limit(const std::string station_limit, const std::string station,
  92                  const std::string structure_limit, const std::string structure,
  93                  const double offset)
  94 {
  95   if (!GetCheckSupportFlag())
  96     StartNewStep("[Support] "+station_limit+" limit of "+station+" is "+structure_limit+" limit of "+structure);
  97
  98   // Check
  99   if ((station_limit != "superior") && (station_limit != "inferior")) {
 100     clitkExceptionMacro("Error station_limit must be 'inferior' or 'superior', not '"<< station_limit);
 101   }
 102   if ((structure_limit != "superior") && (structure_limit != "inferior")) {
 103     clitkExceptionMacro("Error structure_limit must be 'inferior' or 'superior', not '"<< structure_limit);
 104   }
 105
 106   // Get current support
 107   if (m_ListOfSupports.find(station) == m_ListOfSupports.end()) {
 108     // std::cerr << "Warning: support " << station << " not initialized" << std::endl;
 109     m_ListOfSupports[station] = m_Mediastinum;
 110   }
 111   m_Working_Support = m_ListOfSupports[station];
 112
 113   // Get structure or structureZ
 114   double z;
 115   int found=0;
 116   std::string file;
 117
 118   // Try to load structure and compute extrema point
 119   if (this->GetAFDB()->TagExist(structure)) {
 120     MaskImagePointer Structure = this->GetAFDB()->template GetImage <MaskImageType>(structure);
 121     file = this->GetAFDB()->GetTagValue(structure);
 122     MaskImagePointType p;
 123     p[0] = p[1] = p[2] =  0.0; // to avoid warning
 124     if (structure_limit == "superior")
 125       clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Structure, GetBackgroundValue(), 2, false, p);
 126     else
 127       clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Structure, GetBackgroundValue(), 2, true, p);
 128     z = p[2];
 129     found=1;
 130   }
 131
 132   // Try to load structureZ
 133   if ((found==0) && (this->GetAFDB()->TagExist(structure+"Z"))) {
 134     z = this->GetAFDB()->GetDouble(structure+"Z");
 135     found=2;
 136   }
 137
 138   // Try to load structurePoint
 139   if ((found==0) && (this->GetAFDB()->TagExist(structure+"Point"))) {
 140     MaskImagePointType p;
 141     this->GetAFDB()->GetPoint3D(structure+"Point", p);
 142     z = p[2];
 143     found=3;
 144   }
 145
 146   // Try to see if it is an already computed support
 147   if (found==0) {
 148     if (m_ListOfSupports.find(structure) != m_ListOfSupports.end()) {
 149       MaskImagePointer Structure = m_ListOfSupports[structure];
 150       MaskImagePointType p;
 151       if (structure_limit == "superior")
 152         clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Structure, GetBackgroundValue(), 2, false, p);
 153       else
 154         clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Structure, GetBackgroundValue(), 2, true, p);
 155       z = p[2];
 156       found=4;
 157     }
 158   }
 159
 160   // Try special case : "FindApexOfTheChest"
 161   if (structure == "FindApexOfTheChest") {
 162     z = FindApexOfTheChest();
 163     found=5;
 164   }
 165   if (structure == "FindInferiorBorderOfAorticArch") {
 166     z = FindInferiorBorderOfAorticArch();
 167     found=6;
 168   }
 169   if (structure == "FindSuperiorBorderOfAorticArch") {
 170     z = FindSuperiorBorderOfAorticArch();
 171     found=6;
 172   }
 173
 174   // If we find anything
 175   if (found == 0) {
 176     std::cerr << "ERROR : I could not find " << structure << " nor " << structure << "Z nor "
 177               << structure << "Point" << std::endl;
 178     exit(EXIT_FAILURE);
 179   }
 180
 181   // Apply offset
 182   z += offset;
 183
 184   // Remove Lower or greater
 185   if (station_limit == "inferior") {
 186     m_Working_Support =
 187       clitk::CropImageRemoveLowerThan<MaskImageType>(m_Working_Support, 2, z, true, GetBackgroundValue());
 188   }
 189   else {
 190     m_Working_Support =
 191       clitk::CropImageRemoveGreaterThan<MaskImageType>(m_Working_Support, 2, z, true, GetBackgroundValue());
 192   }
 193
 194   // Check: if reference station is given, display information
 195   if (GetCheckSupportFlag())  {
 196     try {
 197       MaskImagePointer Ref = this->GetAFDB()->template GetImage <MaskImageType>(station+"_Ref");
 198       MaskImagePointType p_support;
 199       MaskImagePointType p_ref;
 200       if (station_limit == "superior") {
 201         clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Ref, GetBackgroundValue(), 2, false, p_ref);
 202         clitk::FindExtremaPointInAGivenDirection<MaskImageType>(m_Working_Support, GetBackgroundValue(), 2, false, p_support);
 203       }
 204       else {
 205         clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Ref, GetBackgroundValue(), 2, true, p_ref);
 206         clitk::FindExtremaPointInAGivenDirection<MaskImageType>(m_Working_Support, GetBackgroundValue(), 2, true, p_support);
 207       }
 208       std::ostringstream os;
 209       os << "[Support] \t" << station << "\t" << station_limit << " "
 210          << "Z = " << z << std::setprecision(2) << std::fixed
 211          << "\tSupport = " << p_support[2]
 212          << "\tRef = " << p_ref[2]
 213          << "\tdiff = " << p_support[2]-p_ref[2] << "\t"
 214          << structure << " " << structure_limit;
 215       if (found==1) os << " (S "+file+")";
 216       if (found==2) os << " (Z)";
 217       if (found==3) os << " (P)";
 218       if (found==4) os << " (p)";
 219       if (found==5) os << " (Apex)";
 220       if (found==6) os << " (AorticArch)";
 221       StartNewStep(os.str());
 222     } catch(clitk::ExceptionObject e) { }
 223   }
 224
 225   // Set support
 226   m_ListOfSupports[station] = m_Working_Support;
 227   StopCurrentStep<MaskImageType>(m_Working_Support);
 228 }
 229 //--------------------------------------------------------------------
 230
 231
 232 //--------------------------------------------------------------------
 233 template <class ImageType>
 234 void
 235 clitk::ExtractLymphStationsFilter<ImageType>::
 236 Support_LeftRight_S1R_S1L()
 237 {
 238   // Step S1RL : Left-Right
 239   StartNewStep("[Support] Left-Right S1R S1L");
 240   std::vector<ImagePointType> A;
 241   std::vector<ImagePointType> B;
 242   // Search for centroid positions of trachea
 243   MaskImagePointer Trachea = this->GetAFDB()->template GetImage <MaskImageType>("Trachea");
 244   MaskImagePointer S1RL = m_ListOfSupports["S1R"];
 245   Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, S1RL, GetBackgroundValue());
 246   std::vector<MaskSlicePointer> slices;
 247   clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices);
 248   for(uint i=0; i<slices.size(); i++) {
 249     slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
 250     std::vector<typename MaskSliceType::PointType> c;
 251     clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
 252     ImagePointType a,b;
 253     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(c[1], Trachea, i, a);
 254     A.push_back(a);
 255     b = a;
 256     b[1] += 50;
 257     B.push_back(b);
 258   }
 259   clitk::WriteListOfLandmarks<MaskImageType>(A, "S1LR_A.txt");
 260   clitk::WriteListOfLandmarks<MaskImageType>(B, "S1LR_B.txt");
 261
 262   // Clone support
 263   MaskImagePointer S1R = clitk::Clone<MaskImageType>(S1RL);
 264   MaskImagePointer S1L = clitk::Clone<MaskImageType>(S1RL);
 265
 266   // Right part
 267   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(S1R, A, B,
 268                                                                     GetBackgroundValue(), 0, 10);
 269   S1R = clitk::AutoCrop<MaskImageType>(S1R, GetBackgroundValue());
 270   m_ListOfSupports["S1R"] = S1R;
 271
 272   // Left part
 273   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(S1L, A, B,
 274                                                                     GetBackgroundValue(), 0, -10);
 275   S1L = clitk::AutoCrop<MaskImageType>(S1L, GetBackgroundValue());
 276   m_ListOfSupports["S1L"] = S1L;
 277   StopCurrentStep<MaskImageType>(m_ListOfSupports["S1L"]);
 278 }
 279 //--------------------------------------------------------------------
 280
 281
 282 //--------------------------------------------------------------------
 283 template <class ImageType>
 284 void
 285 clitk::ExtractLymphStationsFilter<ImageType>::
 286 Support_LeftRight_S2R_S2L()
 287 {
 288   // ---------------------------------------------------------------------------
 289   /* Step : S2RL LeftRight
 290      As for lymph node station 4R, 2R includes nodes extending to the
 291      left lateral border of the trachea
 292      Rod says:  "For station 2 there is a shift, dividing 2R from 2L, from midline
 293      to the left paratracheal border."
 294   */
 295   StartNewStep("[Support] Separate 2R/2L according to Trachea");
 296   MaskImagePointer S2R = m_ListOfSupports["S2R"];
 297   MaskImagePointer S2L = m_ListOfSupports["S2L"];
 298   S2R = LimitsWithTrachea(S2R, 0, 1, -10);
 299   S2L = LimitsWithTrachea(S2L, 0, 1, 10);
 300   S2R = clitk::AutoCrop<MaskImageType>(S2R, GetBackgroundValue());
 301   S2L = clitk::AutoCrop<MaskImageType>(S2L, GetBackgroundValue());
 302   m_ListOfSupports["S2R"] = S2R;
 303   m_ListOfSupports["S2L"] = S2L;
 304   this->GetAFDB()->template ReleaseImage<MaskImageType>("Trachea");
 305 }
 306 //--------------------------------------------------------------------
 307
 308
 309 //--------------------------------------------------------------------
 310 template <class ImageType>
 311 void
 312 clitk::ExtractLymphStationsFilter<ImageType>::
 313 Support_LeftRight_S4R_S4L()
 314 {
 315   // ---------------------------------------------------------------------------
 316   /* Step : S4RL LeftRight
 317
 318      - 4R: includes right paratracheal nodes, and pretracheal nodes
 319      extending to the left lateral border of trachea
 320
 321      - 4L: includes nodes to the left of the left lateral border of
 322      the trachea, medial to the ligamentum arteriosum
 323
 324      => same than 2RL
 325   */
 326   StartNewStep("[Support] Left Right separation of 4R/4L");
 327
 328   MaskImagePointer S4R = m_ListOfSupports["S4R"];
 329   MaskImagePointer S4L = m_ListOfSupports["S4L"];
 330   S4R = LimitsWithTrachea(S4R, 0, 1, -10);
 331   S4L = LimitsWithTrachea(S4L, 0, 1, 10);
 332   m_ListOfSupports["S4R"] = S4R;
 333   m_ListOfSupports["S4L"] = S4L;
 334 }
 335 //--------------------------------------------------------------------
 336
 337
 338 //--------------------------------------------------------------------
 339 template <class ImageType>
 340 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
 341 clitk::ExtractLymphStationsFilter<ImageType>::
 342 LimitsWithTrachea(MaskImageType * input, int extremaDirection, int lineDirection,
 343                   double offset)
 344 {
 345   MaskImagePointType min, max;
 346   GetMinMaxBoundary<MaskImageType>(input, min, max);
 347   return LimitsWithTrachea(input, extremaDirection, lineDirection, offset, max[2]);
 348 }
 349 template <class ImageType>
 350 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
 351 clitk::ExtractLymphStationsFilter<ImageType>::
 352 LimitsWithTrachea(MaskImageType * input, int extremaDirection, int lineDirection,
 353                   double offset, double maxSupPosition)
 354 {
 355   /*
 356     Take the input mask, consider the trachea and limit according to
 357     Left border of the trachea. Keep at Left or at Right according to
 358     the offset
 359   */
 360   // Read the trachea
 361   MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
 362
 363   // Find extrema post positions
 364   std::vector<MaskImagePointType> tracheaLeftPositionsA;
 365   std::vector<MaskImagePointType> tracheaLeftPositionsB;
 366
 367   // Crop Trachea only on the Sup-Inf axes, without autocrop
 368   //  Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, input, GetBackgroundValue());
 369   MaskImagePointType min, max;
 370   GetMinMaxBoundary<MaskImageType>(input, min, max);
 371   Trachea = clitk::CropImageAlongOneAxis<MaskImageType>(Trachea, 2, min[2], max[2],
 372                                                         false, GetBackgroundValue());
 373
 374   // Select the main CCL (because of bronchus)
 375   Trachea = SliceBySliceKeepMainCCL<MaskImageType>(Trachea, GetBackgroundValue(), GetForegroundValue());
 376
 377   // Slice by slice, build the separation line
 378   clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition<MaskImageType>(Trachea,
 379                                                                                GetBackgroundValue(), 2,
 380                                                                                extremaDirection, false, // Left
 381                                                                                lineDirection, // Vertical line
 382                                                                                -1, // margins
 383                                                                                tracheaLeftPositionsA,
 384                                                                                tracheaLeftPositionsB);
 385   // Do not consider trachea above the limit
 386   int indexMax=tracheaLeftPositionsA.size();
 387   for(uint i=0; i<tracheaLeftPositionsA.size(); i++) {
 388     if (tracheaLeftPositionsA[i][2] > maxSupPosition) {
 389       indexMax = i;
 390       i = tracheaLeftPositionsA.size(); // stop loop
 391     }
 392   }
 393   tracheaLeftPositionsA.erase(tracheaLeftPositionsA.begin()+indexMax, tracheaLeftPositionsA.end());
 394   tracheaLeftPositionsB.erase(tracheaLeftPositionsB.begin()+indexMax, tracheaLeftPositionsB.end());
 395
 396   // Cut post to this line
 397   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(input,
 398                                                                     tracheaLeftPositionsA,
 399                                                                     tracheaLeftPositionsB,
 400                                                                     GetBackgroundValue(),
 401                                                                     extremaDirection, offset);
 402   MaskImagePointer output = clitk::AutoCrop<MaskImageType>(input, GetBackgroundValue());
 403   return output;
 404 }
 405 //--------------------------------------------------------------------
 406
 407
 408 //--------------------------------------------------------------------
 409 template <class ImageType>
 410 void
 411 clitk::ExtractLymphStationsFilter<ImageType>::
 412 Support_Post_S1S2S4()
 413 {
 414   StartNewStep("[Support] Post limits of S1RL, S2RL, S4RL");
 415
 416   double m_ApexOfTheChest = FindApexOfTheChest();
 417
 418   // Post limits with Trachea
 419   MaskImagePointer S1R = m_ListOfSupports["S1R"];
 420   MaskImagePointer S1L = m_ListOfSupports["S1L"];
 421   MaskImagePointer S2R = m_ListOfSupports["S2R"];
 422   MaskImagePointer S2L = m_ListOfSupports["S2L"];
 423   MaskImagePointer S4R = m_ListOfSupports["S4R"];
 424   MaskImagePointer S4L = m_ListOfSupports["S4L"];
 425   m_ListOfSupports["S1R"] = LimitsWithTrachea(S1L, 1, 0, -10, m_ApexOfTheChest);
 426   m_ListOfSupports["S1L"] = LimitsWithTrachea(S1R, 1, 0, -10, m_ApexOfTheChest);
 427   m_ListOfSupports["S2R"] = LimitsWithTrachea(S2R, 1, 0, -10, m_ApexOfTheChest);
 428   m_ListOfSupports["S2L"] = LimitsWithTrachea(S2L, 1, 0, -10, m_ApexOfTheChest);
 429   m_ListOfSupports["S4R"] = LimitsWithTrachea(S4R, 1, 0, -10, m_ApexOfTheChest);
 430   m_ListOfSupports["S4L"] = LimitsWithTrachea(S4L, 1, 0, -10, m_ApexOfTheChest);
 431 }
 432 //--------------------------------------------------------------------
 433
 434