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