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