itk/clitkSliceBySliceRelativePositionFilter.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 // clitk
  20 #include "clitkCropLikeImageFilter.h"
  21 #include "clitkSegmentationUtils.h"
  22 #include "clitkExtractSliceFilter.h"
  23 #include "clitkResampleImageWithOptionsFilter.h"
  24
  25 // itk
  26 #include <itkJoinSeriesImageFilter.h>
  27
  28 //--------------------------------------------------------------------
  29 template <class ImageType>
  30 clitk::SliceBySliceRelativePositionFilter<ImageType>::
  31 SliceBySliceRelativePositionFilter():
  32   clitk::AddRelativePositionConstraintToLabelImageFilter<ImageType>()
  33 {
  34   SetDirection(2);
  35   UniqueConnectedComponentBySliceFlagOff();
  36   SetIgnoreEmptySliceObjectFlag(false);
  37   UseTheLargestObjectCCLFlagOff();
  38   this->VerboseStepFlagOff();
  39   this->WriteStepFlagOff();
  40   this->SetCombineWithOrFlag(false);
  41   ObjectCCLSelectionFlagOff();
  42   SetObjectCCLSelectionDimension(0);
  43   SetObjectCCLSelectionDirection(1);
  44   ObjectCCLSelectionIgnoreSingleCCLFlagOff();
  45 }
  46 //--------------------------------------------------------------------
  47
  48
  49 //--------------------------------------------------------------------
  50 template <class ImageType>
  51 void
  52 clitk::SliceBySliceRelativePositionFilter<ImageType>::
  53 SetInput(const ImageType * image)
  54 {
  55   // Process object is not const-correct so the const casting is required.
  56   this->SetNthInput(0, const_cast<ImageType *>(image));
  57 }
  58 //--------------------------------------------------------------------
  59
  60
  61 //--------------------------------------------------------------------
  62 template <class ImageType>
  63 void
  64 clitk::SliceBySliceRelativePositionFilter<ImageType>::
  65 SetInputObject(const ImageType * image)
  66 {
  67   // Process object is not const-correct so the const casting is required.
  68   this->SetNthInput(1, const_cast<ImageType *>(image));
  69 }
  70 //--------------------------------------------------------------------
  71
  72
  73 //--------------------------------------------------------------------
  74 template <class ImageType>
  75 void
  76 clitk::SliceBySliceRelativePositionFilter<ImageType>::
  77 PrintOptions(std::ostream & os)
  78 {
  79   os << "Slice direction = " << this->GetDirection() << std::endl
  80      << "BG value        = " << this->GetBackgroundValue() << std::endl;
  81   for(int i=0; i<this->GetNumberOfAngles(); i++) {
  82     os << "Orientation     = " << this->GetOrientationTypeString()[i] << std::endl;
  83     os << "Angles     = " << clitk::rad2deg(this->GetAngle1InRad(i))
  84        << " " << clitk::rad2deg(this->GetAngle2InRad(i)) << std::endl;
  85   }
  86   os << "InverseOrientationFlag  = " << this->GetInverseOrientationFlag() << std::endl
  87      << "SpacingFlag     = " << this->GetIntermediateSpacingFlag() << std::endl
  88      << "Spacing         = " << this->GetIntermediateSpacing() << std::endl
  89      << "FuzzyThreshold  = " << this->GetFuzzyThreshold() << std::endl
  90      << "UniqueConnectedComponentBySliceFlag  = " << this->GetUniqueConnectedComponentBySliceFlag() << std::endl
  91      << "AutoCropFlag    = " << this->GetAutoCropFlag() << std::endl
  92      << "RemoveObjectFlag= " << this->GetRemoveObjectFlag() << std::endl
  93      << "CombineWithOrFlag = " << this->GetCombineWithOrFlag() << std::endl
  94      << "UseTheLargestObjectCCLFlag = " << this->GetUseTheLargestObjectCCLFlag() << std::endl
  95      << "ObjectCCLSelectionFlag = " << this->GetObjectCCLSelectionFlag() << std::endl
  96      << "ObjectCCLSelectionDimension = " << this->GetObjectCCLSelectionDimension() << std::endl
  97      << "ObjectCCLSelectionIgnoreSingleCCLFlag = " << this->GetObjectCCLSelectionIgnoreSingleCCLFlag() << std::endl
  98      << "IgnoreEmptySliceObjectFlag = " << this->GetIgnoreEmptySliceObjectFlag() << std::endl;
  99 }
 100 //--------------------------------------------------------------------
 101
 102
 103 //--------------------------------------------------------------------
 104 template <class ImageType>
 105 void
 106 clitk::SliceBySliceRelativePositionFilter<ImageType>::
 107 GenerateInputRequestedRegion()
 108 {
 109   // Call default
 110   itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
 111   // Get input pointers and set requested region to common region
 112   ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 113   ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 114   input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
 115   input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
 116 }
 117 //--------------------------------------------------------------------
 118
 119
 120 //--------------------------------------------------------------------
 121 template <class ImageType>
 122 void
 123 clitk::SliceBySliceRelativePositionFilter<ImageType>::
 124 GenerateOutputInformation()
 125 {
 126   if (this->GetVerboseOptionFlag()) {
 127     PrintOptions();
 128   }
 129
 130   // Get input pointer
 131   input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 132   object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 133
 134   //--------------------------------------------------------------------
 135   // Resample object to the same spacing than input
 136   if (!clitk::HaveSameSpacing<ImageType, ImageType>(object, input)) {
 137     this->StartNewStep("Resample object to the same spacing than input");
 138     m_working_object = clitk::ResampleImageSpacing<ImageType>(object, input->GetSpacing());
 139     this->template StopCurrentStep<ImageType>(m_working_object);
 140   }
 141   else {
 142     m_working_object = object;
 143   }
 144
 145   //--------------------------------------------------------------------
 146   // Pad object to the same size than input
 147   if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(m_working_object, input)) {
 148     this->StartNewStep("Pad object to the same size than input");
 149     m_working_object = clitk::ResizeImageLike<ImageType>(m_working_object,
 150                                                          input,
 151                                                          this->GetObjectBackgroundValue());
 152     this->template StopCurrentStep<ImageType>(m_working_object);
 153   }
 154   else {
 155   }
 156
 157   /*
 158     - extract vector of slices in input, in object
 159     - slice by slice rel position
 160     - joint result
 161     - post process
 162   */
 163
 164
 165   //--------------------------------------------------------------------
 166   // Extract input slices
 167   this->StartNewStep("Extract input slices");
 168   typedef clitk::ExtractSliceFilter<ImageType> ExtractSliceFilterType;
 169   typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 170   extractSliceFilter->SetInput(input);
 171   extractSliceFilter->SetDirection(GetDirection());
 172   extractSliceFilter->Update();
 173   typedef typename ExtractSliceFilterType::SliceType SliceType;
 174   std::vector<typename SliceType::Pointer> mInputSlices;
 175   extractSliceFilter->GetOutputSlices(mInputSlices);
 176   this->template StopCurrentStep<SliceType>(mInputSlices[0]);
 177
 178   //--------------------------------------------------------------------
 179   // Extract object slices
 180   this->StartNewStep("Extract object slices");
 181   extractSliceFilter = ExtractSliceFilterType::New();
 182   extractSliceFilter->SetInput(m_working_object);//object);
 183   extractSliceFilter->SetDirection(GetDirection());
 184   extractSliceFilter->Update();
 185   std::vector<typename SliceType::Pointer> mObjectSlices;
 186   extractSliceFilter->GetOutputSlices(mObjectSlices);
 187   this->template StopCurrentStep<SliceType>(mObjectSlices[0]);
 188
 189   //--------------------------------------------------------------------
 190   // Prepare fuzzy slices (if needed)
 191   std::vector<typename FloatSliceType::Pointer> mFuzzyMapSlices;
 192   mFuzzyMapSlices.resize(mInputSlices.size());
 193
 194   //--------------------------------------------------------------------
 195   // Perform slice by slice relative position
 196   this->StartNewStep("Perform slice by slice relative position");
 197   for(unsigned int i=0; i<mInputSlices.size(); i++) {
 198
 199     // Count the number of CCL (allow to ignore empty slice)
 200     int nb=0;
 201     mObjectSlices[i] = LabelizeAndCountNumberOfObjects<SliceType>(mObjectSlices[i], 0, true, 1, nb);
 202
 203     // If no object and empty slices :
 204     if ((nb==0) && (this->GetFuzzyMapOnlyFlag())) {
 205       typename FloatSliceType::Pointer one = FloatSliceType::New();
 206       one->CopyInformation(mObjectSlices[0]);
 207       one->SetRegions(mObjectSlices[0]->GetLargestPossibleRegion());
 208       one->Allocate();
 209       one->FillBuffer(2.0);
 210       mFuzzyMapSlices[i] = one;
 211     }
 212     else {
 213       if ((!GetIgnoreEmptySliceObjectFlag()) || (nb!=0)) {
 214
 215         // Select or not a single CCL ?
 216         if (GetUseTheLargestObjectCCLFlag()) {
 217           mObjectSlices[i] = KeepLabels<SliceType>(mObjectSlices[i], 0, 1, 1, 1, true);
 218         }
 219
 220         // Select a single according to a position if more than one CCL
 221         if (GetObjectCCLSelectionFlag()) {
 222           // if several CCL, choose the most extrema according a direction,
 223           // if not -> should we consider this slice ?
 224           if (nb<2) {
 225             if (GetObjectCCLSelectionIgnoreSingleCCLFlag()) {
 226               mObjectSlices[i] = SetBackground<SliceType, SliceType>(mObjectSlices[i], mObjectSlices[i],
 227                                                                      1, this->GetBackgroundValue(),
 228                                                                      true);
 229             }
 230           }
 231           int dim = GetObjectCCLSelectionDimension();
 232           int direction = GetObjectCCLSelectionDirection();
 233           std::vector<typename SliceType::PointType> centroids;
 234           ComputeCentroids<SliceType>(mObjectSlices[i], this->GetBackgroundValue(), centroids);
 235           uint index=1;
 236           for(uint j=1; j<centroids.size(); j++) {
 237             if (direction == 1) {
 238               if (centroids[j][dim] > centroids[index][dim]) index = j;
 239             }
 240             else {
 241               if (centroids[j][dim] < centroids[index][dim]) index = j;
 242             }
 243           }
 244           for(uint v=1; v<centroids.size(); v++) {
 245             if (v != index) {
 246               mObjectSlices[i] = SetBackground<SliceType, SliceType>(mObjectSlices[i], mObjectSlices[i],
 247                                                                      (char)v, this->GetBackgroundValue(),
 248                                                                      true);
 249             }
 250           }
 251         } // end GetbjectCCLSelectionFlag = true
 252
 253         // Relative position
 254         typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> RelPosFilterType;
 255         typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 256
 257         relPosFilter->VerboseStepFlagOff();
 258         relPosFilter->WriteStepFlagOff();
 259         relPosFilter->SetBackgroundValue(this->GetBackgroundValue());
 260         relPosFilter->SetInput(mInputSlices[i]);
 261         relPosFilter->SetInputObject(mObjectSlices[i]);
 262         relPosFilter->SetRemoveObjectFlag(this->GetRemoveObjectFlag());
 263         // This flag (InverseOrientation) *must* be set before
 264         // AddOrientation because AddOrientation can change it.
 265         relPosFilter->SetInverseOrientationFlag(this->GetInverseOrientationFlag());
 266         for(int j=0; j<this->GetNumberOfAngles(); j++) {
 267           //          relPosFilter->AddOrientationTypeString(this->GetOrientationTypeString(j));
 268           relPosFilter->AddAnglesInRad(this->GetAngle1InRad(j), this->GetAngle2InRad(j));
 269           // DD(this->GetOrientationTypeString(j));
 270         }
 271         // DD(this->GetInverseOrientationFlag());
 272         //relPosFilter->SetOrientationType(this->GetOrientationType());
 273         relPosFilter->SetIntermediateSpacing(this->GetIntermediateSpacing());
 274         relPosFilter->SetIntermediateSpacingFlag(this->GetIntermediateSpacingFlag());
 275         relPosFilter->SetFuzzyThreshold(this->GetFuzzyThreshold());
 276         relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
 277         relPosFilter->SetCombineWithOrFlag(this->GetCombineWithOrFlag());
 278
 279         // should we stop after fuzzy map ?
 280         relPosFilter->SetFuzzyMapOnlyFlag(this->GetFuzzyMapOnlyFlag());
 281
 282         // Go !
 283         relPosFilter->Update();
 284
 285         // If we stop after the fuzzy map, store the fuzzy slices
 286         if (this->GetFuzzyMapOnlyFlag()) {
 287           mFuzzyMapSlices[i] = relPosFilter->GetFuzzyMap();
 288           // writeImage<FloatSliceType>(mFuzzyMapSlices[i], "slice_"+toString(i)+".mha");
 289         }
 290         else  {
 291           mInputSlices[i] = relPosFilter->GetOutput();
 292           // Select main CC if needed
 293           if (GetUniqueConnectedComponentBySliceFlag()) {
 294             mInputSlices[i] = Labelize<SliceType>(mInputSlices[i], 0, true, 1);
 295             mInputSlices[i] = KeepLabels<SliceType>(mInputSlices[i], 0, 1, 1, 1, true);
 296           }
 297
 298         }
 299
 300       }
 301       /*
 302       // Select unique CC according to the most in a given direction
 303       if (GetUniqueConnectedComponentBySliceAccordingToADirection()) {
 304       int nb;
 305       mInputSlices[i] = LabelizeAndCountNumberOfObjects<SliceType>(mInputSlices[i], 0, true, 1, nb);
 306       std::vector<typename ImageType::PointType> & centroids;
 307       ComputeCentroids
 308       }
 309       */
 310     }
 311   }
 312
 313   // Join the fuzzy map if needed
 314   if (this->GetFuzzyMapOnlyFlag()) {
 315     this->m_FuzzyMap = clitk::JoinSlices<FloatImageType>(mFuzzyMapSlices, input, GetDirection());
 316     this->template StopCurrentStep<FloatImageType>(this->m_FuzzyMap);
 317     return;
 318   }
 319
 320   // Join the slices
 321   m_working_input = clitk::JoinSlices<ImageType>(mInputSlices, input, GetDirection());
 322   this->template StopCurrentStep<ImageType>(m_working_input);
 323
 324   //--------------------------------------------------------------------
 325   // Step 7: autocrop
 326   if (this->GetAutoCropFlag()) {
 327     this->StartNewStep("Final AutoCrop");
 328     typedef clitk::AutoCropFilter<ImageType> CropFilterType;
 329     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 330     cropFilter->SetInput(m_working_input);
 331     cropFilter->ReleaseDataFlagOff();
 332     cropFilter->Update();
 333     m_working_input = cropFilter->GetOutput();
 334     this->template StopCurrentStep<ImageType>(m_working_input);
 335   }
 336
 337   // Update output info
 338   this->GetOutput(0)->SetRegions(m_working_input->GetLargestPossibleRegion());
 339 }
 340 //--------------------------------------------------------------------
 341
 342
 343 //--------------------------------------------------------------------
 344 template <class ImageType>
 345 void
 346 clitk::SliceBySliceRelativePositionFilter<ImageType>::
 347 GenerateData()
 348 {
 349   // Get input pointer
 350   //--------------------------------------------------------------------
 351   //--------------------------------------------------------------------
 352   // Final Step -> set output
 353   //this->SetNthOutput(0, m_working_input);
 354   if (this->GetFuzzyMapOnlyFlag()) return; // no output in this case
 355   this->GraftOutput(m_working_input);
 356   return;
 357 }
 358 //--------------------------------------------------------------------
 359