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 <QApplication>
26 #include <itkWarpImageFilter.h>
27 #include <itkVectorResampleImageFilter.h>
29 #include "vvMidPosition.h"
30 #include "clitkCommon.h"
31 #include "vvFromITK.h"
33 #include "clitkInvertVFFilter.h"
35 vvMidPosition::vvMidPosition() :
36 slicer_manager(0), error(false),
37 reference_image_index(0),
38 p_bar("Computing mid-position...","Cancel", 0,6),
43 void vvMidPosition::Update()
46 while (this->isRunning())
49 this->update_progress();
50 qApp->processEvents();
55 typedef itk::Vector<float,3> VFPixelType;
56 typedef itk::Image<VFPixelType,4> VFType;
57 typedef itk::Image<VFPixelType,3> OutputVFType;
61 ///This function averages a vector field along the temporal dimension
62 // the progress parameter is a reference to a counter used to update the progress dialog
63 itk::Image<itk::Vector<float,3>,3>::Pointer AverageField(itk::Image<itk::Vector<float,3>,4>::ConstPointer vf,int &progress);
64 ///Warps an image frame using the vf passed in parameter
65 template<class ImagePixelType> vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,int reference_image_index);
67 void vvMidPosition::run()
70 if(slicer_manager->GetImage()->GetNumberOfDimensions() != 4)
71 error_message="Computation of midposition is only supported for 4D image sequences.";
72 else if (!slicer_manager->GetVF())
73 error_message="A VF is required for midposition computation";
74 else if (slicer_manager->GetVF()->GetScalarTypeAsString() != "float")
75 error_message="Please use a vector field of type float.";
78 VFType::ConstPointer vf = vvImageToITK<VFType>(slicer_manager->GetVF());
79 OutputVFType::Pointer avg_vf=AverageField(vf,this->progress);
80 clitk::InvertVFFilter<OutputVFType,OutputVFType>::Pointer inv_filter=
81 clitk::InvertVFFilter<OutputVFType,OutputVFType>::New();
82 inv_filter->SetInput(avg_vf);
85 if (slicer_manager->GetImage()->GetScalarTypeAsString() == "short")
86 this->output=WarpRefImage<short>(inv_filter->GetOutput(),slicer_manager->GetImage(),reference_image_index);
89 error_message="Unsupported image pixel type.";
97 template<class ImagePixelType>
98 vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,int reference_image_index)
100 typedef itk::Image<ImagePixelType,3> ImageType;
101 typedef itk::WarpImageFilter<ImageType,ImageType,itk::Image<VFPixelType,3> > FilterType;
103 typename ImageType::ConstPointer input = vvSingleFrameToITK<3,ImagePixelType>(image,reference_image_index);
105 //We resample the VF because itk's warp filter doesn't like it when the vf and the image have
107 typename itk::VectorResampleImageFilter<OutputVFType, OutputVFType >::Pointer
108 resampler =itk::VectorResampleImageFilter<OutputVFType, OutputVFType >::New();
109 resampler->SetInput(vf);
110 resampler->SetOutputSpacing(input->GetSpacing());
111 resampler->SetOutputOrigin(vf->GetOrigin());
112 //Calculate the new size so that it contains the vf
113 typename ImageType::SizeType newSize;
114 for (unsigned int i=0 ; i <3; i++)
115 newSize[i]=(unsigned int) (vf->GetLargestPossibleRegion().GetSize()[i]*vf->GetSpacing()[i]/input->GetSpacing()[i]);
116 resampler->SetSize( newSize);
118 typename FilterType::Pointer warp_filter = FilterType::New();
119 warp_filter->SetInput(input);
120 warp_filter->SetDeformationField(resampler->GetOutput());
121 warp_filter->SetOutputSpacing(input->GetSpacing());
122 warp_filter->SetOutputOrigin(input->GetOrigin());
123 warp_filter->SetOutputSize(input->GetLargestPossibleRegion().GetSize());
124 warp_filter->SetEdgePaddingValue(-1000);
125 warp_filter->Update();
126 return vvImageFromITK<3,ImagePixelType>(warp_filter->GetOutput());
129 itk::Image<itk::Vector<float,3>,3>::Pointer AverageField(itk::Image<itk::Vector<float,3>,4>::ConstPointer vf, int& progress)
133 VFType::RegionType region4D=vf->GetLargestPossibleRegion();
134 VFType::RegionType::SizeType size4D=region4D.GetSize();
135 VFType::IndexType index4D=region4D.GetIndex();
136 VFType::SpacingType spacing4D=vf->GetSpacing();
137 VFType::PointType origin4D=vf->GetOrigin();
139 OutputVFType::RegionType region;
140 OutputVFType::RegionType::SizeType size;
141 OutputVFType::IndexType index;
142 OutputVFType::SpacingType spacing;
143 OutputVFType::PointType origin;
144 for (unsigned int i=0; i< 3; i++)
148 spacing[i]=spacing4D[i];
149 origin[i]=origin4D[i];
151 region.SetSize(size);
152 region.SetIndex(index);
153 OutputVFType::Pointer output= OutputVFType::New();
154 output->SetRegions(region);
155 output->SetSpacing(spacing);
156 output->SetOrigin(origin);
162 typedef itk::ImageRegionConstIterator<VFType> IteratorType;
163 std::vector<IteratorType> iterators(size4D[3]);
164 for (unsigned int i=0; i< size4D[3]; i++)
166 VFType::RegionType regionIt=region4D;
167 VFType::RegionType::SizeType sizeIt=regionIt.GetSize();
169 regionIt.SetSize(sizeIt);
170 VFType::IndexType indexIt=regionIt.GetIndex();
172 regionIt.SetIndex(indexIt);
173 iterators[i]=IteratorType(vf, regionIt);
177 typedef itk::ImageRegionIterator<OutputVFType> OutputIteratorType;
178 OutputIteratorType avIt(output, output->GetLargestPossibleRegion());
182 VFPixelType zeroVector=itk::NumericTraits<VFPixelType>::Zero;
184 while (!(iterators[0]).IsAtEnd())
187 for (unsigned int i=0; i<size4D[3]; i++)
189 vector+=iterators[i].Get();
201 void vvMidPosition::update_progress()
203 p_bar.setValue(progress);