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 "clitkAutoCropFilter.h"
  26 #include "clitkSegmentationUtils.h"
  27 #include "clitkMorphoMathFilter.h"
  28
  29 // itk
  30 #include <itkBinaryThresholdImageFilter.h>
  31 #include <itkGrayscaleDilateImageFilter.h>
  32 #include <itkMinimumMaximumImageCalculator.h>
  33
  34 //--------------------------------------------------------------------
  35 template <class TImageType>
  36 clitk::ExtractMediastinalVesselsFilter<TImageType>::
  37 ExtractMediastinalVesselsFilter():
  38   clitk::FilterBase(),
  39   clitk::FilterWithAnatomicalFeatureDatabaseManagement(),
  40   itk::ImageToImageFilter<TImageType, MaskImageType>()
  41 {
  42   this->SetNumberOfRequiredInputs(1);
  43   SetBackgroundValue(0);
  44   SetForegroundValue(1);
  45   SetThreshold(140);
  46   SetTemporaryForegroundValue(1);
  47 }
  48 //--------------------------------------------------------------------
  49
  50
  51 //--------------------------------------------------------------------
  52 template <class TImageType>
  53 void
  54 clitk::ExtractMediastinalVesselsFilter<TImageType>::
  55 GenerateOutputInformation() {
  56   // Get inputs
  57   LoadAFDB();
  58   m_Input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
  59
  60   // ------------------------------------------------------------------
  61   // Sup-Inf crop -> Carina
  62   CropSupInf();
  63
  64   // ------------------------------------------------------------------
  65   // Binarize
  66   StartNewStep("Binarize with treshold = "+toString(GetThreshold()));
  67   typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> BinarizeFilterType;
  68   typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
  69   binarizeFilter->SetInput(m_Input);
  70   binarizeFilter->SetLowerThreshold(GetThreshold());
  71   binarizeFilter->SetInsideValue(GetTemporaryForegroundValue());
  72   binarizeFilter->SetOutsideValue(GetBackgroundValue());
  73   binarizeFilter->Update();
  74   m_Mask = binarizeFilter->GetOutput();
  75   clitk::writeImage<MaskImageType>(m_Mask, "m.mhd");
  76   StopCurrentStep<MaskImageType>(m_Mask);
  77
  78   // Keep main CCL ?
  79   m_Mask = clitk::Labelize<MaskImageType>(m_Mask, GetBackgroundValue(), false, 10);
  80   m_Mask = KeepLabels<MaskImageType>(m_Mask, GetBackgroundValue(), GetTemporaryForegroundValue(), 1, 1, true);
  81   clitk::writeImage<MaskImageType>(m_Mask, "m2.mhd");
  82
  83   // ------------------------------------------------------------------
  84   // Extract slices
  85   StartNewStep("Detect vessels (slice by slice reconstruction)");
  86   std::vector<MaskSlicePointer> slices_mask;
  87   clitk::ExtractSlices<MaskImageType>(m_Mask, 2, slices_mask);
  88
  89   DD(slices_mask.size());
  90
  91   std::vector<MaskSlicePointer> debug_eroded;
  92   std::vector<MaskSlicePointer> debug_labeled;
  93   std::vector<MaskSlicePointer> debug_slabeled;
  94
  95   int radius = 3;
  96   DD(radius); // TO PUT IN OPTION
  97
  98   // ------------------------------------------------------------------
  99   // Loop Slice by Slice -> erode find CCL and reconstruct
 100   clitk::MorphoMathFilter<MaskSliceType>::Pointer f= clitk::MorphoMathFilter<MaskSliceType>::New();
 101   for(uint i=0; i<slices_mask.size(); i++) {
 102     // Erosion
 103     f->SetInput(slices_mask[i]);
 104     f->SetBackgroundValue(GetBackgroundValue());
 105     f->SetForegroundValue(GetTemporaryForegroundValue());
 106     f->SetRadius(radius);
 107     f->SetOperationType(0); // Erode
 108     f->VerboseFlagOff();
 109     f->Update();
 110     MaskSlicePointer eroded = f->GetOutput();
 111     //    writeImage<MaskSliceType>(eroded, "er-"+toString(i)+".mhd");
 112     debug_eroded.push_back(eroded);
 113
 114     // CCL
 115     int nb;
 116     MaskSlicePointer labeled =
 117       clitk::LabelizeAndCountNumberOfObjects<MaskSliceType>(eroded, GetBackgroundValue(), false, 1, nb);
 118
 119     // Relabel, large CCL with large label number
 120     for(int n=nb; n>0; n--) {
 121       //        DD(n);
 122       int li = n;
 123       int lo = 2*(nb+1)-li;
 124       labeled = clitk::SetBackground<MaskSliceType, MaskSliceType>(labeled, labeled, li, lo, true);
 125     }
 126     debug_labeled.push_back(labeled);
 127
 128     // Create kernel for GrayscaleDilateImageFilter
 129     typedef itk::BinaryBallStructuringElement<MaskSliceType::PixelType,MaskSliceType::ImageDimension > KernelType;
 130     KernelType k;
 131     k.SetRadius(radius+1);
 132     k.CreateStructuringElement();
 133
 134     // Keep the MAX -> we prefer the opposite su change the label
 135     typedef itk::GrayscaleDilateImageFilter<MaskSliceType, MaskSliceType, KernelType> FilterType;
 136     FilterType::Pointer m = FilterType::New();
 137     m->SetKernel(k);
 138     m->SetInput(labeled);
 139     // DD(m->GetAlgorithm());
 140     // m->SetAlgorithm(3);
 141     m->Update();
 142     MaskSlicePointer s = m->GetOutput();
 143
 144
 145     // Remove Initial BG
 146     s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, slices_mask[i],
 147                                                    GetBackgroundValue(), GetBackgroundValue(), true);
 148
 149     m_slice_recon.push_back(s);
 150   } // end loop
 151   DD("end loop");
 152
 153   MaskImageType::Pointer eroded = clitk::JoinSlices<MaskImageType>(debug_eroded, m_Mask, 2);
 154   clitk::writeImage<MaskImageType>(eroded, "eroded.mhd");
 155
 156   DD("l");
 157   MaskImageType::Pointer l = clitk::JoinSlices<MaskImageType>(debug_labeled, m_Mask, 2);
 158   clitk::writeImage<MaskImageType>(l, "labeled.mhd");
 159
 160   DD("r");
 161   MaskImageType::Pointer r = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
 162   clitk::writeImage<MaskImageType>(r, "recon.mhd");
 163
 164   // ------------------------------------------------------------------
 165   // Loop Slice by Slice -> BCA not found yet
 166   /*  MaskImagePointType BCA_p;
 167   GetAFDB()->GetPoint3D("BrachioCephalicArteryFirstInferiorPoint", BCA_p);
 168   DD(BCA_p);
 169   MaskImagePointType bif1;
 170   MaskImagePointType bif2;
 171   TrackBifurcationFromPoint(r, BCA_p, bif1, bif2);
 172   DD(bif1);
 173   DD(bif2);
 174   */
 175   // Find max label
 176   typedef itk::MinimumMaximumImageCalculator<MaskImageType> MinMaxFilterType;
 177   MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
 178   ff->SetImage(r);
 179   ff->ComputeMaximum();
 180   LabelType newLabel = ff->GetMaximum()+1;
 181   DD(newLabel);
 182
 183   // Get all centroids of the first slice
 184   std::vector<MaskSlicePointType> centroids2D;
 185   clitk::ComputeCentroids<MaskSliceType>(m_slice_recon[0], GetBackgroundValue(), centroids2D);
 186   DD(centroids2D.size());
 187   std::vector<MaskImagePointType> bifurcations;
 188   clitk::PointsUtils<MaskImageType>::Convert2DListTo3DList(centroids2D, 0, r, bifurcations);
 189   DD(bifurcations.size());
 190   for(uint i=1; i<bifurcations.size()+1; i++) {
 191     DD(i);
 192     DD(bifurcations.size());
 193     TrackBifurcationFromPoint(r, m_slice_recon, bifurcations[i], newLabel+i, bifurcations);
 194     DD("end track");
 195     DD(bifurcations.size());
 196     MaskImageType::Pointer rr = clitk::JoinSlices<MaskImageType>(m_slice_recon, m_Mask, 2);
 197     clitk::writeImage<MaskImageType>(rr, "recon"+toString(i)+".mhd");
 198   }
 199
 200 }
 201 //--------------------------------------------------------------------
 202
 203
 204 //--------------------------------------------------------------------
 205 template <class TImageType>
 206 void
 207 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 208 GenerateInputRequestedRegion() {
 209   //DD("GenerateInputRequestedRegion (nothing?)");
 210 }
 211 //--------------------------------------------------------------------
 212
 213
 214 //--------------------------------------------------------------------
 215 template <class TImageType>
 216 void
 217 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 218 GenerateData() {
 219   DD("GenerateData");
 220   // Final Step -> graft output (if SetNthOutput => redo)
 221   MaskImagePointer BrachioCephalicArtery =
 222     GetAFDB()->template GetImage<MaskImageType>("BrachioCephalicArtery");
 223   this->GraftNthOutput(0, BrachioCephalicArtery);
 224 }
 225 //--------------------------------------------------------------------
 226
 227
 228 //--------------------------------------------------------------------
 229 template <class TImageType>
 230 void
 231 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 232 CropSupInf() {
 233   StartNewStep("Inf/Sup limits (carina) and crop with mediastinum");
 234   // Get Trachea and Carina
 235   MaskImagePointer Trachea = GetAFDB()->template GetImage <MaskImageType>("Trachea");
 236
 237   // Get or compute Carina
 238   double m_CarinaZ;
 239   // Get Carina Z position
 240   MaskImagePointer Carina = GetAFDB()->template GetImage<MaskImageType>("Carina");
 241
 242   std::vector<MaskImagePointType> centroids;
 243   clitk::ComputeCentroids<MaskImageType>(Carina, GetBackgroundValue(), centroids);
 244   m_CarinaZ = centroids[1][2];
 245   // DD(m_CarinaZ);
 246   // add one slice to include carina ?
 247   m_CarinaZ += Carina->GetSpacing()[2];
 248   // We dont need Carina structure from now
 249   Carina->Delete();
 250   GetAFDB()->SetDouble("CarinaZ", m_CarinaZ);
 251
 252   // Crop
 253   m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 2,
 254                                                        m_CarinaZ, false, GetBackgroundValue());
 255
 256   // Crop not post to centroid
 257   double m_CarinaY = centroids[1][1];
 258   DD(m_CarinaY);
 259   m_Input = clitk::CropImageRemoveGreaterThan<ImageType>(m_Input, 1, // OLD ABOVE
 260                                                          m_CarinaY, false, GetBackgroundValue());
 261   // Crop not ant to centroid
 262   m_Input = clitk::CropImageRemoveLowerThan<ImageType>(m_Input, 1,
 263                                                          m_CarinaY-80, false, GetBackgroundValue());
 264
 265   // AutoCrop with Mediastinum
 266   m_Mediastinum  = GetAFDB()->template GetImage<MaskImageType>("Mediastinum");
 267   // Resize like input (sup to carina)
 268   m_Mediastinum = clitk::ResizeImageLike<MaskImageType>(m_Mediastinum, m_Input, GetBackgroundValue());
 269   // Auto crop
 270   m_Mediastinum = clitk::AutoCrop<MaskImageType>(m_Mediastinum, GetBackgroundValue());
 271   // Resize input
 272   m_Input = clitk::ResizeImageLike<ImageType>(m_Input, m_Mediastinum, GetBackgroundValue());
 273
 274   // End
 275   StopCurrentStep<ImageType>(m_Input);
 276 }
 277 //--------------------------------------------------------------------
 278
 279
 280 //--------------------------------------------------------------------
 281 template <class TImageType>
 282 void
 283 clitk::ExtractMediastinalVesselsFilter<TImageType>::
 284 //SearchBrachioCephalicArtery(int & BCA_first_slice, LabelType & BCA_first_label) {
 285 TrackBifurcationFromPoint(MaskImagePointer & recon,
 286                           std::vector<MaskSlicePointer> & slices_recon,
 287                           MaskImagePointType point3D,
 288                           LabelType newLabel,
 289                           std::vector<MaskImagePointType> & bifurcations) {
 290   StartNewStep("Search for BCA first slice and label");
 291   DD(newLabel);
 292
 293   // Extract slices
 294   //  std::vector<MaskSlicePointer> slices_recon;
 295   //clitk::ExtractSlices<MaskImageType>(recon, 2, slices_recon);
 296
 297   // Find first slice
 298   MaskImageIndexType index;
 299   recon->TransformPhysicalPointToIndex(point3D, index);
 300   DD(point3D);
 301   DD(index);
 302
 303   uint i=index[2];
 304   bool found = false;
 305   LabelType previous_largest_label=recon->GetPixel(index);
 306   DD(previous_largest_label);
 307   do {
 308     DD(i);
 309     // Consider current reconstructed slice
 310     MaskSlicePointer s = slices_recon[i];
 311     MaskSlicePointer previous;
 312     if (i==index[2]) previous = s;
 313     else previous = slices_recon[i-1];
 314
 315     // Get centroids of the labels in the current slice
 316     typedef typename MaskSliceType::PointType SlicePointType;
 317     static const unsigned int Dim = MaskSliceType::ImageDimension;
 318     typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType;
 319     typedef itk::LabelMap< LabelObjectType > LabelMapType;
 320     typedef itk::LabelImageToLabelMapFilter<MaskSliceType, LabelMapType> ImageToMapFilterType;
 321     typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New();
 322     typedef itk::ShapeLabelMapFilter<LabelMapType, MaskSliceType> ShapeFilterType;
 323     typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New();
 324     imageToLabelFilter->SetBackgroundValue(GetBackgroundValue());
 325     imageToLabelFilter->SetInput(s);
 326     statFilter->SetInput(imageToLabelFilter->GetOutput());
 327     statFilter->Update();
 328     typename LabelMapType::Pointer labelMap = statFilter->GetOutput();
 329
 330     // Look what centroid inside the previous largest one
 331     std::vector<SlicePointType> centroids;
 332     std::vector<LabelType> centroids_label;
 333     for(uint c=0; c<labelMap->GetNumberOfLabelObjects(); c++) {
 334       int label = labelMap->GetLabels()[c];
 335       DD(label);
 336       SlicePointType center = labelMap->GetLabelObject(label)->GetCentroid();
 337       SlicePointType center_previous = center;
 338       center_previous[2] -= m_Input->GetSpacing()[2];
 339       // Get label into previous slice
 340       typename MaskSliceType::IndexType index;
 341       previous->TransformPhysicalPointToIndex(center_previous, index);
 342       LabelType l = previous->GetPixel(index);
 343       DD(l);
 344       if (l == previous_largest_label) {
 345         centroids.push_back(center);
 346         centroids_label.push_back(label);
 347       }
 348     }
 349     DD(centroids.size());
 350
 351     // If several centroids, we found a bifurcation
 352     if (centroids.size() > 1) {
 353       found = true;
 354       for(uint c=0; c<centroids.size(); c++) {
 355         ImagePointType bif;
 356         clitk::PointsUtils<MaskImageType>::Convert2DTo3D(centroids[c], m_Mask, i, bif);
 357         bifurcations.push_back(bif);
 358         s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[c], newLabel+c+1, true);
 359         //        slices_recon[i] = s; // (useful ?)
 360       }
 361       DD("FOUND");
 362     }
 363
 364     // if only one centroids, we change the current image with the current label
 365     if (centroids.size() == 1) {
 366       s = clitk::SetBackground<MaskSliceType, MaskSliceType>(s, s, centroids_label[0], newLabel, true);
 367       previous_largest_label = newLabel;
 368       /*typedef itk::BinaryThresholdImageFilter<MaskSliceType, MaskSliceType> BinarizeFilterType;
 369       typename BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New();
 370       binarizeFilter->SetInput(s);
 371       binarizeFilter->SetLowerThreshold(centroids_label[0]);
 372       binarizeFilter->SetUpperThreshold(centroids_label[0]+1);
 373       binarizeFilter->SetInsideValue(previous_largest_label);
 374       binarizeFilter->SetOutsideValue(GetBackgroundValue());
 375       binarizeFilter->Update();
 376       s = binarizeFilter->GetOutput();*/
 377       slices_recon[i] = s; // (not useful ?)
 378     }
 379
 380     if (centroids.size() == 0) {
 381       DD("no centroid, I stop");
 382       found = true;
 383     }
 384
 385     if (i == slices_recon.size()-1) found = true;
 386
 387     // iterate
 388     ++i;
 389   } while (!found);
 390
 391
 392   //MaskImageType::Pointer rr = clitk::JoinSlices<MaskImageType>(slices_recon, m_Mask, 2);
 393   //clitk::writeImage<MaskImageType>(rr, "recon2.mhd");
 394
 395   // End
 396   StopCurrentStep();
 397 }
 398 //--------------------------------------------------------------------
 399
 400
 401
 402 #endif //#define CLITKEXTRACTMEDIASTINALVESSELSFILTER_TXX