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