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
  23 //--------------------------------------------------------------------
  24 template <class ImageType>
  25 clitk::SliceBySliceRelativePositionFilter<ImageType>::
  26 SliceBySliceRelativePositionFilter():
  27   clitk::FilterBase(),
  28   itk::ImageToImageFilter<ImageType, ImageType>()
  29 {
  30   this->SetNumberOfRequiredInputs(2);
  31   SetDirection(2);
  32   SetObjectBackgroundValue(0);
  33   SetFuzzyThreshold(0.6);
  34   SetOrientationType(RelPosFilterType::LeftTo);
  35   SetIntermediateSpacing(10);
  36   ResampleBeforeRelativePositionFilterOff();
  37 }
  38 //--------------------------------------------------------------------
  39
  40
  41 //--------------------------------------------------------------------
  42 template <class ImageType>
  43 void
  44 clitk::SliceBySliceRelativePositionFilter<ImageType>::
  45 SetInput(const ImageType * image)
  46 {
  47   // Process object is not const-correct so the const casting is required.
  48   this->SetNthInput(0, const_cast<ImageType *>(image));
  49 }
  50 //--------------------------------------------------------------------
  51
  52
  53 //--------------------------------------------------------------------
  54 template <class ImageType>
  55 void
  56 clitk::SliceBySliceRelativePositionFilter<ImageType>::
  57 SetInputObject(const ImageType * image)
  58 {
  59   // Process object is not const-correct so the const casting is required.
  60   this->SetNthInput(1, const_cast<ImageType *>(image));
  61 }
  62 //--------------------------------------------------------------------
  63
  64
  65 //--------------------------------------------------------------------
  66 template <class ImageType>
  67 void
  68 clitk::SliceBySliceRelativePositionFilter<ImageType>::
  69 GenerateOutputInformation()
  70 {
  71   ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
  72   ImagePointer outputImage = this->GetOutput(0);
  73   outputImage->SetRegions(input->GetLargestPossibleRegion());
  74 }
  75 //--------------------------------------------------------------------
  76
  77
  78 //--------------------------------------------------------------------
  79 template <class ImageType>
  80 void
  81 clitk::SliceBySliceRelativePositionFilter<ImageType>::
  82 GenerateInputRequestedRegion()
  83 {
  84   // Call default
  85   itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
  86   // Get input pointers and set requested region to common region
  87   ImagePointer input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
  88   ImagePointer input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
  89   input1->SetRequestedRegion(input1->GetLargestPossibleRegion());
  90   input2->SetRequestedRegion(input2->GetLargestPossibleRegion());
  91 }
  92 //--------------------------------------------------------------------
  93
  94
  95 //--------------------------------------------------------------------
  96 template <class ImageType>
  97 void
  98 clitk::SliceBySliceRelativePositionFilter<ImageType>::
  99 GenerateData()
 100 {
 101   // Get input pointer
 102   input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
 103   object = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
 104
 105   //--------------------------------------------------------------------
 106   // Resample object to the same spacing than input
 107   if (!clitk::HaveSameSpacing<ImageType, ImageType>(object, input)) {
 108     StartNewStep("Resample object to the same spacing than input");
 109     m_working_object = clitk::ResampleImageSpacing<ImageType>(object, input->GetSpacing());
 110     StopCurrentStep<ImageType>(m_working_object);
 111   }
 112   else {
 113     DD("no resampling");
 114     m_working_object = object;
 115   }
 116
 117   //--------------------------------------------------------------------
 118   // Pad object to the same size than input
 119   if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(m_working_object, input)) {
 120     StartNewStep("Pad object to the same size than input");
 121     m_working_object = clitk::EnlargeImageLike<ImageType>(m_working_object,
 122                                                           input,
 123                                                           GetObjectBackgroundValue());
 124     StopCurrentStep<ImageType>(m_working_object);
 125   }
 126   else {
 127     DD("no pad");
 128   }
 129
 130   /*
 131     - extract vector of slices in input, in object
 132     - slice by slice rel position
 133     - joint result
 134     - post process
 135   */
 136
 137
 138   //--------------------------------------------------------------------
 139   // Extract input slices
 140   StartNewStep("Extract input slices");
 141   typedef clitk::ExtractSliceFilter<ImageType> ExtractSliceFilterType;
 142   typename ExtractSliceFilterType::Pointer extractSliceFilter = ExtractSliceFilterType::New();
 143   extractSliceFilter->SetInput(input);
 144   extractSliceFilter->SetDirection(GetDirection());
 145   extractSliceFilter->Update();
 146   typedef typename ExtractSliceFilterType::SliceType SliceType;
 147   std::vector<typename SliceType::Pointer> mInputSlices;
 148   extractSliceFilter->GetOutputSlices(mInputSlices);
 149   StopCurrentStep<SliceType>(mInputSlices[0]);
 150
 151   //--------------------------------------------------------------------
 152   // Extract object slices
 153   StartNewStep("Extract object slices");
 154   extractSliceFilter = ExtractSliceFilterType::New();
 155   extractSliceFilter->SetInput(object);
 156   extractSliceFilter->SetDirection(GetDirection());
 157   extractSliceFilter->Update();
 158   std::vector<typename SliceType::Pointer> mObjectSlices;
 159   extractSliceFilter->GetOutputSlices(mObjectSlices);
 160   StopCurrentStep<SliceType>(mObjectSlices[0]);
 161
 162   //--------------------------------------------------------------------
 163   // Perform slice by slice relative position
 164   StartNewStep("Perform slice by slice relative position");
 165   for(unsigned int i=0; i<mInputSlices.size(); i++) {
 166     // DD(i);
 167     //     DD(mInputSlices[i]->GetOrigin());
 168     //     writeImage<SliceType>(mInputSlices[i], "inp"+clitk::toString(i)+".mhd");
 169
 170     // Select main CC in each object slice : this should be the main bronchus
 171     mObjectSlices[i] = Labelize<SliceType>(mObjectSlices[i], 0, true, 1);
 172     mObjectSlices[i] = KeepLabels<SliceType>(mObjectSlices[i], 0, 1, 1, 1, true);
 173
 174     // Relative position
 175     typedef clitk::AddRelativePositionConstraintToLabelImageFilter<SliceType> RelPosFilterType;
 176     typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New();
 177     relPosFilter->VerboseStepOff();
 178     relPosFilter->WriteStepOff();
 179     relPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId());
 180     relPosFilter->SetInput(mInputSlices[i]);
 181     relPosFilter->SetInputObject(mObjectSlices[i]);
 182     relPosFilter->SetOrientationType(this->GetOrientationType());
 183     relPosFilter->SetIntermediateSpacing(this->GetIntermediateSpacing());
 184     relPosFilter->SetResampleBeforeRelativePositionFilter(this->GetResampleBeforeRelativePositionFilter());
 185     relPosFilter->SetFuzzyThreshold(this->GetFuzzyThreshold());
 186     relPosFilter->AutoCropOff(); // important ! because we join the slices after this loop
 187     relPosFilter->Update();
 188     // writeImage<SliceType>(relPosFilter->GetOutput(), "inp-after"+clitk::toString(i)+".mhd");
 189     mInputSlices[i] = relPosFilter->GetOutput();
 190   }
 191   DD(this->GetIntermediateSpacing());
 192   DD(this->GetResampleBeforeRelativePositionFilter());
 193   DD("End slice");
 194
 195   typedef itk::JoinSeriesImageFilter<SliceType, ImageType> JoinSeriesFilterType;
 196   typename JoinSeriesFilterType::Pointer joinFilter = JoinSeriesFilterType::New();
 197   joinFilter->SetOrigin(input->GetOrigin()[GetDirection()]);
 198   joinFilter->SetSpacing(input->GetSpacing()[GetDirection()]);
 199   for(unsigned int i=0; i<mInputSlices.size(); i++) {
 200   // DD(mInputSlices[i]->GetLargestPossibleRegion().GetIndex());
 201 //   DD(mInputSlices[i]->GetLargestPossibleRegion().GetSize());
 202 //   DD(mInputSlices[i]->GetRequestedRegion().GetIndex());
 203 //   DD(mInputSlices[i]->GetRequestedRegion().GetSize());
 204     joinFilter->PushBackInput(mInputSlices[i]);
 205     //SetInput(i, mInputSlices[i]);
 206   }
 207   DD("before update");
 208   joinFilter->Update();
 209   DD("after update");
 210   m_working_input = joinFilter->GetOutput();
 211
 212   // Update the origin
 213   DD(input->GetSpacing());
 214   DD(input->GetOrigin());
 215   DD(mInputSlices[0]->GetSpacing());
 216   DD(mInputSlices[0]->GetOrigin());
 217   DD(m_working_input->GetSpacing());
 218   DD(m_working_input->GetOrigin());
 219   // typename ImageType::PointType origin = m_working_input->GetOrigin();
 220 //   origin[GetDirection()] = input->GetOrigin()[GetDirection()];
 221 //   m_working_input->SetOrigin(origin);
 222 //   DD(m_working_input->GetOrigin());
 223   StopCurrentStep<ImageType>(m_working_input);
 224
 225   //--------------------------------------------------------------------
 226   //--------------------------------------------------------------------
 227   // Final Step -> set output
 228   this->SetNthOutput(0, m_working_input);
 229   return;
 230 }
 231 //--------------------------------------------------------------------
 232