segmentation/clitkExtractLymphStationsFilter.txx

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