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