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      << "(RP) FastFlag              = " << this->GetFastFlag() << std::endl
 100      << "(RP) Radius                = " << this->GetRadius() << std::endl;
 101 }
 102 //--------------------------------------------------------------------
 103
 104
 105 //--------------------------------------------------------------------
 106 template <class ImageType>
 107 void
 108 clitk::SliceBySliceRelativePositionFilter<ImageType>::
 109 GenerateInputRequestedRegion()
 110 {
 111   // Call default
 112   itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
 113   // Get input pointers and set requested region to common region
 114   ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 115   ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 116   input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
 117   input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
 118 }
 119 //--------------------------------------------------------------------
 120
 121
 122 //--------------------------------------------------------------------
 123 template <class ImageType>
 124 void
 125 clitk::SliceBySliceRelativePositionFilter<ImageType>::
 126 GenerateOutputInformation()
 127 {
 128   if (this->GetVerboseOptionFlag()) {
 129     PrintOptions();
 130   }
 131
 132   //  if (this->GetFuzzyMapOnlyFlag()) this->ComputeFuzzyMapFlagOn();
 133
 134   // Get input pointer
 135   input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 136   object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 137   m_working_object = object;
 138   m_working_input = input;
 139
 140   //--------------------------------------------------------------------
 141   // Resample object to the same spacing than input
 142   if (!clitk::HaveSameSpacing<ImageType, ImageType>(object, input)) {
 143     this->StartNewStep("Resample object to the same spacing than input");
 144     m_working_object = clitk::ResampleImageSpacing<ImageType>(object, input->GetSpacing());
 145     this->template StopCurrentStep<ImageType>(m_working_object);
 146   }
 147
 148   //--------------------------------------------------------------------
 149   // Resize image according to common area (except in Z)
 150   if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(m_working_object, input)) {
 151     this->StartNewStep("Resize images (union in XY and like input in Z)");
 152
 153     /* OLD STUFF
 154     this->StartNewStep("Pad object to the same size than input");
 155     m_working_object = clitk::ResizeImageLike<ImageType>(m_working_object,
 156     input,
 157     this->GetObjectBackgroundValue());
 158     this->template StopCurrentStep<ImageType>(m_working_object);
 159     */
 160
 161     // Compute union of bounding boxes in X and Y
 162     static const unsigned int dim = ImageType::ImageDimension;
 163     typedef itk::BoundingBox<unsigned long, dim> BBType;
 164     typename BBType::Pointer bb1 = BBType::New();
 165     ComputeBBFromImageRegion<ImageType>(m_working_object, m_working_object->GetLargestPossibleRegion(), bb1);
 166     typename BBType::Pointer bb2 = BBType::New();
 167     ComputeBBFromImageRegion<ImageType>(input, input->GetLargestPossibleRegion(), bb2);
 168     typename BBType::Pointer bbo = BBType::New();
 169     ComputeBBUnion<dim>(bbo, bb1, bb2);
 170
 171     //We set Z BB like input
 172     typename ImageType::PointType maxs = bbo->GetMaximum();
 173     typename ImageType::PointType mins = bbo->GetMinimum();
 174     maxs[2] = bb2->GetMaximum()[2];
 175     mins[2] = bb2->GetMinimum()[2];
 176     bbo->SetMaximum(maxs);
 177     bbo->SetMinimum(mins);
 178
 179     // Crop
 180     m_working_input = clitk::ResizeImageLike<ImageType>(input, bbo, this->GetBackgroundValue());
 181     m_working_object = clitk::ResizeImageLike<ImageType>(m_working_object,
 182                                                          m_working_input,
 183                                                          this->GetObjectBackgroundValue());
 184
 185     // End
 186     this->template StopCurrentStep<ImageType>(m_working_input);
 187   }
 188
 189   //--------------------------------------------------------------------
 190   /* Steps :
 191     - extract vector of slices in input, in object
 192     - slice by slice rel position
 193     - joint result
 194     - post process
 195   */
 196
 197   //--------------------------------------------------------------------
 198   // Extract input slices
 199   this->StartNewStep("Extract input slices");
 200   typedef clitk::ExtractSliceFilter<ImageType> ExtractSliceFilterType;
 201   typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 202   extractSliceFilter->SetInput(m_working_input);
 203   extractSliceFilter->SetDirection(GetDirection());
 204   extractSliceFilter->Update();
 205   typedef typename ExtractSliceFilterType::SliceType SliceType;
 206   std::vector<typename SliceType::Pointer> mInputSlices;
 207   extractSliceFilter->GetOutputSlices(mInputSlices);
 208   this->template StopCurrentStep<SliceType>(mInputSlices[0]);
 209
 210   //--------------------------------------------------------------------
 211   // Extract object slices
 212   this->StartNewStep("Extract object slices");
 213   extractSliceFilter = ExtractSliceFilterType::New();
 214   extractSliceFilter->SetInput(m_working_object);
 215   extractSliceFilter->SetDirection(GetDirection());
 216   extractSliceFilter->Update();
 217   std::vector<typename SliceType::Pointer> mObjectSlices;
 218   extractSliceFilter->GetOutputSlices(mObjectSlices);
 219   this->template StopCurrentStep<SliceType>(mObjectSlices[0]);
 220
 221   //--------------------------------------------------------------------
 222   // Prepare fuzzy slices (if needed)
 223   std::vector<typename FloatSliceType::Pointer> mFuzzyMapSlices;
 224   mFuzzyMapSlices.resize(mInputSlices.size());
 225
 226   //--------------------------------------------------------------------
 227   // Perform slice by slice relative position
 228   this->StartNewStep("Perform slice by slice relative position ("+toString(mInputSlices.size())+")");
 229   for(unsigned int i=0; i<mInputSlices.size(); i++) {
 230
 231     // Count the number of CCL (allow to ignore empty slice)
 232     int nb=0;
 233     mObjectSlices[i] = LabelizeAndCountNumberOfObjects<SliceType>(mObjectSlices[i], 0, true, 1, nb);
 234
 235     // If no object and empty slices and if we need the full fuzzy map, create a dummy one.
 236     if ((nb==0) && (this->GetFuzzyMapOnlyFlag())) {
 237       typename FloatSliceType::Pointer one = FloatSliceType::New();
 238       one->CopyInformation(mObjectSlices[0]);
 239       one->SetRegions(mObjectSlices[0]->GetLargestPossibleRegion());
 240       one->Allocate();
 241       one->FillBuffer(2.0);
 242       mFuzzyMapSlices[i] = one;
 243     } // End nb==0 && GetComputeFuzzyMapFlag
 244     else {
 245       if ((!GetIgnoreEmptySliceObjectFlag()) || (nb!=0)) {
 246
 247         // Select or not a single CCL ?
 248         if (GetUseTheLargestObjectCCLFlag()) {
 249           mObjectSlices[i] = KeepLabels<SliceType>(mObjectSlices[i], 0, 1, 1, 1, true);
 250         }
 251
 252         // Select a single according to a position if more than one CCL
 253         if (GetObjectCCLSelectionFlag()) {
 254           // if several CCL, choose the most extrema according a direction,
 255           // if not -> should we consider this slice ?
 256           if (nb<2) {
 257             if (GetObjectCCLSelectionIgnoreSingleCCLFlag()) {
 258               mObjectSlices[i] = SetBackground<SliceType, SliceType>(mObjectSlices[i], mObjectSlices[i],
 259                                                                      1, this->GetBackgroundValue(),
 260                                                                      true);
 261             }
 262           }
 263           int dim = GetObjectCCLSelectionDimension();
 264           int direction = GetObjectCCLSelectionDirection();
 265           std::vector<typename SliceType::PointType> centroids;
 266           ComputeCentroids<SliceType>(mObjectSlices[i], this->GetBackgroundValue(), centroids);
 267           uint index=1;
 268           for(uint j=1; j<centroids.size(); j++) {
 269             if (direction == 1) {
 270               if (centroids[j][dim] > centroids[index][dim]) index = j;
 271             }
 272             else {
 273               if (centroids[j][dim] < centroids[index][dim]) index = j;
 274             }
 275           }
 276           for(uint v=1; v<centroids.size(); v++) {
 277             if (v != index) {
 278               mObjectSlices[i] = SetBackground<SliceType, SliceType>(mObjectSlices[i], mObjectSlices[i],
 279                                                                      (char)v, this->GetBackgroundValue(),
 280                                                                      true);
 281             }
 282           }
 283         } // end GetbjectCCLSelectionFlag = true
 284
 285         // Relative position
 286         typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> RelPosFilterType;
 287         typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 288         relPosFilter->VerboseStepFlagOff();
 289         relPosFilter->WriteStepFlagOff();
 290         // relPosFilter->VerboseMemoryFlagOn();
 291         relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()+"-"+toString(i));
 292         relPosFilter->SetBackgroundValue(this->GetBackgroundValue());
 293         relPosFilter->SetInput(mInputSlices[i]);
 294         relPosFilter->SetInputObject(mObjectSlices[i]);
 295         relPosFilter->SetRemoveObjectFlag(this->GetRemoveObjectFlag());
 296         // This flag (InverseOrientation) *must* be set before
 297         // AddOrientation because AddOrientation can change it.
 298         relPosFilter->SetInverseOrientationFlag(this->GetInverseOrientationFlag());
 299         for(int j=0; j<this->GetNumberOfAngles(); j++) {
 300           relPosFilter->AddAnglesInRad(this->GetAngle1InRad(j), this->GetAngle2InRad(j));
 301         }
 302         relPosFilter->SetIntermediateSpacing(this->GetIntermediateSpacing());
 303         relPosFilter->SetIntermediateSpacingFlag(this->GetIntermediateSpacingFlag());
 304         relPosFilter->SetFuzzyThreshold(this->GetFuzzyThreshold());
 305         relPosFilter->AutoCropFlagOff(); // important ! because we join the slices after this loop
 306         relPosFilter->SetCombineWithOrFlag(this->GetCombineWithOrFlag());
 307         // should we stop after fuzzy map ?
 308         relPosFilter->SetFuzzyMapOnlyFlag(this->GetFuzzyMapOnlyFlag());
 309         //        relPosFilter->SetComputeFuzzyMapFlag(this->GetComputeFuzzyMapFlag());
 310         relPosFilter->SetFastFlag(this->GetFastFlag());
 311         relPosFilter->SetRadius(this->GetRadius());
 312
 313         // Go !
 314         relPosFilter->Update();
 315
 316         // If we stop after the fuzzy map, store the fuzzy slices
 317         if (this->GetFuzzyMapOnlyFlag()) {
 318           mFuzzyMapSlices[i] = relPosFilter->GetFuzzyMap();
 319           // writeImage<FloatSliceType>(mFuzzyMapSlices[i], "slice_"+toString(i)+".mha");
 320         }
 321
 322         // Set input slices
 323         if (!this->GetFuzzyMapOnlyFlag())  {
 324           mInputSlices[i] = relPosFilter->GetOutput();
 325           // Select main CC if needed
 326           if (GetUniqueConnectedComponentBySliceFlag()) {
 327             mInputSlices[i] = Labelize<SliceType>(mInputSlices[i], 0, true, 1);
 328             mInputSlices[i] = KeepLabels<SliceType>(mInputSlices[i], 0, 1, 1, 1, true);
 329           }
 330         }
 331
 332       }
 333
 334       /*
 335       // Select unique CC according to the most in a given direction
 336       if (GetUniqueConnectedComponentBySliceAccordingToADirection()) {
 337       int nb;
 338       mInputSlices[i] = LabelizeAndCountNumberOfObjects<SliceType>(mInputSlices[i], 0, true, 1, nb);
 339       std::vector<typename ImageType::PointType> & centroids;
 340       ComputeCentroids
 341       }
 342       */
 343
 344     } // End nb!=0 || GetComputeFuzzyMapFlagOFF
 345
 346   } // end for i mInputSlices
 347
 348   // Join the slices
 349   m_working_input = clitk::JoinSlices<ImageType>(mInputSlices, m_working_input, GetDirection());
 350   this->template StopCurrentStep<ImageType>(m_working_input);
 351
 352   // Join the fuzzy map if needed
 353   if (this->GetFuzzyMapOnlyFlag()) {
 354     this->m_FuzzyMap = clitk::JoinSlices<FloatImageType>(mFuzzyMapSlices, m_working_input, GetDirection());
 355     this->template StopCurrentStep<FloatImageType>(this->m_FuzzyMap);
 356     if (this->GetFuzzyMapOnlyFlag()) return;
 357   }
 358
 359   //--------------------------------------------------------------------
 360   // Step 7: autocrop
 361   if (this->GetAutoCropFlag()) {
 362     this->StartNewStep("Final AutoCrop");
 363     typedef clitk::AutoCropFilter<ImageType> CropFilterType;
 364     typename CropFilterType::Pointer cropFilter = CropFilterType::New();
 365     cropFilter->SetInput(m_working_input);
 366     cropFilter->ReleaseDataFlagOff();
 367     cropFilter->Update();
 368     m_working_input = cropFilter->GetOutput();
 369     this->template StopCurrentStep<ImageType>(m_working_input);
 370   }
 371
 372   // Update output info
 373   this->GetOutput(0)->SetRegions(m_working_input->GetLargestPossibleRegion());
 374 }
 375 //--------------------------------------------------------------------
 376
 377
 378 //--------------------------------------------------------------------
 379 template <class ImageType>
 380 void
 381 clitk::SliceBySliceRelativePositionFilter<ImageType>::
 382 GenerateData()
 383 {
 384   // Get input pointer
 385   //--------------------------------------------------------------------
 386   //--------------------------------------------------------------------
 387   // Final Step -> set output
 388   //this->SetNthOutput(0, m_working_input);
 389   if (this->GetFuzzyMapOnlyFlag()) return; // no output in this case
 390   this->GraftOutput(m_working_input);
 391   return;
 392 }
 393 //--------------------------------------------------------------------
 394