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