segmentation/clitkExtractMediastinalVesselsFilter.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 "clitkExtractMediastinalVesselsFilter.h"
  25 #include "clitkSegmentationUtils.h"
  26 #include "clitkReconstructWithConditionalGrayscaleDilateImageFilter.h"
  27
  28 // itk
  29 #include <itkBinaryThresholdImageFilter.h>
  30 #include <itkMinimumMaximumImageCalculator.h>
  31
  32 template<class T> struct index_cmp {
  33   index_cmp(const T varr) : arr(varr) {}
  34   bool operator()(const size_t a, const size_t b) const
  35   { return arr[a] < arr[b]; }
  36   const T arr;
  37 };
  38
  39
  40 //--------------------------------------------------------------------
  41 template <class TImageType>
  42 clitk::ExtractMediastinalVesselsFilter<TImageType>::
  43 ExtractMediastinalVesselsFilter():
  44   clitk::FilterBase(),
  45   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  46   itk::ImageToImageFilter<TImageType, MaskImageType>()
  47 {
  48   this->SetNumberOfRequiredInputs(1);
  49   SetBackgroundValue(0);
  50   SetForegroundValue(1);
  51   SetThresholdHigh(140);
  52   SetThresholdLow(55);
  53   SetErosionRadius(2);
  54   SetDilatationRadius(9);
  55   SetMaxDistancePostToCarina(10);
  56   SetMaxDistanceAntToCarina(40);
  57   SetMaxDistanceLeftToCarina(35);
  58   SetMaxDistanceRightToCarina(35);
  59   SetSoughtVesselSeedName("NoSeedNameGiven");
  60   SetFinalOpeningRadius(1);
  61   VerboseTrackingFlagOff();
  62 }
  63 //--------------------------------------------------------------------
  64
  65
  66 //--------------------------------------------------------------------
  67 template <class TImageType>
  68 void
  69 clitk::ExtractMediastinalVesselsFilter<TImageType>::
  70 GenerateOutputInformation() {
  71   // Get inputs
  72   LoadAFDB();
  73   m_Input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
  74
  75   // ------------------------------------------------------------------
  76   // Crop the initial image superiorly and inferiorly.
  77   // TODO : add options for x cm above/below
  78   CropInputImage();
  79
  80   // ------------------------------------------------------------------
  81   // Binarize the image. Need two thresholds, one high to select
  82   // structures (CCL) that are almost not connected (after erosion),
  83   // and one low thresholds to select the real contours. Will be
  84   // reconstructed later.
  85   StartNewStep("Binarize with high threshold = "+toString(GetThresholdHigh()));
  86   typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> BinarizeFilterType;
  87   typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
  88   binarizeFilter->SetInput(m_Input);
  89   binarizeFilter->SetLowerThreshold(GetThresholdHigh());
  90   binarizeFilter->SetInsideValue(GetForegroundValue());
  91   binarizeFilter->SetOutsideValue(GetBackgroundValue());
  92   binarizeFilter->Update();
  93   m_Mask = binarizeFilter->GetOutput();
  94   StopCurrentStep<MaskImageType>(m_Mask);
  95
  96   // ------------------------------------------------------------------
  97   StartNewStep("Binarize with low threshold = "+toString(GetThresholdLow()));
  98   binarizeFilter = BinarizeFilterType::New();
  99   binarizeFilter->SetInput(m_Input);
 100   binarizeFilter->SetLowerThreshold(GetThresholdLow());
 101   binarizeFilter->SetInsideValue(GetForegroundValue());
 102   binarizeFilter->SetOutsideValue(GetBackgroundValue());
 103   binarizeFilter->Update();
 104   MaskImagePointer m_Mask2 = binarizeFilter->GetOutput();
 105   StopCurrentStep<MaskImageType>(m_Mask2);
 106
 107   // ------------------------------------------------------------------
 108   // Extract slices
 109   StartNewStep("Detect objects : erosion, then slice by slice reconstruction");
 110   std::vector<MaskSlicePointer> slices_mask;
 111   clitk::ExtractSlices<MaskImageType>(m_Mask, 2, slices_mask);
 112   std::vector<MaskSlicePointer> slices_mask2;
 113   clitk::ExtractSlices<MaskImageType>(m_Mask2, 2, slices_mask2);
 114   int radius = GetErosionRadius();
 115
 116   // List of working slices (debug only)
 117   std::vector<MaskSlicePointer> debug_eroded;
 118   // std::vector<MaskSlicePointer> debug_labeled;
 119
 120   // ------------------------------------------------------------------
 121   // Loop Slice by Slice in order to break connectivity between
 122   // CCL. Erode and reconstruct all labels at the same time without
 123   // merging them.
 124   for(uint i=0; i<slices_mask.size(); i++) {
 125
 126     /*// Erosion kernel
 127       typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,2> KernelType;
 128       KernelType structuringElement;
 129       structuringElement.SetRadius(radius);
 130       structuringElement.CreateStructuringElement();
 131     */
 132
 133     // Erosion -> we break the connectivity between structure
 134     MaskSliceType::SizeType r;
 135     r[0] = r[1] = radius;
 136     MaskSlicePointer eroded = clitk::Opening<MaskSliceType>(slices_mask[i],
 137                                                             r,
 138                                                             GetBackgroundValue(),
 139                                                             GetForegroundValue());
 140     /*
 141     //typedef itk::BinaryErodeImageFilter<MaskSliceType, MaskSliceType, KernelType> ErodeFilterType;
 142     typedef itk::BinaryMorphologicalOpeningImageFilter<MaskSliceType, MaskSliceType, KernelType> ErodeFilterType;
 143     typename ErodeFilterType::Pointer eroder = ErodeFilterType::New();
 144     eroder->SetInput(slices_mask[i]);
 145     eroder->SetBackgroundValue(GetBackgroundValue());
 146     eroder->SetForegroundValue(GetForegroundValue());
 147     //eroder->SetBoundaryToForeground(true); // ?? for BinaryErodeImageFilter
 148     eroder->SetKernel(structuringElement);
 149     eroder->Update();
 150     MaskSlicePointer eroded = eroder->GetOutput();
 151     */
 152
 153     // Keep slice for debug
 154     debug_eroded.push_back(eroded);
 155
 156     // Labelize (CCL)
 157     MaskSlicePointer labeled =
 158       clitk::Labelize<MaskSliceType>(eroded, GetBackgroundValue(), true, 1); // Fully connected !
 159     // debug_labeled.push_back(labeled);
 160
 161     // Make Reconstruction filter : dilation all labels at the same
 162     // time, prevent to merge them.
 163     typedef clitk::ReconstructWithConditionalGrayscaleDilateImageFilter<MaskSliceType> ReconstructFilterType;
 164     typename ReconstructFilterType::Pointer reconstructor = ReconstructFilterType::New();
 165     reconstructor->SetInput(labeled);
 166     reconstructor->SetIterationNumber(radius+GetDilatationRadius());
 167     reconstructor->Update();
 168     MaskSlicePointer s = reconstructor->GetOutput();
 169
 170     // Remove Initial BG of the second tresholded image
 171     s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, slices_mask2[i],
 172                                                            GetBackgroundValue(), GetBackgroundValue(), true);
 173     m_slice_recon.push_back(s);
 174
 175   } // end loop
 176
 177   // Build 3D images from the slice by slice processing
 178   MaskImageType::Pointer eroded = clitk::JoinSlices<MaskImageType>(debug_eroded, m_Mask, 2);
 179   writeImage<MaskImageType>(eroded, "erode.mhd");
 180   //MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);
 181   MaskImageType::Pointer r = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
 182   writeImage<MaskImageType>(r, "recon1.mhd");
 183
 184   // ------------------------------------------------------------------
 185   // Track the SoughtVessel from the given first point
 186   // superiorly. This is done by TrackBifurcationFromPoint
 187   MaskImagePointType SoughtVesselSeedPoint;
 188   GetAFDB()->GetPoint3D(m_SoughtVesselSeedName, SoughtVesselSeedPoint);
 189   MaskImagePointType MaxSlicePoint;
 190   if (GetAFDB()->TagExist(m_SoughtVesselSeedName+"Max")) {
 191     GetAFDB()->GetPoint3D(m_SoughtVesselSeedName+"Max", MaxSlicePoint);
 192   }
 193   else {
 194     MaxSlicePoint = SoughtVesselSeedPoint;
 195     MaxSlicePoint[2] += 1000;
 196   }
 197
 198   // Find the label with the maximum value to set the result
 199   typedef itk::MinimumMaximumImageCalculator<MaskImageType> MinMaxFilterType;
 200   MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
 201   ff->SetImage(r);
 202   ff->ComputeMaximum();
 203   LabelType newLabel = ff->GetMaximum()+1;
 204
 205   // the following bifurcations point will the centroids of the
 206   // components obtain when (hopefully!) the SoughtVessel
 207   // split into CommonArtery and SubclavianArtery.
 208   std::vector<MaskImagePointType> bifurcations;
 209   //  TrackBifurcationFromPoint(r, m_slice_recon, SoughtVesselSeedPoint,
 210   //                         MaxSlicePoint, newLabel, bifurcations);
 211
 212   TrackVesselsFromPoint(r, m_slice_recon, SoughtVesselSeedPoint,
 213                         MaxSlicePoint, newLabel);
 214
 215   // Build the final 3D image from the previous slice by slice processing
 216   m_SoughtVessel = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
 217   // writeImage<MaskImageType>(m_SoughtVessel, "recon2.mhd");
 218
 219   // Set binary image, (remove other labels).
 220   m_SoughtVessel =
 221     clitk::Binarize<MaskImageType>(m_SoughtVessel, newLabel, newLabel,
 222                                    GetBackgroundValue(), GetForegroundValue());
 223
 224   //  writeImage<MaskImageType>(m_SoughtVessel, "afterbinarize.mhd");
 225   m_SoughtVessel = clitk::AutoCrop<MaskImageType>(m_SoughtVessel, GetBackgroundValue());
 226
 227   //  writeImage<MaskImageType>(m_SoughtVessel, "afterautocrop.mhd");
 228
 229   // Clean the image : Opening (not in Z direction)
 230   typename MaskImageType::SizeType rad;
 231   rad[0] = rad[1] = GetFinalOpeningRadius();
 232   rad[2] = 0;
 233   m_SoughtVessel = clitk::Opening<MaskImageType>(m_SoughtVessel, rad,
 234                                                  GetBackgroundValue(), GetForegroundValue());
 235
 236   //  writeImage<MaskImageType>(m_SoughtVessel, "afteropen.mhd");
 237
 238   // Clean the image : keep main CCL slice by slice
 239   m_SoughtVessel = clitk::SliceBySliceKeepMainCCL<MaskImageType>(m_SoughtVessel,
 240                                                                  GetBackgroundValue(),
 241                                                                  GetForegroundValue());
 242 }
 243 //--------------------------------------------------------------------
 244
 245
 246 //--------------------------------------------------------------------
 247 template <class TImageType>
 248 void
 249 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 250 GenerateInputRequestedRegion() {
 251   //DD("GenerateInputRequestedRegion (nothing?)");
 252 }
 253 //--------------------------------------------------------------------
 254
 255
 256 //--------------------------------------------------------------------
 257 template <class TImageType>
 258 void
 259 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 260 GenerateData() {
 261   // Save in the AFDB (not write on the disk here)
 262   GetAFDB()->SetImageFilename(GetSoughtVesselName(), GetOutputFilename());
 263   WriteAFDB();
 264   // Final Step -> graft output
 265   this->GraftNthOutput(0, m_SoughtVessel);
 266 }
 267 //--------------------------------------------------------------------
 268
 269
 270 //--------------------------------------------------------------------
 271 template <class TImageType>
 272 void
 273 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 274 CropInputImage() {
 275   StartNewStep("Crop the input image: SI,AP limits with carina and crop with mediastinum");
 276   /*
 277     Need : Trachea, Carina (roi not point),
 278   */
 279   // Get Trachea and Carina
 280   MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
 281
 282   // Compute Carina position
 283   double m_CarinaZ;
 284   MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
 285   std::vector<MaskImagePointType> centroids;
 286   clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
 287   m_CarinaZ = centroids[1][2];
 288   // add one slice to include carina
 289   m_CarinaZ += Carina->GetSpacing()[2];
 290   // We dont need Carina structure from now
 291   Carina->Delete();
 292   GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
 293
 294   // Crop Inf, remove below Carina
 295   m_Input =
 296     clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 2, m_CarinaZ, false, GetBackgroundValue());
 297
 298   // // Get seed, crop around
 299   // MaskImagePointType SoughtVesselSeedPoint;
 300   // GetAFDB()->GetPoint3D(m_SoughtVesselSeedName, SoughtVesselSeedPoint);
 301
 302   // Crop post
 303   double m_CarinaY = centroids[1][1];
 304   m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 1,
 305                                                          m_CarinaY+GetMaxDistancePostToCarina(),
 306                                                          false, GetBackgroundValue());
 307   // Crop ant
 308   m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 1,
 309                                                        m_CarinaY-GetMaxDistanceAntToCarina(),
 310                                                        false, GetBackgroundValue());
 311   // Crop Right
 312   double m_CarinaX = centroids[1][0];
 313   m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 0,
 314                                                        m_CarinaX-GetMaxDistanceRightToCarina(),
 315                                                        false, GetBackgroundValue());
 316   // Crop Left
 317   m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 0,
 318                                                          m_CarinaX+GetMaxDistanceLeftToCarina(),
 319                                                          false, GetBackgroundValue());
 320
 321   /*
 322   // AutoCrop with Mediastinum, generaly only allow to remove few slices (superiorly)
 323   m_Mediastinum  = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
 324   // Resize like input (sup to carina)
 325   m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Input, GetBackgroundValue());
 326   // Auto crop
 327   m_Mediastinum = clitk::AutoCrop<MaskImageType>(m_Mediastinum, GetBackgroundValue());
 328   // Resize input
 329   m_Input = clitk::ResizeImageLike<ImageType>(m_Input, m_Mediastinum, GetBackgroundValue());
 330   */
 331
 332   //  writeImage<ImageType>(m_Input, "crop.mhd");
 333   // End
 334   StopCurrentStep<ImageType>(m_Input);
 335 }
 336 //--------------------------------------------------------------------
 337
 338
 339 //--------------------------------------------------------------------
 340 template <class TImageType>
 341 void
 342 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 343 TrackBifurcationFromPoint(MaskImagePointer & recon,
 344                           std::vector<MaskSlicePointer> & slices_recon,
 345                           MaskImagePointType point3D,
 346                           MaskImagePointType pointMaxSlice,
 347                           LabelType newLabel,
 348                           std::vector<MaskImagePointType> & bifurcations) {
 349   StartNewStep("Track the SoughtVessel from the seed point");
 350
 351   // Find first slice index
 352   MaskImageIndexType index;
 353   recon->TransformPhysicalPointToIndex(point3D, index);
 354   int numberOfBifurcation = 0;
 355   typedef typename MaskSliceType::PointType SlicePointType;
 356   SlicePointType previousCenter;
 357
 358   // Max slice
 359   MaskImageIndexType indexMaxSlice;
 360   recon->TransformPhysicalPointToIndex(pointMaxSlice, indexMaxSlice);
 361   uint maxSlice = indexMaxSlice[2];
 362
 363   // Get current label at the point3D of interest
 364   uint currentSlice=index[2];
 365   bool found = false;
 366   LabelType previous_slice_label=recon->GetPixel(index);
 367
 368
 369   if (GetVerboseTrackingFlag()) {
 370     std::cout << "TrackBifurcationFromPoint " << std::endl;
 371     std::cout << "\t point3D = " << point3D << std::endl;
 372     std::cout << "\t pointMaxSlice = " << pointMaxSlice << std::endl;
 373     std::cout << "\t newLabel = " << newLabel << std::endl;
 374   }
 375
 376   //  DD(slices_recon.size());
 377   do {
 378     if (GetVerboseTrackingFlag()) {
 379       std::cout << "currentSlice = " << currentSlice << std::endl;
 380     }
 381     // Consider current reconstructed slice
 382     MaskSlicePointer s = slices_recon[currentSlice];
 383     MaskSlicePointer previous;
 384     if (currentSlice == index[2]) previous = s;
 385     else {
 386       previous = slices_recon[currentSlice-1];
 387     }
 388
 389     // Get centroids of the labels in the current slice
 390     static const unsigned int Dim = MaskSliceType::ImageDimension;
 391     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 392     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 393     typedef itk::LabelImageToLabelMapFilter<MaskSliceType, LabelMapType> ImageToMapFilterType;
 394     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 395     typedef itk::ShapeLabelMapFilter<LabelMapType, MaskSliceType> ShapeFilterType;
 396     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 397     imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
 398     imageToLabelFilter->SetInput(s);
 399     statFilter->SetInput(imageToLabelFilter->GetOutput());
 400     // statFilter->SetComputeFeretDiameter( true );
 401     statFilter->ComputePerimeterOn(); // To be able to get proper Roundness value
 402     statFilter->Update();
 403     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 404
 405     // Look what centroid inside the previous largest one
 406     std::vector<SlicePointType> centroids;
 407     std::vector<LabelType> centroids_label;
 408     std::vector<double> labels_size;
 409     for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
 410       int label = labelMap->GetLabels()[c];
 411       //      DD(label);
 412       SlicePointType center = labelMap->GetLabelObject(label)->GetCentroid();
 413       //      DD(center);
 414       // Get label into previous slice
 415       typename MaskSliceType::IndexType centerIndex;
 416       previous->TransformPhysicalPointToIndex(center, centerIndex);
 417       LabelType labelInPreviousSlice = previous->GetPixel(centerIndex);
 418       // if this current centroid was in the current label, add it to bifurcations
 419       if (labelInPreviousSlice == previous_slice_label) {
 420         centroids.push_back(center);
 421         centroids_label.push_back(label);
 422         labels_size.push_back(labelMap->GetLabelObject(label)->GetPhysicalSize());
 423         //DD(labels_size.back());
 424         //DD(labelMap->GetLabelObject(label)->GetRoundness());
 425         // previousCenter = centroids.back();
 426       }
 427     }
 428
 429     // -------------------------
 430     // If no centroid were found
 431     if (centroids.size() == 0) {
 432       if (GetVerboseTrackingFlag()) {
 433         std::cout << "no centroid" << std::endl;
 434       }
 435        // Last attempt to find -> check if previous centroid is inside a CCL
 436       //      if in s -> get value, getcentroid add.
 437       //      DD(currentSlice);
 438       //DD("Last change to find");
 439       //DD(previous_slice_label);
 440       //DD(newLabel);
 441       typename MaskSliceType::IndexType previousCenterIndex;
 442       s->TransformPhysicalPointToIndex(previousCenter, previousCenterIndex);
 443       //DD(previousCenter);
 444       LabelType labelInSlice = s->GetPixel(previousCenterIndex);
 445       //DD(labelInSlice);
 446       if (labelInSlice != GetBackgroundValue()) {
 447         centroids.push_back(labelMap->GetLabelObject(labelInSlice)->GetCentroid());
 448         centroids_label.push_back(labelInSlice);
 449         labels_size.push_back(labelMap->GetLabelObject(labelInSlice)->GetPhysicalSize());
 450       }
 451     }
 452
 453     // Some centroid were found
 454     // If several centroids, we found a bifurcation
 455     if (centroids.size() > 1) {
 456       if (GetVerboseTrackingFlag()) {
 457         std::cout << "Centroid" << centroids.size() << std::endl;
 458       }
 459       //      int n = centroids.size();
 460       // Debug point
 461       std::vector<ImagePointType> d;
 462       clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(centroids, currentSlice, m_Mask, d);
 463       //      DDV(d, d.size());
 464
 465       /*
 466       // try one or all centroids
 467       std::vector<
 468       for(uint a<=0; a<nb++; a++) {
 469       DD(a);
 470       // Create the list of candidates
 471       std::vector<int> c;
 472       if (a==nb) { for(uint x=0; x<nb; x++) c.push_back(x); }
 473       else c.push_back(a);
 474       DD(a.size());
 475
 476       // Test size
 477       double size=0.0;
 478       for(uint x=0; x<c.size(); c++) { size += labels_size[c[x]]; }
 479       DD(size);
 480       */
 481
 482
 483       // new potential bifurcation found
 484       numberOfBifurcation++;
 485       // If the number of bifurcation is greater than the required one, we stop
 486       if (numberOfBifurcation > GetMaxNumberOfFoundBifurcation()) {
 487         found = true;
 488         //DD("max bif reach");
 489         for(uint c=0; c<centroids.size(); c++) {
 490           ImagePointType bif;
 491           clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, currentSlice, bif);
 492           bifurcations.push_back(bif);
 493         }
 494       }
 495       // Else we continue along the main (largest) connected component
 496       else {
 497         int indexOfLargest = 0;
 498         for(uint b=0; b<centroids.size(); b++) {
 499           if (labels_size[b] > labels_size[indexOfLargest]) {
 500             indexOfLargest = b;
 501           }
 502         }
 503         //DD(indexOfLargest);
 504         //DD(labels_size[indexOfLargest]);
 505         SlicePointType c = centroids[indexOfLargest];
 506         LabelType l = centroids_label[indexOfLargest];
 507         //DD(l);
 508         //DD(c);
 509         centroids.clear();
 510         centroids.push_back(c);
 511         centroids_label.clear();
 512         centroids_label.push_back(l);
 513       }
 514     }
 515
 516
 517     /*    ==> here all centroid are considered as ok.*/
 518
 519     // REMOVE IF CENTROID=1, REPLACE BY >0
 520
 521     // if only one centroids, we change the current image with the current label
 522     if (centroids.size() == 1) {
 523       if (GetVerboseTrackingFlag()) {
 524         std::cout << "Centroid" << centroids.size() << std::endl;
 525       }
 526
 527       //DD(centroids_label[0]);
 528       s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
 529       slices_recon[currentSlice] = s;
 530       previous_slice_label = newLabel;
 531       // It can happend that several CCL share this same label. To
 532       // prevent this case, we only consider the one that contains
 533       // the centroid.
 534       MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
 535       //      writeImage<MaskSliceType>(temp, "bin-"+toString(currentSlice)+".mhd");
 536       temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
 537       //writeImage<MaskSliceType>(temp, "label-"+toString(currentSlice)+".mhd");
 538       typename MaskSliceType::IndexType centroids_index;
 539       temp->TransformPhysicalPointToIndex(centroids[0], centroids_index);
 540       typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
 541
 542       // It can happend that the centroid is inside the BG, so we keep
 543       // the largest CCL (the first);
 544       if (v == GetBackgroundValue()) {
 545         //        DD(currentSlice);
 546         //        DD("inside BG");
 547         //        DD(centroids[0]);
 548         v = 1; // largest one
 549       }
 550
 551       //DD(v);
 552       temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
 553       //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
 554       s = temp;
 555       slices_recon[currentSlice] = s;
 556
 557       // I need to recompute the centroid if we have removed some
 558       // connected component.
 559       clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), centroids);
 560       previousCenter = centroids[1];
 561     }
 562
 563     if (GetVerboseTrackingFlag()) {
 564         std::cout << "End iteration c=" << centroids.size() << std::endl;
 565     }
 566
 567     if (centroids.size() == 0) {
 568       //      DD("ZERO");
 569       found = true;
 570     }
 571
 572     if (currentSlice == slices_recon.size()-1) {
 573       //      DD("end of slices");
 574       found = true;
 575     }
 576
 577     if (currentSlice == maxSlice) {
 578       //      DD("end max slice");
 579       found = true;
 580     }
 581
 582     // iterate
 583     ++currentSlice;
 584   } while (!found);
 585
 586   // End
 587   StopCurrentStep();
 588 }
 589 //--------------------------------------------------------------------
 590
 591
 592 //--------------------------------------------------------------------
 593 template <class TImageType>
 594 void
 595 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 596 TrackVesselsFromPoint(MaskImagePointer & recon,
 597                       std::vector<MaskSlicePointer> & slices,
 598                       MaskImagePointType seedPoint,
 599                       MaskImagePointType pointMaxSlice,
 600                       LabelType newLabel) {
 601   StartNewStep("Track the SoughtVessel from the seed point");
 602
 603   // Find first slice index
 604   MaskImageIndexType seedIndex;
 605   recon->TransformPhysicalPointToIndex(seedPoint, seedIndex);
 606   //  int numberOfBifurcation = 0;
 607   typedef typename MaskSliceType::PointType SlicePointType;
 608   SlicePointType previousCentroid;
 609   previousCentroid[0] = seedPoint[0];
 610   previousCentroid[1] = seedPoint[1];
 611
 612   // Max slice
 613   MaskImageIndexType indexMaxSlice;
 614   recon->TransformPhysicalPointToIndex(pointMaxSlice, indexMaxSlice);
 615   uint maxSlice = std::min((uint)indexMaxSlice[2], (uint)slices.size());
 616   //  DD(maxSlice);
 617
 618   // Get current label at the seedPoint of interest
 619   uint currentSlice=seedIndex[2];
 620   bool found = false;
 621   // LabelType previous_slice_label=recon->GetPixel(seedIndex);
 622
 623   // Currrent label map variable
 624   typedef itk::ShapeLabelObject< LabelType, 2> LabelObjectType;
 625   typedef itk::LabelMap< LabelObjectType > LabelMapType;
 626   typename LabelMapType::Pointer labelMap;
 627   std::vector<typename LabelObjectType::Pointer> shapeObjectsList;
 628   std::vector<MaskSlicePointer> shapeObjectsSliceList; // keep slice, useful for 'union'
 629   typename LabelObjectType::Pointer previousShapeObject;
 630
 631  if (GetVerboseTrackingFlag()) {
 632     std::cout << "TrackBifurcationFromPoint " << std::endl;
 633     std::cout << "\t seedPoint = " << seedPoint << std::endl;
 634     std::cout << "\t pointMaxSlice = " << pointMaxSlice << std::endl;
 635     std::cout << "\t newLabel = " << newLabel << std::endl;
 636   }
 637
 638   do {
 639     // Debug
 640     //std::cout << std::endl;
 641     //DD(currentSlice);
 642     ImagePointType c;
 643     clitk::PointsUtils<MaskImageType>::Convert2DTo3D(previousCentroid, m_Mask, currentSlice-1, c);
 644     //DD(c);
 645
 646     if (GetVerboseTrackingFlag()) {
 647       std::cout << "Loop slice = " << currentSlice << " c = " << c << std::endl;
 648     }
 649
 650     // Consider current reconstructed slice
 651     MaskSlicePointer s = slices[currentSlice];
 652     MaskSlicePointer previous;
 653     shapeObjectsList.clear();
 654     shapeObjectsSliceList.clear();
 655
 656     // Get shape of all labels in the current slice (it is already labelized)
 657
 658     // Normal -> same CCL with different label
 659     // PB -> sometimes same label in different CCL ! car growing
 660     //ADD les deux l+s ? mais avec max entre chaque ccl number  (bof)
 661     /*
 662       for each label in s -> map avec label in l; if different -> change
 663      */
 664     MaskSlicePointer ll = clitk::Labelize<MaskSliceType>(s, GetBackgroundValue(), true, 1);
 665     //    writeImage<MaskSliceType>(s, "slice-"+toString(currentSlice)+".mhd");
 666     //writeImage<MaskSliceType>(ll, "slice-label-"+toString(currentSlice)+".mhd");
 667     typedef itk::ImageRegionIteratorWithIndex<MaskSliceType> IteratorType;
 668     IteratorType its(s, s->GetLargestPossibleRegion());
 669     IteratorType itl(ll, ll->GetLargestPossibleRegion());
 670     std::map<LabelType, LabelType> labelInL;
 671     std::map<LabelType, std::map<LabelType, LabelType> > labelToChange;
 672     its.GoToBegin();
 673     itl.GoToBegin();
 674     int currentLabel = newLabel+10;
 675     while (!its.IsAtEnd()) {
 676       LabelType labs = its.Get();
 677       if (labs != GetBackgroundValue()) {
 678         LabelType labl = itl.Get();
 679         if (labelInL.find(labs) == labelInL.end()) { // Not found in map, first time
 680           //          DD("first");
 681           labelInL[labs] = labl;
 682           //DD(labs);
 683           //DD(labl);
 684         }
 685         else {
 686           if (labelInL[labs] != labl) { // I found a labs with a different labl. Store it.
 687             if (labelToChange[labs].find(labl) == labelToChange[labs].end()) { // if not already found
 688               //DD("found");
 689               //DD(labs);
 690               //DD(labl);
 691               //DD(labelInL[labs]);
 692               //DD(currentLabel);
 693               labelToChange[labs][labl] = currentLabel;
 694               ++currentLabel;
 695             }
 696           }
 697         }
 698       }
 699       ++its;
 700       ++itl;
 701     }
 702
 703     its.GoToBegin();
 704     itl.GoToBegin();
 705     while (!its.IsAtEnd()) {
 706       LabelType labs = its.Get();
 707       if (labs != GetBackgroundValue()) { // if not BG
 708         LabelType labl = itl.Get();
 709         if (labelToChange[labs].find(labl) != labelToChange[labs].end()) { // if some labs can change their label
 710           its.Set(labelToChange[labs][labl]); // change with the label for <labs-labl>
 711         }
 712       }
 713       ++its;
 714       ++itl;
 715     } // end while
 716
 717     //    writeImage<MaskSliceType>(s, "slice-final"+toString(currentSlice)+".mhd");
 718
 719
 720     labelMap = clitk::ComputeLabelMap<MaskSliceType, LabelType>(s, GetBackgroundValue(), true);
 721     //    DD(labelMap->GetNumberOfLabelObjects());
 722
 723     // If this is the first slice, get the object that contains the seed
 724     if (currentSlice == seedIndex[2]) {
 725       // DD("First slice");
 726       LabelType l = recon->GetPixel(seedIndex);
 727       // DD(l);
 728       shapeObjectsList.push_back(labelMap->GetLabelObject(l));
 729       shapeObjectsSliceList.push_back(s);
 730       previous = s;
 731       previousCentroid = shapeObjectsList[0]->GetCentroid();
 732       previousShapeObject = shapeObjectsList[0];
 733     }
 734     else {
 735       previous = slices[currentSlice-1];
 736       // Loop on labels to check if centroid is on the previous contour
 737       for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
 738         // Get the current label number
 739         int label = labelMap->GetLabels()[c];
 740         //DD(label);
 741         // Get the centroids
 742         SlicePointType centroid = labelMap->GetLabelObject(label)->GetCentroid();
 743         // Convert centroid into index in previous slice (same coordinate)
 744         typename MaskSliceType::IndexType centroidIndex;
 745         previous->TransformPhysicalPointToIndex(centroid, centroidIndex);
 746         LabelType labelInPreviousSlice = previous->GetPixel(centroidIndex);
 747         // if this current centroid was in the current label, we keep it
 748         //DD(labelInPreviousSlice);
 749         if (labelInPreviousSlice == newLabel) {
 750           shapeObjectsList.push_back(labelMap->GetLabelObject(label));
 751           shapeObjectsSliceList.push_back(s);
 752         }
 753       }
 754     }
 755
 756
 757     // Potentially the previous centroid could be inside another
 758     // labels, we consider i
 759     typename MaskSliceType::IndexType previousCentroidIndex;
 760     s->TransformPhysicalPointToIndex(previousCentroid, previousCentroidIndex);
 761     LabelType l = s->GetPixel(previousCentroidIndex);
 762     //DD(l);
 763     if (l != 0) { // if is not the background label
 764       int index = -1;
 765       for(uint c=0; c<shapeObjectsList.size(); c++) {
 766         if (shapeObjectsList[c]->GetLabel() == l) {
 767           index = c;
 768         }
 769       }
 770       if (index == -1) {
 771         //DD("not inside");
 772         shapeObjectsList.push_back(labelMap->GetLabelObject(l));
 773         shapeObjectsSliceList.push_back(s);
 774       }
 775       else {
 776         // DD("already inside");
 777       }
 778     }
 779
 780     // for(uint c=0; c<shapeObjectsList.size(); c++) {
 781     //   std::cout << c << " " << shapeObjectsList[c]->GetLabel() << " "
 782     //             << shapeObjectsList[c]->GetCentroid() << std::endl;
 783     // }
 784
 785
 786     // If several candidates, add one more with the union of all candidates
 787     MaskSlicePointer temp;
 788     if (shapeObjectsList.size() > 1) {
 789       //DD("add union");
 790       // Copy the slice
 791       temp = clitk::Clone<MaskSliceType>(s);
 792       // change label to a single label
 793       LabelType l = newLabel+1;
 794       for(uint c=0; c<shapeObjectsList.size(); c++) {
 795         LabelType ll = shapeObjectsList[c]->GetLabel();
 796         temp = clitk::SetBackground<MaskSliceType, MaskSliceType>(temp, temp, ll, l, true);
 797       }
 798       // Compute Label object properties
 799       labelMap = clitk::ComputeLabelMap<MaskSliceType, LabelType>(temp, GetBackgroundValue(), true);
 800       shapeObjectsList.push_back(labelMap->GetLabelObject(l));
 801       shapeObjectsSliceList.push_back(temp);
 802     }
 803
 804     /*
 805     for(uint c=0; c<shapeObjectsList.size(); c++) {
 806       std::cout << c << " " << shapeObjectsList[c]->GetLabel() << " "
 807                 << shapeObjectsList[c]->GetCentroid() << std::endl;
 808     }
 809     */
 810
 811
 812     for(uint c=0; c<shapeObjectsList.size(); c++) {
 813       ImagePointType cc;
 814       clitk::PointsUtils<MaskImageType>::Convert2DTo3D(shapeObjectsList[c]->GetCentroid(), m_Mask, currentSlice, cc);
 815       // std::cout << c << " " << shapeObjectsList[c]->GetLabel() << " "
 816       //   //                << shapeObjectsList[c]->GetCentroid() << " "
 817       //           << cc << " "
 818       //           << shapeObjectsList[c]->GetPhysicalSize() << " "
 819       //           << shapeObjectsList[c]->GetRoundness() << std::endl;
 820     }
 821
 822
 823     if (shapeObjectsList.size() == 0) {
 824       found = true;
 825     }
 826     else {
 827       // Heuristic to select the good one. For each candidate, we consider the size
 828       std::vector<double> sizes;
 829       std::vector<double> roundness;
 830       std::vector<size_t> index_sizes;
 831       std::vector<size_t> index_roundness;
 832       double previousSize = previousShapeObject->GetPhysicalSize();
 833       //DD(previousSize);
 834       for(uint c=0; c<shapeObjectsList.size(); c++) {
 835         double s = shapeObjectsList[c]->GetPhysicalSize();
 836         sizes.push_back(fabs(previousSize-s)/previousSize);
 837         roundness.push_back(fabs(1.0-shapeObjectsList[c]->GetRoundness()));
 838         index_sizes.push_back(c);
 839         index_roundness.push_back(c);
 840       }
 841       //DDV(sizes, sizes.size());
 842       //DDV(roundness, roundness.size());
 843       // DDV(index_sizes, index_sizes.size());
 844       // DDV(index_roundness, index_roundness.size());
 845       sort(index_sizes.begin(), index_sizes.end(), index_cmp<std::vector<double>&>(sizes));
 846       sort(index_roundness.begin(), index_roundness.end(), index_cmp<std::vector<double>&>(roundness));
 847       //DDV(index_sizes, index_sizes.size());
 848       // DDV(index_roundness, index_roundness.size());
 849
 850       // TEMPORARY GET THE FIRST
 851       int best = index_sizes[0];
 852       // if (currentSlice == seedIndex[2]) { // first contour => idiot, first = single contour
 853       //   best = index_roundness[0]; // best is more round
 854       // }
 855       LabelType label = shapeObjectsList[best]->GetLabel();
 856       //      DD(label);
 857       s = shapeObjectsSliceList[best];
 858       s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, label, newLabel, true);
 859
 860       // HERE
 861
 862       // It can happend that several CCL share this same label. To
 863       // prevent this case, we only consider the one that contains
 864       // the centroid.
 865
 866       // TODO -> PREVIOUS CENTROID ???
 867
 868       MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
 869       temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
 870       typename MaskSliceType::IndexType centroids_index;
 871       temp->TransformPhysicalPointToIndex(shapeObjectsList[best]->GetCentroid(), centroids_index);
 872       typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
 873         if (v == GetBackgroundValue()) {
 874           // DD(currentSlice);
 875           // DD("inside BG");
 876           //DD(centroids[0]);
 877           v = 1; // largest one
 878         }
 879
 880         //DD(v);
 881         temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
 882         //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
 883         s = temp;
 884
 885
 886       // end
 887       slices[currentSlice] = s;
 888       previousCentroid = shapeObjectsList[best]->GetCentroid();
 889       previousShapeObject = shapeObjectsList[best];
 890     }
 891
 892     ++currentSlice;
 893
 894     if (currentSlice == maxSlice) {
 895       // DD("end max slice");
 896       found = true;
 897     }
 898
 899   } while (!found);
 900
 901   /*
 902   // -------------------------
 903   // If no centroid were found
 904   if (shapeObjectsList.size() == 0) {
 905   // Last attempt to find -> check if previous centroid is inside a CCL
 906   //      if in s -> get value, getcentroid add.
 907   DD(currentSlice);
 908   DD("Last change to find");
 909   DD(previous_slice_label);
 910   DD(newLabel);
 911   typename MaskSliceType::IndexType previousCentroidIndex;
 912   s->TransformPhysicalPointToIndex(previousCentroid, previousCentroidIndex);
 913   DD(previousCentroid);
 914   LabelType labelInSlice = s->GetPixel(previousCentroidIndex);
 915   DD(labelInSlice);
 916   if (labelInSlice != GetBackgroundValue()) {
 917   centroids.push_back(labelMap->GetLabelObject(labelInSlice)->GetCentroid());
 918   centroids_label.push_back(labelInSlice);
 919   labels_size.push_back(labelMap->GetLabelObject(labelInSlice)->GetPhysicalSize());
 920   }
 921   }
 922
 923   // Some centroid were found
 924   // If several centroids, we found a bifurcation
 925   if (centroids.size() > 1) {
 926   //      int n = centroids.size();
 927   // Debug point
 928   std::vector<ImagePointType> d;
 929   clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(centroids, currentSlice, m_Mask, d);
 930   DDV(d, d.size());
 931
 932   // new potential bifurcation found
 933   numberOfBifurcation++;
 934   // If the number of bifurcation is greater than the required one, we stop
 935   if (numberOfBifurcation > GetMaxNumberOfFoundBifurcation()) {
 936   found = true;
 937   DD("max bif reach");
 938   for(uint c=0; c<centroids.size(); c++) {
 939   ImagePointType bif;
 940   clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, currentSlice, bif);
 941   bifurcations.push_back(bif);
 942   }
 943   }
 944   // Else we continue along the main (largest) connected component
 945   else {
 946   int indexOfLargest = 0;
 947   for(uint b=0; b<centroids.size(); b++) {
 948   if (labels_size[b] > labels_size[indexOfLargest]) {
 949   indexOfLargest = b;
 950   }
 951   }
 952   DD(indexOfLargest);
 953   DD(labels_size[indexOfLargest]);
 954   SlicePointType c = centroids[indexOfLargest];
 955   LabelType l = centroids_label[indexOfLargest];
 956   DD(l);
 957   DD(c);
 958   centroids.clear();
 959   centroids.push_back(c);
 960   centroids_label.clear();
 961   centroids_label.push_back(l);
 962   }
 963   }
 964   */
 965
 966   /*    ==> here all centroid are considered as ok.*/
 967
 968   /*
 969   // REMOVE IF CENTROID=1, REPLACE BY >0
 970
 971   // if only one centroids, we change the current image with the current label
 972   if (centroids.size() == 1) {
 973   DD(centroids_label[0]);
 974   s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
 975   slices[currentSlice] = s;
 976   previous_slice_label = newLabel;
 977   // It can happend that several CCL share this same label. To
 978   // prevent this case, we only consider the one that contains
 979   // the centroid.
 980   MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
 981   //      writeImage<MaskSliceType>(temp, "bin-"+toString(currentSlice)+".mhd");
 982   temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
 983   //writeImage<MaskSliceType>(temp, "label-"+toString(currentSlice)+".mhd");
 984   typename MaskSliceType::IndexType centroids_index;
 985   temp->TransformPhysicalPointToIndex(centroids[0], centroids_index);
 986   typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
 987
 988   // It can happend that the centroid is inside the BG, so we keep
 989   // the largest CCL (the first);
 990   if (v == GetBackgroundValue()) {
 991   DD(currentSlice);
 992   DD("inside BG");
 993   DD(centroids[0]);
 994   v = 1; // largest one
 995   }
 996
 997   //DD(v);
 998   temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
 999   //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
1000   s = temp;
1001   slices[currentSlice] = s;
1002
1003   // I need to recompute the centroid if we have removed some
1004   // connected component.
1005   clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), centroids);
1006   previousCentroid = centroids[1];
1007   }
1008
1009   if (centroids.size() == 0) {
1010   DD("ZERO");
1011   found = true;
1012   }
1013
1014   if (currentSlice == slices.size()-1) {
1015   DD("end of slices");
1016   found = true;
1017   }
1018
1019   if (currentSlice == maxSlice) {
1020   DD("end max slice");
1021   found = true;
1022   }
1023
1024   // iterate
1025   ++currentSlice;
1026   } while (!found);
1027   */
1028
1029   // End
1030   StopCurrentStep();
1031 }
1032 //--------------------------------------------------------------------
1033
1034
1035 #endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX