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