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