segmentation/clitkExtractLymphStation_3P.txx

   1
   2 #include <itkBinaryDilateImageFilter.h>
   3 #include <itkMirrorPadImageFilter.h>
   4
   5 //--------------------------------------------------------------------
   6 template <class ImageType>
   7 void
   8 clitk::ExtractLymphStationsFilter<ImageType>::
   9 ExtractStation_3P_SetDefaultValues()
  10 {
  11 }
  12 //--------------------------------------------------------------------
  13
  14
  15 //--------------------------------------------------------------------
  16 template <class ImageType>
  17 void
  18 clitk::ExtractLymphStationsFilter<ImageType>::
  19 ExtractStation_3P_SI_Limits()
  20 {
  21   /*
  22     Apex of the chest & Carina.
  23   */
  24   StartNewStep("[Station 3P] Inf/Sup limits with apex of the chest and carina");
  25
  26   // Get Carina position (has been determined in Station8)
  27   m_CarinaZ = GetAFDB()->GetDouble("CarinaZ");
  28
  29   // Get Apex of the Chest. The "lungs" structure is autocroped so we
  30   // consider the most superior point.
  31   MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
  32   MaskImageIndexType index = Lungs->GetLargestPossibleRegion().GetIndex();
  33   index += Lungs->GetLargestPossibleRegion().GetSize();
  34   MaskImagePointType p;
  35   Lungs->TransformIndexToPhysicalPoint(index, p);
  36   p[2] -= 5; // 5 mm below because the autocrop is slightly above the apex
  37   double m_ApexOfTheChest = p[2];
  38
  39   /* Crop support :
  40      Superior limit = carina
  41      Inferior limit = Apex of the chest */
  42   m_Working_Support =
  43     clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
  44                                                 m_CarinaZ,
  45                                                 m_ApexOfTheChest, true,
  46                                                 GetBackgroundValue());
  47
  48   StopCurrentStep<MaskImageType>(m_Working_Support);
  49   m_ListOfStations["3P"] = m_Working_Support;
  50 }
  51 //--------------------------------------------------------------------
  52
  53
  54 //--------------------------------------------------------------------
  55 template <class ImageType>
  56 void
  57 clitk::ExtractLymphStationsFilter<ImageType>::
  58 ExtractStation_3P_Remove_Structures()
  59 {
  60   /*
  61     3 - (sup) remove Aorta, VB, subcl
  62     not LR subcl ! -> a séparer LR ?
  63     (inf) remove Eso, Aorta, Azygousvein
  64   */
  65
  66   StartNewStep("[Station 3P] Remove some structures.");
  67
  68   Remove_Structures("Aorta");
  69   Remove_Structures("VertebralBody");
  70   Remove_Structures("SubclavianArtery");
  71   Remove_Structures("Esophagus");
  72   Remove_Structures("Azygousvein");
  73   Remove_Structures("Thyroid");
  74   Remove_Structures("VertebralArtery");
  75
  76   StopCurrentStep<MaskImageType>(m_Working_Support);
  77   m_ListOfStations["3P"] = m_Working_Support;
  78 }
  79 //--------------------------------------------------------------------
  80
  81
  82 //--------------------------------------------------------------------
  83 template <class ImageType>
  84 void
  85 clitk::ExtractLymphStationsFilter<ImageType>::
  86 ExtractStation_3P_Ant_Limits()
  87 {
  88   /*
  89     Ant Post limit :
  90
  91     Anteriorly, it is in contact with the posterior aspect of Stations
  92     1–2 superiorly (Fig. 2A–C) and with Stations 4R and 4L inferiorly
  93     (Fig. 2D–I and 3A–C). The anterior limit of Station 3P is kept
  94     posterior to the trachea, which is defined by an imaginary
  95     horizontal line running along the posterior wall of the trachea
  96     (Fig. 2B,E, red line). Posteriorly, it is delineated along the
  97     anterior and lateral borders of the vertebral body until an
  98     imaginary horizontal line running 1 cm posteriorly to the
  99     anterior border of the vertebral body (Fig. 2D).
 100
 101     1 - post to the trachea : define an imaginary line just post the
 102     Trachea and remove what is anterior.
 103
 104   */
 105   StartNewStep("[Station 3P] Ant limits with Trachea ");
 106
 107   // Load Trachea
 108   MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
 109
 110   // Crop like current support (need by SliceBySliceSetBackgroundFromLineSeparation after)
 111   Trachea =
 112     clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
 113
 114   // Slice by slice, determine the most post point of the trachea (A)
 115   // and consider a second point on the same horizontal line (B)
 116   std::vector<MaskImagePointType> p_A;
 117   std::vector<MaskImagePointType> p_B;
 118   std::vector<typename MaskSliceType::Pointer> slices;
 119   clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices);
 120   MaskImagePointType p;
 121   typename MaskSliceType::PointType sp;
 122   for(uint i=0; i<slices.size(); i++) {
 123     // First only consider the main CCL (trachea, not bronchus)
 124     slices[i] = Labelize<MaskSliceType>(slices[i], 0, true, 100);
 125     slices[i] = KeepLabels<MaskSliceType>(slices[i], GetBackgroundValue(),
 126                                           GetForegroundValue(), 1, 1, true);
 127     // Find most posterior point
 128     clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i], GetBackgroundValue(),
 129                                                             1, false, sp);
 130     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp, Trachea, i, p);
 131     p_A.push_back(p);
 132     p[0] -= 20;
 133     p_B.push_back(p);
 134   }
 135   clitk::WriteListOfLandmarks<MaskImageType>(p_A, "trachea-post.txt");
 136
 137   // Remove Ant part above those lines
 138   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 139                                                                     p_A, p_B,
 140                                                                     GetBackgroundValue(),
 141                                                                     1, 10);
 142   // End, release memory
 143   GetAFDB()->template ReleaseImage<MaskImageType>("Trachea");
 144   StopCurrentStep<MaskImageType>(m_Working_Support);
 145   m_ListOfStations["3P"] = m_Working_Support;
 146 }
 147 //--------------------------------------------------------------------
 148
 149
 150 //--------------------------------------------------------------------
 151 template <class ImageType>
 152 void
 153 clitk::ExtractLymphStationsFilter<ImageType>::
 154 ExtractStation_3P_Post_Limits()
 155 {
 156   /*
 157     Ant Post limit :
 158
 159     Anteriorly, it is in contact with the posterior aspect of Stations
 160     1–2 superiorly (Fig. 2A–C) and with Stations 4R and 4L inferiorly
 161     (Fig. 2D–I and 3A–C). The anterior limit of Station 3P is kept
 162     posterior to the trachea, which is defined by an imaginary
 163     horizontal line running along the posterior wall of the trachea
 164     (Fig. 2B,E, red line). Posteriorly, it is delineated along the
 165     anterior and lateral borders of the vertebral body until an
 166     imaginary horizontal line running 1 cm posteriorly to the
 167     anterior border of the vertebral body (Fig. 2D).
 168
 169     2 - post to the trachea : define an imaginary line just post the
 170     Trachea and remove what is anterior.
 171
 172   */
 173   StartNewStep("[Station 3P] Post limits with VertebralBody ");
 174
 175   // Load VertebralBody
 176   MaskImagePointer VertebralBody = GetAFDB()->template GetImage<MaskImageType>("VertebralBody");
 177
 178   // Crop like current support (need by SliceBySliceSetBackgroundFromLineSeparation after)
 179   VertebralBody = clitk::ResizeImageLike<MaskImageType>(VertebralBody, m_Working_Support, GetBackgroundValue());
 180
 181   writeImage<MaskImageType>(VertebralBody, "vb.mhd");
 182
 183   // Compute VertebralBody most Ant position (again because slices
 184   // changes). Slice by slice, determine the most post point of the
 185   // trachea (A) and consider a second point on the same horizontal
 186   // line (B)
 187   std::vector<MaskImagePointType> p_A;
 188   std::vector<MaskImagePointType> p_B;
 189   std::vector<typename MaskSliceType::Pointer> slices;
 190   clitk::ExtractSlices<MaskImageType>(VertebralBody, 2, slices);
 191   MaskImagePointType p;
 192   typename MaskSliceType::PointType sp;
 193   for(uint i=0; i<slices.size(); i++) {
 194     // Find most anterior point
 195     bool found = clitk::FindExtremaPointInAGivenDirection<MaskSliceType>(slices[i], GetBackgroundValue(),
 196                                                                          1, true, sp);
 197
 198     // If the VertebralBody stop superiorly before the end of
 199     // m_Working_Support, we consider the same previous point.
 200     if (!found) {
 201       p = p_A.back();
 202       p[2] += VertebralBody->GetSpacing()[2];
 203       p_A.push_back(p);
 204       p = p_B.back();
 205       p[2] += VertebralBody->GetSpacing()[2];
 206       p_B.push_back(p);
 207     }
 208     else {
 209       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(sp, VertebralBody, i, p);
 210       p[1] += 10; // Consider 10 mm more post
 211       p_A.push_back(p);
 212       p[0] -= 20;
 213       p_B.push_back(p);
 214     }
 215   }
 216   clitk::WriteListOfLandmarks<MaskImageType>(p_A, "vb-ant.txt");
 217
 218
 219   // Remove Ant part above those lines
 220   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 221                                                                     p_A, p_B,
 222                                                                     GetBackgroundValue(),
 223                                                                     1, -10);
 224   writeImage<MaskImageType>(m_Working_Support, "a.mhd");
 225
 226   StopCurrentStep<MaskImageType>(m_Working_Support);
 227   m_ListOfStations["3P"] = m_Working_Support;
 228 }
 229 //--------------------------------------------------------------------
 230
 231
 232 //--------------------------------------------------------------------
 233 template <class ImageType>
 234 void
 235 clitk::ExtractLymphStationsFilter<ImageType>::
 236 ExtractStation_3P_LR_sup_Limits()
 237 {
 238   /*
 239     "On the right side, the limit is deﬁned by the air–soft-tissue
 240     interface. On the left side, it is deﬁned by the air–tissue
 241     interface superiorly (Fig. 2A–C) and the aorta inferiorly
 242     (Figs. 2D–I and 3A–C)."
 243
 244     sup :
 245     Resizelike support : Trachea, SubclavianArtery
 246     Trachea, slice by slice, get centroid trachea
 247     SubclavianArtery, slice by slice, CCL
 248     prendre la 1ère à L et R, not at Left
 249
 250   */
 251   StartNewStep("[Station 3P] Left/Right limits (superior part) ");
 252
 253   // Load structures
 254   MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
 255   MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
 256
 257   // Crop like current support
 258   Trachea = clitk::ResizeImageLike<MaskImageType>(Trachea, m_Working_Support, GetBackgroundValue());
 259   SubclavianArtery = clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, m_Working_Support, GetBackgroundValue());
 260
 261   writeImage<MaskImageType>(Trachea, "tr.mhd");
 262   writeImage<MaskImageType>(SubclavianArtery, "sca.mhd");
 263
 264   // Get list of slices
 265   std::vector<MaskSlicePointer> slices_support;
 266   std::vector<MaskSlicePointer> slices_trachea;
 267   std::vector<MaskSlicePointer> slices_subclavianartery;
 268   clitk::ExtractSlices<MaskImageType>(m_Working_Support, 2, slices_support);
 269   clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_trachea);
 270   clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_subclavianartery);
 271
 272   // Loop on slices
 273   std::vector<MaskImagePointType> points;
 274   MaskImagePointType p;
 275   for(uint i=0; i<slices_support.size(); i++) {
 276     // Get Trachea centroid
 277     std::vector<typename MaskSliceType::PointType> centroids;
 278     typename MaskSliceType::PointType c;
 279     ComputeCentroids<MaskSliceType>(slices_trachea[i], GetBackgroundValue(), centroids);
 280     c = centroids[1];
 281
 282     // [debug] Store point
 283     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[1], Trachea, i, p);
 284     points.push_back(p);
 285
 286     // Get Right and Left CCL in SubclavianArtery
 287     slices_subclavianartery[i] = Labelize<MaskSliceType>(slices_subclavianartery[i], 0, true, 10);
 288     ComputeCentroids<MaskSliceType>(slices_subclavianartery[i], GetBackgroundValue(), centroids);
 289
 290     if (centroids.size() > 1) {
 291       // Determine the one at Right/Left -> first after Trachea
 292       // centroid
 293       typename MaskSliceType::PointType right;
 294       typename MaskSliceType::PointType left;
 295       int label_right=-1;
 296       int label_left=-1;
 297       right[0] = c[0]-100;
 298       left[0] = c[0]+100;
 299       for(uint j=1; j<centroids.size(); j++) {
 300         if (centroids[j][0] < c[0]) { // At Right of Trachea centroid
 301           if (centroids[j][0] >= right[0]) {
 302             right = centroids[j];
 303             label_right = j;
 304           }
 305         }
 306         if (centroids[j][0] > c[0]) { // At Left of Trachea centroid
 307           if (centroids[j][0] <= left[0]) {
 308             left = centroids[j];
 309             label_left = j;
 310           }
 311         }
 312       }
 313
 314       if (label_right != -1) {
 315
 316         // Debug points
 317         clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[label_right], SubclavianArtery, i, p);
 318         points.push_back(p);
 319
 320         // Set Background and ForegroundValue according to label_right
 321         MaskSlicePointer object =
 322           clitk::Binarize<MaskSliceType>(slices_subclavianartery[i], label_right, label_right,
 323                                          GetBackgroundValue(), GetForegroundValue());
 324
 325         // Relative Position : not at Right
 326         typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskSliceType> RelPosFilterType;
 327         typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 328         relPosFilter->VerboseStepFlagOff();
 329         relPosFilter->WriteStepFlagOff();
 330         relPosFilter->SetBackgroundValue(GetBackgroundValue());
 331         relPosFilter->SetInput(slices_support[i]);
 332         relPosFilter->SetInputObject(object);
 333         relPosFilter->AddOrientationTypeString("R");
 334         relPosFilter->SetInverseOrientationFlag(true);
 335         //      relPosFilter->SetIntermediateSpacing(3);
 336         relPosFilter->SetIntermediateSpacingFlag(false);
 337         relPosFilter->SetFuzzyThreshold(0.7);
 338         relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
 339         relPosFilter->Update();
 340         slices_support[i] = relPosFilter->GetOutput();
 341
 342         // Relative Position : not Anterior
 343         relPosFilter = RelPosFilterType::New();
 344         relPosFilter->VerboseStepFlagOff();
 345         relPosFilter->WriteStepFlagOff();
 346         relPosFilter->SetBackgroundValue(GetBackgroundValue());
 347         relPosFilter->SetInput(slices_support[i]);
 348         relPosFilter->SetInputObject(object);
 349         relPosFilter->AddOrientationTypeString("A");
 350         relPosFilter->SetInverseOrientationFlag(true);
 351         //      relPosFilter->SetIntermediateSpacing(3);
 352         relPosFilter->SetIntermediateSpacingFlag(false);
 353         relPosFilter->SetFuzzyThreshold(0.7);
 354         relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
 355         relPosFilter->Update();
 356         slices_support[i] = relPosFilter->GetOutput();
 357
 358       } // End RelativePosition for Right
 359
 360
 361       if (label_left != -1) {
 362
 363         // Debug points
 364         clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[label_left], SubclavianArtery, i, p);
 365         points.push_back(p);
 366
 367         // Set Background and ForegroundValue according to label_right
 368         MaskSlicePointer object =
 369           clitk::Binarize<MaskSliceType>(slices_subclavianartery[i], label_left, label_left,
 370                                          GetBackgroundValue(), GetForegroundValue());
 371
 372         // Relative Position : not at Right
 373         typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskSliceType> RelPosFilterType;
 374         typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 375         relPosFilter->VerboseStepFlagOff();
 376         relPosFilter->WriteStepFlagOff();
 377         relPosFilter->SetBackgroundValue(GetBackgroundValue());
 378         relPosFilter->SetInput(slices_support[i]);
 379         relPosFilter->SetInputObject(object);
 380         relPosFilter->AddOrientationTypeString("L");
 381         relPosFilter->SetInverseOrientationFlag(true);
 382         //      relPosFilter->SetIntermediateSpacing(3);
 383         relPosFilter->SetIntermediateSpacingFlag(false);
 384         relPosFilter->SetFuzzyThreshold(0.7);
 385         relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
 386         relPosFilter->Update();
 387         slices_support[i] = relPosFilter->GetOutput();
 388
 389         // Relative Position : not Anterior
 390         relPosFilter = RelPosFilterType::New();
 391         relPosFilter->VerboseStepFlagOff();
 392         relPosFilter->WriteStepFlagOff();
 393         relPosFilter->SetBackgroundValue(GetBackgroundValue());
 394         relPosFilter->SetInput(slices_support[i]);
 395         relPosFilter->SetInputObject(object);
 396         relPosFilter->AddOrientationTypeString("A");
 397         relPosFilter->SetInverseOrientationFlag(true);
 398         //      relPosFilter->SetIntermediateSpacing(3);
 399         relPosFilter->SetIntermediateSpacingFlag(false);
 400         relPosFilter->SetFuzzyThreshold(0.7);
 401         relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
 402         relPosFilter->Update();
 403         slices_support[i] = relPosFilter->GetOutput();
 404
 405       }
 406
 407
 408     } // if centroids.size > 1
 409   } // End loop slices
 410
 411   // Joint slices
 412   m_Working_Support = clitk::JoinSlices<MaskImageType>(slices_support, m_Working_Support, 2);
 413
 414   // Save list points
 415   clitk::WriteListOfLandmarks<MaskImageType>(points, "subcl-lr.txt");
 416
 417
 418   StopCurrentStep<MaskImageType>(m_Working_Support);
 419   m_ListOfStations["3P"] = m_Working_Support;
 420 }
 421 //--------------------------------------------------------------------
 422
 423 //--------------------------------------------------------------------
 424 template <class ImageType>
 425 void
 426 clitk::ExtractLymphStationsFilter<ImageType>::
 427 ExtractStation_3P_LR_sup_Limits_2()
 428 {
 429   /*
 430     Use VertebralArtery to limit.
 431
 432
 433   StartNewStep("[Station 3P] Left/Right limits with VertebralArtery");
 434
 435   // Load structures
 436   MaskImagePointer VertebralArtery = GetAFDB()->template GetImage<MaskImageType>("VertebralArtery");
 437
 438   clitk::AndNot<MaskImageType>(m_Working_Support, VertebralArtery);
 439
 440   StopCurrentStep<MaskImageType>(m_Working_Support);
 441   m_ListOfStations["3P"] = m_Working_Support;
 442   */
 443 }
 444 //--------------------------------------------------------------------
 445
 446 //--------------------------------------------------------------------
 447 template <class ImageType>
 448 void
 449 clitk::ExtractLymphStationsFilter<ImageType>::
 450 ExtractStation_3P_LR_inf_Limits()
 451 {
 452   /*
 453     "On the right side, the limit is deﬁned by the air–soft-tissue
 454     interface. On the left side, it is deﬁned by the air–tissue
 455     interface superiorly (Fig. 2A–C) and the aorta inferiorly
 456     (Figs. 2D–I and 3A–C)."
 457
 458     inf : not Right to Azygousvein
 459   */
 460   StartNewStep("[Station 3P] Left/Right limits (inferior part) with Azygousvein and Aorta");
 461
 462   // Load structures
 463   MaskImagePointer AzygousVein = GetAFDB()->template GetImage<MaskImageType>("AzygousVein");
 464   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 465
 466   typedef clitk::AddRelativePositionConstraintToLabelImageFilter<MaskImageType> RelPosFilterType;
 467   typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 468   relPosFilter->VerboseStepFlagOff();
 469   relPosFilter->WriteStepFlagOff();
 470   relPosFilter->SetBackgroundValue(GetBackgroundValue());
 471   relPosFilter->SetInput(m_Working_Support);
 472   relPosFilter->SetInputObject(AzygousVein);
 473   relPosFilter->AddOrientationTypeString("R");
 474   relPosFilter->SetInverseOrientationFlag(true);
 475   //      relPosFilter->SetIntermediateSpacing(3);
 476   relPosFilter->SetIntermediateSpacingFlag(false);
 477   relPosFilter->SetFuzzyThreshold(0.7);
 478   relPosFilter->AutoCropFlagOn();
 479   relPosFilter->Update();
 480   m_Working_Support = relPosFilter->GetOutput();
 481
 482   writeImage<MaskImageType>(m_Working_Support, "before-L-aorta.mhd");
 483   relPosFilter->SetInput(m_Working_Support);
 484   relPosFilter->SetInputObject(Aorta);
 485   relPosFilter->AddOrientationTypeString("L");
 486   relPosFilter->SetInverseOrientationFlag(true);
 487   //      relPosFilter->SetIntermediateSpacing(3);
 488   relPosFilter->SetIntermediateSpacingFlag(false);
 489   relPosFilter->SetFuzzyThreshold(0.7);
 490   relPosFilter->AutoCropFlagOn();
 491   relPosFilter->Update();
 492   m_Working_Support = relPosFilter->GetOutput();
 493   writeImage<MaskImageType>(m_Working_Support, "after-L-aorta.mhd");
 494
 495   StopCurrentStep<MaskImageType>(m_Working_Support);
 496   m_ListOfStations["3P"] = m_Working_Support;
 497 }
 498 //--------------------------------------------------------------------