1 /*=========================================================================
5 Author : Joel Schaerer (joel.schaerer@insa-lyon.fr)
8 Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
9 CREATIS-LRMN http://www.creatis.insa-lyon.fr
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.
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.
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/>.
23 =========================================================================*/
24 #include <QProgressDialog>
27 #include <itkWarpImageFilter.h>
28 #include <itkJoinSeriesImageFilter.h>
29 #include <itkImageFileReader.h>
30 #include <itkSubtractImageFilter.h>
31 #include <itkDisplacementFieldJacobianDeterminantFilter.h>
33 #include "vvFromITK.h"
35 #include "vvImageWarp.h"
36 #include "clitkCommon.h"
38 //====================================================================
39 vvImageWarp::vvImageWarp(vvImage::Pointer input,vvImage::Pointer vf,unsigned int ref_image,QWidget* parent):
41 parent_window(parent),
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);
56 // TRY_TYPE(int); // no uint ...
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;
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...");
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;
81 typename itk::VTKImageToImageFilter<VectorImageType>::Pointer vf_connector=
82 itk::VTKImageToImageFilter<VectorImageType>::New();
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();
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;
102 for (unsigned int num = 0; num < input.size(); num++)
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);
121 for (typename std::vector<typename ImageType::Pointer>::const_iterator i = warped_images.begin(); i!=warped_images.end();i++)
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);
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());
138 //====================================================================
141 bool vvImageWarp::ComputeWarpedImage()
143 for (int i=0;i<mInputImage->GetNumberOfDimensions();i++)
145 if (mInputImage->GetSpacing()[i] != mVF->GetSpacing()[i])
148 switch (mInputImage->GetNumberOfDimensions())
150 // case 2: this->Update_WithDim<2>(); break;;
151 // case 3: this->Update_WithDim<3>(); break;;
153 this->Update_WithDim<3>();
156 DD("Error: dimension not handled.");