]> Creatis software - clitk.git/blob - vv/vvMidPosition.cxx
Fusion windows level is now 4 decimals
[clitk.git] / vv / vvMidPosition.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 <QApplication>
19
20 #include <itkWarpImageFilter.h>
21 #include <itkVectorResampleImageFilter.h>
22
23 #include "vvMidPosition.h"
24 #include "clitkCommon.h"
25 #include "vvFromITK.h"
26 #include "vvToITK.h"
27 #include "clitkInvertVFFilter.h"
28
29 vvMidPosition::vvMidPosition() :
30   slicer_manager(0), error(false),
31   reference_image_index(0),
32   p_bar("Computing mid-position...","Cancel", 0,6),
33   progress(0)
34 {
35 }
36
37 void vvMidPosition::Update()
38 {
39   this->start();
40   while (this->isRunning()) {
41     this->wait(50);
42     this->update_progress();
43     qApp->processEvents();
44   }
45 }
46
47 //Common typedefs
48 typedef itk::Vector<float,3> VFPixelType;
49 typedef itk::Image<VFPixelType,4> VFType;
50 typedef itk::Image<VFPixelType,3> OutputVFType;
51
52 ///Internal functions
53
54 ///This function averages a vector field along the temporal dimension
55 // the progress parameter is a reference to a counter used to update the progress dialog
56 itk::Image<itk::Vector<float,3>,3>::Pointer AverageField(itk::Image<itk::Vector<float,3>,4>::ConstPointer vf,int &progress);
57 ///Warps an image frame using the vf passed in parameter
58 template<class ImagePixelType> vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,int reference_image_index);
59
60 void vvMidPosition::run()
61 {
62   error=true;
63   if(slicer_manager->GetImage()->GetNumberOfDimensions() != 4)
64     error_message="Computation of midposition is only supported for 4D image sequences.";
65   else if (!slicer_manager->GetVF())
66     error_message="A VF is required for midposition computation";
67   else if (slicer_manager->GetVF()->GetScalarTypeAsITKString() != "float")
68     error_message="Please use a vector field of type float.";
69   else {
70     VFType::ConstPointer vf = vvImageToITK<VFType>(slicer_manager->GetVF());
71     OutputVFType::Pointer avg_vf=AverageField(vf,this->progress);
72     clitk::InvertVFFilter<OutputVFType,OutputVFType>::Pointer inv_filter=
73       clitk::InvertVFFilter<OutputVFType,OutputVFType>::New();
74     inv_filter->SetInput(avg_vf);
75     inv_filter->Update();
76     progress++;
77     if (slicer_manager->GetImage()->GetScalarTypeAsITKString() == "short")
78       this->output=WarpRefImage<short>(inv_filter->GetOutput(),slicer_manager->GetImage(),reference_image_index);
79     else {
80       error_message="Unsupported image pixel type.";
81       return;
82     }
83     progress++;
84     error=false;
85   }
86 }
87
88 template<class ImagePixelType>
89 vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,int reference_image_index)
90 {
91   typedef itk::Image<ImagePixelType,3> ImageType;
92   typedef itk::WarpImageFilter<ImageType,ImageType,itk::Image<VFPixelType,3> > FilterType;
93
94   typename ImageType::ConstPointer input = vvSingleFrameToITK<3,ImagePixelType>(image,reference_image_index);
95
96   //We resample the VF because itk's warp filter doesn't like it when the vf and the image have
97   //different spacings
98   typename itk::VectorResampleImageFilter<OutputVFType, OutputVFType >::Pointer
99   resampler =itk::VectorResampleImageFilter<OutputVFType, OutputVFType >::New();
100   resampler->SetInput(vf);
101   resampler->SetOutputSpacing(input->GetSpacing());
102   resampler->SetOutputOrigin(vf->GetOrigin());
103   //Calculate the new size so that it contains the vf
104   typename ImageType::SizeType newSize;
105   for (unsigned int i=0 ; i <3; i++)
106     newSize[i]=(unsigned int) (vf->GetLargestPossibleRegion().GetSize()[i]*vf->GetSpacing()[i]/input->GetSpacing()[i]);
107   resampler->SetSize( newSize);
108
109   typename FilterType::Pointer warp_filter = FilterType::New();
110   warp_filter->SetInput(input);
111   warp_filter->SetDeformationField(resampler->GetOutput());
112   warp_filter->SetOutputSpacing(input->GetSpacing());
113   warp_filter->SetOutputOrigin(input->GetOrigin());
114   warp_filter->SetOutputSize(input->GetLargestPossibleRegion().GetSize());
115   warp_filter->SetEdgePaddingValue(-1000);
116   warp_filter->Update();
117   return vvImageFromITK<3,ImagePixelType>(warp_filter->GetOutput());
118 }
119
120 itk::Image<itk::Vector<float,3>,3>::Pointer AverageField(itk::Image<itk::Vector<float,3>,4>::ConstPointer vf, int& progress)
121 {
122   progress++;
123
124   VFType::RegionType region4D=vf->GetLargestPossibleRegion();
125   VFType::RegionType::SizeType size4D=region4D.GetSize();
126   VFType::IndexType index4D=region4D.GetIndex();
127   VFType::SpacingType spacing4D=vf->GetSpacing();
128   VFType::PointType origin4D=vf->GetOrigin();
129
130   OutputVFType::RegionType region;
131   OutputVFType::RegionType::SizeType size;
132   OutputVFType::IndexType index;
133   OutputVFType::SpacingType spacing;
134   OutputVFType::PointType origin;
135   for (unsigned int i=0; i< 3; i++) {
136     size[i]=size4D[i];
137     index[i]=index4D[i];
138     spacing[i]=spacing4D[i];
139     origin[i]=origin4D[i];
140   }
141   region.SetSize(size);
142   region.SetIndex(index);
143   OutputVFType::Pointer output= OutputVFType::New();
144   output->SetRegions(region);
145   output->SetSpacing(spacing);
146   output->SetOrigin(origin);
147   output->Allocate();
148   progress++;
149
150
151   // Region iterators
152   typedef itk::ImageRegionConstIterator<VFType> IteratorType;
153   std::vector<IteratorType> iterators(size4D[3]);
154   for (unsigned int i=0; i< size4D[3]; i++) {
155     VFType::RegionType regionIt=region4D;
156     VFType::RegionType::SizeType sizeIt=regionIt.GetSize();
157     sizeIt[3]=1;
158     regionIt.SetSize(sizeIt);
159     VFType::IndexType indexIt=regionIt.GetIndex();
160     indexIt[3]=i;
161     regionIt.SetIndex(indexIt);
162     iterators[i]=IteratorType(vf, regionIt);
163   }
164   progress++;
165
166   typedef itk::ImageRegionIterator<OutputVFType> OutputIteratorType;
167   OutputIteratorType avIt(output, output->GetLargestPossibleRegion());
168
169   // Average
170   VFPixelType vector;
171   VFPixelType zeroVector=itk::NumericTraits<VFPixelType>::Zero;
172
173   while (!(iterators[0]).IsAtEnd()) {
174     vector=zeroVector;
175     for (unsigned int i=0; i<size4D[3]; i++) {
176       vector+=iterators[i].Get();
177       ++(iterators[i]);
178     }
179     vector/=size4D[3];
180     avIt.Set(vector);
181     ++avIt;
182   }
183   progress++;
184   return output;
185 }
186
187
188 void vvMidPosition::update_progress()
189 {
190   p_bar.setValue(progress);
191   p_bar.show();
192 }
193