segmentation/clitkExtractLymphStation_3P.txx

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