segmentation/clitkExtractLymphStationsFilter.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://oncora1.lyon.fnclcc.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 #ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
  20 #define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX
  21
  22 // clitk
  23 #include "clitkCommon.h"
  24 #include "clitkExtractLymphStationsFilter.h"
  25 #include "clitkAddRelativePositionConstraintToLabelImageFilter.h"
  26 #include "clitkSegmentationUtils.h"
  27 #include "clitkAutoCropFilter.h"
  28 #include "clitkSegmentationUtils.h"
  29 #include "clitkSliceBySliceRelativePositionFilter.h"
  30
  31 // itk
  32 #include <itkStatisticsLabelMapFilter.h>
  33 #include <itkLabelImageToStatisticsLabelMapFilter.h>
  34 #include <itkRegionOfInterestImageFilter.h>
  35 #include <itkBinaryThresholdImageFilter.h>
  36 #include <itkImageSliceConstIteratorWithIndex.h>
  37 #include <itkImageSliceIteratorWithIndex.h>
  38 #include <itkBinaryThinningImageFilter.h>
  39
  40 // itk ENST
  41 #include "RelativePositionPropImageFilter.h"
  42
  43 //--------------------------------------------------------------------
  44 template <class TImageType>
  45 clitk::ExtractLymphStationsFilter<TImageType>::
  46 ExtractLymphStationsFilter():
  47   clitk::FilterBase(),
  48   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  49   itk::ImageToImageFilter<TImageType, MaskImageType>()
  50 {
  51   this->SetNumberOfRequiredInputs(1);
  52   SetBackgroundValue(0);
  53   SetForegroundValue(1);
  54
  55   // Default values
  56   ExtractStation_8_SetDefaultValues();
  57   ExtractStation_3P_SetDefaultValues();
  58   ExtractStation_2RL_SetDefaultValues();
  59   ExtractStation_3A_SetDefaultValues();
  60   ExtractStation_7_SetDefaultValues();
  61 }
  62 //--------------------------------------------------------------------
  63
  64
  65 //--------------------------------------------------------------------
  66 template <class TImageType>
  67 void
  68 clitk::ExtractLymphStationsFilter<TImageType>::
  69 GenerateOutputInformation() {
  70   // Get inputs
  71   LoadAFDB();
  72   m_Input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
  73   m_Mediastinum = GetAFDB()->template GetImage <MaskImageType>("Mediastinum");
  74
  75   // Global supports for stations
  76   StartNewStep("Supports for stations");
  77   StartSubStep();
  78   ExtractStationSupports();
  79   StopSubStep();
  80
  81   // Extract Station8
  82   StartNewStep("Station 8");
  83   StartSubStep();
  84   ExtractStation_8();
  85   StopSubStep();
  86
  87   // Extract Station3P
  88   StartNewStep("Station 3P");
  89   StartSubStep();
  90   ExtractStation_3P();
  91   StopSubStep();
  92
  93   // Extract Station2RL
  94   /*
  95     StartNewStep("Station 2RL");
  96     StartSubStep();
  97     ExtractStation_2RL();
  98     StopSubStep();
  99   */
 100
 101   // Extract Station3A
 102   StartNewStep("Station 3A");
 103   StartSubStep();
 104   ExtractStation_3A();
 105   StopSubStep();
 106
 107   // Extract Station7
 108   StartNewStep("Station 7");
 109   StartSubStep();
 110   ExtractStation_7();
 111   StopSubStep();
 112
 113   if (0) { // temporary suppress
 114     // Extract Station4RL
 115     StartNewStep("Station 4RL");
 116     StartSubStep();
 117     //ExtractStation_4RL();
 118     StopSubStep();
 119   }
 120 }
 121 //--------------------------------------------------------------------
 122
 123
 124 //--------------------------------------------------------------------
 125 template <class TImageType>
 126 void
 127 clitk::ExtractLymphStationsFilter<TImageType>::
 128 GenerateInputRequestedRegion() {
 129   //DD("GenerateInputRequestedRegion (nothing?)");
 130 }
 131 //--------------------------------------------------------------------
 132
 133
 134 //--------------------------------------------------------------------
 135 template <class TImageType>
 136 void
 137 clitk::ExtractLymphStationsFilter<TImageType>::
 138 GenerateData() {
 139   // Final Step -> graft output (if SetNthOutput => redo)
 140   this->GraftOutput(m_ListOfStations["8"]);
 141 }
 142 //--------------------------------------------------------------------
 143
 144
 145 //--------------------------------------------------------------------
 146 template <class TImageType>
 147 bool
 148 clitk::ExtractLymphStationsFilter<TImageType>::
 149 CheckForStation(std::string station)
 150 {
 151   // Compute Station name
 152   std::string s = "Station"+station;
 153
 154
 155   // Check if station already exist in DB
 156   bool found = false;
 157   if (GetAFDB()->TagExist(s)) {
 158     m_ListOfStations[station] = GetAFDB()->template GetImage<MaskImageType>(s);
 159     found = true;
 160   }
 161
 162   // Define the starting support
 163   if (found && GetComputeStation(station)) {
 164     std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl;
 165   }
 166   if (!found || GetComputeStation(station)) {
 167     m_Working_Support = m_Mediastinum = GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
 168     return true;
 169   }
 170   else {
 171     std::cout << "Station " << station << " found. I used it" << std::endl;
 172     return false;
 173   }
 174 }
 175 //--------------------------------------------------------------------
 176
 177
 178 //--------------------------------------------------------------------
 179 template <class TImageType>
 180 bool
 181 clitk::ExtractLymphStationsFilter<TImageType>::
 182 GetComputeStation(std::string station)
 183 {
 184   return (m_ComputeStationMap.find(station) != m_ComputeStationMap.end());
 185 }
 186 //--------------------------------------------------------------------
 187
 188
 189 //--------------------------------------------------------------------
 190 template <class TImageType>
 191 void
 192 clitk::ExtractLymphStationsFilter<TImageType>::
 193 AddComputeStation(std::string station)
 194 {
 195   m_ComputeStationMap[station] = true;
 196 }
 197 //--------------------------------------------------------------------
 198
 199
 200 //--------------------------------------------------------------------
 201 template <class TImageType>
 202 void
 203 clitk::ExtractLymphStationsFilter<TImageType>::
 204 ExtractStation_3P()
 205 {
 206   if (CheckForStation("3P")) {
 207     ExtractStation_3P_SI_Limits();
 208     ExtractStation_3P_Ant_Limits();
 209     ExtractStation_3P_Post_Limits();
 210     ExtractStation_3P_LR_sup_Limits();
 211     //    ExtractStation_3P_LR_sup_Limits_2();
 212     ExtractStation_3P_LR_inf_Limits();
 213     ExtractStation_8_Single_CCL_Limits(); // YES 8 !
 214     ExtractStation_3P_Remove_Structures(); // after CCL
 215     // Store image filenames into AFDB
 216     writeImage<MaskImageType>(m_ListOfStations["3P"], "seg/Station3P.mhd");
 217     GetAFDB()->SetImageFilename("Station3P", "seg/Station3P.mhd");
 218     WriteAFDB();
 219   }
 220 }
 221 //--------------------------------------------------------------------
 222
 223
 224 //--------------------------------------------------------------------
 225 template <class TImageType>
 226 void
 227 clitk::ExtractLymphStationsFilter<TImageType>::
 228 ExtractStation_3A()
 229 {
 230   if (CheckForStation("3A")) {
 231     ExtractStation_3A_SI_Limits();
 232     ExtractStation_3A_Ant_Limits();
 233     ExtractStation_3A_Post_Limits();
 234     // Store image filenames into AFDB
 235     writeImage<MaskImageType>(m_ListOfStations["3A"], "seg/Station3A.mhd");
 236     GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd");
 237     WriteAFDB();
 238   }
 239 }
 240 //--------------------------------------------------------------------
 241
 242
 243 //--------------------------------------------------------------------
 244 template <class TImageType>
 245 void
 246 clitk::ExtractLymphStationsFilter<TImageType>::
 247 ExtractStation_2RL()
 248 {
 249   if (CheckForStation("2RL")) {
 250     ExtractStation_2RL_SI_Limits();
 251     ExtractStation_2RL_Post_Limits();
 252     ExtractStation_2RL_Ant_Limits2();
 253     ExtractStation_2RL_Ant_Limits();
 254     ExtractStation_2RL_LR_Limits();
 255     ExtractStation_2RL_Remove_Structures();
 256     ExtractStation_2RL_SeparateRL();
 257
 258     // Store image filenames into AFDB
 259     writeImage<MaskImageType>(m_ListOfStations["2R"], "seg/Station2R.mhd");
 260     writeImage<MaskImageType>(m_ListOfStations["2L"], "seg/Station2L.mhd");
 261     GetAFDB()->SetImageFilename("Station2R", "seg/Station2R.mhd");
 262     GetAFDB()->SetImageFilename("Station2L", "seg/Station2L.mhd");
 263     WriteAFDB();
 264   }
 265 }
 266 //--------------------------------------------------------------------
 267
 268
 269 //--------------------------------------------------------------------
 270 template <class TImageType>
 271 void
 272 clitk::ExtractLymphStationsFilter<TImageType>::
 273 ExtractStation_4RL() {
 274   DD("TODO");
 275   exit(0);
 276   /*
 277     WARNING ONLY 4R FIRST !!! (not same inf limits)
 278   */
 279   ExtractStation_4RL_SI_Limits();
 280   ExtractStation_4RL_LR_Limits();
 281   ExtractStation_4RL_AP_Limits();
 282 }
 283 //--------------------------------------------------------------------
 284
 285
 286 //--------------------------------------------------------------------
 287 template <class ImageType>
 288 void
 289 clitk::ExtractLymphStationsFilter<ImageType>::
 290 Remove_Structures(std::string station, std::string s)
 291 {
 292   try {
 293     StartNewStep("[Station"+station+"] Remove "+s);
 294     MaskImagePointer Structure = GetAFDB()->template GetImage<MaskImageType>(s);
 295     clitk::AndNot<MaskImageType>(m_Working_Support, Structure, GetBackgroundValue());
 296   }
 297   catch(clitk::ExceptionObject e) {
 298     std::cout << s << " not found, skip." << std::endl;
 299   }
 300 }
 301 //--------------------------------------------------------------------
 302
 303
 304 //--------------------------------------------------------------------
 305 template <class TImageType>
 306 void
 307 clitk::ExtractLymphStationsFilter<TImageType>::
 308 SetFuzzyThreshold(std::string station, std::string tag, double value)
 309 {
 310   m_FuzzyThreshold[station][tag] = value;
 311 }
 312 //--------------------------------------------------------------------
 313
 314
 315 //--------------------------------------------------------------------
 316 template <class TImageType>
 317 double
 318 clitk::ExtractLymphStationsFilter<TImageType>::
 319 GetFuzzyThreshold(std::string station, std::string tag)
 320 {
 321   if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
 322     clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
 323     return 0.0;
 324   }
 325
 326   if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
 327     clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
 328     return 0.0;
 329   }
 330
 331   return m_FuzzyThreshold[station][tag];
 332 }
 333 //--------------------------------------------------------------------
 334
 335
 336 //--------------------------------------------------------------------
 337 template <class TImageType>
 338 void
 339 clitk::ExtractLymphStationsFilter<TImageType>::
 340 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
 341 {
 342   // Create line from A to B with
 343   // A = upper border of LLL at left
 344   // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
 345
 346   try {
 347     GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
 348     GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
 349   }
 350   catch(clitk::ExceptionObject & o) {
 351
 352     DD("FindLineForS7S8Separation");
 353     // Load LeftLowerLobeBronchus and get centroid point
 354     MaskImagePointer LeftLowerLobeBronchus =
 355       GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
 356     std::vector<MaskImagePointType> c;
 357     clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
 358     A = c[1];
 359
 360     // Load RightMiddleLobeBronchus and get superior point (not centroid here)
 361     MaskImagePointer RightMiddleLobeBronchus =
 362       GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
 363     bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
 364                                                               GetBackgroundValue(),
 365                                                               2, false, B);
 366     if (!b) {
 367       clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
 368     }
 369
 370     // Insert into the DB
 371     GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
 372     GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
 373   }
 374 }
 375 //--------------------------------------------------------------------
 376
 377
 378 //--------------------------------------------------------------------
 379 template <class TImageType>
 380 double
 381 clitk::ExtractLymphStationsFilter<TImageType>::
 382 FindCarinaSlicePosition()
 383 {
 384   double z;
 385   try {
 386     z = GetAFDB()->GetDouble("CarinaZ");
 387   }
 388   catch(clitk::ExceptionObject e) {
 389     DD("FindCarinaSlicePosition");
 390     // Get Carina
 391     MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
 392
 393     // Get Centroid and Z value
 394     std::vector<MaskImagePointType> centroids;
 395     clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
 396
 397     // We dont need Carina structure from now
 398     Carina->Delete();
 399
 400     // Put inside the AFDB
 401     GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
 402     GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
 403     WriteAFDB();
 404     z = centroids[1][2];
 405   }
 406   return z;
 407 }
 408 //--------------------------------------------------------------------
 409
 410
 411 //--------------------------------------------------------------------
 412 template <class TImageType>
 413 void
 414 clitk::ExtractLymphStationsFilter<TImageType>::
 415 FindLeftAndRightBronchi()
 416 {
 417   try {
 418     m_RightBronchus = GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
 419     m_LeftBronchus = GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
 420   }
 421   catch(clitk::ExceptionObject & o) {
 422
 423     DD("FindLeftAndRightBronchi");
 424     // The goal is to separate the trachea inferiorly to the carina into
 425     // a Left and Right bronchus.
 426
 427     // Get the trachea
 428     MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
 429
 430     // Get the Carina position
 431     m_CarinaZ = FindCarinaSlicePosition();
 432
 433     // Consider only inferiorly to the Carina
 434     MaskImagePointer m_Working_Trachea =
 435       clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
 436                                                        GetBackgroundValue());
 437
 438     // Labelize the trachea
 439     m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
 440
 441     // Carina position must at the first slice that separate the two
 442     // main bronchus (not superiorly). We thus first check that the
 443     // upper slice is composed of at least two labels
 444     MaskImagePointer RightBronchus;
 445     MaskImagePointer LeftBronchus;
 446     typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
 447     SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
 448     iter.SetFirstDirection(0);
 449     iter.SetSecondDirection(1);
 450     iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
 451     int maxLabel=0;
 452     while (!iter.IsAtReverseEndOfSlice()) {
 453       while (!iter.IsAtReverseEndOfLine()) {
 454         if (iter.Get() > maxLabel) maxLabel = iter.Get();
 455         --iter;
 456       }
 457       iter.PreviousLine();
 458     }
 459     if (maxLabel < 2) {
 460       clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
 461     }
 462
 463     // Compute 3D centroids of both parts to identify the left from the
 464     // right bronchus
 465     std::vector<ImagePointType> c;
 466     clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
 467     ImagePointType C1 = c[1];
 468     ImagePointType C2 = c[2];
 469
 470     ImagePixelType rightLabel;
 471     ImagePixelType leftLabel;
 472     if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
 473     else { rightLabel = 2; leftLabel = 1; }
 474
 475     // Select LeftLabel (set one label to Backgroundvalue)
 476     RightBronchus =
 477       clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
 478                                      GetBackgroundValue(), GetForegroundValue());
 479     /*
 480       SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
 481       leftLabel, GetBackgroundValue(), false);
 482     */
 483     LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
 484                                                   GetBackgroundValue(), GetForegroundValue());
 485     /*
 486       SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
 487       rightLabel, GetBackgroundValue(), false);
 488     */
 489
 490     // Crop images
 491     RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
 492     LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
 493
 494     // Insert int AFDB if need after
 495     GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
 496                                                  RightBronchus, true);
 497     GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
 498                                                  LeftBronchus, true);
 499   }
 500 }
 501 //--------------------------------------------------------------------
 502
 503 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX