vv/vvMidPosition.cxx

   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://www.centreleonberard.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 #include <QApplication>
  19
  20 #include <itkWarpImageFilter.h>
  21 #if ( ITK_VERSION_MAJOR < 5 )
  22 #include <itkVectorResampleImageFilter.h>
  23 #else
  24 #include <itkResampleImageFilter.h>
  25 #endif
  26
  27 #include "vvMidPosition.h"
  28 #include "clitkCommon.h"
  29 #include "vvFromITK.h"
  30 #include "vvToITK.h"
  31 #include "clitkInvertVFFilter.h"
  32
  33 vvMidPosition::vvMidPosition() :
  34   slicer_manager(0), error(false),
  35   reference_image_index(0),
  36   p_bar("Computing mid-position...","Cancel", 0,6),
  37   progress(0)
  38 {
  39 }
  40
  41 void vvMidPosition::Update()
  42 {
  43   this->start();
  44   while (this->isRunning()) {
  45     this->wait(50);
  46     this->update_progress();
  47     qApp->processEvents();
  48   }
  49 }
  50
  51 //Common typedefs
  52 typedef itk::Vector<float,3> VFPixelType;
  53 typedef itk::Image<VFPixelType,4> VFType;
  54 typedef itk::Image<VFPixelType,3> OutputVFType;
  55
  56 ///Internal functions
  57
  58 ///This function averages a vector field along the temporal dimension
  59 // the progress parameter is a reference to a counter used to update the progress dialog
  60 itk::Image<itk::Vector<float,3>,3>::Pointer AverageField(itk::Image<itk::Vector<float,3>,4>::ConstPointer vf,int &progress);
  61 ///Warps an image frame using the vf passed in parameter
  62 template<class ImagePixelType> vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,int reference_image_index);
  63
  64 void vvMidPosition::run()
  65 {
  66   error=true;
  67   if(slicer_manager->GetImage()->GetNumberOfDimensions() != 4)
  68     error_message="Computation of midposition is only supported for 4D image sequences.";
  69   else if (!slicer_manager->GetVF())
  70     error_message="A VF is required for midposition computation";
  71   else if (slicer_manager->GetVF()->GetScalarTypeAsITKString() != "float")
  72     error_message="Please use a vector field of type float.";
  73   else {
  74     VFType::ConstPointer vf = vvImageToITK<VFType>(slicer_manager->GetVF());
  75     OutputVFType::Pointer avg_vf=AverageField(vf,this->progress);
  76     clitk::InvertVFFilter<OutputVFType,OutputVFType>::Pointer inv_filter=
  77       clitk::InvertVFFilter<OutputVFType,OutputVFType>::New();
  78     inv_filter->SetInput(avg_vf);
  79     inv_filter->Update();
  80     progress++;
  81     if (slicer_manager->GetImage()->GetScalarTypeAsITKString() == "short")
  82       this->output=WarpRefImage<short>(inv_filter->GetOutput(),slicer_manager->GetImage(),reference_image_index);
  83     else {
  84       error_message="Unsupported image pixel type.";
  85       return;
  86     }
  87     progress++;
  88     error=false;
  89   }
  90 }
  91
  92 template<class ImagePixelType>
  93 vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,int reference_image_index)
  94 {
  95   typedef itk::Image<ImagePixelType,3> ImageType;
  96   typedef itk::WarpImageFilter<ImageType,ImageType,itk::Image<VFPixelType,3> > FilterType;
  97
  98   typename ImageType::ConstPointer input = vvSingleFrameToITK<3,ImagePixelType>(image,reference_image_index);
  99
 100   //We resample the VF because itk's warp filter doesn't like it when the vf and the image have
 101   //different spacings
 102 #if ( ITK_VERSION_MAJOR < 5 )
 103   typename itk::VectorResampleImageFilter<OutputVFType, OutputVFType >::Pointer
 104   resampler =itk::VectorResampleImageFilter<OutputVFType, OutputVFType >::New();
 105 #else
 106   typename itk::ResampleImageFilter<OutputVFType, OutputVFType >::Pointer
 107   resampler =itk::ResampleImageFilter<OutputVFType, OutputVFType >::New();
 108 #endif
 109   resampler->SetInput(vf);
 110   resampler->SetOutputSpacing(input->GetSpacing());
 111   resampler->SetOutputOrigin(vf->GetOrigin());
 112   //Calculate the new size so that it contains the vf
 113   typename ImageType::SizeType newSize;
 114   for (unsigned int i=0 ; i <3; i++)
 115     newSize[i]=(unsigned int) (vf->GetLargestPossibleRegion().GetSize()[i]*vf->GetSpacing()[i]/input->GetSpacing()[i]);
 116   resampler->SetSize( newSize);
 117
 118   typename FilterType::Pointer warp_filter = FilterType::New();
 119   warp_filter->SetInput(input);
 120   warp_filter->SetDisplacementField(resampler->GetOutput());
 121   warp_filter->SetOutputSpacing(input->GetSpacing());
 122   warp_filter->SetOutputOrigin(input->GetOrigin());
 123   warp_filter->SetOutputSize(input->GetLargestPossibleRegion().GetSize());
 124   warp_filter->SetEdgePaddingValue(-1000);
 125   warp_filter->Update();
 126   return vvImageFromITK<3,ImagePixelType>(warp_filter->GetOutput());
 127 }
 128
 129 itk::Image<itk::Vector<float,3>,3>::Pointer AverageField(itk::Image<itk::Vector<float,3>,4>::ConstPointer vf, int& progress)
 130 {
 131   progress++;
 132
 133   VFType::RegionType region4D=vf->GetLargestPossibleRegion();
 134   VFType::RegionType::SizeType size4D=region4D.GetSize();
 135   VFType::IndexType index4D=region4D.GetIndex();
 136   VFType::SpacingType spacing4D=vf->GetSpacing();
 137   VFType::PointType origin4D=vf->GetOrigin();
 138
 139   OutputVFType::RegionType region;
 140   OutputVFType::RegionType::SizeType size;
 141   OutputVFType::IndexType index;
 142   OutputVFType::SpacingType spacing;
 143   OutputVFType::PointType origin;
 144   for (unsigned int i=0; i< 3; i++) {
 145     size[i]=size4D[i];
 146     index[i]=index4D[i];
 147     spacing[i]=spacing4D[i];
 148     origin[i]=origin4D[i];
 149   }
 150   region.SetSize(size);
 151   region.SetIndex(index);
 152   OutputVFType::Pointer output= OutputVFType::New();
 153   output->SetRegions(region);
 154   output->SetSpacing(spacing);
 155   output->SetOrigin(origin);
 156   output->Allocate();
 157   progress++;
 158
 159
 160   // Region iterators
 161   typedef itk::ImageRegionConstIterator<VFType> IteratorType;
 162   std::vector<IteratorType> iterators(size4D[3]);
 163   for (unsigned int i=0; i< size4D[3]; i++) {
 164     VFType::RegionType regionIt=region4D;
 165     VFType::RegionType::SizeType sizeIt=regionIt.GetSize();
 166     sizeIt[3]=1;
 167     regionIt.SetSize(sizeIt);
 168     VFType::IndexType indexIt=regionIt.GetIndex();
 169     indexIt[3]=i;
 170     regionIt.SetIndex(indexIt);
 171     iterators[i]=IteratorType(vf, regionIt);
 172   }
 173   progress++;
 174
 175   typedef itk::ImageRegionIterator<OutputVFType> OutputIteratorType;
 176   OutputIteratorType avIt(output, output->GetLargestPossibleRegion());
 177
 178   // Average
 179   VFPixelType vector;
 180   VFPixelType zeroVector;//=itk::NumericTraits<VFPixelType>::Zero;
 181   for(unsigned int i=0;i <VFPixelType::Dimension; i++) zeroVector[i] = 0.0;
 182
 183   while (!(iterators[0]).IsAtEnd()) {
 184     vector=zeroVector;
 185     for (unsigned int i=0; i<size4D[3]; i++) {
 186       vector+=iterators[i].Get();
 187       ++(iterators[i]);
 188     }
 189     vector/=size4D[3];
 190     avIt.Set(vector);
 191     ++avIt;
 192   }
 193   progress++;
 194   return output;
 195 }
 196
 197
 198 void vvMidPosition::update_progress()
 199 {
 200   p_bar.setValue(progress);
 201   p_bar.show();
 202 }