vv/vvMidPosition.cxx

   1 /*=========================================================================
   2
   3  Program:   vv
   4  Language:  C++
   5  Author :   Joel Schaerer (joel.schaerer@insa-lyon.fr)
   6
   7 Copyright (C) 2008
   8 Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
   9 CREATIS-LRMN http://www.creatis.insa-lyon.fr
  10
  11 This program is free software: you can redistribute it and/or modify
  12 it under the terms of the GNU General Public License as published by
  13 the Free Software Foundation, version 3 of the License.
  14
  15 This program is distributed in the hope that it will be useful,
  16 but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18 GNU General Public License for more details.
  19
  20 You should have received a copy of the GNU General Public License
  21 along with this program.  If not, see <http://www.gnu.org/licenses/>.
  22
  23 =========================================================================*/
  24 #include <QApplication>
  25
  26 #include <itkWarpImageFilter.h>
  27 #include <itkVectorResampleImageFilter.h>
  28
  29 #include "vvMidPosition.h"
  30 #include "clitkCommon.h"
  31 #include "vvFromITK.h"
  32 #include "vvToITK.h"
  33 #include "clitkInvertVFFilter.h"
  34
  35 vvMidPosition::vvMidPosition() :
  36     slicer_manager(0), error(false),
  37     reference_image_index(0),
  38     p_bar("Computing mid-position...","Cancel", 0,6),
  39     progress(0)
  40 {
  41 }
  42
  43 void vvMidPosition::Update()
  44 {
  45     this->start();
  46     while (this->isRunning())
  47     {
  48         this->wait(50);
  49         this->update_progress();
  50         qApp->processEvents();
  51     }
  52 }
  53
  54 //Common typedefs
  55 typedef itk::Vector<float,3> VFPixelType;
  56 typedef itk::Image<VFPixelType,4> VFType;
  57 typedef itk::Image<VFPixelType,3> OutputVFType;
  58
  59 ///Internal functions
  60
  61 ///This function averages a vector field along the temporal dimension
  62 // the progress parameter is a reference to a counter used to update the progress dialog
  63 itk::Image<itk::Vector<float,3>,3>::Pointer AverageField(itk::Image<itk::Vector<float,3>,4>::ConstPointer vf,int &progress);
  64 ///Warps an image frame using the vf passed in parameter
  65 template<class ImagePixelType> vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,int reference_image_index);
  66
  67 void vvMidPosition::run()
  68 {
  69     error=true;
  70     if(slicer_manager->GetImage()->GetNumberOfDimensions() != 4)
  71         error_message="Computation of midposition is only supported for 4D image sequences.";
  72     else if (not slicer_manager->GetVF())
  73         error_message="A VF is required for midposition computation";
  74     else if (slicer_manager->GetVF()->GetScalarTypeAsString() != "float")
  75         error_message="Please use a vector field of type float.";
  76     else
  77     {
  78         VFType::ConstPointer vf = vvImageToITK<VFType>(slicer_manager->GetVF());
  79         OutputVFType::Pointer avg_vf=AverageField(vf,this->progress);
  80         clitk::InvertVFFilter<OutputVFType,OutputVFType>::Pointer inv_filter=
  81             clitk::InvertVFFilter<OutputVFType,OutputVFType>::New();
  82         inv_filter->SetInput(avg_vf);
  83         inv_filter->Update();
  84         progress++;
  85         if (slicer_manager->GetImage()->GetScalarTypeAsString() == "short")
  86             this->output=WarpRefImage<short>(inv_filter->GetOutput(),slicer_manager->GetImage(),reference_image_index);
  87         else
  88         {
  89             error_message="Unsupported image pixel type.";
  90             return;
  91         }
  92         progress++;
  93         error=false;
  94     }
  95 }
  96
  97 template<class ImagePixelType>
  98 vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,int reference_image_index)
  99 {
 100     typedef itk::Image<ImagePixelType,3> ImageType;
 101     typedef itk::WarpImageFilter<ImageType,ImageType,itk::Image<VFPixelType,3> > FilterType;
 102
 103     typename ImageType::ConstPointer input = vvSingleFrameToITK<3,ImagePixelType>(image,reference_image_index);
 104
 105     //We resample the VF because itk's warp filter doesn't like it when the vf and the image have
 106     //different spacings
 107         typename itk::VectorResampleImageFilter<OutputVFType, OutputVFType >::Pointer
 108         resampler =itk::VectorResampleImageFilter<OutputVFType, OutputVFType >::New();
 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->SetDeformationField(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     {
 146         size[i]=size4D[i];
 147         index[i]=index4D[i];
 148         spacing[i]=spacing4D[i];
 149         origin[i]=origin4D[i];
 150     }
 151     region.SetSize(size);
 152     region.SetIndex(index);
 153     OutputVFType::Pointer output= OutputVFType::New();
 154     output->SetRegions(region);
 155     output->SetSpacing(spacing);
 156     output->SetOrigin(origin);
 157     output->Allocate();
 158     progress++;
 159
 160
 161     // Region iterators
 162     typedef itk::ImageRegionConstIterator<VFType> IteratorType;
 163     std::vector<IteratorType> iterators(size4D[3]);
 164     for (unsigned int i=0; i< size4D[3]; i++)
 165     {
 166         VFType::RegionType regionIt=region4D;
 167         VFType::RegionType::SizeType sizeIt=regionIt.GetSize();
 168         sizeIt[3]=1;
 169         regionIt.SetSize(sizeIt);
 170         VFType::IndexType indexIt=regionIt.GetIndex();
 171         indexIt[3]=i;
 172         regionIt.SetIndex(indexIt);
 173         iterators[i]=IteratorType(vf, regionIt);
 174     }
 175     progress++;
 176
 177     typedef itk::ImageRegionIterator<OutputVFType> OutputIteratorType;
 178     OutputIteratorType avIt(output, output->GetLargestPossibleRegion());
 179
 180     // Average
 181     VFPixelType vector;
 182     VFPixelType zeroVector=itk::NumericTraits<VFPixelType>::Zero;
 183
 184     while (!(iterators[0]).IsAtEnd())
 185     {
 186         vector=zeroVector;
 187         for (unsigned int i=0; i<size4D[3]; i++)
 188         {
 189             vector+=iterators[i].Get();
 190             ++(iterators[i]);
 191         }
 192         vector/=size4D[3];
 193         avIt.Set(vector);
 194         ++avIt;
 195     }
 196     progress++;
 197     return output;
 198 }
 199
 200
 201 void vvMidPosition::update_progress()
 202 {
 203     p_bar.setValue(progress);
 204     p_bar.show();
 205 }
 206