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 #define TRY_TYPE(TYPE)                                                  \
  45   if (clitk::IsSameType<TYPE>(mInputImage->GetScalarTypeAsString())) { this->Update_WithDimAndPixelType<Dim, TYPE>(); return; }
  46 //  TRY_TYPE(signed char);
  47 //  TRY_TYPE(uchar);
  48     TRY_TYPE(short);
  49 //  TRY_TYPE(ushort);
  50 //  TRY_TYPE(int); // no uint ...
  51 //  TRY_TYPE(float);
  52     //TRY_TYPE(double);
  53 #undef TRY_TYPE
  54
  55     std::string list = clitk::CreateListOfTypes<char, clitk::uchar, short, ushort, int, float, double>();
  56     std::cerr << "Error, I don't know the type '" << mInputImage->GetScalarTypeAsString() << "' for the input image. "
  57               << std::endl << "Known types are " << list << std::endl;
  58     exit(0);
  59 }
  60
  61 //====================================================================
  62 template<unsigned int Dim, class PixelType>
  63 void vvImageWarp::Update_WithDimAndPixelType() {
  64     QProgressDialog progress(parent_window);
  65     progress.setCancelButton(0);
  66     progress.setLabelText("Computing warped and ventilation images...");
  67     progress.show();
  68     typedef itk::Image<PixelType,Dim> ImageType;
  69     typedef itk::Image<float,Dim> JacobianImageType;
  70     typedef itk::Image<itk::Vector<float,Dim>,Dim> VectorImageType;
  71     typedef std::vector<typename ImageType::ConstPointer> ImageSeriesType;
  72     ImageSeriesType input = vvImageToITKImageVector<Dim,PixelType>(mInputImage);
  73     ImageSeriesType output;
  74
  75     typename itk::VTKImageToImageFilter<VectorImageType>::Pointer vf_connector=
  76         itk::VTKImageToImageFilter<VectorImageType>::New();
  77
  78     //Warp, then join, then convert to vv
  79     typedef itk::WarpImageFilter<ImageType,ImageType,VectorImageType> WarpFilterType;
  80     typedef itk::DisplacementFieldJacobianDeterminantFilter<VectorImageType> JacobianFilterType;
  81     vvImage::Pointer result=vvImage::New();
  82     typedef itk::JoinSeriesImageFilter< ImageType,itk::Image<PixelType,Dim+1> > JoinFilterType;
  83     typedef itk::JoinSeriesImageFilter< JacobianImageType,itk::Image<float,Dim+1> > JacobianJoinFilterType;
  84     typename JoinFilterType::Pointer join=JoinFilterType::New();
  85     typename JoinFilterType::Pointer diff_join=JoinFilterType::New();
  86     typename JacobianJoinFilterType::Pointer jacobian_join=JacobianJoinFilterType::New();
  87     join->SetSpacing(1);
  88     join->SetOrigin(0); //Set the temporal origin
  89     diff_join->SetSpacing(1);
  90     diff_join->SetOrigin(0);
  91     jacobian_join->SetSpacing(1);
  92     jacobian_join->SetOrigin(0);
  93     typedef itk::SubtractImageFilter<ImageType,ImageType,ImageType> DiffFilter;
  94     std::vector<typename ImageType::Pointer> warped_images;
  95
  96     for (unsigned int num = 0; num < input.size(); num++)
  97     {
  98         typename WarpFilterType::Pointer warp_filter=WarpFilterType::New();
  99         typename JacobianFilterType::Pointer jacobian_filter=JacobianFilterType::New();
 100         jacobian_filter->SetUseImageSpacingOn();
 101         vf_connector->SetInput(mVF->GetVTKImages()[num]);
 102         warp_filter->SetInput(input[num]);
 103         warp_filter->SetDeformationField(vf_connector->GetOutput());
 104         jacobian_filter->SetInput(vf_connector->GetOutput());
 105         warp_filter->SetOutputSpacing(input[num]->GetSpacing());
 106         warp_filter->SetOutputOrigin(input[num]->GetOrigin());
 107         warp_filter->SetEdgePaddingValue(-1000);
 108         warp_filter->Update();
 109         jacobian_filter->Update();
 110         warped_images.push_back(warp_filter->GetOutput());
 111         jacobian_join->PushBackInput(jacobian_filter->GetOutput());
 112         join->PushBackInput(warp_filter->GetOutput());
 113         progress.setValue(progress.value()+1);
 114     }
 115     for (typename std::vector<typename ImageType::Pointer>::const_iterator i = warped_images.begin(); i!=warped_images.end();i++)
 116     {
 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     {
 139         if (mInputImage->GetSpacing()[i] != mVF->GetSpacing()[i])
 140             return false;
 141     }
 142     switch (mInputImage->GetNumberOfDimensions())
 143     {
 144 //      case 2: this->Update_WithDim<2>(); break;;
 145 //      case 3: this->Update_WithDim<3>(); break;;
 146     case 4:
 147         this->Update_WithDim<3>();
 148         break;;
 149     default:
 150         DD("Error: dimension not handled.");
 151     }
 152     return true;
 153 }