]> Creatis software - clitk.git/blob - vv/vvImageWarp.cxx
Add clitkImage2DicomSeries tool to write a Dicom Series from an image without a corre...
[clitk.git] / 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     warp_filter->SetDisplacementField(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 }