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