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://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 <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 #if ITK_VERSION_MAJOR >= 4
 105     warp_filter->SetDisplacementField(vf_connector->GetOutput());
 106 #else
 107     warp_filter->SetDeformationField(vf_connector->GetOutput());
 108 #endif
 109     jacobian_filter->SetInput(vf_connector->GetOutput());
 110     warp_filter->SetOutputSpacing(input[num]->GetSpacing());
 111     warp_filter->SetOutputOrigin(input[num]->GetOrigin());
 112     warp_filter->SetEdgePaddingValue(-1000);
 113     warp_filter->Update();
 114     jacobian_filter->Update();
 115     warped_images.push_back(warp_filter->GetOutput());
 116     jacobian_join->PushBackInput(jacobian_filter->GetOutput());
 117     join->PushBackInput(warp_filter->GetOutput());
 118     progress.setValue(progress.value()+1);
 119   }
 120   for (typename std::vector<typename ImageType::Pointer>::const_iterator i = warped_images.begin(); i!=warped_images.end(); i++) {
 121     typename DiffFilter::Pointer diff_filter = DiffFilter::New();
 122     diff_filter->SetInput2(*i);
 123     diff_filter->SetInput1(*(warped_images.begin()+mRefImage));
 124     diff_filter->Update();
 125     diff_join->PushBackInput(diff_filter->GetOutput());
 126     progress.setValue(progress.value()+1);
 127   }
 128   join->Update();
 129   diff_join->Update();
 130   jacobian_join->Update();
 131   mWarpedImage = vvImageFromITK<Dim+1,PixelType>(join->GetOutput());
 132   mDiffImage = vvImageFromITK<Dim+1,PixelType>(diff_join->GetOutput());
 133   mJacobianImage = vvImageFromITK<Dim+1,float>(jacobian_join->GetOutput());
 134   //mJacobianImage = vvImageFromITK<Dim,float>(temporal_filter->GetOutput());
 135 }
 136 //====================================================================
 137 //
 138
 139 bool vvImageWarp::ComputeWarpedImage()
 140 {
 141   for (int i=0; i<mInputImage->GetNumberOfDimensions(); i++) {
 142     if (mInputImage->GetSpacing()[i] != mVF->GetSpacing()[i])
 143       return false;
 144   }
 145   switch (mInputImage->GetNumberOfDimensions()) {
 146 //      case 2: this->Update_WithDim<2>(); break;;
 147 //      case 3: this->Update_WithDim<3>(); break;;
 148   case 4:
 149     this->Update_WithDim<3>();
 150     break;;
 151   default:
 152     DD("Error: dimension not handled.");
 153   }
 154   return true;
 155 }