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