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 #include "clitkMeshToBinaryImageFilter.h"
  31
  32 // itk
  33 #include <itkStatisticsLabelMapFilter.h>
  34 #include <itkLabelImageToStatisticsLabelMapFilter.h>
  35 #include <itkRegionOfInterestImageFilter.h>
  36 #include <itkBinaryThresholdImageFilter.h>
  37 #include <itkImageSliceConstIteratorWithIndex.h>
  38 #include <itkImageSliceIteratorWithIndex.h>
  39 #include <itkBinaryThinningImageFilter.h>
  40
  41 // itk ENST
  42 #include "RelativePositionPropImageFilter.h"
  43
  44 // vtk
  45 #include <vtkAppendPolyData.h>
  46 #include <vtkPolyDataWriter.h>
  47 #include <vtkCellArray.h>
  48
  49 //--------------------------------------------------------------------
  50 template <class TImageType>
  51 clitk::ExtractLymphStationsFilter<TImageType>::
  52 ExtractLymphStationsFilter():
  53   clitk::StructuresExtractionFilter<ImageType>()
  54 {
  55   this->SetNumberOfRequiredInputs(1);
  56   SetBackgroundValue(0);
  57   SetForegroundValue(1);
  58   ForceSupportsFlagOn();
  59   SetSupportLimitsFilename("none");
  60   CheckSupportFlagOff();
  61
  62   // Default values
  63   ExtractStation_3P_SetDefaultValues();
  64   ExtractStation_2RL_SetDefaultValues();
  65   ExtractStation_3A_SetDefaultValues();
  66   ExtractStation_1RL_SetDefaultValues();
  67   ExtractStation_4RL_SetDefaultValues();
  68   ExtractStation_5_SetDefaultValues();
  69   ExtractStation_6_SetDefaultValues();
  70
  71   // TODO
  72   ExtractStation_7_SetDefaultValues();
  73   ExtractStation_8_SetDefaultValues();
  74 }
  75 //--------------------------------------------------------------------
  76
  77
  78 //--------------------------------------------------------------------
  79 template <class TImageType>
  80 void
  81 clitk::ExtractLymphStationsFilter<TImageType>::
  82 GenerateOutputInformation() {
  83   // Get inputs
  84   this->LoadAFDB();
  85   m_Input = dynamic_cast<const ImageType*>(itk::ProcessObject::GetInput(0));
  86   m_Mediastinum = this->GetAFDB()->template GetImage <MaskImageType>("Mediastinum");
  87
  88   // Clean some computer landmarks to force the recomputation
  89   // FIXME -> to put elsewhere ?
  90   this->GetAFDB()->RemoveTag("AntPostVesselsSeparation");
  91
  92   // Must I compute the supports ?
  93   bool supportsExist = true;
  94   if (!GetForceSupportsFlag()) {
  95     try {
  96       m_ListOfSupports["S1R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1R");
  97       m_ListOfSupports["S1L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S1L");
  98       m_ListOfSupports["S2R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2R");
  99       m_ListOfSupports["S2L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L");
 100       m_ListOfSupports["S4R"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4R");
 101       m_ListOfSupports["S4L"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S4L");
 102
 103       m_ListOfSupports["S3A"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A");
 104       m_ListOfSupports["S3P"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S3P");
 105       m_ListOfSupports["S5"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S5");
 106       m_ListOfSupports["S6"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S6");
 107       m_ListOfSupports["S7"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S7");
 108       m_ListOfSupports["S8"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S8");
 109       m_ListOfSupports["S9"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S9");
 110       m_ListOfSupports["S10"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S10");
 111       m_ListOfSupports["S11"] = this->GetAFDB()->template GetImage<MaskImageType>("Support_S11");
 112     } catch(clitk::ExceptionObject o) {
 113       supportsExist = false;
 114     }
 115   }
 116
 117   if (!supportsExist || GetForceSupportsFlag()) {
 118     this->StartNewStep("Supports for stations");
 119     this->StartSubStep();
 120
 121     // FIXME : why should I remove theses tags ???
 122     this->GetAFDB()->RemoveTag("CarinaZ");
 123     this->GetAFDB()->RemoveTag("ApexOfTheChestZ");
 124     this->GetAFDB()->RemoveTag("ApexOfTheChest");
 125     this->GetAFDB()->RemoveTag("RightBronchus");
 126     this->GetAFDB()->RemoveTag("LeftBronchus");
 127     this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArchZ");
 128     this->GetAFDB()->RemoveTag("SuperiorBorderOfAorticArch");
 129     this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArchZ");
 130     this->GetAFDB()->RemoveTag("InferiorBorderOfAorticArch");
 131     ExtractStationSupports();
 132     this->StopSubStep();
 133   }
 134
 135   // Extract Stations
 136   ExtractStation_1RL();
 137   ExtractStation_2RL();
 138   ExtractStation_3P();
 139   ExtractStation_3A();
 140   ExtractStation_4R();
 141   ExtractStation_4L();
 142   ExtractStation_5();
 143   ExtractStation_6();
 144
 145   // ---------- todo -----------------------
 146
 147   // Extract Station8
 148   //  ExtractStation_8();
 149
 150   // Extract Station7
 151   //this->StartNewStep("Station 7");
 152   //this->StartSubStep();
 153   //ExtractStation_7();
 154   //this->StopSubStep();
 155
 156 }
 157 //--------------------------------------------------------------------
 158
 159
 160 //--------------------------------------------------------------------
 161 template <class TImageType>
 162 void
 163 clitk::ExtractLymphStationsFilter<TImageType>::
 164 GenerateInputRequestedRegion() {
 165   //DD("GenerateInputRequestedRegion (nothing?)");
 166 }
 167 //--------------------------------------------------------------------
 168
 169
 170 //--------------------------------------------------------------------
 171 template <class TImageType>
 172 void
 173 clitk::ExtractLymphStationsFilter<TImageType>::
 174 GenerateData() {
 175   // Final Step -> graft output (if SetNthOutput => redo)
 176   // this->GraftOutput(m_ListOfStations["8"]);
 177 }
 178 //--------------------------------------------------------------------
 179
 180
 181 //--------------------------------------------------------------------
 182 template <class TImageType>
 183 bool
 184 clitk::ExtractLymphStationsFilter<TImageType>::
 185 CheckForStation(std::string station)
 186 {
 187   // Compute Station name
 188   std::string s = "Station"+station;
 189
 190
 191   // Define the starting support
 192   // if (GetComputeStation(station)) {
 193   //   std::cout << "Station " << station << " already exists, but re-computation forced." << std::endl;
 194   // }
 195   if (GetComputeStation(station)) {
 196     m_Working_Support = m_Mediastinum = this->GetAFDB()->template GetImage<MaskImageType>("Mediastinum", true);
 197     return true;
 198   }
 199   else return false;
 200   // else {
 201   //   std::cout << "Station " << station << " found. I used it" << std::endl;
 202   //   return false;
 203   // }
 204
 205   // Check if station already exist in DB
 206
 207   // FIXME -> do nothing if not on the command line. Is it what I want ?
 208   //bool found = false;
 209   if (this->GetAFDB()->TagExist(s)) {
 210     m_ListOfStations[station] = this->GetAFDB()->template GetImage<MaskImageType>(s);
 211     //found = true;
 212   }
 213
 214 }
 215 //--------------------------------------------------------------------
 216
 217
 218 //--------------------------------------------------------------------
 219 template <class TImageType>
 220 bool
 221 clitk::ExtractLymphStationsFilter<TImageType>::
 222 GetComputeStation(std::string station)
 223 {
 224   return (m_ComputeStationMap.find(station) != m_ComputeStationMap.end());
 225 }
 226 //--------------------------------------------------------------------
 227
 228
 229 //--------------------------------------------------------------------
 230 template <class TImageType>
 231 void
 232 clitk::ExtractLymphStationsFilter<TImageType>::
 233 AddComputeStation(std::string station)
 234 {
 235   m_ComputeStationMap[station] = true;
 236 }
 237 //--------------------------------------------------------------------
 238
 239
 240
 241 //--------------------------------------------------------------------
 242 template<class PointType>
 243 class comparePointsX {
 244 public:
 245   bool operator() (PointType i, PointType j) { return (i[0]<j[0]); }
 246 };
 247 //--------------------------------------------------------------------
 248
 249
 250 //--------------------------------------------------------------------
 251 template<class PairType>
 252 class comparePointsWithAngle {
 253 public:
 254   bool operator() (PairType i, PairType j) { return (i.second < j.second); }
 255 };
 256 //--------------------------------------------------------------------
 257
 258
 259 //--------------------------------------------------------------------
 260 template<int Dim>
 261 void HypercubeCorners(std::vector<itk::Point<double, Dim> > & out) {
 262   std::vector<itk::Point<double, Dim-1> > previous;
 263   HypercubeCorners<Dim-1>(previous);
 264   out.resize(previous.size()*2);
 265   for(unsigned int i=0; i<out.size(); i++) {
 266     itk::Point<double, Dim> p;
 267     if (i<previous.size()) p[0] = 0;
 268     else p[0] = 1;
 269     for(int j=0; j<Dim-1; j++)
 270       {
 271         p[j+1] = previous[i%previous.size()][j];
 272       }
 273     out[i] = p;
 274   }
 275 }
 276 //--------------------------------------------------------------------
 277
 278
 279 //--------------------------------------------------------------------
 280 template<>
 281 void HypercubeCorners<1>(std::vector<itk::Point<double, 1> > & out) {
 282   out.resize(2);
 283   out[0][0] = 0;
 284   out[1][0] = 1;
 285 }
 286 //--------------------------------------------------------------------
 287
 288
 289 //--------------------------------------------------------------------
 290 template<class ImageType>
 291 void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image,
 292                                        std::vector<typename ImageType::PointType> & bounds)
 293 {
 294   // Get image max/min coordinates
 295   const unsigned int dim=ImageType::ImageDimension;
 296   typedef typename ImageType::PointType PointType;
 297   typedef typename ImageType::IndexType IndexType;
 298   PointType min_c, max_c;
 299   IndexType min_i, max_i;
 300   min_i = image->GetLargestPossibleRegion().GetIndex();
 301   for(unsigned int i=0; i<dim; i++)
 302     max_i[i] = image->GetLargestPossibleRegion().GetSize()[i] + min_i[i];
 303   image->TransformIndexToPhysicalPoint(min_i, min_c);
 304   image->TransformIndexToPhysicalPoint(max_i, max_c);
 305
 306   // Get corners coordinates
 307   HypercubeCorners<ImageType::ImageDimension>(bounds);
 308   for(unsigned int i=0; i<bounds.size(); i++) {
 309     for(unsigned int j=0; j<dim; j++) {
 310       if (bounds[i][j] == 0) bounds[i][j] = min_c[j];
 311       if (bounds[i][j] == 1) bounds[i][j] = max_c[j];
 312     }
 313   }
 314 }
 315 //--------------------------------------------------------------------
 316
 317
 318 //--------------------------------------------------------------------
 319 template <class ImageType>
 320 void
 321 clitk::ExtractLymphStationsFilter<ImageType>::
 322 Remove_Structures(std::string station, std::string s)
 323 {
 324   try {
 325     this->StartNewStep("[Station"+station+"] Remove "+s);
 326     MaskImagePointer Structure = this->GetAFDB()->template GetImage<MaskImageType>(s);
 327     clitk::AndNot<MaskImageType>(m_Working_Support, Structure, GetBackgroundValue());
 328   }
 329   catch(clitk::ExceptionObject e) {
 330     std::cout << s << " not found, skip." << std::endl;
 331   }
 332 }
 333 //--------------------------------------------------------------------
 334
 335
 336 //--------------------------------------------------------------------
 337 template <class TImageType>
 338 void
 339 clitk::ExtractLymphStationsFilter<TImageType>::
 340 SetFuzzyThreshold(std::string station, std::string tag, double value)
 341 {
 342   m_FuzzyThreshold[station][tag] = value;
 343 }
 344 //--------------------------------------------------------------------
 345
 346
 347 //--------------------------------------------------------------------
 348 template <class TImageType>
 349 void
 350 clitk::ExtractLymphStationsFilter<TImageType>::
 351 SetThreshold(std::string station, std::string tag, double value)
 352 {
 353   m_Threshold[station][tag] = value;
 354 }
 355 //--------------------------------------------------------------------
 356
 357
 358 //--------------------------------------------------------------------
 359 template <class TImageType>
 360 double
 361 clitk::ExtractLymphStationsFilter<TImageType>::
 362 GetFuzzyThreshold(std::string station, std::string tag)
 363 {
 364   if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
 365     clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
 366     return 0.0;
 367   }
 368
 369   if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
 370     clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
 371     return 0.0;
 372   }
 373
 374   return m_FuzzyThreshold[station][tag];
 375 }
 376 //--------------------------------------------------------------------
 377
 378
 379 //--------------------------------------------------------------------
 380 template <class TImageType>
 381 double
 382 clitk::ExtractLymphStationsFilter<TImageType>::
 383 GetThreshold(std::string station, std::string tag)
 384 {
 385   if (m_Threshold.find(station) == m_Threshold.end()) {
 386     clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
 387     return 0.0;
 388   }
 389
 390   if (m_Threshold[station].find(tag) == m_Threshold[station].end()) {
 391     clitkExceptionMacro("Could not find options "+tag+" in the list of Threshold for station "+station+".");
 392     return 0.0;
 393   }
 394
 395   return m_Threshold[station][tag];
 396 }
 397 //--------------------------------------------------------------------
 398
 399
 400 //--------------------------------------------------------------------
 401 template <class TImageType>
 402 void
 403 clitk::ExtractLymphStationsFilter<TImageType>::
 404 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
 405 {
 406   // Create line from A to B with
 407   // A = upper border of LLL at left
 408   // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
 409
 410   try {
 411     this->GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
 412     this->GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
 413   }
 414   catch(clitk::ExceptionObject & o) {
 415
 416     DD("FindLineForS7S8Separation");
 417     // Load LeftLowerLobeBronchus and get centroid point
 418     MaskImagePointer LeftLowerLobeBronchus =
 419       this->GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
 420     std::vector<MaskImagePointType> c;
 421     clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
 422     A = c[1];
 423
 424     // Load RightMiddleLobeBronchus and get superior point (not centroid here)
 425     MaskImagePointer RightMiddleLobeBronchus =
 426       this->GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
 427     bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
 428                                                               GetBackgroundValue(),
 429                                                               2, false, B);
 430     if (!b) {
 431       clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
 432     }
 433
 434     // Insert into the DB
 435     this->GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
 436     this->GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
 437   }
 438 }
 439 //--------------------------------------------------------------------
 440
 441
 442 //--------------------------------------------------------------------
 443 template <class TImageType>
 444 double
 445 clitk::ExtractLymphStationsFilter<TImageType>::
 446 FindCarina()
 447 {
 448   double z;
 449   try {
 450     z = this->GetAFDB()->GetDouble("CarinaZ");
 451   }
 452   catch(clitk::ExceptionObject e) {
 453     //DD("FindCarinaSlicePosition");
 454     // Get Carina
 455     MaskImagePointer Carina = this->GetAFDB()->template GetImage<MaskImageType>("Carina");
 456
 457     // Get Centroid and Z value
 458     std::vector<MaskImagePointType> centroids;
 459     clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
 460
 461     // We dont need Carina structure from now
 462     this->GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
 463
 464     // Put inside the AFDB
 465     this->GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
 466     this->GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
 467     this->WriteAFDB();
 468     z = centroids[1][2];
 469   }
 470   return z;
 471 }
 472 //--------------------------------------------------------------------
 473
 474
 475 //--------------------------------------------------------------------
 476 template <class TImageType>
 477 double
 478 clitk::ExtractLymphStationsFilter<TImageType>::
 479 FindApexOfTheChest()
 480 {
 481   double z;
 482   try {
 483     z = this->GetAFDB()->GetDouble("ApexOfTheChestZ");
 484   }
 485   catch(clitk::ExceptionObject e) {
 486
 487     /*
 488     //DD("FindApexOfTheChestPosition");
 489     MaskImagePointer Lungs = this->GetAFDB()->template GetImage<MaskImageType>("Lungs");
 490     MaskImagePointType p;
 491     p[0] = p[1] = p[2] =  0.0; // to avoid warning
 492     clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
 493
 494     // We dont need Lungs structure from now
 495     this->GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
 496
 497     // Put inside the AFDB
 498     this->GetAFDB()->SetPoint3D("ApexOfTheChest", p);
 499     p[2] -= 5; // We consider 5 mm lower
 500     this->GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
 501     this->WriteAFDB();
 502     z = p[2];
 503     */
 504
 505     // the superior border becomes the more inferior of the two apices
 506     MaskImagePointer RightLung = this->GetAFDB()->template GetImage<MaskImageType>("RightLung");
 507     MaskImagePointer LeftLung = this->GetAFDB()->template GetImage<MaskImageType>("LeftLung");
 508     MaskImagePointType pr;
 509     MaskImagePointType pl;
 510     clitk::FindExtremaPointInAGivenDirection<MaskImageType>(RightLung, GetBackgroundValue(), 2, false, pr);
 511     clitk::FindExtremaPointInAGivenDirection<MaskImageType>(LeftLung, GetBackgroundValue(), 2, false, pl);
 512     // We dont need Lungs structure from now
 513     this->GetAFDB()->template ReleaseImage<MaskImageType>("RightLung");
 514     this->GetAFDB()->template ReleaseImage<MaskImageType>("LeftLung");
 515     // Put inside the AFDB
 516     if (pr[2] < pl[2]) z = pr[2];
 517     else z = pl[2];
 518     this->GetAFDB()->SetDouble("ApexOfTheChestZ", z);
 519     this->WriteAFDB();
 520   }
 521   return z;
 522 }
 523 //--------------------------------------------------------------------
 524
 525
 526 //--------------------------------------------------------------------
 527 template <class TImageType>
 528 void
 529 clitk::ExtractLymphStationsFilter<TImageType>::
 530 FindLeftAndRightBronchi()
 531 {
 532   try {
 533     m_RightBronchus = this->GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
 534     m_LeftBronchus = this->GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
 535   }
 536   catch(clitk::ExceptionObject & o) {
 537
 538     DD("FindLeftAndRightBronchi");
 539     // The goal is to separate the trachea inferiorly to the carina into
 540     // a Left and Right bronchus.
 541
 542     // Get the trachea
 543     MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
 544
 545     // Get the Carina position
 546     double m_CarinaZ = FindCarina();
 547
 548     // Consider only inferiorly to the Carina
 549     MaskImagePointer m_Working_Trachea =
 550       clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
 551                                                        GetBackgroundValue());
 552
 553     // Labelize the trachea
 554     m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
 555
 556     // Carina position must at the first slice that separate the two
 557     // main bronchus (not superiorly). We thus first check that the
 558     // upper slice is composed of at least two labels
 559     MaskImagePointer RightBronchus;
 560     MaskImagePointer LeftBronchus;
 561     typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
 562     SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
 563     iter.SetFirstDirection(0);
 564     iter.SetSecondDirection(1);
 565     iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
 566     int maxLabel=0;
 567     while (!iter.IsAtReverseEndOfSlice()) {
 568       while (!iter.IsAtReverseEndOfLine()) {
 569         if (iter.Get() > maxLabel) maxLabel = iter.Get();
 570         --iter;
 571       }
 572       iter.PreviousLine();
 573     }
 574     if (maxLabel < 2) {
 575       clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
 576     }
 577
 578     // Compute 3D centroids of both parts to identify the left from the
 579     // right bronchus
 580     std::vector<ImagePointType> c;
 581     clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
 582     ImagePointType C1 = c[1];
 583     ImagePointType C2 = c[2];
 584
 585     ImagePixelType rightLabel;
 586     ImagePixelType leftLabel;
 587     if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
 588     else { rightLabel = 2; leftLabel = 1; }
 589
 590     // Select LeftLabel (set one label to Backgroundvalue)
 591     RightBronchus =
 592       clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
 593                                      GetBackgroundValue(), GetForegroundValue());
 594     /*
 595       SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
 596       leftLabel, GetBackgroundValue(), false);
 597     */
 598     LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
 599                                                   GetBackgroundValue(), GetForegroundValue());
 600     /*
 601       SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
 602       rightLabel, GetBackgroundValue(), false);
 603     */
 604
 605     // Crop images
 606     RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
 607     LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
 608
 609     // Insert int AFDB if need after
 610     this->GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
 611                                                  RightBronchus, true);
 612     this->GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
 613                                                  LeftBronchus, true);
 614   }
 615 }
 616 //--------------------------------------------------------------------
 617
 618
 619 //--------------------------------------------------------------------
 620 template <class TImageType>
 621 double
 622 clitk::ExtractLymphStationsFilter<TImageType>::
 623 FindSuperiorBorderOfAorticArch()
 624 {
 625   double z;
 626   try {
 627     z = this->GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
 628   }
 629   catch(clitk::ExceptionObject e) {
 630     //    DD("FindSuperiorBorderOfAorticArch");
 631     MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
 632     MaskImagePointType p;
 633     p[0] = p[1] = p[2] =  0.0; // to avoid warning
 634     clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
 635     p[2] += Aorta->GetSpacing()[2]; // the slice above
 636
 637     // We dont need Lungs structure from now
 638     this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
 639
 640     // Put inside the AFDB
 641     this->GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
 642     this->GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
 643     this->WriteAFDB();
 644     z = p[2];
 645   }
 646   return z;
 647 }
 648 //--------------------------------------------------------------------
 649
 650
 651 //--------------------------------------------------------------------
 652 template <class TImageType>
 653 double
 654 clitk::ExtractLymphStationsFilter<TImageType>::
 655 FindInferiorBorderOfAorticArch()
 656 {
 657   double z;
 658   try {
 659     z = this->GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
 660   }
 661   catch(clitk::ExceptionObject e) {
 662     //DD("FindInferiorBorderOfAorticArch");
 663     MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
 664     std::vector<MaskSlicePointer> slices;
 665     clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
 666     bool found=false;
 667     uint i = slices.size()-1;
 668     while (!found) {
 669       slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
 670       std::vector<typename MaskSliceType::PointType> c;
 671       clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
 672       if (c.size()>2) {
 673         found = true;
 674       }
 675       else {
 676         i--;
 677       }
 678     }
 679     MaskImageIndexType index;
 680     index[0] = index[1] = 0.0;
 681     index[2] = i+1;
 682     MaskImagePointType lower;
 683     Aorta->TransformIndexToPhysicalPoint(index, lower);
 684
 685     // We dont need Lungs structure from now
 686     this->GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
 687
 688     // Put inside the AFDB
 689     this->GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
 690     this->GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
 691     this->WriteAFDB();
 692     z = lower[2];
 693   }
 694   return z;
 695 }
 696 //--------------------------------------------------------------------
 697
 698
 699 //--------------------------------------------------------------------
 700 template <class ImageType>
 701 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
 702 clitk::ExtractLymphStationsFilter<ImageType>::
 703 FindAntPostVesselsOLD()
 704 {
 705   // -----------------------------------------------------
 706   /* Rod says: "The anterior border, as with the Atlas – UM, is
 707     posterior to the vessels (right subclavian vein, left
 708     brachiocephalic vein, right brachiocephalic vein, left subclavian
 709     artery, left common carotid artery and brachiocephalic trunk).
 710     These vessels are not included in the nodal station.  The anterior
 711     border is drawn to the midpoint of the vessel and an imaginary
 712     line joins the middle of these vessels.  Between the vessels,
 713     station 2 is in contact with station 3a." */
 714
 715   // Check if no already done
 716   bool found = true;
 717   MaskImagePointer binarizedContour;
 718   try {
 719     binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
 720   }
 721   catch(clitk::ExceptionObject e) {
 722     found = false;
 723   }
 724   if (found) {
 725     return binarizedContour;
 726   }
 727
 728   /* Here, we consider the vessels form a kind of anterior barrier. We
 729      link all vessels centroids and remove what is post to it.
 730     - select the list of structure
 731             vessel1 = BrachioCephalicArtery
 732             vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
 733             vessel3 = CommonCarotidArtery
 734             vessel4 = SubclavianArtery
 735             other   = Thyroid
 736             other   = Aorta
 737      - crop images as needed
 738      - slice by slice, choose the CCL and extract centroids
 739      - slice by slice, sort according to polar angle wrt Trachea centroid.
 740      - slice by slice, link centroids and close contour
 741      - remove outside this contour
 742      - merge with support
 743   */
 744
 745   // Read structures
 746   MaskImagePointer BrachioCephalicArtery = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
 747   MaskImagePointer BrachioCephalicVein = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
 748   MaskImagePointer CommonCarotidArtery = this->GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
 749   MaskImagePointer SubclavianArtery = this->GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
 750   MaskImagePointer Thyroid = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
 751   MaskImagePointer Aorta = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
 752   MaskImagePointer Trachea = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
 753
 754   // Create a temporay support
 755   // From first slice of BrachioCephalicVein to end of 3A
 756   MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
 757   MaskImagePointType p;
 758   p[0] = p[1] = p[2] =  0.0; // to avoid warning
 759   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
 760   double inf = p [2];
 761   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
 762                                                           GetBackgroundValue(), 2, false, p);
 763   double sup = p [2];
 764   support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
 765                                                         false, GetBackgroundValue());
 766
 767   // Resize all structures like support
 768   BrachioCephalicArtery =
 769     clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
 770   CommonCarotidArtery =
 771     clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
 772   SubclavianArtery =
 773     clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
 774   Thyroid =
 775     clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
 776   Aorta =
 777     clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
 778   BrachioCephalicVein =
 779     clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
 780   Trachea =
 781     clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
 782
 783   // Extract slices
 784   std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
 785   clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
 786   std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
 787   clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
 788   std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
 789   clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
 790   std::vector<MaskSlicePointer> slices_SubclavianArtery;
 791   clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
 792   std::vector<MaskSlicePointer> slices_Thyroid;
 793   clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
 794   std::vector<MaskSlicePointer> slices_Aorta;
 795   clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
 796   std::vector<MaskSlicePointer> slices_Trachea;
 797   clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
 798   unsigned int n= slices_BrachioCephalicArtery.size();
 799
 800   // Get the boundaries of one slice
 801   std::vector<MaskSlicePointType> bounds;
 802   ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
 803
 804   // For all slices, for all structures, find the centroid and build the contour
 805   // List of 3D points (for debug)
 806   std::vector<MaskImagePointType> p3D;
 807
 808   vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
 809   for(unsigned int i=0; i<n; i++) {
 810     // Labelize the slices
 811     slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
 812                                                             GetBackgroundValue(), true, 1);
 813     slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
 814                                                          GetBackgroundValue(), true, 1);
 815     slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
 816                                                              GetBackgroundValue(), true, 1);
 817     slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
 818                                                             GetBackgroundValue(), true, 1);
 819     slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
 820                                                 GetBackgroundValue(), true, 1);
 821     slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
 822                                               GetBackgroundValue(), true, 1);
 823
 824     // Search centroids
 825     std::vector<MaskSlicePointType> points2D;
 826     std::vector<MaskSlicePointType> centroids1;
 827     std::vector<MaskSlicePointType> centroids2;
 828     std::vector<MaskSlicePointType> centroids3;
 829     std::vector<MaskSlicePointType> centroids4;
 830     std::vector<MaskSlicePointType> centroids5;
 831     std::vector<MaskSlicePointType> centroids6;
 832     ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
 833     ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
 834     ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
 835     ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
 836     ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
 837     ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
 838
 839     // BrachioCephalicVein -> when it is separated into two CCL, we
 840     // only consider the most at Right one
 841     if (centroids6.size() > 2) {
 842       if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
 843       else centroids6.erase(centroids6.begin()+1);
 844     }
 845
 846     // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
 847     // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
 848     if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
 849       centroids6.clear();
 850     }
 851
 852     for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
 853     for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
 854     for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
 855     for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
 856     for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
 857     for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
 858
 859     // Sort by angle according to trachea centroid and vertical line,
 860     // in polar coordinates :
 861     // http://en.wikipedia.org/wiki/Polar_coordinate_system
 862     std::vector<MaskSlicePointType> centroids_trachea;
 863     ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
 864     typedef std::pair<MaskSlicePointType, double> PointAngleType;
 865     std::vector<PointAngleType> angles;
 866     for(unsigned int j=0; j<points2D.size(); j++) {
 867       //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
 868       double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
 869       double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
 870       double angle = 0;
 871       if (x>0) angle = atan(y/x);
 872       if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
 873       if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
 874       if (x==0) {
 875         if (y>0) angle = M_PI/2.0;
 876         if (y<0) angle = -M_PI/2.0;
 877         if (y==0) angle = 0;
 878       }
 879       angle = clitk::rad2deg(angle);
 880       // Angle is [-180;180] wrt the X axis. We change the X axis to
 881       // be the vertical line, because we want to sort from Right to
 882       // Left from Post to Ant.
 883       if (angle>0)
 884         angle = (270-angle);
 885       if (angle<0) {
 886         angle = -angle-90;
 887         if (angle<0) angle = 360-angle;
 888       }
 889       PointAngleType a;
 890       a.first = points2D[j];
 891       a.second = angle;
 892       angles.push_back(a);
 893     }
 894
 895     // Do nothing if less than 2 points --> n
 896     if (points2D.size() < 3) { //continue;
 897       continue;
 898     }
 899
 900     // Sort points2D according to polar angles
 901     std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
 902     for(unsigned int j=0; j<angles.size(); j++) {
 903       points2D[j] = angles[j].first;
 904     }
 905     //    DDV(points2D, points2D.size());
 906
 907     /* When vessels are far away, we try to replace the line segment
 908        with a curved line that join the two vessels but stay
 909        approximately at the same distance from the trachea centroids
 910        than the vessels.
 911
 912        For that:
 913        - let a and c be two successive vessels centroids
 914        - id distance(a,c) < threshold, next point
 915
 916        TODO HERE
 917
 918        - compute mid position between two successive points -
 919        compute dist to trachea centroid for the 3 pts - if middle too
 920        low, add one point
 921     */
 922     std::vector<MaskSlicePointType> toadd;
 923     unsigned int index = 0;
 924     double dmax = 5;
 925     while (index<points2D.size()-1) {
 926       MaskSlicePointType a = points2D[index];
 927       MaskSlicePointType c = points2D[index+1];
 928
 929       double dac = a.EuclideanDistanceTo(c);
 930       if (dac>dmax) {
 931
 932         MaskSlicePointType b;
 933         b[0] = a[0]+(c[0]-a[0])/2.0;
 934         b[1] = a[1]+(c[1]-a[1])/2.0;
 935
 936         // Compute distance to trachea centroid
 937         MaskSlicePointType m = centroids_trachea[1];
 938         double da = m.EuclideanDistanceTo(a);
 939         double db = m.EuclideanDistanceTo(b);
 940         //double dc = m.EuclideanDistanceTo(c);
 941
 942         // Mean distance, find point on the line from trachea centroid
 943         // to b
 944         double alpha = (da+db)/2.0;
 945         MaskSlicePointType v;
 946         double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
 947         v[0] = (b[0]-m[0])/n;
 948         v[1] = (b[1]-m[1])/n;
 949         MaskSlicePointType r;
 950         r[0] = m[0]+alpha*v[0];
 951         r[1] = m[1]+alpha*v[1];
 952         points2D.insert(points2D.begin()+index+1, r);
 953       }
 954       else {
 955         index++;
 956       }
 957     }
 958     //    DDV(points2D, points2D.size());
 959
 960     // Add some points to close the contour
 961     // - H line towards Right
 962     MaskSlicePointType p = points2D[0];
 963     p[0] = bounds[0][0];
 964     points2D.insert(points2D.begin(), p);
 965     // - corner Right/Post
 966     p = bounds[0];
 967     points2D.insert(points2D.begin(), p);
 968     // - H line towards Left
 969     p = points2D.back();
 970     p[0] = bounds[2][0];
 971     points2D.push_back(p);
 972     // - corner Right/Post
 973     p = bounds[2];
 974     points2D.push_back(p);
 975     // Close contour with the first point
 976     points2D.push_back(points2D[0]);
 977     //    DDV(points2D, points2D.size());
 978
 979     // Build 3D points from the 2D points
 980     std::vector<ImagePointType> points3D;
 981     clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
 982     for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
 983
 984     // Build the mesh from the contour's points
 985     vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
 986     append->AddInput(mesh);
 987   }
 988
 989   // DEBUG: write points3D
 990   clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
 991
 992   // Build the final 3D mesh form the list 2D mesh
 993   append->Update();
 994   vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
 995
 996   // Debug, write the mesh
 997   /*
 998     vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
 999     w->SetInput(mesh);
1000     w->SetFileName("bidon.vtk");
1001     w->Write();
1002   */
1003
1004   // Compute a single binary 3D image from the list of contours
1005   clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1006     clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1007   filter->SetMesh(mesh);
1008   filter->SetLikeImage(support);
1009   filter->Update();
1010   binarizedContour = filter->GetOutput();
1011
1012   // Crop
1013   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
1014   inf = p[2];
1015   DD(p);
1016   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
1017   sup = p[2];
1018   DD(p);
1019   binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1020                                                         false, GetBackgroundValue());
1021   // Update the AFDB
1022   writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mha");
1023   this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mha");
1024   this->WriteAFDB();
1025   return binarizedContour;
1026
1027   /*
1028   // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1029   ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1030   p[1] += 50; // 50 mm Post from this point must be kept
1031   ImageIndexType index;
1032   binarizedContour->TransformPhysicalPointToIndex(p, index);
1033   bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1034
1035   // remove from support
1036   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1037   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1038   boolFilter->InPlaceOn();
1039   boolFilter->SetInput1(m_Working_Support);
1040   boolFilter->SetInput2(binarizedContour);
1041   boolFilter->SetBackgroundValue1(GetBackgroundValue());
1042   boolFilter->SetBackgroundValue2(GetBackgroundValue());
1043   if (isInside)
1044     boolFilter->SetOperationType(BoolFilterType::And);
1045   else
1046     boolFilter->SetOperationType(BoolFilterType::AndNot);
1047   boolFilter->Update();
1048   m_Working_Support = boolFilter->GetOutput();
1049   */
1050
1051   // End
1052   //StopCurrentStep<MaskImageType>(m_Working_Support);
1053   //m_ListOfStations["2R"] = m_Working_Support;
1054   //m_ListOfStations["2L"] = m_Working_Support;
1055 }
1056 //--------------------------------------------------------------------
1057
1058
1059 //--------------------------------------------------------------------
1060 template <class ImageType>
1061 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1062 clitk::ExtractLymphStationsFilter<ImageType>::
1063 FindAntPostVessels2()
1064 {
1065   // -----------------------------------------------------
1066   /* Rod says: "The anterior border, as with the Atlas – UM, is
1067     posterior to the vessels (right subclavian vein, left
1068     brachiocephalic vein, right brachiocephalic vein, left subclavian
1069     artery, left common carotid artery and brachiocephalic trunk).
1070     These vessels are not included in the nodal station.  The anterior
1071     border is drawn to the midpoint of the vessel and an imaginary
1072     line joins the middle of these vessels.  Between the vessels,
1073     station 2 is in contact with station 3a." */
1074
1075   // Check if no already done
1076   bool found = true;
1077   MaskImagePointer binarizedContour;
1078   try {
1079     binarizedContour = this->GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
1080   }
1081   catch(clitk::ExceptionObject e) {
1082     found = false;
1083   }
1084   if (found) {
1085     return binarizedContour;
1086   }
1087
1088   /* Here, we consider the vessels form a kind of anterior barrier. We
1089      link all vessels centroids and remove what is post to it.
1090     - select the list of structure
1091             vessel1 = BrachioCephalicArtery
1092             vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
1093             vessel3 = CommonCarotidArtery
1094             vessel4 = SubclavianArtery
1095             other   = Thyroid
1096             other   = Aorta
1097      - crop images as needed
1098      - slice by slice, choose the CCL and extract centroids
1099      - slice by slice, sort according to polar angle wrt Trachea centroid.
1100      - slice by slice, link centroids and close contour
1101      - remove outside this contour
1102      - merge with support
1103   */
1104
1105   // Read structures
1106   std::map<std::string, MaskImagePointer> MapOfStructures;
1107   typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
1108   MapOfStructures["BrachioCephalicArtery"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
1109   MapOfStructures["BrachioCephalicVein"] = this->GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
1110   MapOfStructures["CommonCarotidArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("LeftCommonCarotidArtery");
1111   MapOfStructures["CommonCarotidArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("RightCommonCarotidArtery");
1112   MapOfStructures["SubclavianArteryLeft"] = this->GetAFDB()->template GetImage<MaskImageType>("LeftSubclavianArtery");
1113   MapOfStructures["SubclavianArteryRight"] = this->GetAFDB()->template GetImage<MaskImageType>("RightSubclavianArtery");
1114   MapOfStructures["Thyroid"] = this->GetAFDB()->template GetImage<MaskImageType>("Thyroid");
1115   MapOfStructures["Aorta"] = this->GetAFDB()->template GetImage<MaskImageType>("Aorta");
1116   MapOfStructures["Trachea"] = this->GetAFDB()->template GetImage<MaskImageType>("Trachea");
1117
1118   std::vector<std::string> ListOfStructuresNames;
1119
1120   // Create a temporay support
1121   // From first slice of BrachioCephalicVein to end of 3A or end of 2RL
1122   MaskImagePointer support = this->GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
1123   MaskImagePointType p;
1124   p[0] = p[1] = p[2] =  0.0; // to avoid warning
1125   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
1126                                                           GetBackgroundValue(), 2, true, p);
1127   double inf = p[2];
1128   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
1129                                                           GetBackgroundValue(), 2, false, p);
1130   MaskImagePointType p2;
1131   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(this->GetAFDB()->template GetImage<MaskImageType>("Support_S2L"),
1132                                                           GetBackgroundValue(), 2, false, p2);
1133   if (p2[2] > p[2]) p = p2;
1134
1135   double sup = p[2]+support->GetSpacing()[2];//one slice more ?
1136   support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
1137                                                         false, GetBackgroundValue());
1138
1139   // Resize all structures like support
1140   for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1141     it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
1142   }
1143
1144   // Extract slices
1145   typedef std::vector<MaskSlicePointer> SlicesType;
1146   std::map<std::string, SlicesType> MapOfSlices;
1147   for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1148     SlicesType s;
1149     clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
1150     MapOfSlices[it->first] = s;
1151   }
1152
1153   unsigned int n= MapOfSlices["Trachea"].size();
1154
1155   // Get the boundaries of one slice
1156   std::vector<MaskSlicePointType> bounds;
1157   ComputeImageBoundariesCoordinates<MaskSliceType>(MapOfSlices["Trachea"][0], bounds);
1158
1159   // LOOP ON SLICES
1160   // For all slices, for all structures, find the centroid and build the contour
1161   // List of 3D points (for debug)
1162   std::vector<MaskImagePointType> p3D;
1163   vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
1164   for(unsigned int i=0; i<n; i++) {
1165
1166     // Labelize the slices
1167     for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1168       MaskSlicePointer & s = MapOfSlices[it->first][i];
1169       s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
1170     }
1171
1172     // Search centroids
1173     std::vector<MaskSlicePointType> points2D;
1174     typedef std::vector<MaskSlicePointType> CentroidsType;
1175     std::map<std::string, CentroidsType> MapOfCentroids;
1176     for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1177       std::string structure = it->first;
1178       MaskSlicePointer & s = MapOfSlices[structure][i];
1179       CentroidsType c;
1180       clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), c);
1181       MapOfCentroids[structure] = c;
1182     }
1183
1184
1185     // BrachioCephalicVein -> when it is separated into two CCL, we
1186     // only consider the most at Right one
1187     CentroidsType & c = MapOfCentroids["BrachioCephalicVein"];
1188     if (c.size() > 2) {
1189       if (c[1][0] < c[2][0]) c.erase(c.begin()+2);
1190       else c.erase(c.begin()+1);
1191     }
1192
1193     /*
1194     // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
1195     // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
1196     if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1)
1197         && (MapOfCentroids["SubclavianArtery"].size() > 2)) {
1198       MapOfCentroids["BrachioCephalicVein"].clear();
1199     }
1200     */
1201
1202     // Add all 2D points
1203     for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1204       std::string structure = it->first;
1205       if (structure != "Trachea") { // do not add centroid of trachea
1206         CentroidsType & c = MapOfCentroids[structure];
1207         for(unsigned int j=1; j<c.size(); j++) points2D.push_back(c[j]);
1208       }
1209     }
1210
1211     // Sort by angle according to trachea centroid and vertical line,
1212     // in polar coordinates :
1213     // http://en.wikipedia.org/wiki/Polar_coordinate_system
1214     //    std::vector<MaskSlicePointType> centroids_trachea;
1215     //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
1216     CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
1217     typedef std::pair<MaskSlicePointType, double> PointAngleType;
1218     std::vector<PointAngleType> angles;
1219     for(unsigned int j=0; j<points2D.size(); j++) {
1220       //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
1221       double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
1222       double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
1223       double angle = 0;
1224       if (x>0) angle = atan(y/x);
1225       if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
1226       if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
1227       if (x==0) {
1228         if (y>0) angle = M_PI/2.0;
1229         if (y<0) angle = -M_PI/2.0;
1230         if (y==0) angle = 0;
1231       }
1232       angle = clitk::rad2deg(angle);
1233       // Angle is [-180;180] wrt the X axis. We change the X axis to
1234       // be the vertical line, because we want to sort from Right to
1235       // Left from Post to Ant.
1236       if (angle>0)
1237         angle = (270-angle);
1238       if (angle<0) {
1239         angle = -angle-90;
1240         if (angle<0) angle = 360-angle;
1241       }
1242       PointAngleType a;
1243       a.first = points2D[j];
1244       a.second = angle;
1245       angles.push_back(a);
1246     }
1247
1248     // Do nothing if less than 2 points --> n
1249     if (points2D.size() < 3) { //continue;
1250       continue;
1251     }
1252
1253     // Sort points2D according to polar angles
1254     std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
1255     for(unsigned int j=0; j<angles.size(); j++) {
1256       points2D[j] = angles[j].first;
1257     }
1258     //    DDV(points2D, points2D.size());
1259
1260     /* When vessels are far away, we try to replace the line segment
1261        with a curved line that join the two vessels but stay
1262        approximately at the same distance from the trachea centroids
1263        than the vessels.
1264
1265        For that:
1266        - let a and c be two successive vessels centroids
1267        - id distance(a,c) < threshold, next point
1268
1269        TODO HERE
1270
1271        - compute mid position between two successive points -
1272        compute dist to trachea centroid for the 3 pts - if middle too
1273        low, add one point
1274     */
1275     std::vector<MaskSlicePointType> toadd;
1276     unsigned int index = 0;
1277     double dmax = 5;
1278     while (index<points2D.size()-1) {
1279       MaskSlicePointType a = points2D[index];
1280       MaskSlicePointType c = points2D[index+1];
1281
1282       double dac = a.EuclideanDistanceTo(c);
1283       if (dac>dmax) {
1284
1285         MaskSlicePointType b;
1286         b[0] = a[0]+(c[0]-a[0])/2.0;
1287         b[1] = a[1]+(c[1]-a[1])/2.0;
1288
1289         // Compute distance to trachea centroid
1290         MaskSlicePointType m = centroids_trachea[1];
1291         double da = m.EuclideanDistanceTo(a);
1292         double db = m.EuclideanDistanceTo(b);
1293         //double dc = m.EuclideanDistanceTo(c);
1294
1295         // Mean distance, find point on the line from trachea centroid
1296         // to b
1297         double alpha = (da+db)/2.0;
1298         MaskSlicePointType v;
1299         double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
1300         v[0] = (b[0]-m[0])/n;
1301         v[1] = (b[1]-m[1])/n;
1302         MaskSlicePointType r;
1303         r[0] = m[0]+alpha*v[0];
1304         r[1] = m[1]+alpha*v[1];
1305         points2D.insert(points2D.begin()+index+1, r);
1306       }
1307       else {
1308         index++;
1309       }
1310     }
1311     //    DDV(points2D, points2D.size());
1312
1313     // Add some points to close the contour
1314     // - H line towards Right
1315     MaskSlicePointType p = points2D[0];
1316     p[0] = bounds[0][0];
1317     points2D.insert(points2D.begin(), p);
1318     // - corner Right/Post
1319     p = bounds[0];
1320     points2D.insert(points2D.begin(), p);
1321     // - H line towards Left
1322     p = points2D.back();
1323     p[0] = bounds[2][0];
1324     points2D.push_back(p);
1325     // - corner Right/Post
1326     p = bounds[2];
1327     points2D.push_back(p);
1328     // Close contour with the first point
1329     points2D.push_back(points2D[0]);
1330     //    DDV(points2D, points2D.size());
1331
1332     // Build 3D points from the 2D points
1333     std::vector<ImagePointType> points3D;
1334     clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
1335     for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
1336
1337     // Build the mesh from the contour's points
1338     vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1339     append->AddInput(mesh);
1340     // if (i ==n-1) { // last slice
1341     //   clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
1342     //   vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1343     //   append->AddInput(mesh);
1344     // }
1345   }
1346
1347   // DEBUG: write points3D
1348   clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
1349
1350   // Build the final 3D mesh form the list 2D mesh
1351   append->Update();
1352   vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
1353
1354   // Debug, write the mesh
1355   /*
1356   vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1357   w->SetInput(mesh);
1358   w->SetFileName("bidon.vtk");
1359   w->Write();
1360   */
1361
1362   // Compute a single binary 3D image from the list of contours
1363   clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1364     clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1365   filter->SetMesh(mesh);
1366   filter->SetLikeImage(support);
1367   filter->Update();
1368   binarizedContour = filter->GetOutput();
1369
1370   // Crop
1371   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1372                                                           GetForegroundValue(), 2, true, p);
1373   inf = p[2];
1374   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1375                                                           GetForegroundValue(), 2, false, p);
1376   sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
1377   binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1378                                                         false, GetBackgroundValue());
1379
1380   // Update the AFDB
1381   writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mha");
1382   this->GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mha");
1383   this->WriteAFDB();
1384   return binarizedContour;
1385
1386   /*
1387   // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1388   ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1389   p[1] += 50; // 50 mm Post from this point must be kept
1390   ImageIndexType index;
1391   binarizedContour->TransformPhysicalPointToIndex(p, index);
1392   bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1393
1394   // remove from support
1395   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1396   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1397   boolFilter->InPlaceOn();
1398   boolFilter->SetInput1(m_Working_Support);
1399   boolFilter->SetInput2(binarizedContour);
1400   boolFilter->SetBackgroundValue1(GetBackgroundValue());
1401   boolFilter->SetBackgroundValue2(GetBackgroundValue());
1402   if (isInside)
1403     boolFilter->SetOperationType(BoolFilterType::And);
1404   else
1405     boolFilter->SetOperationType(BoolFilterType::AndNot);
1406   boolFilter->Update();
1407   m_Working_Support = boolFilter->GetOutput();
1408   */
1409
1410   // End
1411   //StopCurrentStep<MaskImageType>(m_Working_Support);
1412   //m_ListOfStations["2R"] = m_Working_Support;
1413   //m_ListOfStations["2L"] = m_Working_Support;
1414 }
1415 //--------------------------------------------------------------------
1416
1417 //--------------------------------------------------------------------
1418 template <class ImageType>
1419 void
1420 clitk::ExtractLymphStationsFilter<ImageType>::
1421 WriteImageSupport(std::string support)
1422 {
1423   writeImage<MaskImageType>(m_ListOfSupports[support], this->GetAFDBPath()+"/"+"seg/Support_"+support+".mha");
1424   this->GetAFDB()->SetImageFilename("Support_"+support, "seg/Support_"+support+".mha");
1425 }
1426 //--------------------------------------------------------------------
1427
1428
1429 //--------------------------------------------------------------------
1430 template <class ImageType>
1431 void
1432 clitk::ExtractLymphStationsFilter<ImageType>::
1433 WriteImageStation(std::string station)
1434 {
1435   writeImage<MaskImageType>(m_ListOfStations[station], GetAFDB()->GetPath()+"/seg/Station"+station+".mha");
1436   GetAFDB()->SetImageFilename("Station"+station, "seg/Station"+station+".mha");
1437   WriteAFDB();
1438 }
1439 //--------------------------------------------------------------------
1440
1441
1442 //--------------------------------------------------------------------
1443 template <class ImageType>
1444 void
1445 clitk::ExtractLymphStationsFilter<ImageType>::
1446 ComputeOverlapWithRef(std::string station)
1447 {
1448   if (GetComputeStation(station)) {
1449     MaskImagePointer ref = this->GetAFDB()->template GetImage <MaskImageType>("Station"+station+"_Ref");
1450     typedef clitk::LabelImageOverlapMeasureFilter<MaskImageType> FilterType;
1451     typename FilterType::Pointer filter = FilterType::New();
1452     filter->SetInput(0, m_ListOfStations[station]);
1453     filter->SetInput(1, ref);
1454     filter->Update();
1455   }
1456 }
1457 //--------------------------------------------------------------------
1458
1459 //--------------------------------------------------------------------
1460 template <class ImageType>
1461 void
1462 clitk::ExtractLymphStationsFilter<ImageType>::
1463 ReadSupportLimits(std::string filename)
1464 {
1465   m_ListOfSupportLimits.clear();
1466   ifstream is;
1467   openFileForReading(is, filename);
1468   while (is) {
1469     skipComment(is);
1470     SupportLimitsType s;
1471     s.Read(is);
1472     if (is) m_ListOfSupportLimits.push_back(s);
1473   }
1474 }
1475 //--------------------------------------------------------------------
1476
1477 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX
1478