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