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