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 #ifndef CLITKIMAGEARITHMGENERICFILTER_TXX
19 #define CLITKIMAGEARITHMGENERICFILTER_TXX
23 //--------------------------------------------------------------------
24 template<class args_info_type>
25 ImageArithmGenericFilter<args_info_type>::ImageArithmGenericFilter()
26 :ImageToImageGenericFilter<Self>("ImageArithmGenericFilter"),mTypeOfOperation(0) {
27 InitializeImageType<2>();
28 InitializeImageType<3>();
29 InitializeImageType<4>();
30 mIsOperationUseASecondImage = false;
31 mOverwriteInputImage = true;
33 //--------------------------------------------------------------------
36 //--------------------------------------------------------------------
37 template<class args_info_type>
38 template<unsigned int Dim>
39 void ImageArithmGenericFilter<args_info_type>::InitializeImageType() {
40 ADD_IMAGE_TYPE(Dim, char);
41 ADD_IMAGE_TYPE(Dim, uchar);
42 ADD_IMAGE_TYPE(Dim, short);
43 ADD_IMAGE_TYPE(Dim, ushort);
44 ADD_IMAGE_TYPE(Dim, int);
45 ADD_IMAGE_TYPE(Dim, float);
46 ADD_IMAGE_TYPE(Dim, double);
48 //--------------------------------------------------------------------
51 //--------------------------------------------------------------------
52 template<class args_info_type>
53 void ImageArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b) {
54 mOverwriteInputImage = b;
56 //--------------------------------------------------------------------
59 //--------------------------------------------------------------------
60 template<class args_info_type>
61 void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a) {
65 SetIOVerbose(mArgsInfo.verbose_flag);
66 mTypeOfOperation = mArgsInfo.operation_arg;
67 mDefaultPixelValue = mArgsInfo.pixelValue_arg;
68 mScalar = mArgsInfo.scalar_arg;
69 mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
71 if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
73 if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
74 if (mArgsInfo.input2_given) {
75 mIsOperationUseASecondImage = true;
76 AddInputFilename(mArgsInfo.input2_arg);
79 if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
81 // Check type of operation (with scalar or with other image)
82 if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
83 std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
86 if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
87 if (mArgsInfo.operation_arg < 5) {
88 std::cerr << "Such operation need the --scalar option." << std::endl;
93 //--------------------------------------------------------------------
96 //--------------------------------------------------------------------
97 template<class args_info_type>
98 template<class ImageType>
99 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType() {
101 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
103 // Set input image iterator
104 typedef itk::ImageRegionIterator<ImageType> IteratorType;
105 IteratorType it(input1, input1->GetLargestPossibleRegion());
108 typename ImageType::Pointer input2 = this->template GetInput<ImageType>(1);
111 if (mIsOperationUseASecondImage) {
113 input2 = this->template GetInput<ImageType>(1);
114 // Set input image iterator
115 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
118 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
119 if (mOverwriteInputImage && mOutputIsFloat && (typeid(ImageType::PixelType) != typeid(float))) {
120 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
121 // << typeid(PixelType).name()
123 mOverwriteInputImage = false;
126 // ---------------- Overwrite input Image ---------------------
127 if (mOverwriteInputImage) {
128 // Set output iterator (to input1)
129 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
130 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
131 else ComputeImage(it, ito);
132 this->template SetNextOutput<ImageType>(input1);
135 // ---------------- Create new output Image ---------------------
137 if (mOutputIsFloat) {
138 // Create output image
139 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
140 typename OutputImageType::Pointer output = OutputImageType::New();
141 output->SetRegions(input1->GetLargestPossibleRegion());
142 output->SetOrigin(input1->GetOrigin());
143 output->SetSpacing(input1->GetSpacing());
145 // Set output iterator
146 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
147 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
148 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
149 else ComputeImage(it, ito);
150 this->template SetNextOutput<OutputImageType>(output);
153 // Create output image
154 typedef ImageType OutputImageType;
155 typename OutputImageType::Pointer output = OutputImageType::New();
156 output->SetRegions(input1->GetLargestPossibleRegion());
157 output->SetOrigin(input1->GetOrigin());
158 output->SetSpacing(input1->GetSpacing());
160 // Set output iterator
161 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
162 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
163 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
164 else ComputeImage(it, ito);
165 this->template SetNextOutput<OutputImageType>(output);
169 //--------------------------------------------------------------------
171 //--------------------------------------------------------------------
172 template<class args_info_type>
173 template<class Iter1, class Iter2, class Iter3>
174 void ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito) {
178 typedef typename Iter3::PixelType PixelType;
180 switch (mTypeOfOperation) {
182 while (!ito.IsAtEnd()) {
183 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
188 while (!ito.IsAtEnd()) {
189 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
194 while (!ito.IsAtEnd()) {
196 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
197 else ito.Set(mDefaultPixelValue);
202 while (!ito.IsAtEnd()) {
203 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
208 while (!ito.IsAtEnd()) {
209 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
213 case 5: // Absolute difference
214 while (!ito.IsAtEnd()) {
215 ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
219 case 6: // Squared differences
220 while (!ito.IsAtEnd()) {
221 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
225 case 7: // Difference
226 while (!ito.IsAtEnd()) {
227 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
231 case 8: // Relative Difference
232 while (!ito.IsAtEnd()) {
233 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
239 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
243 //--------------------------------------------------------------------
246 //--------------------------------------------------------------------
247 template<class args_info_type>
248 template<class Iter1, class Iter2>
249 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito) {
252 typedef typename Iter2::PixelType PixelType;
255 switch (mTypeOfOperation) {
257 while (!it.IsAtEnd()) {
258 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
263 while (!it.IsAtEnd()) {
264 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
269 while (!it.IsAtEnd()) {
271 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
272 else ito.Set(mDefaultPixelValue);
277 while (!it.IsAtEnd()) {
278 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
279 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
284 while (!it.IsAtEnd()) {
285 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
286 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
290 case 5: // Absolute value
291 while (!it.IsAtEnd()) {
292 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
293 // <= zero to avoid warning for unsigned types
294 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
298 case 6: // Squared value
299 while (!it.IsAtEnd()) {
300 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
305 while (!it.IsAtEnd()) {
307 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
308 else ito.Set(mDefaultPixelValue);
313 while (!it.IsAtEnd()) {
314 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
319 while (!it.IsAtEnd()) {
321 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
323 if (it.Get() ==0) ito.Set(0);
324 else ito.Set(mDefaultPixelValue);
330 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
334 //--------------------------------------------------------------------
338 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX