segmentation/clitkExtractLymphStation_7.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 //--------------------------------------------------------------------
  20 template <class TImageType>
  21 void
  22 clitk::ExtractLymphStationsFilter<TImageType>::
  23 ExtractStation_7_SetDefaultValues()
  24 {
  25   SetFuzzyThresholdForS7("Bronchi", 0.1);
  26   SetFuzzyThresholdForS7("LeftSuperiorPulmonaryVein", 0.3);
  27   SetFuzzyThresholdForS7("RightSuperiorPulmonaryVein", 0.2);
  28   SetFuzzyThresholdForS7("RightPulmonaryArtery", 0.3);
  29   SetFuzzyThresholdForS7("LeftPulmonaryArtery", 0.5);
  30   SetFuzzyThresholdForS7("SVC", 0.2);
  31 }
  32 //--------------------------------------------------------------------
  33
  34 //--------------------------------------------------------------------
  35 template <class TImageType>
  36 void
  37 clitk::ExtractLymphStationsFilter<TImageType>::
  38 SetFuzzyThresholdForS7(std::string tag, double value)
  39 {
  40   m_FuzzyThresholdForS7[tag] = value;
  41 }
  42 //--------------------------------------------------------------------
  43
  44
  45 //--------------------------------------------------------------------
  46 template <class TImageType>
  47 double
  48 clitk::ExtractLymphStationsFilter<TImageType>::
  49 GetFuzzyThresholdForS7(std::string tag)
  50 {
  51   if (m_FuzzyThresholdForS7.find(tag) != m_FuzzyThresholdForS7.end()) {
  52     return m_FuzzyThresholdForS7[tag];
  53   }
  54   else {
  55     clitkExceptionMacro("Could not find options "+tag+" in the m_FuzzyThresholdForS7 list");
  56   }
  57 }
  58 //--------------------------------------------------------------------
  59
  60
  61 //--------------------------------------------------------------------
  62 template <class TImageType>
  63 void
  64 clitk::ExtractLymphStationsFilter<TImageType>::
  65 ExtractStation_7_SI_Limits()
  66 {
  67   StartNewStep("[Station7] Inf/Sup mediastinum limits with carina/LLLBronchus");
  68   // Get Inputs
  69   MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
  70
  71   // We suppoe that CarinaZ was already computed (S8)
  72   double m_CarinaZ = GetAFDB()->GetDouble("CarinaZ");
  73
  74   //  double m_OriginOfRightMiddleLobeBronchusZ = GetAFDB()->GetPoint3D("OriginOfRightMiddleLobeBronchus", 2);
  75   // DD(m_OriginOfRightMiddleLobeBronchusZ);
  76   MaskImagePointer UpperBorderOfLLLBronchus = GetAFDB()->template GetImage<MaskImageType>("UpperBorderOfLLLBronchus");
  77
  78   // Search most inf point (WHY ? IS IT THE RIGHT STRUCTURE ??)
  79   MaskImagePointType ps = UpperBorderOfLLLBronchus->GetOrigin(); // initialise to avoid warning
  80   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(UpperBorderOfLLLBronchus, GetBackgroundValue(), 2, true, ps);
  81   double m_UpperBorderOfLLLBronchusZ = ps[2];
  82
  83   /*
  84   std::vector<MaskImagePointType> centroids;
  85   clitk::ComputeCentroids<MaskImageType>(UpperBorderOfLLLBronchus, GetBackgroundValue(), centroids);
  86   double m_UpperBorderOfLLLBronchusZ = centroids[1][2];
  87   DD(m_UpperBorderOfLLLBronchusZ)
  88   */
  89
  90   /* Crop support */
  91   m_Working_Support =
  92     clitk::CropImageAlongOneAxis<MaskImageType>(m_Working_Support, 2,
  93                                                 m_UpperBorderOfLLLBronchusZ,
  94                                                 m_CarinaZ, true,
  95                                                 GetBackgroundValue());
  96   /* Crop trachea */
  97   m_Working_Trachea =
  98     clitk::CropImageAlongOneAxis<MaskImageType>(Trachea, 2,
  99                                                 m_UpperBorderOfLLLBronchusZ,
 100                                                 m_CarinaZ, true,
 101                                                 GetBackgroundValue());
 102
 103   StopCurrentStep<MaskImageType>(m_Working_Support);
 104   m_ListOfStations["7"] = m_Working_Support;
 105 }
 106 //--------------------------------------------------------------------
 107
 108
 109 //--------------------------------------------------------------------
 110 template <class TImageType>
 111 void
 112 clitk::ExtractLymphStationsFilter<TImageType>::
 113 ExtractStation_7_RL_Limits()
 114 {
 115   // ----------------------------------------------------------------
 116   StartNewStep("[Station7] Limits with bronchus : RightTo the left bronchus");
 117
 118   // First consider bronchus and keep the main CCL by slice
 119   m_RightBronchus = GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
 120   m_LeftBronchus = GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
 121
 122   // Extract slices, Label, compute centroid, keep most central connected component
 123   std::vector<MaskSlicePointer> slices_leftbronchus;
 124   std::vector<MaskSlicePointer> slices_rightbronchus;
 125   clitk::ExtractSlices<MaskImageType>(m_LeftBronchus, 2, slices_leftbronchus);
 126   clitk::ExtractSlices<MaskImageType>(m_RightBronchus, 2, slices_rightbronchus);
 127
 128   // Loop on slices
 129   for(uint i=0; i<slices_leftbronchus.size(); i++) {
 130     slices_leftbronchus[i] = Labelize<MaskSliceType>(slices_leftbronchus[i], 0, false, 10);
 131     std::vector<typename MaskSliceType::PointType> c;
 132     clitk::ComputeCentroids<MaskSliceType>(slices_leftbronchus[i], GetBackgroundValue(), c);
 133     if (c.size() > 1) {
 134       double most_at_left = c[1][0];
 135       int most_at_left_index=1;
 136       for(uint j=1; j<c.size(); j++) {
 137         if (c[j][0] < most_at_left) {
 138           most_at_left = c[j][0];
 139           most_at_left_index = j;
 140         }
 141       }
 142       // Put all other CCL to Background
 143       slices_leftbronchus[i] =
 144         clitk::Binarize<MaskSliceType>(slices_leftbronchus[i], most_at_left_index,
 145                                        most_at_left_index, GetBackgroundValue(), GetForegroundValue());
 146     } // end c.size
 147   }
 148
 149   for(uint i=0; i<slices_rightbronchus.size(); i++) {
 150     slices_rightbronchus[i] = Labelize<MaskSliceType>(slices_rightbronchus[i], 0, false, 10);
 151     std::vector<typename MaskSliceType::PointType> c;
 152     clitk::ComputeCentroids<MaskSliceType>(slices_rightbronchus[i], GetBackgroundValue(), c);
 153     if (c.size() > 1) {
 154       double most_at_right = c[1][0];
 155       int most_at_right_index=1;
 156       for(uint j=1; j<c.size(); j++) {
 157         if (c[j][0] > most_at_right) {
 158           most_at_right = c[j][0];
 159           most_at_right_index = j;
 160         }
 161       }
 162       // Put all other CCL to Background
 163       slices_rightbronchus[i] =
 164         clitk::Binarize<MaskSliceType>(slices_rightbronchus[i], most_at_right_index,
 165                                        most_at_right_index, GetBackgroundValue(), GetForegroundValue());
 166     } // end c.size
 167   }
 168
 169   // Joint slices
 170   m_LeftBronchus = clitk::JoinSlices<MaskImageType>(slices_leftbronchus, m_LeftBronchus, 2);
 171   m_RightBronchus = clitk::JoinSlices<MaskImageType>(slices_rightbronchus, m_RightBronchus, 2);
 172
 173   writeImage<MaskImageType>(m_LeftBronchus, "step-left.mhd");
 174   writeImage<MaskImageType>(m_RightBronchus, "step-right.mhd");
 175
 176   m_Working_Support =
 177     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, m_LeftBronchus, 2,
 178                                                        GetFuzzyThresholdForS7("Bronchi"), "RightTo",
 179                                                        false, 3, false);
 180   StopCurrentStep<MaskImageType>(m_Working_Support);
 181
 182
 183   // ----------------------------------------------------------------
 184   StartNewStep("[Station7] Limits with bronchus : LeftTo the right bronchus");
 185   m_Working_Support =
 186     clitk::SliceBySliceRelativePosition<MaskImageType>(m_Working_Support, m_RightBronchus, 2,
 187                                                        GetFuzzyThresholdForS7("Bronchi"), "LeftTo",
 188                                                        false, 3, false);
 189   StopCurrentStep<MaskImageType>(m_Working_Support);
 190
 191
 192   // ----------------------------------------------------------------
 193   StartNewStep("[Station7] Limits with LeftSuperiorPulmonaryVein");
 194   try {
 195     MaskImagePointer LeftSuperiorPulmonaryVein = GetAFDB()->template GetImage<MaskImageType>("LeftSuperiorPulmonaryVein");
 196     typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
 197     typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 198     sliceRelPosFilter->SetInput(m_Working_Support);
 199     sliceRelPosFilter->SetInputObject(LeftSuperiorPulmonaryVein);
 200     sliceRelPosFilter->SetDirection(2);
 201     sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS7("LeftSuperiorPulmonaryVein"));
 202     sliceRelPosFilter->AddOrientationTypeString("NotLeftTo");
 203     sliceRelPosFilter->AddOrientationTypeString("NotAntTo");
 204     sliceRelPosFilter->SetIntermediateSpacingFlag(true);
 205     sliceRelPosFilter->SetIntermediateSpacing(3);
 206     sliceRelPosFilter->SetUniqueConnectedComponentBySlice(false);
 207     sliceRelPosFilter->SetAutoCropFlag(false);
 208     sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 209     sliceRelPosFilter->Update();
 210     m_Working_Support = sliceRelPosFilter->GetOutput();
 211     StopCurrentStep<MaskImageType>(m_Working_Support);
 212   }
 213   catch (clitk::ExceptionObject e) {
 214     std::cout << "Not LeftSuperiorPulmonaryVein, skip" << std::endl;
 215   }
 216
 217   // ----------------------------------------------------------------
 218   StartNewStep("[Station7] Limits with RightSuperiorPulmonaryVein");
 219   try {
 220     MaskImagePointer RightSuperiorPulmonaryVein = GetAFDB()->template GetImage<MaskImageType>("RightSuperiorPulmonaryVein");
 221     typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
 222     typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 223     sliceRelPosFilter->SetInput(m_Working_Support);
 224     sliceRelPosFilter->SetInputObject(RightSuperiorPulmonaryVein);
 225     sliceRelPosFilter->SetDirection(2);
 226     sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS7("RightSuperiorPulmonaryVein"));
 227     sliceRelPosFilter->AddOrientationTypeString("NotRightTo");
 228     sliceRelPosFilter->AddOrientationTypeString("NotAntTo");
 229     sliceRelPosFilter->AddOrientationTypeString("NotPostTo");
 230     sliceRelPosFilter->SetIntermediateSpacingFlag(true);
 231     sliceRelPosFilter->SetIntermediateSpacing(3);
 232     sliceRelPosFilter->SetUniqueConnectedComponentBySlice(false);
 233     sliceRelPosFilter->SetAutoCropFlag(false);
 234     sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 235     sliceRelPosFilter->Update();
 236     m_Working_Support = sliceRelPosFilter->GetOutput();
 237     StopCurrentStep<MaskImageType>(m_Working_Support);
 238   }
 239   catch (clitk::ExceptionObject e) {
 240     std::cout << "Not RightSuperiorPulmonaryVein, skip" << std::endl;
 241   }
 242
 243   // ----------------------------------------------------------------
 244   StartNewStep("[Station7] Limits with RightPulmonaryArtery");
 245   MaskImagePointer RightPulmonaryArtery = GetAFDB()->template GetImage<MaskImageType>("RightPulmonaryArtery");
 246   typedef SliceBySliceRelativePositionFilter<MaskImageType> SliceRelPosFilterType;
 247   typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New();
 248   sliceRelPosFilter->SetInput(m_Working_Support);
 249   sliceRelPosFilter->SetInputObject(RightPulmonaryArtery);
 250   sliceRelPosFilter->SetDirection(2);
 251   sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS7("RightPulmonaryArtery"));
 252   sliceRelPosFilter->AddOrientationTypeString("NotAntTo");
 253   sliceRelPosFilter->SetIntermediateSpacingFlag(true);
 254   sliceRelPosFilter->SetIntermediateSpacing(3);
 255   sliceRelPosFilter->SetUniqueConnectedComponentBySlice(false);
 256   sliceRelPosFilter->SetAutoCropFlag(false);
 257   sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 258   sliceRelPosFilter->Update();
 259   m_Working_Support = sliceRelPosFilter->GetOutput();
 260   StopCurrentStep<MaskImageType>(m_Working_Support);
 261
 262   // ----------------------------------------------------------------
 263   StartNewStep("[Station7] Limits with LeftPulmonaryArtery");
 264   MaskImagePointer LeftPulmonaryArtery = GetAFDB()->template GetImage<MaskImageType>("LeftPulmonaryArtery");
 265   sliceRelPosFilter = SliceRelPosFilterType::New();
 266   sliceRelPosFilter->SetInput(m_Working_Support);
 267   sliceRelPosFilter->SetInputObject(LeftPulmonaryArtery);
 268   sliceRelPosFilter->SetDirection(2);
 269   sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS7("LeftPulmonaryArtery"));
 270   sliceRelPosFilter->AddOrientationTypeString("NotAntTo");
 271   sliceRelPosFilter->SetIntermediateSpacingFlag(true);
 272   sliceRelPosFilter->SetIntermediateSpacing(3);
 273   sliceRelPosFilter->SetUniqueConnectedComponentBySlice(false);
 274   sliceRelPosFilter->SetAutoCropFlag(false);
 275   sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 276   sliceRelPosFilter->Update();
 277   m_Working_Support = sliceRelPosFilter->GetOutput();
 278   StopCurrentStep<MaskImageType>(m_Working_Support);
 279
 280   StartNewStep("[Station7] Limits with SVC");
 281   MaskImagePointer SVC = GetAFDB()->template GetImage<MaskImageType>("SVC");
 282   sliceRelPosFilter = SliceRelPosFilterType::New();
 283   sliceRelPosFilter->SetInput(m_Working_Support);
 284   sliceRelPosFilter->SetInputObject(SVC);
 285   sliceRelPosFilter->SetDirection(2);
 286   sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS7("SVC"));
 287   sliceRelPosFilter->AddOrientationTypeString("NotRightTo");
 288   sliceRelPosFilter->AddOrientationTypeString("NotAntTo");
 289   sliceRelPosFilter->SetIntermediateSpacingFlag(true);
 290   sliceRelPosFilter->SetIntermediateSpacing(3);
 291   sliceRelPosFilter->SetUniqueConnectedComponentBySlice(false);
 292   sliceRelPosFilter->SetAutoCropFlag(true);
 293   sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn();
 294   sliceRelPosFilter->Update();
 295   m_Working_Support = sliceRelPosFilter->GetOutput();
 296   StopCurrentStep<MaskImageType>(m_Working_Support);
 297
 298   // End
 299   m_ListOfStations["7"] = m_Working_Support;
 300 }
 301 //--------------------------------------------------------------------
 302
 303
 304 //--------------------------------------------------------------------
 305 template <class TImageType>
 306 void
 307 clitk::ExtractLymphStationsFilter<TImageType>::
 308 ExtractStation_7_Posterior_Limits()
 309 {
 310   StartNewStep("[Station7] Posterior limits -> must be AntTo post wall of the bronchi (OLD CLASSIF)");
 311
 312   // Search for points that are the most left/post/ant and most
 313   // right/post/ant of the left and right bronchus
 314
 315   // extract, loop slices, label/keep, find extrema x 3
 316   /*  FindExtremaPointsInBronchus(m_LeftBronchus, 0, 15, m_RightMostInLeftBronchus,
 317                               m_AntMostInLeftBronchus, m_PostMostInLeftBronchus);
 318   FindExtremaPointsInBronchus(m_RightBronchus, 1, 15, m_LeftMostInRightBronchus,
 319                               m_AntMostInRightBronchus, m_PostMostInRightBronchus);
 320   */
 321
 322   // First cut bronchus to the correct sup/inf support
 323   MaskImagePointer RightBronchus = clitk::ResizeImageLike<MaskImageType>(m_RightBronchus, m_Working_Support, GetBackgroundValue());
 324   MaskImagePointer LeftBronchus = clitk::ResizeImageLike<MaskImageType>(m_LeftBronchus, m_Working_Support, GetBackgroundValue());
 325
 326   // Find extrema points
 327   FindExtremaPointsInBronchus(RightBronchus, 0, 10, m_LeftMostInRightBronchus,
 328                               m_AntMostInRightBronchus, m_PostMostInRightBronchus);
 329
 330   FindExtremaPointsInBronchus(LeftBronchus, 1, 10, m_RightMostInLeftBronchus,
 331                               m_AntMostInLeftBronchus, m_PostMostInLeftBronchus);
 332
 333
 334
 335   // DEBUG
 336   std::ofstream osrl; openFileForWriting(osrl, "osrl.txt"); osrl << "LANDMARKS1" << std::endl;
 337   std::ofstream osal; openFileForWriting(osal, "osal.txt"); osal << "LANDMARKS1" << std::endl;
 338   std::ofstream ospl; openFileForWriting(ospl, "ospl.txt"); ospl << "LANDMARKS1" << std::endl;
 339   std::ofstream osrr; openFileForWriting(osrr, "osrr.txt"); osrr << "LANDMARKS1" << std::endl;
 340   std::ofstream osar; openFileForWriting(osar, "osar.txt"); osar << "LANDMARKS1" << std::endl;
 341   std::ofstream ospr; openFileForWriting(ospr, "ospr.txt"); ospr << "LANDMARKS1" << std::endl;
 342
 343   for(uint i=0; i<m_RightMostInLeftBronchus.size(); i++) {
 344     osrl << i << " " << m_RightMostInLeftBronchus[i][0] << " " << m_RightMostInLeftBronchus[i][1]
 345          << " " << m_RightMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
 346     osal << i << " " << m_AntMostInLeftBronchus[i][0] << " " << m_AntMostInLeftBronchus[i][1]
 347          << " " << m_AntMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
 348     ospl << i << " " << m_PostMostInLeftBronchus[i][0] << " " << m_PostMostInLeftBronchus[i][1]
 349          << " " << m_PostMostInLeftBronchus[i][2] << " 0 0 " << std::endl;
 350   }
 351
 352   for(uint i=0; i<m_LeftMostInRightBronchus.size(); i++) {
 353     osrr << i << " " << m_LeftMostInRightBronchus[i][0] << " " << m_LeftMostInRightBronchus[i][1]
 354          << " " << m_LeftMostInRightBronchus[i][2] << " 0 0 " << std::endl;
 355     osar << i << " " << m_AntMostInRightBronchus[i][0] << " " << m_AntMostInRightBronchus[i][1]
 356          << " " << m_AntMostInRightBronchus[i][2] << " 0 0 " << std::endl;
 357     ospr << i << " " << m_PostMostInRightBronchus[i][0] << " " << m_PostMostInRightBronchus[i][1]
 358          << " " << m_PostMostInRightBronchus[i][2] << " 0 0 " << std::endl;
 359   }
 360   osrl.close();
 361   osal.close();
 362   ospl.close();
 363   osrr.close();
 364   osar.close();
 365   ospr.close();
 366
 367   clitk::SliceBySliceSetBackgroundFromLineSeparation<MaskImageType>(m_Working_Support,
 368                                                                     m_PostMostInRightBronchus,
 369                                                                     m_PostMostInLeftBronchus,
 370                                                                     GetBackgroundValue(), 1, -10);
 371   // If needed -> can do the same with AntMost.
 372
 373   // End
 374   StopCurrentStep<MaskImageType>(m_Working_Support);
 375   m_ListOfStations["7"] = m_Working_Support;
 376 }
 377 //--------------------------------------------------------------------
 378
 379 //--------------------------------------------------------------------
 380 template <class ImageType>
 381 void
 382 clitk::ExtractLymphStationsFilter<ImageType>::
 383 ExtractStation_7_Remove_Structures()
 384 {
 385
 386   //--------------------------------------------------------------------
 387   StartNewStep("[Station7] remove some structures");
 388
 389   Remove_Structures("AzygousVein");
 390   Remove_Structures("Aorta");
 391   Remove_Structures("Esophagus");
 392   Remove_Structures("RightPulmonaryArtery");
 393   Remove_Structures("LeftPulmonaryArtery");
 394   Remove_Structures("LeftSuperiorPulmonaryVein");
 395   Remove_Structures("PulmonaryTrunk");
 396   Remove_Structures("VertebralBody");
 397
 398   // END
 399   StopCurrentStep<MaskImageType>(m_Working_Support);
 400   m_ListOfStations["7"] = m_Working_Support;
 401   return;
 402 }
 403 //--------------------------------------------------------------------
 404
 405
 406