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 //--------------------------------------------------------------------
  33 template <class TImageType>
  34 clitk::ExtractMediastinalVesselsFilter<TImageType>::
  35 ExtractMediastinalVesselsFilter():
  36   clitk::FilterBase(),
  37   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  38   itk::ImageToImageFilter<TImageType, MaskImageType>()
  39 {
  40   this->SetNumberOfRequiredInputs(1);
  41   SetBackgroundValue(0);
  42   SetForegroundValue(1);
  43   SetThresholdHigh(140);
  44   SetThresholdLow(55);
  45   SetErosionRadius(2);
  46   SetDilatationRadius(9);
  47   SetMaxDistancePostToCarina(10);
  48   SetMaxDistanceAntToCarina(40);
  49   SetMaxDistanceLeftToCarina(35);
  50   SetMaxDistanceRightToCarina(35);
  51   SetSoughtVesselSeedName("NoSeedNameGiven");
  52 }
  53 //--------------------------------------------------------------------
  54
  55
  56 //--------------------------------------------------------------------
  57 template <class TImageType>
  58 void
  59 clitk::ExtractMediastinalVesselsFilter<TImageType>::
  60 GenerateOutputInformation() {
  61   // Get inputs
  62   LoadAFDB();
  63   m_Input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
  64
  65   // ------------------------------------------------------------------
  66   // Crop the initial image superiorly and inferiorly.
  67   // TODO : add options for x cm above/below
  68   CropInputImage();
  69
  70   // ------------------------------------------------------------------
  71   // Binarize the image. Need two thresholds, one high to select
  72   // structures (CCL) that are almost not connected (after erosion),
  73   // and one low thresholds to select the real contours. Will be
  74   // reconstructed later.
  75   StartNewStep("Binarize with high threshold = "+toString(GetThresholdHigh()));
  76   typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> BinarizeFilterType;
  77   typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
  78   binarizeFilter->SetInput(m_Input);
  79   binarizeFilter->SetLowerThreshold(GetThresholdHigh());
  80   binarizeFilter->SetInsideValue(GetForegroundValue());
  81   binarizeFilter->SetOutsideValue(GetBackgroundValue());
  82   binarizeFilter->Update();
  83   m_Mask = binarizeFilter->GetOutput();
  84   StopCurrentStep<MaskImageType>(m_Mask);
  85
  86   // ------------------------------------------------------------------
  87   StartNewStep("Binarize with low threshold = "+toString(GetThresholdLow()));
  88   binarizeFilter = BinarizeFilterType::New();
  89   binarizeFilter->SetInput(m_Input);
  90   binarizeFilter->SetLowerThreshold(GetThresholdLow());
  91   binarizeFilter->SetInsideValue(GetForegroundValue());
  92   binarizeFilter->SetOutsideValue(GetBackgroundValue());
  93   binarizeFilter->Update();
  94   MaskImagePointer m_Mask2 = binarizeFilter->GetOutput();
  95   StopCurrentStep<MaskImageType>(m_Mask2);
  96
  97   // ------------------------------------------------------------------
  98   // Extract slices
  99   StartNewStep("Detect objects : erosion, then slice by slice reconstruction");
 100   std::vector<MaskSlicePointer> slices_mask;
 101   clitk::ExtractSlices<MaskImageType>(m_Mask, 2, slices_mask);
 102   std::vector<MaskSlicePointer> slices_mask2;
 103   clitk::ExtractSlices<MaskImageType>(m_Mask2, 2, slices_mask2);
 104   int radius = GetErosionRadius();
 105
 106   // List of working slices (debug only)
 107   std::vector<MaskSlicePointer> debug_eroded;
 108   // std::vector<MaskSlicePointer> debug_labeled;
 109
 110   // ------------------------------------------------------------------
 111   // Loop Slice by Slice in order to break connectivity between
 112   // CCL. Erode and reconstruct all labels at the same time without
 113   // merging them.
 114   for(uint i=0; i<slices_mask.size(); i++) {
 115     // Erosion kernel
 116     typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,2> KernelType;
 117     KernelType structuringElement;
 118     structuringElement.SetRadius(radius);
 119     structuringElement.CreateStructuringElement();
 120
 121     // Erosion -> we break the connectivity between structure
 122     typedef itk::BinaryErodeImageFilter<MaskSliceType, MaskSliceType, KernelType> ErodeFilterType;
 123     typename ErodeFilterType::Pointer eroder = ErodeFilterType::New();
 124     eroder->SetInput(slices_mask[i]);
 125     eroder->SetBackgroundValue(GetBackgroundValue());
 126     eroder->SetForegroundValue(GetForegroundValue());
 127     eroder->SetBoundaryToForeground(true); // ??
 128     eroder->SetKernel(structuringElement);
 129     eroder->Update();
 130     MaskSlicePointer eroded = eroder->GetOutput();
 131     debug_eroded.push_back(eroded);
 132
 133     // Labelize (CCL)
 134     MaskSlicePointer labeled =
 135       clitk::Labelize<MaskSliceType>(eroded, GetBackgroundValue(), true, 1); // Fully connected !
 136     // debug_labeled.push_back(labeled);
 137
 138     // Make Reconstruction filter : dilation all labels at the same
 139     // time, prevent to merge them.
 140     typedef clitk::ReconstructWithConditionalGrayscaleDilateImageFilter<MaskSliceType> ReconstructFilterType;
 141     typename ReconstructFilterType::Pointer reconstructor = ReconstructFilterType::New();
 142     reconstructor->SetInput(labeled);
 143     reconstructor->SetIterationNumber(radius+GetDilatationRadius());
 144     reconstructor->Update();
 145     MaskSlicePointer s = reconstructor->GetOutput();
 146
 147     // Remove Initial BG of the second tresholded image
 148     s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, slices_mask2[i],
 149                                                            GetBackgroundValue(), GetBackgroundValue(), true);
 150     m_slice_recon.push_back(s);
 151
 152   } // end loop
 153
 154   // Build 3D images from the slice by slice processing
 155   MaskImageType::Pointer eroded = clitk::JoinSlices<MaskImageType>(debug_eroded, m_Mask, 2);
 156   writeImage<MaskImageType>(eroded, "erode.mhd");
 157   //MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);
 158   MaskImageType::Pointer r = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
 159   writeImage<MaskImageType>(r, "recon1.mhd");
 160
 161   // ------------------------------------------------------------------
 162   // Track the SoughtVessel from the given first point
 163   // superiorly. This is done by TrackBifurcationFromPoint
 164   MaskImagePointType SoughtVesselSeedPoint;
 165   GetAFDB()->GetPoint3D(m_SoughtVesselSeedName, SoughtVesselSeedPoint);
 166
 167   // Find the label with the maximum value to set the result
 168   typedef itk::MinimumMaximumImageCalculator<MaskImageType> MinMaxFilterType;
 169   MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
 170   ff->SetImage(r);
 171   ff->ComputeMaximum();
 172   LabelType newLabel = ff->GetMaximum()+1;
 173
 174   // the following bifurcations point will the centroids of the
 175   // components obtain when (hopefully!) the SoughtVessel
 176   // split into CommonArtery and SubclavianArtery.
 177   std::vector<MaskImagePointType> bifurcations;
 178   TrackBifurcationFromPoint(r, m_slice_recon, SoughtVesselSeedPoint, newLabel, bifurcations);
 179
 180   // Build the final 3D image from the previous slice by slice processing
 181   m_SoughtVessel = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
 182   writeImage<MaskImageType>(m_SoughtVessel, "recon2.mhd");
 183
 184   // Set binary image, (remove other labels).
 185   // TODO: keep labeled image to track SubclavianArtery and CommonArtery
 186   m_SoughtVessel =
 187     clitk::Binarize<MaskImageType>(m_SoughtVessel, newLabel, newLabel,
 188                                    GetBackgroundValue(), GetForegroundValue());
 189
 190   writeImage<MaskImageType>(m_SoughtVessel, "afterbinarize.mhd");
 191
 192   m_SoughtVessel = clitk::AutoCrop<MaskImageType>(m_SoughtVessel, GetBackgroundValue());
 193
 194   writeImage<MaskImageType>(m_SoughtVessel, "afterautocrop.mhd");
 195
 196   // Clean the image : Opening (not in Z direction)
 197   typename MaskImageType::SizeType rad;
 198   rad[0] = rad[1] = 2;
 199   rad[2] = 0;
 200   m_SoughtVessel = clitk::Opening<MaskImageType>(m_SoughtVessel, rad,
 201                                                           GetBackgroundValue(), GetForegroundValue());
 202
 203   writeImage<MaskImageType>(m_SoughtVessel, "afteropen.mhd");
 204
 205   // Clean the image : keep main CCL slice by slice
 206   m_SoughtVessel = clitk::SliceBySliceKeepMainCCL<MaskImageType>(m_SoughtVessel,
 207                                                                           GetBackgroundValue(),
 208                                                                           GetForegroundValue());
 209 }
 210 //--------------------------------------------------------------------
 211
 212
 213 //--------------------------------------------------------------------
 214 template <class TImageType>
 215 void
 216 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 217 GenerateInputRequestedRegion() {
 218   //DD("GenerateInputRequestedRegion (nothing?)");
 219 }
 220 //--------------------------------------------------------------------
 221
 222
 223 //--------------------------------------------------------------------
 224 template <class TImageType>
 225 void
 226 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 227 GenerateData() {
 228   // Save in the AFDB (not write on the disk here)
 229   GetAFDB()->SetImageFilename(GetSoughtVesselName(), GetOutputFilename());
 230   WriteAFDB();
 231   // Final Step -> graft output
 232   this->GraftNthOutput(0, m_SoughtVessel);
 233 }
 234 //--------------------------------------------------------------------
 235
 236
 237 //--------------------------------------------------------------------
 238 template <class TImageType>
 239 void
 240 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 241 CropInputImage() {
 242   StartNewStep("Crop the input image: SI,AP limits with carina and crop with mediastinum");
 243   /*
 244     Need : Trachea, Carina (roi not point),
 245    */
 246   // Get Trachea and Carina
 247   MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
 248
 249   // Compute Carina position
 250   double m_CarinaZ;
 251   MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
 252   std::vector<MaskImagePointType> centroids;
 253   clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
 254   m_CarinaZ = centroids[1][2];
 255   // add one slice to include carina
 256   m_CarinaZ += Carina->GetSpacing()[2];
 257   // We dont need Carina structure from now
 258   Carina->Delete();
 259   GetAFDB()->SetPoint3D("CarinaPoint", centroids[1]);
 260
 261   // Crop Inf, remove below Carina
 262   m_Input =
 263     clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 2, m_CarinaZ, false, GetBackgroundValue());
 264
 265   // Crop post
 266   double m_CarinaY = centroids[1][1];
 267   m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 1,
 268                                                          m_CarinaY+GetMaxDistancePostToCarina(),
 269                                                          false, GetBackgroundValue());
 270   // Crop ant
 271   m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 1,
 272                                                        m_CarinaY-GetMaxDistanceAntToCarina(),
 273                                                        false, GetBackgroundValue());
 274
 275   // Crop Right
 276   double m_CarinaX = centroids[1][0];
 277   m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 0,
 278                                                        m_CarinaX-GetMaxDistanceRightToCarina(),
 279                                                        false, GetBackgroundValue());
 280   // Crop Left
 281   m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 0,
 282                                                          m_CarinaX+GetMaxDistanceLeftToCarina(),
 283                                                          false, GetBackgroundValue());
 284
 285   /*
 286   // AutoCrop with Mediastinum, generaly only allow to remove few slices (superiorly)
 287   m_Mediastinum  = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
 288   // Resize like input (sup to carina)
 289   m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Input, GetBackgroundValue());
 290   // Auto crop
 291   m_Mediastinum = clitk::AutoCrop<MaskImageType>(m_Mediastinum, GetBackgroundValue());
 292   // Resize input
 293   m_Input = clitk::ResizeImageLike<ImageType>(m_Input, m_Mediastinum, GetBackgroundValue());
 294   */
 295
 296   // End
 297   StopCurrentStep<ImageType>(m_Input);
 298 }
 299 //--------------------------------------------------------------------
 300
 301
 302 //--------------------------------------------------------------------
 303 template <class TImageType>
 304 void
 305 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 306 TrackBifurcationFromPoint(MaskImagePointer & recon,
 307                           std::vector<MaskSlicePointer> & slices_recon,
 308                           MaskImagePointType point3D,
 309                           LabelType newLabel,
 310                           std::vector<MaskImagePointType> & bifurcations) {
 311   StartNewStep("Track the SoughtVessel from the seed point");
 312
 313   // Find first slice index
 314   MaskImageIndexType index;
 315   recon->TransformPhysicalPointToIndex(point3D, index);
 316   int numberOfBifurcation = 0;
 317   typedef typename MaskSliceType::PointType SlicePointType;
 318   SlicePointType previousCenter;
 319
 320   // Get current label at the point3D of interest
 321   uint currentSlice=index[2];
 322   bool found = false;
 323   LabelType previous_slice_label=recon->GetPixel(index);
 324   //  DD(slices_recon.size());
 325   do {
 326     //    DD(currentSlice);
 327     // Consider current reconstructed slice
 328     MaskSlicePointer s = slices_recon[currentSlice];
 329     MaskSlicePointer previous;
 330     if (currentSlice == index[2]) previous = s;
 331     else {
 332       previous = slices_recon[currentSlice-1];
 333     }
 334
 335     // Get centroids of the labels in the current slice
 336     static const unsigned int Dim = MaskSliceType::ImageDimension;
 337     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 338     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 339     typedef itk::LabelImageToLabelMapFilter<MaskSliceType, LabelMapType> ImageToMapFilterType;
 340     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 341     typedef itk::ShapeLabelMapFilter<LabelMapType, MaskSliceType> ShapeFilterType;
 342     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 343     imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
 344     imageToLabelFilter->SetInput(s);
 345     statFilter->SetInput(imageToLabelFilter->GetOutput());
 346     statFilter->Update();
 347     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 348
 349     // Look what centroid inside the previous largest one
 350     std::vector<SlicePointType> centroids;
 351     std::vector<LabelType> centroids_label;
 352     std::vector<double> labels_size;
 353     for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
 354       int label = labelMap->GetLabels()[c];
 355       //      DD(label);
 356       SlicePointType center = labelMap->GetLabelObject(label)->GetCentroid();
 357       //      DD(center);
 358       // Get label into previous slice
 359       typename MaskSliceType::IndexType centerIndex;
 360       previous->TransformPhysicalPointToIndex(center, centerIndex);
 361       LabelType labelInPreviousSlice = previous->GetPixel(centerIndex);
 362       // if this current centroid was in the current label, add it to bifurcations
 363       if (labelInPreviousSlice == previous_slice_label) {
 364         centroids.push_back(center);
 365         centroids_label.push_back(label);
 366         labels_size.push_back(labelMap->GetLabelObject(label)->GetPhysicalSize());
 367       }
 368     }
 369
 370     if (centroids.size() == 0) {
 371       // Last attempt to find -> check if previous centroid is inside a CCL
 372       //      if in s -> get value, getcentroid add.
 373       DD(currentSlice);
 374       DD("Last change to find");
 375       typename MaskSliceType::IndexType previousCenterIndex;
 376       s->TransformPhysicalPointToIndex(previousCenter, previousCenterIndex);
 377       DD(previousCenter);
 378       LabelType labelInSlice = s->GetPixel(previousCenterIndex);
 379       DD(labelInSlice);
 380       if (labelInSlice != GetBackgroundValue()) {
 381         centroids.push_back(labelMap->GetLabelObject(labelInSlice)->GetCentroid());
 382         centroids_label.push_back(labelInSlice);
 383         labels_size.push_back(labelMap->GetLabelObject(labelInSlice)->GetPhysicalSize());
 384       }
 385     }
 386
 387
 388     //    DD(centroids.size());
 389
 390     // If several centroids, we found a bifurcation
 391     if (centroids.size() > 1) {
 392       numberOfBifurcation++;
 393       // If the number of bifurcation is greater than the required one, we stop
 394       if (numberOfBifurcation > GetMaxNumberOfFoundBifurcation()) {
 395         found = true;
 396         DD("max bif reach");
 397         for(uint c=0; c<centroids.size(); c++) {
 398           ImagePointType bif;
 399           clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, currentSlice, bif);
 400           bifurcations.push_back(bif);
 401         }
 402       }
 403       // Else we continue along the main (largest) connected component
 404       else {
 405         int indexOfLargest = 0;
 406         for(uint b=0; b<centroids.size(); b++) {
 407           if (labels_size[b] > labels_size[indexOfLargest]) {
 408             indexOfLargest = b;
 409           }
 410         }
 411         SlicePointType c = centroids[indexOfLargest];
 412         LabelType l = centroids_label[indexOfLargest];
 413         centroids.clear();
 414         centroids.push_back(c);
 415         centroids_label.push_back(l);
 416       }
 417     }
 418
 419     // if only one centroids, we change the current image with the current label
 420     if (centroids.size() == 1) {
 421       s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
 422       slices_recon[currentSlice] = s;
 423       previous_slice_label = newLabel;
 424       // It can happend that several CCL share this same label. To
 425       // prevent this case, we only consider the one that contains
 426       // the centroid.
 427       MaskSlicePointer temp = clitk::Binarize<MaskSliceType>(s, newLabel, newLabel, GetBackgroundValue(), GetForegroundValue());
 428       //      writeImage<MaskSliceType>(temp, "bin-"+toString(currentSlice)+".mhd");
 429       temp = clitk::Labelize<MaskSliceType>(temp, GetBackgroundValue(), true, 1);
 430       //writeImage<MaskSliceType>(temp, "label-"+toString(currentSlice)+".mhd");
 431       typename MaskSliceType::IndexType centroids_index;
 432       temp->TransformPhysicalPointToIndex(centroids[0], centroids_index);
 433       typename MaskSliceType::PixelType v = temp->GetPixel(centroids_index);
 434
 435       // It can happend that the centroid is inside the BG, so we keep
 436       // the largest CCL (the first);
 437       if (v == GetBackgroundValue()) {
 438         DD(currentSlice);
 439         DD("inside BG");
 440         DD(centroids[0]);
 441         v = 1; // largest one
 442       }
 443
 444       //DD(v);
 445       temp = clitk::Binarize<MaskSliceType>(temp, v, v, GetBackgroundValue(), newLabel);
 446       //writeImage<MaskSliceType>(temp, "relabel-"+toString(currentSlice)+".mhd");
 447       s = temp;
 448       slices_recon[currentSlice] = s;
 449
 450       // I need to recompute the centroid if we have removed some
 451       // connected component.
 452       clitk::ComputeCentroids<MaskSliceType>(s, GetBackgroundValue(), centroids);
 453       previousCenter = centroids[1];
 454     }
 455
 456     if (centroids.size() == 0) {
 457       DD("ZERO");
 458       found = true;
 459     }
 460
 461     if (currentSlice == slices_recon.size()-1) {
 462       DD("end of slices");
 463       found = true;
 464     }
 465
 466     // iterate
 467     ++currentSlice;
 468   } while (!found);
 469
 470   // End
 471   StopCurrentStep();
 472 }
 473 //--------------------------------------------------------------------
 474
 475
 476 #endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX