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