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://www.centreleonberard.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 #ifndef clitkMedianTemporalDimensionGenericFilter_txx
19 #define clitkMedianTemporalDimensionGenericFilter_txx
21 /* =================================================
22 * @file clitkMedianTemporalDimensionGenericFilter.txx
28 ===================================================*/
35 //-----------------------------------------------------------
37 //-----------------------------------------------------------
38 template<class args_info_type>
39 MedianTemporalDimensionGenericFilter<args_info_type>::MedianTemporalDimensionGenericFilter()
46 //-----------------------------------------------------------
48 //-----------------------------------------------------------
49 template<class args_info_type>
50 void MedianTemporalDimensionGenericFilter<args_info_type>::Update()
52 // Read the Dimension and PixelType
53 int Dimension, Components;
54 std::string PixelType;
55 ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components);
59 if (m_ArgsInfo.input_given>1) Dimension+=1;
60 if(Dimension==2) UpdateWithDim<2>(PixelType, Components);
61 else if(Dimension==3)UpdateWithDim<3>(PixelType, Components);
62 else if (Dimension==4)UpdateWithDim<4>(PixelType, Components);
64 std::cout<<"Error, Only for 2, 3 or 4D!!!"<<std::endl ;
69 //-------------------------------------------------------------------
70 // Update with the number of dimensions
71 //-------------------------------------------------------------------
72 template<class args_info_type>
73 template<unsigned int Dimension>
75 MedianTemporalDimensionGenericFilter<args_info_type>::UpdateWithDim(const std::string PixelType, const int Components)
77 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D with "<<Components<<" component(s) of "<< PixelType<<"..."<<std::endl;
80 if(PixelType == "short") {
81 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
82 UpdateWithDimAndPixelType<Dimension, signed short>();
84 else if (PixelType == "unsigned_char") {
85 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
86 UpdateWithDimAndPixelType<Dimension, unsigned char>();
88 else if (PixelType == "float") {
89 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
90 UpdateWithDimAndPixelType<Dimension, float>();
92 else if (PixelType == "double") {
93 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
94 UpdateWithDimAndPixelType<Dimension, double>();
96 else std::cerr<<"Number of components is "<<Components<<", not supported!"<<std::endl;
101 //-------------------------------------------------------------------
102 // Update with the number of dimensions and the pixeltype
103 //-------------------------------------------------------------------
104 template<class args_info_type>
105 template <unsigned int Dimension, class PixelType>
107 MedianTemporalDimensionGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
111 typedef itk::Image<PixelType, Dimension> InputImageType;
112 typedef itk::Image<PixelType, Dimension-1> OutputImageType;
113 typename InputImageType::Pointer input;
115 if (m_ArgsInfo.input_given ==1 ) {
117 typedef itk::ImageFileReader<InputImageType> InputReaderType;
118 typename InputReaderType::Pointer reader = InputReaderType::New();
119 reader->SetFileName( m_InputFileName);
121 input= reader->GetOutput();
125 // Read and join multiple inputs
126 if (m_Verbose) std::cout<<m_ArgsInfo.input_given<<" inputs given..."<<std::endl;
127 std::vector<std::string> filenames;
128 for(unsigned int i=0; i<m_ArgsInfo.input_given; i++) {
129 if (m_Verbose) std::cout<<m_ArgsInfo.input_arg[i]<<std::endl;
130 filenames.push_back(m_ArgsInfo.input_arg[i]);
133 typedef itk::ImageSeriesReader<InputImageType> ImageReaderType;
134 typename ImageReaderType::Pointer reader= ImageReaderType::New();
135 reader->SetFileNames(filenames);
137 input =reader->GetOutput();
142 typename OutputImageType::RegionType region;
143 typename OutputImageType::RegionType::SizeType size;
144 typename OutputImageType::IndexType index;
145 typename OutputImageType::SpacingType spacing;
146 typename OutputImageType::PointType origin;
147 typename InputImageType::RegionType region4D=input->GetLargestPossibleRegion();
148 typename InputImageType::RegionType::SizeType size4D=region4D.GetSize();
149 typename InputImageType::IndexType index4D=region4D.GetIndex();
150 typename InputImageType::SpacingType spacing4D=input->GetSpacing();
151 typename InputImageType::PointType origin4D=input->GetOrigin();
153 for (unsigned int i=0; i< Dimension-1; i++) {
156 spacing[i]=spacing4D[i];
157 origin[i]=origin4D[i];
159 region.SetSize(size);
160 region.SetIndex(index);
161 typename OutputImageType::Pointer output= OutputImageType::New();
162 output->SetRegions(region);
163 output->SetSpacing(spacing);
164 output->SetOrigin(origin);
169 typedef itk::ImageRegionIterator<InputImageType> IteratorType;
170 std::vector<IteratorType> iterators(size4D[Dimension-1]);
171 for (unsigned int i=0; i< size4D[Dimension-1]; i++) {
172 typename InputImageType::RegionType regionIt=region4D;
173 typename InputImageType::RegionType::SizeType sizeIt=regionIt.GetSize();
174 sizeIt[Dimension-1]=1;
175 regionIt.SetSize(sizeIt);
176 typename InputImageType::IndexType indexIt=regionIt.GetIndex();
177 indexIt[Dimension-1]=i;
178 regionIt.SetIndex(indexIt);
179 iterators[i]=IteratorType(input, regionIt);
182 typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
183 OutputIteratorType avIt(output, output->GetLargestPossibleRegion());
187 typename InputImageType::IndexValueType temporal_dimension = size4D[Dimension-1];
188 std::vector<PixelType> temp(temporal_dimension);
189 while (!(iterators[0]).IsAtEnd()) {
191 for (unsigned int i=0; i<temporal_dimension; i++) {
192 temp[i] = iterators[i].Get();
195 if (temporal_dimension & 1) {
196 nth_element(temp.begin(), temp.begin() + temporal_dimension/2,temp.end());
197 value = temp[temporal_dimension/2];
199 nth_element(temp.begin(), temp.begin() + temporal_dimension/2 - 1, temp.end());
200 value = temp[temporal_dimension/2 - 1];
201 nth_element(temp.begin(), temp.begin() + temporal_dimension/2, temp.end());
202 value += temp[temporal_dimension/2];
210 typedef itk::ImageFileWriter<OutputImageType> WriterType;
211 typename WriterType::Pointer writer = WriterType::New();
212 writer->SetFileName(m_ArgsInfo.output_arg);
213 writer->SetInput(output);
221 #endif //#define clitkMedianTemporalDimensionGenericFilter_txx