1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
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>
20 #include <itkWarpImageFilter.h>
21 #include <itkVectorResampleImageFilter.h>
23 #include "vvMidPosition.h"
24 #include "clitkCommon.h"
25 #include "vvFromITK.h"
27 #include "clitkInvertVFFilter.h"
29 vvMidPosition::vvMidPosition() :
30 slicer_manager(0), error(false),
31 reference_image_index(0),
32 p_bar("Computing mid-position...","Cancel", 0,6),
37 void vvMidPosition::Update()
40 while (this->isRunning())
43 this->update_progress();
44 qApp->processEvents();
49 typedef itk::Vector<float,3> VFPixelType;
50 typedef itk::Image<VFPixelType,4> VFType;
51 typedef itk::Image<VFPixelType,3> OutputVFType;
55 ///This function averages a vector field along the temporal dimension
56 // the progress parameter is a reference to a counter used to update the progress dialog
57 itk::Image<itk::Vector<float,3>,3>::Pointer AverageField(itk::Image<itk::Vector<float,3>,4>::ConstPointer vf,int &progress);
58 ///Warps an image frame using the vf passed in parameter
59 template<class ImagePixelType> vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,int reference_image_index);
61 void vvMidPosition::run()
64 if(slicer_manager->GetImage()->GetNumberOfDimensions() != 4)
65 error_message="Computation of midposition is only supported for 4D image sequences.";
66 else if (!slicer_manager->GetVF())
67 error_message="A VF is required for midposition computation";
68 else if (slicer_manager->GetVF()->GetScalarTypeAsString() != "float")
69 error_message="Please use a vector field of type float.";
72 VFType::ConstPointer vf = vvImageToITK<VFType>(slicer_manager->GetVF());
73 OutputVFType::Pointer avg_vf=AverageField(vf,this->progress);
74 clitk::InvertVFFilter<OutputVFType,OutputVFType>::Pointer inv_filter=
75 clitk::InvertVFFilter<OutputVFType,OutputVFType>::New();
76 inv_filter->SetInput(avg_vf);
79 if (slicer_manager->GetImage()->GetScalarTypeAsString() == "short")
80 this->output=WarpRefImage<short>(inv_filter->GetOutput(),slicer_manager->GetImage(),reference_image_index);
83 error_message="Unsupported image pixel type.";
91 template<class ImagePixelType>
92 vvImage::Pointer WarpRefImage(OutputVFType::Pointer vf,vvImage::Pointer image,int reference_image_index)
94 typedef itk::Image<ImagePixelType,3> ImageType;
95 typedef itk::WarpImageFilter<ImageType,ImageType,itk::Image<VFPixelType,3> > FilterType;
97 typename ImageType::ConstPointer input = vvSingleFrameToITK<3,ImagePixelType>(image,reference_image_index);
99 //We resample the VF because itk's warp filter doesn't like it when the vf and the image have
101 typename itk::VectorResampleImageFilter<OutputVFType, OutputVFType >::Pointer
102 resampler =itk::VectorResampleImageFilter<OutputVFType, OutputVFType >::New();
103 resampler->SetInput(vf);
104 resampler->SetOutputSpacing(input->GetSpacing());
105 resampler->SetOutputOrigin(vf->GetOrigin());
106 //Calculate the new size so that it contains the vf
107 typename ImageType::SizeType newSize;
108 for (unsigned int i=0 ; i <3; i++)
109 newSize[i]=(unsigned int) (vf->GetLargestPossibleRegion().GetSize()[i]*vf->GetSpacing()[i]/input->GetSpacing()[i]);
110 resampler->SetSize( newSize);
112 typename FilterType::Pointer warp_filter = FilterType::New();
113 warp_filter->SetInput(input);
114 warp_filter->SetDeformationField(resampler->GetOutput());
115 warp_filter->SetOutputSpacing(input->GetSpacing());
116 warp_filter->SetOutputOrigin(input->GetOrigin());
117 warp_filter->SetOutputSize(input->GetLargestPossibleRegion().GetSize());
118 warp_filter->SetEdgePaddingValue(-1000);
119 warp_filter->Update();
120 return vvImageFromITK<3,ImagePixelType>(warp_filter->GetOutput());
123 itk::Image<itk::Vector<float,3>,3>::Pointer AverageField(itk::Image<itk::Vector<float,3>,4>::ConstPointer vf, int& progress)
127 VFType::RegionType region4D=vf->GetLargestPossibleRegion();
128 VFType::RegionType::SizeType size4D=region4D.GetSize();
129 VFType::IndexType index4D=region4D.GetIndex();
130 VFType::SpacingType spacing4D=vf->GetSpacing();
131 VFType::PointType origin4D=vf->GetOrigin();
133 OutputVFType::RegionType region;
134 OutputVFType::RegionType::SizeType size;
135 OutputVFType::IndexType index;
136 OutputVFType::SpacingType spacing;
137 OutputVFType::PointType origin;
138 for (unsigned int i=0; i< 3; i++)
142 spacing[i]=spacing4D[i];
143 origin[i]=origin4D[i];
145 region.SetSize(size);
146 region.SetIndex(index);
147 OutputVFType::Pointer output= OutputVFType::New();
148 output->SetRegions(region);
149 output->SetSpacing(spacing);
150 output->SetOrigin(origin);
156 typedef itk::ImageRegionConstIterator<VFType> IteratorType;
157 std::vector<IteratorType> iterators(size4D[3]);
158 for (unsigned int i=0; i< size4D[3]; i++)
160 VFType::RegionType regionIt=region4D;
161 VFType::RegionType::SizeType sizeIt=regionIt.GetSize();
163 regionIt.SetSize(sizeIt);
164 VFType::IndexType indexIt=regionIt.GetIndex();
166 regionIt.SetIndex(indexIt);
167 iterators[i]=IteratorType(vf, regionIt);
171 typedef itk::ImageRegionIterator<OutputVFType> OutputIteratorType;
172 OutputIteratorType avIt(output, output->GetLargestPossibleRegion());
176 VFPixelType zeroVector=itk::NumericTraits<VFPixelType>::Zero;
178 while (!(iterators[0]).IsAtEnd())
181 for (unsigned int i=0; i<size4D[3]; i++)
183 vector+=iterators[i].Get();
195 void vvMidPosition::update_progress()
197 p_bar.setValue(progress);