vv/vvMidPosition.cxx

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