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 void
 357 clitk::ExtractLymphStationsFilter<TImageType>::
 358 SetThreshold(std::string station, std::string tag, double value)
 359 {
 360   m_Threshold[station][tag] = value;
 361 }
 362 //--------------------------------------------------------------------
 363
 364
 365 //--------------------------------------------------------------------
 366 template <class TImageType>
 367 double
 368 clitk::ExtractLymphStationsFilter<TImageType>::
 369 GetFuzzyThreshold(std::string station, std::string tag)
 370 {
 371   if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) {
 372     clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
 373     return 0.0;
 374   }
 375
 376   if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) {
 377     clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+".");
 378     return 0.0;
 379   }
 380
 381   return m_FuzzyThreshold[station][tag];
 382 }
 383 //--------------------------------------------------------------------
 384
 385
 386 //--------------------------------------------------------------------
 387 template <class TImageType>
 388 double
 389 clitk::ExtractLymphStationsFilter<TImageType>::
 390 GetThreshold(std::string station, std::string tag)
 391 {
 392   if (m_Threshold.find(station) == m_Threshold.end()) {
 393     clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+").");
 394     return 0.0;
 395   }
 396
 397   if (m_Threshold[station].find(tag) == m_Threshold[station].end()) {
 398     clitkExceptionMacro("Could not find options "+tag+" in the list of Threshold for station "+station+".");
 399     return 0.0;
 400   }
 401
 402   return m_Threshold[station][tag];
 403 }
 404 //--------------------------------------------------------------------
 405
 406
 407 //--------------------------------------------------------------------
 408 template <class TImageType>
 409 void
 410 clitk::ExtractLymphStationsFilter<TImageType>::
 411 FindLineForS7S8Separation(MaskImagePointType & A, MaskImagePointType & B)
 412 {
 413   // Create line from A to B with
 414   // A = upper border of LLL at left
 415   // B = lower border of bronchus intermedius (BI) or RightMiddleLobeBronchus
 416
 417   try {
 418     GetAFDB()->GetPoint3D("LineForS7S8Separation_Begin", A);
 419     GetAFDB()->GetPoint3D("LineForS7S8Separation_End", B);
 420   }
 421   catch(clitk::ExceptionObject & o) {
 422
 423     DD("FindLineForS7S8Separation");
 424     // Load LeftLowerLobeBronchus and get centroid point
 425     MaskImagePointer LeftLowerLobeBronchus =
 426       GetAFDB()->template GetImage <MaskImageType>("LeftLowerLobeBronchus");
 427     std::vector<MaskImagePointType> c;
 428     clitk::ComputeCentroids<MaskImageType>(LeftLowerLobeBronchus, GetBackgroundValue(), c);
 429     A = c[1];
 430
 431     // Load RightMiddleLobeBronchus and get superior point (not centroid here)
 432     MaskImagePointer RightMiddleLobeBronchus =
 433       GetAFDB()->template GetImage <MaskImageType>("RightMiddleLobeBronchus");
 434     bool b = FindExtremaPointInAGivenDirection<MaskImageType>(RightMiddleLobeBronchus,
 435                                                               GetBackgroundValue(),
 436                                                               2, false, B);
 437     if (!b) {
 438       clitkExceptionMacro("Error while searching most superior point in RightMiddleLobeBronchus. Abort");
 439     }
 440
 441     // Insert into the DB
 442     GetAFDB()->SetPoint3D("LineForS7S8Separation_Begin", A);
 443     GetAFDB()->SetPoint3D("LineForS7S8Separation_End", B);
 444   }
 445 }
 446 //--------------------------------------------------------------------
 447
 448
 449 //--------------------------------------------------------------------
 450 template <class TImageType>
 451 double
 452 clitk::ExtractLymphStationsFilter<TImageType>::
 453 FindCarina()
 454 {
 455   double z;
 456   try {
 457     z = GetAFDB()->GetDouble("CarinaZ");
 458   }
 459   catch(clitk::ExceptionObject e) {
 460     DD("FindCarinaSlicePosition");
 461     // Get Carina
 462     MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
 463
 464     // Get Centroid and Z value
 465     std::vector<MaskImagePointType> centroids;
 466     clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
 467
 468     // We dont need Carina structure from now
 469     GetAFDB()->template ReleaseImage<MaskImageType>("Carina");
 470
 471     // Put inside the AFDB
 472     GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
 473     GetAFDB()->SetDouble("CarinaZ", centroids[1][2]);
 474     WriteAFDB();
 475     z = centroids[1][2];
 476   }
 477   return z;
 478 }
 479 //--------------------------------------------------------------------
 480
 481
 482 //--------------------------------------------------------------------
 483 template <class TImageType>
 484 double
 485 clitk::ExtractLymphStationsFilter<TImageType>::
 486 FindApexOfTheChest()
 487 {
 488   double z;
 489   try {
 490     z = GetAFDB()->GetDouble("ApexOfTheChestZ");
 491   }
 492   catch(clitk::ExceptionObject e) {
 493     DD("FindApexOfTheChestPosition");
 494     MaskImagePointer Lungs = GetAFDB()->template GetImage<MaskImageType>("Lungs");
 495     MaskImagePointType p;
 496     p[0] = p[1] = p[2] =  0.0; // to avoid warning
 497     clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Lungs, GetBackgroundValue(), 2, false, p);
 498
 499     // We dont need Lungs structure from now
 500     GetAFDB()->template ReleaseImage<MaskImageType>("Lungs");
 501
 502     // Put inside the AFDB
 503     GetAFDB()->SetPoint3D("ApexOfTheChest", p);
 504     p[2] -= 5; // We consider 5 mm lower
 505     GetAFDB()->SetDouble("ApexOfTheChestZ", p[2]);
 506     WriteAFDB();
 507     z = p[2];
 508   }
 509   return z;
 510 }
 511 //--------------------------------------------------------------------
 512
 513
 514 //--------------------------------------------------------------------
 515 template <class TImageType>
 516 void
 517 clitk::ExtractLymphStationsFilter<TImageType>::
 518 FindLeftAndRightBronchi()
 519 {
 520   try {
 521     m_RightBronchus = GetAFDB()->template GetImage <MaskImageType>("RightBronchus");
 522     m_LeftBronchus = GetAFDB()->template GetImage <MaskImageType>("LeftBronchus");
 523   }
 524   catch(clitk::ExceptionObject & o) {
 525
 526     DD("FindLeftAndRightBronchi");
 527     // The goal is to separate the trachea inferiorly to the carina into
 528     // a Left and Right bronchus.
 529
 530     // Get the trachea
 531     MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
 532
 533     // Get the Carina position
 534     double m_CarinaZ = FindCarina();
 535
 536     // Consider only inferiorly to the Carina
 537     MaskImagePointer m_Working_Trachea =
 538       clitk::CropImageRemoveGreaterThan<MaskImageType>(Trachea, 2, m_CarinaZ, true, // AutoCrop
 539                                                        GetBackgroundValue());
 540
 541     // Labelize the trachea
 542     m_Working_Trachea = Labelize<MaskImageType>(m_Working_Trachea, 0, true, 1);
 543
 544     // Carina position must at the first slice that separate the two
 545     // main bronchus (not superiorly). We thus first check that the
 546     // upper slice is composed of at least two labels
 547     MaskImagePointer RightBronchus;
 548     MaskImagePointer LeftBronchus;
 549     typedef itk::ImageSliceIteratorWithIndex<MaskImageType> SliceIteratorType;
 550     SliceIteratorType iter(m_Working_Trachea, m_Working_Trachea->GetLargestPossibleRegion());
 551     iter.SetFirstDirection(0);
 552     iter.SetSecondDirection(1);
 553     iter.GoToReverseBegin(); // Start from the end (because image is IS not SI)
 554     int maxLabel=0;
 555     while (!iter.IsAtReverseEndOfSlice()) {
 556       while (!iter.IsAtReverseEndOfLine()) {
 557         if (iter.Get() > maxLabel) maxLabel = iter.Get();
 558         --iter;
 559       }
 560       iter.PreviousLine();
 561     }
 562     if (maxLabel < 2) {
 563       clitkExceptionMacro("First slice from Carina does not seems to seperate the two main bronchus. Abort");
 564     }
 565
 566     // Compute 3D centroids of both parts to identify the left from the
 567     // right bronchus
 568     std::vector<ImagePointType> c;
 569     clitk::ComputeCentroids<MaskImageType>(m_Working_Trachea, GetBackgroundValue(), c);
 570     ImagePointType C1 = c[1];
 571     ImagePointType C2 = c[2];
 572
 573     ImagePixelType rightLabel;
 574     ImagePixelType leftLabel;
 575     if (C1[0] < C2[0]) { rightLabel = 1; leftLabel = 2; }
 576     else { rightLabel = 2; leftLabel = 1; }
 577
 578     // Select LeftLabel (set one label to Backgroundvalue)
 579     RightBronchus =
 580       clitk::Binarize<MaskImageType>(m_Working_Trachea, rightLabel, rightLabel,
 581                                      GetBackgroundValue(), GetForegroundValue());
 582     /*
 583       SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
 584       leftLabel, GetBackgroundValue(), false);
 585     */
 586     LeftBronchus = clitk::Binarize<MaskImageType>(m_Working_Trachea, leftLabel, leftLabel,
 587                                                   GetBackgroundValue(), GetForegroundValue());
 588     /*
 589       SetBackground<MaskImageType, MaskImageType>(m_Working_Trachea, m_Working_Trachea,
 590       rightLabel, GetBackgroundValue(), false);
 591     */
 592
 593     // Crop images
 594     RightBronchus = clitk::AutoCrop<MaskImageType>(RightBronchus, GetBackgroundValue());
 595     LeftBronchus = clitk::AutoCrop<MaskImageType>(LeftBronchus, GetBackgroundValue());
 596
 597     // Insert int AFDB if need after
 598     GetAFDB()->template SetImage <MaskImageType>("RightBronchus", "seg/rightBronchus.mhd",
 599                                                  RightBronchus, true);
 600     GetAFDB()->template SetImage <MaskImageType>("LeftBronchus", "seg/leftBronchus.mhd",
 601                                                  LeftBronchus, true);
 602   }
 603 }
 604 //--------------------------------------------------------------------
 605
 606
 607 //--------------------------------------------------------------------
 608 template <class TImageType>
 609 double
 610 clitk::ExtractLymphStationsFilter<TImageType>::
 611 FindSuperiorBorderOfAorticArch()
 612 {
 613   double z;
 614   try {
 615     z = GetAFDB()->GetDouble("SuperiorBorderOfAorticArchZ");
 616   }
 617   catch(clitk::ExceptionObject e) {
 618     DD("FindSuperiorBorderOfAorticArch");
 619     MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 620     MaskImagePointType p;
 621     p[0] = p[1] = p[2] =  0.0; // to avoid warning
 622     clitk::FindExtremaPointInAGivenDirection<MaskImageType>(Aorta, GetBackgroundValue(), 2, false, p);
 623     p[2] += Aorta->GetSpacing()[2]; // the slice above
 624
 625     // We dont need Lungs structure from now
 626     GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
 627
 628     // Put inside the AFDB
 629     GetAFDB()->SetPoint3D("SuperiorBorderOfAorticArch", p);
 630     GetAFDB()->SetDouble("SuperiorBorderOfAorticArchZ", p[2]);
 631     WriteAFDB();
 632     z = p[2];
 633   }
 634   return z;
 635 }
 636 //--------------------------------------------------------------------
 637
 638
 639 //--------------------------------------------------------------------
 640 template <class TImageType>
 641 double
 642 clitk::ExtractLymphStationsFilter<TImageType>::
 643 FindInferiorBorderOfAorticArch()
 644 {
 645   double z;
 646   try {
 647     z = GetAFDB()->GetDouble("InferiorBorderOfAorticArchZ");
 648   }
 649   catch(clitk::ExceptionObject e) {
 650     DD("FindInferiorBorderOfAorticArch");
 651     MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 652     std::vector<MaskSlicePointer> slices;
 653     clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices);
 654     bool found=false;
 655     uint i = slices.size()-1;
 656     while (!found) {
 657       slices[i] = Labelize<MaskSliceType>(slices[i], 0, false, 10);
 658       std::vector<typename MaskSliceType::PointType> c;
 659       clitk::ComputeCentroids<MaskSliceType>(slices[i], GetBackgroundValue(), c);
 660       if (c.size()>2) {
 661         found = true;
 662       }
 663       else {
 664         i--;
 665       }
 666     }
 667     MaskImageIndexType index;
 668     index[0] = index[1] = 0.0;
 669     index[2] = i+1;
 670     MaskImagePointType lower;
 671     Aorta->TransformIndexToPhysicalPoint(index, lower);
 672
 673     // We dont need Lungs structure from now
 674     GetAFDB()->template ReleaseImage<MaskImageType>("Aorta");
 675
 676     // Put inside the AFDB
 677     GetAFDB()->SetPoint3D("InferiorBorderOfAorticArch", lower);
 678     GetAFDB()->SetDouble("InferiorBorderOfAorticArchZ", lower[2]);
 679     WriteAFDB();
 680     z = lower[2];
 681   }
 682   return z;
 683 }
 684 //--------------------------------------------------------------------
 685
 686
 687 //--------------------------------------------------------------------
 688 template <class ImageType>
 689 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
 690 clitk::ExtractLymphStationsFilter<ImageType>::
 691 FindAntPostVessels()
 692 {
 693   // -----------------------------------------------------
 694   /* Rod says: "The anterior border, as with the Atlas – UM, is
 695     posterior to the vessels (right subclavian vein, left
 696     brachiocephalic vein, right brachiocephalic vein, left subclavian
 697     artery, left common carotid artery and brachiocephalic trunk).
 698     These vessels are not included in the nodal station.  The anterior
 699     border is drawn to the midpoint of the vessel and an imaginary
 700     line joins the middle of these vessels.  Between the vessels,
 701     station 2 is in contact with station 3a." */
 702
 703   // Check if no already done
 704   bool found = true;
 705   MaskImagePointer binarizedContour;
 706   try {
 707     DD("FindAntPostVessels try to get");
 708     binarizedContour = GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
 709   }
 710   catch(clitk::ExceptionObject e) {
 711     DD("not found");
 712     found = false;
 713   }
 714   if (found) {
 715     return binarizedContour;
 716   }
 717
 718   /* Here, we consider the vessels form a kind of anterior barrier. We
 719      link all vessels centroids and remove what is post to it.
 720     - select the list of structure
 721             vessel1 = BrachioCephalicArtery
 722             vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
 723             vessel3 = CommonCarotidArtery
 724             vessel4 = SubclavianArtery
 725             other   = Thyroid
 726             other   = Aorta
 727      - crop images as needed
 728      - slice by slice, choose the CCL and extract centroids
 729      - slice by slice, sort according to polar angle wrt Trachea centroid.
 730      - slice by slice, link centroids and close contour
 731      - remove outside this contour
 732      - merge with support
 733   */
 734
 735   // Read structures
 736   MaskImagePointer BrachioCephalicArtery = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
 737   MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
 738   MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArtery");
 739   MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage<MaskImageType>("SubclavianArtery");
 740   MaskImagePointer Thyroid = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
 741   MaskImagePointer Aorta = GetAFDB()->template GetImage<MaskImageType>("Aorta");
 742   MaskImagePointer Trachea = GetAFDB()->template GetImage<MaskImageType>("Trachea");
 743
 744   // Create a temporay support
 745   // From first slice of BrachioCephalicVein to end of 3A
 746   MaskImagePointer support = GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
 747   MaskImagePointType p;
 748   p[0] = p[1] = p[2] =  0.0; // to avoid warning
 749   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(BrachioCephalicVein, GetBackgroundValue(), 2, true, p);
 750   double inf = p [2];
 751   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
 752                                                           GetBackgroundValue(), 2, false, p);
 753   double sup = p [2];
 754   support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
 755                                                         false, GetBackgroundValue());
 756
 757   // Resize all structures like support
 758   BrachioCephalicArtery =
 759     clitk::ResizeImageLike<MaskImageType>(BrachioCephalicArtery, support, GetBackgroundValue());
 760   CommonCarotidArtery =
 761     clitk::ResizeImageLike<MaskImageType>(CommonCarotidArtery, support, GetBackgroundValue());
 762   SubclavianArtery =
 763     clitk::ResizeImageLike<MaskImageType>(SubclavianArtery, support, GetBackgroundValue());
 764   Thyroid =
 765     clitk::ResizeImageLike<MaskImageType>(Thyroid, support, GetBackgroundValue());
 766   Aorta =
 767     clitk::ResizeImageLike<MaskImageType>(Aorta, support, GetBackgroundValue());
 768   BrachioCephalicVein =
 769     clitk::ResizeImageLike<MaskImageType>(BrachioCephalicVein, support, GetBackgroundValue());
 770   Trachea =
 771     clitk::ResizeImageLike<MaskImageType>(Trachea, support, GetBackgroundValue());
 772
 773   // Extract slices
 774   std::vector<MaskSlicePointer> slices_BrachioCephalicArtery;
 775   clitk::ExtractSlices<MaskImageType>(BrachioCephalicArtery, 2, slices_BrachioCephalicArtery);
 776   std::vector<MaskSlicePointer> slices_BrachioCephalicVein;
 777   clitk::ExtractSlices<MaskImageType>(BrachioCephalicVein, 2, slices_BrachioCephalicVein);
 778   std::vector<MaskSlicePointer> slices_CommonCarotidArtery;
 779   clitk::ExtractSlices<MaskImageType>(CommonCarotidArtery, 2, slices_CommonCarotidArtery);
 780   std::vector<MaskSlicePointer> slices_SubclavianArtery;
 781   clitk::ExtractSlices<MaskImageType>(SubclavianArtery, 2, slices_SubclavianArtery);
 782   std::vector<MaskSlicePointer> slices_Thyroid;
 783   clitk::ExtractSlices<MaskImageType>(Thyroid, 2, slices_Thyroid);
 784   std::vector<MaskSlicePointer> slices_Aorta;
 785   clitk::ExtractSlices<MaskImageType>(Aorta, 2, slices_Aorta);
 786   std::vector<MaskSlicePointer> slices_Trachea;
 787   clitk::ExtractSlices<MaskImageType>(Trachea, 2, slices_Trachea);
 788   unsigned int n= slices_BrachioCephalicArtery.size();
 789
 790   // Get the boundaries of one slice
 791   std::vector<MaskSlicePointType> bounds;
 792   ComputeImageBoundariesCoordinates<MaskSliceType>(slices_BrachioCephalicArtery[0], bounds);
 793
 794   // For all slices, for all structures, find the centroid and build the contour
 795   // List of 3D points (for debug)
 796   std::vector<MaskImagePointType> p3D;
 797
 798   vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
 799   for(unsigned int i=0; i<n; i++) {
 800     // Labelize the slices
 801     slices_CommonCarotidArtery[i] = Labelize<MaskSliceType>(slices_CommonCarotidArtery[i],
 802                                                             GetBackgroundValue(), true, 1);
 803     slices_SubclavianArtery[i] = Labelize<MaskSliceType>(slices_SubclavianArtery[i],
 804                                                          GetBackgroundValue(), true, 1);
 805     slices_BrachioCephalicArtery[i] = Labelize<MaskSliceType>(slices_BrachioCephalicArtery[i],
 806                                                              GetBackgroundValue(), true, 1);
 807     slices_BrachioCephalicVein[i] = Labelize<MaskSliceType>(slices_BrachioCephalicVein[i],
 808                                                             GetBackgroundValue(), true, 1);
 809     slices_Thyroid[i] = Labelize<MaskSliceType>(slices_Thyroid[i],
 810                                                 GetBackgroundValue(), true, 1);
 811     slices_Aorta[i] = Labelize<MaskSliceType>(slices_Aorta[i],
 812                                               GetBackgroundValue(), true, 1);
 813
 814     // Search centroids
 815     std::vector<MaskSlicePointType> points2D;
 816     std::vector<MaskSlicePointType> centroids1;
 817     std::vector<MaskSlicePointType> centroids2;
 818     std::vector<MaskSlicePointType> centroids3;
 819     std::vector<MaskSlicePointType> centroids4;
 820     std::vector<MaskSlicePointType> centroids5;
 821     std::vector<MaskSlicePointType> centroids6;
 822     ComputeCentroids<MaskSliceType>(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1);
 823     ComputeCentroids<MaskSliceType>(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2);
 824     ComputeCentroids<MaskSliceType>(slices_BrachioCephalicArtery[i], GetBackgroundValue(), centroids3);
 825     ComputeCentroids<MaskSliceType>(slices_Thyroid[i], GetBackgroundValue(), centroids4);
 826     ComputeCentroids<MaskSliceType>(slices_Aorta[i], GetBackgroundValue(), centroids5);
 827     ComputeCentroids<MaskSliceType>(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6);
 828
 829     // BrachioCephalicVein -> when it is separated into two CCL, we
 830     // only consider the most at Right one
 831     if (centroids6.size() > 2) {
 832       if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2);
 833       else centroids6.erase(centroids6.begin()+1);
 834     }
 835
 836     // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
 837     // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
 838     if ((centroids3.size() ==1) && (centroids2.size() > 2)) {
 839       centroids6.clear();
 840     }
 841
 842     for(unsigned int j=1; j<centroids1.size(); j++) points2D.push_back(centroids1[j]);
 843     for(unsigned int j=1; j<centroids2.size(); j++) points2D.push_back(centroids2[j]);
 844     for(unsigned int j=1; j<centroids3.size(); j++) points2D.push_back(centroids3[j]);
 845     for(unsigned int j=1; j<centroids4.size(); j++) points2D.push_back(centroids4[j]);
 846     for(unsigned int j=1; j<centroids5.size(); j++) points2D.push_back(centroids5[j]);
 847     for(unsigned int j=1; j<centroids6.size(); j++) points2D.push_back(centroids6[j]);
 848
 849     // Sort by angle according to trachea centroid and vertical line,
 850     // in polar coordinates :
 851     // http://en.wikipedia.org/wiki/Polar_coordinate_system
 852     std::vector<MaskSlicePointType> centroids_trachea;
 853     ComputeCentroids<MaskSliceType>(slices_Trachea[i], GetBackgroundValue(), centroids_trachea);
 854     typedef std::pair<MaskSlicePointType, double> PointAngleType;
 855     std::vector<PointAngleType> angles;
 856     for(unsigned int j=0; j<points2D.size(); j++) {
 857       //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
 858       double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
 859       double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
 860       double angle = 0;
 861       if (x>0) angle = atan(y/x);
 862       if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
 863       if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
 864       if (x==0) {
 865         if (y>0) angle = M_PI/2.0;
 866         if (y<0) angle = -M_PI/2.0;
 867         if (y==0) angle = 0;
 868       }
 869       angle = clitk::rad2deg(angle);
 870       // Angle is [-180;180] wrt the X axis. We change the X axis to
 871       // be the vertical line, because we want to sort from Right to
 872       // Left from Post to Ant.
 873       if (angle>0)
 874         angle = (270-angle);
 875       if (angle<0) {
 876         angle = -angle-90;
 877         if (angle<0) angle = 360-angle;
 878       }
 879       PointAngleType a;
 880       a.first = points2D[j];
 881       a.second = angle;
 882       angles.push_back(a);
 883     }
 884
 885     // Do nothing if less than 2 points --> n
 886     if (points2D.size() < 3) { //continue;
 887       continue;
 888     }
 889
 890     // Sort points2D according to polar angles
 891     std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
 892     for(unsigned int j=0; j<angles.size(); j++) {
 893       points2D[j] = angles[j].first;
 894     }
 895     //    DDV(points2D, points2D.size());
 896
 897     /* When vessels are far away, we try to replace the line segment
 898        with a curved line that join the two vessels but stay
 899        approximately at the same distance from the trachea centroids
 900        than the vessels.
 901
 902        For that:
 903        - let a and c be two successive vessels centroids
 904        - id distance(a,c) < threshold, next point
 905
 906        TODO HERE
 907
 908        - compute mid position between two successive points -
 909        compute dist to trachea centroid for the 3 pts - if middle too
 910        low, add one point
 911     */
 912     std::vector<MaskSlicePointType> toadd;
 913     unsigned int index = 0;
 914     double dmax = 5;
 915     while (index<points2D.size()-1) {
 916       MaskSlicePointType a = points2D[index];
 917       MaskSlicePointType c = points2D[index+1];
 918
 919       double dac = a.EuclideanDistanceTo(c);
 920       if (dac>dmax) {
 921
 922         MaskSlicePointType b;
 923         b[0] = a[0]+(c[0]-a[0])/2.0;
 924         b[1] = a[1]+(c[1]-a[1])/2.0;
 925
 926         // Compute distance to trachea centroid
 927         MaskSlicePointType m = centroids_trachea[1];
 928         double da = m.EuclideanDistanceTo(a);
 929         double db = m.EuclideanDistanceTo(b);
 930         //double dc = m.EuclideanDistanceTo(c);
 931
 932         // Mean distance, find point on the line from trachea centroid
 933         // to b
 934         double alpha = (da+db)/2.0;
 935         MaskSlicePointType v;
 936         double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
 937         v[0] = (b[0]-m[0])/n;
 938         v[1] = (b[1]-m[1])/n;
 939         MaskSlicePointType r;
 940         r[0] = m[0]+alpha*v[0];
 941         r[1] = m[1]+alpha*v[1];
 942         points2D.insert(points2D.begin()+index+1, r);
 943       }
 944       else {
 945         index++;
 946       }
 947     }
 948     //    DDV(points2D, points2D.size());
 949
 950     // Add some points to close the contour
 951     // - H line towards Right
 952     MaskSlicePointType p = points2D[0];
 953     p[0] = bounds[0][0];
 954     points2D.insert(points2D.begin(), p);
 955     // - corner Right/Post
 956     p = bounds[0];
 957     points2D.insert(points2D.begin(), p);
 958     // - H line towards Left
 959     p = points2D.back();
 960     p[0] = bounds[2][0];
 961     points2D.push_back(p);
 962     // - corner Right/Post
 963     p = bounds[2];
 964     points2D.push_back(p);
 965     // Close contour with the first point
 966     points2D.push_back(points2D[0]);
 967     //    DDV(points2D, points2D.size());
 968
 969     // Build 3D points from the 2D points
 970     std::vector<ImagePointType> points3D;
 971     clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
 972     for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
 973
 974     // Build the mesh from the contour's points
 975     vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
 976     append->AddInput(mesh);
 977   }
 978
 979   // DEBUG: write points3D
 980   clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
 981
 982   // Build the final 3D mesh form the list 2D mesh
 983   append->Update();
 984   vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
 985
 986   // Debug, write the mesh
 987   /*
 988     vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
 989     w->SetInput(mesh);
 990     w->SetFileName("bidon.vtk");
 991     w->Write();
 992   */
 993
 994   // Compute a single binary 3D image from the list of contours
 995   clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
 996     clitk::MeshToBinaryImageFilter<MaskImageType>::New();
 997   filter->SetMesh(mesh);
 998   filter->SetLikeImage(support);
 999   filter->Update();
