vv/vvImageWarp.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 <QProgressDialog>
  19 #include <QWidget>
  20 #include <itkImage.h>
  21 #include <itkWarpImageFilter.h>
  22 #include <itkJoinSeriesImageFilter.h>
  23 #include <itkImageFileReader.h>
  24 #include <itkSubtractImageFilter.h>
  25 #include <itkDisplacementFieldJacobianDeterminantFilter.h>
  26 #include "vvToITK.h"
  27 #include "vvFromITK.h"
  28
  29 #include "vvImageWarp.h"
  30 #include "clitkCommon.h"
  31
  32 //====================================================================
  33 vvImageWarp::vvImageWarp(vvImage::Pointer input,vvImage::Pointer vf,unsigned int ref_image,QWidget* parent):
  34   mRefImage(ref_image),
  35   parent_window(parent),
  36   mInputImage(input),
  37   mVF(vf)
  38 {}
  39
  40
  41 //====================================================================
  42 template<unsigned int Dim>
  43 void vvImageWarp::Update_WithDim()
  44 {
  45 #define TRY_TYPE(TYPE)                                                  \
  46   if (clitk::IsSameType<TYPE>(mInputImage->GetScalarTypeAsITKString())) { this->Update_WithDimAndPixelType<Dim, TYPE>(); return; }
  47 //  TRY_TYPE(signed char);
  48 //  TRY_TYPE(uchar);
  49   TRY_TYPE(short);
  50 //  TRY_TYPE(ushort);
  51 //  TRY_TYPE(int); // no uint ...
  52 //  TRY_TYPE(float);
  53   //TRY_TYPE(double);
  54 #undef TRY_TYPE
  55
  56   std::string list = clitk::CreateListOfTypes<char, clitk::uchar, short, ushort, int, float, double>();
  57   std::cerr << "Error, I don't know the type '" << mInputImage->GetScalarTypeAsITKString() << "' for the input image. "
  58             << std::endl << "Known types are " << list << std::endl;
  59   exit(0);
  60 }
  61
  62 //====================================================================
  63 template<unsigned int Dim, class PixelType>
  64 void vvImageWarp::Update_WithDimAndPixelType()
  65 {
  66   QProgressDialog progress(parent_window);
  67   progress.setCancelButton(0);
  68   progress.setLabelText("Computing warped and ventilation images...");
  69   progress.show();
  70   typedef itk::Image<PixelType,Dim> ImageType;
  71   typedef itk::Image<float,Dim> JacobianImageType;
  72   typedef itk::Image<itk::Vector<float,Dim>,Dim> VectorImageType;
  73   typedef std::vector<typename ImageType::ConstPointer> ImageSeriesType;
  74   ImageSeriesType input = vvImageToITKImageVector<Dim,PixelType>(mInputImage);
  75   ImageSeriesType output;
  76
  77   typename itk::VTKImageToImageFilter<VectorImageType>::Pointer vf_connector=
  78     itk::VTKImageToImageFilter<VectorImageType>::New();
  79
  80   //Warp, then join, then convert to vv
  81   typedef itk::WarpImageFilter<ImageType,ImageType,VectorImageType> WarpFilterType;
  82   typedef itk::DisplacementFieldJacobianDeterminantFilter<VectorImageType> JacobianFilterType;
  83   vvImage::Pointer result=vvImage::New();
  84   typedef itk::JoinSeriesImageFilter< ImageType,itk::Image<PixelType,Dim+1> > JoinFilterType;
  85   typedef itk::JoinSeriesImageFilter< JacobianImageType,itk::Image<float,Dim+1> > JacobianJoinFilterType;
  86   typename JoinFilterType::Pointer join=JoinFilterType::New();
  87   typename JoinFilterType::Pointer diff_join=JoinFilterType::New();
  88   typename JacobianJoinFilterType::Pointer jacobian_join=JacobianJoinFilterType::New();
  89   join->SetSpacing(1);
  90   join->SetOrigin(0); //Set the temporal origin
  91   diff_join->SetSpacing(1);
  92   diff_join->SetOrigin(0);
  93   jacobian_join->SetSpacing(1);
  94   jacobian_join->SetOrigin(0);
  95   typedef itk::SubtractImageFilter<ImageType,ImageType,ImageType> DiffFilter;
  96   std::vector<typename ImageType::Pointer> warped_images;
  97
  98   for (unsigned int num = 0; num < input.size(); num++) {
  99     typename WarpFilterType::Pointer warp_filter=WarpFilterType::New();
 100     typename JacobianFilterType::Pointer jacobian_filter=JacobianFilterType::New();
 101     jacobian_filter->SetUseImageSpacingOn();
 102     vf_connector->SetInput(mVF->GetVTKImages()[num]);
 103     warp_filter->SetInput(input[num]);
 104     warp_filter->SetDeformationField(vf_connector->GetOutput());
 105     jacobian_filter->SetInput(vf_connector->GetOutput());
 106     warp_filter->SetOutputSpacing(input[num]->GetSpacing());
 107     warp_filter->SetOutputOrigin(input[num]->GetOrigin());
 108     warp_filter->SetEdgePaddingValue(-1000);
 109     warp_filter->Update();
 110     jacobian_filter->Update();
 111     warped_images.push_back(warp_filter->GetOutput());
 112     jacobian_join->PushBackInput(jacobian_filter->GetOutput());
 113     join->PushBackInput(warp_filter->GetOutput());
 114     progress.setValue(progress.value()+1);
 115   }
 116   for (typename std::vector<typename ImageType::Pointer>::const_iterator i = warped_images.begin(); i!=warped_images.end(); i++) {
 117     typename DiffFilter::Pointer diff_filter = DiffFilter::New();
 118     diff_filter->SetInput2(*i);
 119     diff_filter->SetInput1(*(warped_images.begin()+mRefImage));
 120     diff_filter->Update();
 121     diff_join->PushBackInput(diff_filter->GetOutput());
 122     progress.setValue(progress.value()+1);
 123   }
 124   join->Update();
 125   diff_join->Update();
 126   jacobian_join->Update();
 127   mWarpedImage = vvImageFromITK<Dim+1,PixelType>(join->GetOutput());
 128   mDiffImage = vvImageFromITK<Dim+1,PixelType>(diff_join->GetOutput());
 129   mJacobianImage = vvImageFromITK<Dim+1,float>(jacobian_join->GetOutput());
 130   //mJacobianImage = vvImageFromITK<Dim,float>(temporal_filter->GetOutput());
 131 }
 132 //====================================================================
 133 //
 134
 135 bool vvImageWarp::ComputeWarpedImage()
 136 {
 137   for (int i=0; i<mInputImage->GetNumberOfDimensions(); i++) {
 138     if (mInputImage->GetSpacing()[i] != mVF->GetSpacing()[i])
 139       return false;
 140   }
 141   switch (mInputImage->GetNumberOfDimensions()) {
 142 //      case 2: this->Update_WithDim<2>(); break;;
 143 //      case 3: this->Update_WithDim<3>(); break;;
 144   case 4:
 145     this->Update_WithDim<3>();
 146     break;;
 147   default:
 148     DD("Error: dimension not handled.");
 149   }
 150   return true;
 151 }