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