1000   binarizedContour = filter->GetOutput();
1001
1002   // Crop
1003   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, true, p);
1004   inf = p[2];
1005   DD(p);
1006   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour, GetForegroundValue(), 2, false, p);
1007   sup = p[2];
1008   DD(p);
1009   binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1010                                                         false, GetBackgroundValue());
1011   // Update the AFDB
1012   writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
1013   GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
1014   WriteAFDB();
1015   return binarizedContour;
1016
1017   /*
1018   // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1019   ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1020   p[1] += 50; // 50 mm Post from this point must be kept
1021   ImageIndexType index;
1022   binarizedContour->TransformPhysicalPointToIndex(p, index);
1023   bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1024
1025   // remove from support
1026   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1027   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1028   boolFilter->InPlaceOn();
1029   boolFilter->SetInput1(m_Working_Support);
1030   boolFilter->SetInput2(binarizedContour);
1031   boolFilter->SetBackgroundValue1(GetBackgroundValue());
1032   boolFilter->SetBackgroundValue2(GetBackgroundValue());
1033   if (isInside)
1034     boolFilter->SetOperationType(BoolFilterType::And);
1035   else
1036     boolFilter->SetOperationType(BoolFilterType::AndNot);
1037   boolFilter->Update();
1038   m_Working_Support = boolFilter->GetOutput();
1039   */
1040
1041   // End
1042   //StopCurrentStep<MaskImageType>(m_Working_Support);
1043   //m_ListOfStations["2R"] = m_Working_Support;
1044   //m_ListOfStations["2L"] = m_Working_Support;
1045 }
1046 //--------------------------------------------------------------------
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064 //--------------------------------------------------------------------
1065 template <class ImageType>
1066 typename clitk::ExtractLymphStationsFilter<ImageType>::MaskImagePointer
1067 clitk::ExtractLymphStationsFilter<ImageType>::
1068 FindAntPostVessels2()
1069 {
1070   // -----------------------------------------------------
1071   /* Rod says: "The anterior border, as with the Atlas – UM, is
1072     posterior to the vessels (right subclavian vein, left
1073     brachiocephalic vein, right brachiocephalic vein, left subclavian
1074     artery, left common carotid artery and brachiocephalic trunk).
1075     These vessels are not included in the nodal station.  The anterior
1076     border is drawn to the midpoint of the vessel and an imaginary
1077     line joins the middle of these vessels.  Between the vessels,
1078     station 2 is in contact with station 3a." */
1079
1080   // Check if no already done
1081   bool found = true;
1082   MaskImagePointer binarizedContour;
1083   try {
1084     DD("FindAntPostVessels try to get");
1085     binarizedContour = GetAFDB()->template GetImage <MaskImageType>("AntPostVesselsSeparation");
1086   }
1087   catch(clitk::ExceptionObject e) {
1088     DD("not found");
1089     found = false;
1090   }
1091   if (found) {
1092     return binarizedContour;
1093   }
1094
1095   /* Here, we consider the vessels form a kind of anterior barrier. We
1096      link all vessels centroids and remove what is post to it.
1097     - select the list of structure
1098             vessel1 = BrachioCephalicArtery
1099             vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right)
1100             vessel3 = CommonCarotidArtery
1101             vessel4 = SubclavianArtery
1102             other   = Thyroid
1103             other   = Aorta
1104      - crop images as needed
1105      - slice by slice, choose the CCL and extract centroids
1106      - slice by slice, sort according to polar angle wrt Trachea centroid.
1107      - slice by slice, link centroids and close contour
1108      - remove outside this contour
1109      - merge with support
1110   */
1111
1112   // Read structures
1113   std::map<std::string, MaskImagePointer> MapOfStructures;
1114   typedef std::map<std::string, MaskImagePointer>::iterator MapIter;
1115   MapOfStructures["BrachioCephalicArtery"] = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
1116   MapOfStructures["BrachioCephalicVein"] = GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicVein");
1117   MapOfStructures["CommonCarotidArteryLeft"] = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryLeft");
1118   MapOfStructures["CommonCarotidArteryRight"] = GetAFDB()->template GetImage<MaskImageType>("CommonCarotidArteryRight");
1119   MapOfStructures["SubclavianArteryLeft"] = GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryLeft");
1120   MapOfStructures["SubclavianArteryRight"] = GetAFDB()->template GetImage<MaskImageType>("SubclavianArteryRight");
1121   MapOfStructures["Thyroid"] = GetAFDB()->template GetImage<MaskImageType>("Thyroid");
1122   MapOfStructures["Aorta"] = GetAFDB()->template GetImage<MaskImageType>("Aorta");
1123   MapOfStructures["Trachea"] = GetAFDB()->template GetImage<MaskImageType>("Trachea");
1124
1125   std::vector<std::string> ListOfStructuresNames;
1126
1127   // Create a temporay support
1128   // From first slice of BrachioCephalicVein to end of 3A
1129   MaskImagePointer support = GetAFDB()->template GetImage<MaskImageType>("Support_Sup_Carina");
1130   MaskImagePointType p;
1131   p[0] = p[1] = p[2] =  0.0; // to avoid warning
1132   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(MapOfStructures["BrachioCephalicVein"],
1133                                                           GetBackgroundValue(), 2, true, p);
1134   double inf = p[2];
1135   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(GetAFDB()->template GetImage<MaskImageType>("Support_S3A"),
1136                                                           GetBackgroundValue(), 2, false, p);
1137   double sup = p[2]+support->GetSpacing()[2];//one slice more ?
1138   support = clitk::CropImageAlongOneAxis<MaskImageType>(support, 2, inf, sup,
1139                                                         false, GetBackgroundValue());
1140
1141   // Resize all structures like support
1142   for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1143     it->second = clitk::ResizeImageLike<MaskImageType>(it->second, support, GetBackgroundValue());
1144   }
1145
1146   // Extract slices
1147   typedef std::vector<MaskSlicePointer> SlicesType;
1148   std::map<std::string, SlicesType> MapOfSlices;
1149   for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1150     SlicesType s;
1151     clitk::ExtractSlices<MaskImageType>(it->second, 2, s);
1152     MapOfSlices[it->first] = s;
1153   }
1154
1155   unsigned int n= MapOfSlices["Trachea"].size();
1156
1157   // Get the boundaries of one slice
1158   std::vector<MaskSlicePointType> bounds;
1159   ComputeImageBoundariesCoordinates<MaskSliceType>(MapOfSlices["Trachea"][0], bounds);
1160
1161   // LOOP ON SLICES
1162   // For all slices, for all structures, find the centroid and build the contour
1163   // List of 3D points (for debug)
1164   std::vector<MaskImagePointType> p3D;
1165   vtkSmartPointer<vtkAppendPolyData> append = vtkSmartPointer<vtkAppendPolyData>::New();
1166   for(unsigned int i=0; i<n; i++) {
1167
1168     // Labelize the slices
1169     for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1170       MaskSlicePointer & s = MapOfSlices[it->first][i];
1171       s = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
1172     }
1173
1174     // Search centroids
1175     std::vector<MaskSlicePointType> points2D;
1176     typedef std::vector<MaskSlicePointType> CentroidsType;
1177     std::map<std::string, CentroidsType> MapOfCentroids;
1178     for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1179       std::string structure = it->first;
1180       MaskSlicePointer & s = MapOfSlices[structure][i];
1181       CentroidsType c;
1182       clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), c);
1183       MapOfCentroids[structure] = c;
1184     }
1185
1186
1187     // BrachioCephalicVein -> when it is separated into two CCL, we
1188     // only consider the most at Right one
1189     CentroidsType & c = MapOfCentroids["BrachioCephalicVein"];
1190     if (c.size() > 2) {
1191       if (c[1][0] < c[2][0]) c.erase(c.begin()+2);
1192       else c.erase(c.begin()+1);
1193     }
1194
1195     /*
1196     // BrachioCephalicVein -> when SubclavianArtery has 2 CCL
1197     // (BrachioCephalicArtery is divided) -> forget BrachioCephalicVein
1198     if ((MapOfCentroids["BrachioCephalicArtery"].size() ==1)
1199         && (MapOfCentroids["SubclavianArtery"].size() > 2)) {
1200       MapOfCentroids["BrachioCephalicVein"].clear();
1201     }
1202     */
1203
1204     // Add all 2D points
1205     for (MapIter it = MapOfStructures.begin(); it != MapOfStructures.end(); ++it) {
1206       std::string structure = it->first;
1207       if (structure != "Trachea") { // do not add centroid of trachea
1208         CentroidsType & c = MapOfCentroids[structure];
1209         for(unsigned int j=1; j<c.size(); j++) points2D.push_back(c[j]);
1210       }
1211     }
1212
1213     // Sort by angle according to trachea centroid and vertical line,
1214     // in polar coordinates :
1215     // http://en.wikipedia.org/wiki/Polar_coordinate_system
1216     //    std::vector<MaskSlicePointType> centroids_trachea;
1217     //ComputeCentroids<MaskSliceType>(MapOfSlices["Trachea"][i], GetBackgroundValue(), centroids_trachea);
1218     CentroidsType centroids_trachea = MapOfCentroids["Trachea"];
1219     typedef std::pair<MaskSlicePointType, double> PointAngleType;
1220     std::vector<PointAngleType> angles;
1221     for(unsigned int j=0; j<points2D.size(); j++) {
1222       //double r = centroids_trachea[1].EuclideanDistanceTo(points2D[j]);
1223       double x = (points2D[j][0]-centroids_trachea[1][0]); // X : Right to Left
1224       double y = (centroids_trachea[1][1]-points2D[j][1]); // Y : Post to Ant
1225       double angle = 0;
1226       if (x>0) angle = atan(y/x);
1227       if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI;
1228       if ((x<0) && (y<0)) angle = atan(y/x)-M_PI;
1229       if (x==0) {
1230         if (y>0) angle = M_PI/2.0;
1231         if (y<0) angle = -M_PI/2.0;
1232         if (y==0) angle = 0;
1233       }
1234       angle = clitk::rad2deg(angle);
1235       // Angle is [-180;180] wrt the X axis. We change the X axis to
1236       // be the vertical line, because we want to sort from Right to
1237       // Left from Post to Ant.
1238       if (angle>0)
1239         angle = (270-angle);
1240       if (angle<0) {
1241         angle = -angle-90;
1242         if (angle<0) angle = 360-angle;
1243       }
1244       PointAngleType a;
1245       a.first = points2D[j];
1246       a.second = angle;
1247       angles.push_back(a);
1248     }
1249
1250     // Do nothing if less than 2 points --> n
1251     if (points2D.size() < 3) { //continue;
1252       continue;
1253     }
1254
1255     // Sort points2D according to polar angles
1256     std::sort(angles.begin(), angles.end(), comparePointsWithAngle<PointAngleType>());
1257     for(unsigned int j=0; j<angles.size(); j++) {
1258       points2D[j] = angles[j].first;
1259     }
1260     //    DDV(points2D, points2D.size());
1261
1262     /* When vessels are far away, we try to replace the line segment
1263        with a curved line that join the two vessels but stay
1264        approximately at the same distance from the trachea centroids
1265        than the vessels.
1266
1267        For that:
1268        - let a and c be two successive vessels centroids
1269        - id distance(a,c) < threshold, next point
1270
1271        TODO HERE
1272
1273        - compute mid position between two successive points -
1274        compute dist to trachea centroid for the 3 pts - if middle too
1275        low, add one point
1276     */
1277     std::vector<MaskSlicePointType> toadd;
1278     unsigned int index = 0;
1279     double dmax = 5;
1280     while (index<points2D.size()-1) {
1281       MaskSlicePointType a = points2D[index];
1282       MaskSlicePointType c = points2D[index+1];
1283
1284       double dac = a.EuclideanDistanceTo(c);
1285       if (dac>dmax) {
1286
1287         MaskSlicePointType b;
1288         b[0] = a[0]+(c[0]-a[0])/2.0;
1289         b[1] = a[1]+(c[1]-a[1])/2.0;
1290
1291         // Compute distance to trachea centroid
1292         MaskSlicePointType m = centroids_trachea[1];
1293         double da = m.EuclideanDistanceTo(a);
1294         double db = m.EuclideanDistanceTo(b);
1295         //double dc = m.EuclideanDistanceTo(c);
1296
1297         // Mean distance, find point on the line from trachea centroid
1298         // to b
1299         double alpha = (da+db)/2.0;
1300         MaskSlicePointType v;
1301         double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2));
1302         v[0] = (b[0]-m[0])/n;
1303         v[1] = (b[1]-m[1])/n;
1304         MaskSlicePointType r;
1305         r[0] = m[0]+alpha*v[0];
1306         r[1] = m[1]+alpha*v[1];
1307         points2D.insert(points2D.begin()+index+1, r);
1308       }
1309       else {
1310         index++;
1311       }
1312     }
1313     //    DDV(points2D, points2D.size());
1314
1315     // Add some points to close the contour
1316     // - H line towards Right
1317     MaskSlicePointType p = points2D[0];
1318     p[0] = bounds[0][0];
1319     points2D.insert(points2D.begin(), p);
1320     // - corner Right/Post
1321     p = bounds[0];
1322     points2D.insert(points2D.begin(), p);
1323     // - H line towards Left
1324     p = points2D.back();
1325     p[0] = bounds[2][0];
1326     points2D.push_back(p);
1327     // - corner Right/Post
1328     p = bounds[2];
1329     points2D.push_back(p);
1330     // Close contour with the first point
1331     points2D.push_back(points2D[0]);
1332     //    DDV(points2D, points2D.size());
1333
1334     // Build 3D points from the 2D points
1335     std::vector<ImagePointType> points3D;
1336     clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i, support, points3D);
1337     for(unsigned int x=0; x<points3D.size(); x++) p3D.push_back(points3D[x]);
1338
1339     // Build the mesh from the contour's points
1340     vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1341     append->AddInput(mesh);
1342     // if (i ==n-1) { // last slice
1343     //   clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(points2D, i+1, support, points3D);
1344     //   vtkSmartPointer<vtkPolyData> mesh = Build3DMeshFrom2DContour(points3D);
1345     //   append->AddInput(mesh);
1346     // }
1347   }
1348
1349   // DEBUG: write points3D
1350   clitk::WriteListOfLandmarks<MaskImageType>(p3D, "vessels-centroids.txt");
1351
1352   // Build the final 3D mesh form the list 2D mesh
1353   append->Update();
1354   vtkSmartPointer<vtkPolyData> mesh = append->GetOutput();
1355
1356   // Debug, write the mesh
1357   /*
1358   vtkSmartPointer<vtkPolyDataWriter> w = vtkSmartPointer<vtkPolyDataWriter>::New();
1359   w->SetInput(mesh);
1360   w->SetFileName("bidon.vtk");
1361   w->Write();
1362   */
1363
1364   // Compute a single binary 3D image from the list of contours
1365   clitk::MeshToBinaryImageFilter<MaskImageType>::Pointer filter =
1366     clitk::MeshToBinaryImageFilter<MaskImageType>::New();
1367   filter->SetMesh(mesh);
1368   filter->SetLikeImage(support);
1369   filter->Update();
1370   binarizedContour = filter->GetOutput();
1371
1372   // Crop
1373   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1374                                                           GetForegroundValue(), 2, true, p);
1375   inf = p[2];
1376   clitk::FindExtremaPointInAGivenDirection<MaskImageType>(binarizedContour,
1377                                                           GetForegroundValue(), 2, false, p);
1378   sup = p[2]+binarizedContour->GetSpacing()[2]; // extend to include the last slice
1379   binarizedContour = clitk::CropImageAlongOneAxis<MaskImageType>(binarizedContour, 2, inf, sup,
1380                                                         false, GetBackgroundValue());
1381
1382   // Update the AFDB
1383   writeImage<MaskImageType>(binarizedContour, "seg/AntPostVesselsSeparation.mhd");
1384   GetAFDB()->SetImageFilename("AntPostVesselsSeparation", "seg/AntPostVesselsSeparation.mhd");
1385   WriteAFDB();
1386   return binarizedContour;
1387
1388   /*
1389   // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse
1390   ImagePointType p = p3D[2]; // This is the first centroid of the first slice
1391   p[1] += 50; // 50 mm Post from this point must be kept
1392   ImageIndexType index;
1393   binarizedContour->TransformPhysicalPointToIndex(p, index);
1394   bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue());
1395
1396   // remove from support
1397   typedef clitk::BooleanOperatorLabelImageFilter<MaskImageType> BoolFilterType;
1398   typename BoolFilterType::Pointer boolFilter = BoolFilterType::New();
1399   boolFilter->InPlaceOn();
1400   boolFilter->SetInput1(m_Working_Support);
1401   boolFilter->SetInput2(binarizedContour);
1402   boolFilter->SetBackgroundValue1(GetBackgroundValue());
1403   boolFilter->SetBackgroundValue2(GetBackgroundValue());
1404   if (isInside)
1405     boolFilter->SetOperationType(BoolFilterType::And);
1406   else
1407     boolFilter->SetOperationType(BoolFilterType::AndNot);
1408   boolFilter->Update();
1409   m_Working_Support = boolFilter->GetOutput();
1410   */
1411
1412   // End
1413   //StopCurrentStep<MaskImageType>(m_Working_Support);
1414   //m_ListOfStations["2R"] = m_Working_Support;
1415   //m_ListOfStations["2L"] = m_Working_Support;
1416 }
1417 //--------------------------------------------------------------------
1418
1419
1420
1421 #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX