]> Creatis software - clitk.git/blob - tools/clitkImageArithmGenericFilter.txx
Add "normalize" and "divide by a scalar"
[clitk.git] / tools / clitkImageArithmGenericFilter.txx
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://www.centreleonberard.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 #ifndef CLITKIMAGEARITHMGENERICFILTER_TXX
19 #define CLITKIMAGEARITHMGENERICFILTER_TXX
20
21 #include "clitkImageCommon.h"
22
23 #include "itkMinimumMaximumImageCalculator.h"
24
25 namespace clitk
26 {
27
28 //--------------------------------------------------------------------
29 template<class args_info_type>
30 ImageArithmGenericFilter<args_info_type>::ImageArithmGenericFilter()
31   :ImageToImageGenericFilter<Self>("ImageArithmGenericFilter"),mTypeOfOperation(0)
32 {
33   InitializeImageType<2>();
34   InitializeImageType<3>();
35   InitializeImageType<4>();
36   mIsOperationUseASecondImage = false;
37   mOverwriteInputImage = true;
38 }
39 //--------------------------------------------------------------------
40
41
42 //--------------------------------------------------------------------
43 template<class args_info_type>
44 template<unsigned int Dim>
45 void ImageArithmGenericFilter<args_info_type>::InitializeImageType()
46 {
47   ADD_DEFAULT_IMAGE_TYPES(Dim);
48 }
49 //--------------------------------------------------------------------
50
51
52 //--------------------------------------------------------------------
53 template<class args_info_type>
54 void ImageArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b)
55 {
56   mOverwriteInputImage = b;
57 }
58 //--------------------------------------------------------------------
59
60
61 //--------------------------------------------------------------------
62 template<class args_info_type>
63 void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
64 {
65   mArgsInfo=a;
66
67   // Set value
68   SetIOVerbose(mArgsInfo.verbose_flag);
69   mTypeOfOperation = mArgsInfo.operation_arg;
70   mDefaultPixelValue = mArgsInfo.pixelValue_arg;
71   mScalar = mArgsInfo.scalar_arg;
72   mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
73
74   if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
75
76   if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
77   if (mArgsInfo.input2_given) {
78     mIsOperationUseASecondImage = true;
79     AddInputFilename(mArgsInfo.input2_arg);
80   }
81
82   if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
83
84   // Check type of operation (with scalar or with other image)
85   if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
86     std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
87     exit(-1);
88   }
89   if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
90     if (mArgsInfo.operation_arg < 5) {
91       std::cerr << "Such operation need the --scalar option." << std::endl;
92       exit(-1);
93     }
94   }
95 }
96 //--------------------------------------------------------------------
97
98
99 //--------------------------------------------------------------------
100 template<class args_info_type>
101 template<class ImageType>
102 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
103 {
104   // Read input1
105   typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
106
107   // Set input image iterator
108   typedef itk::ImageRegionIterator<ImageType> IteratorType;
109   IteratorType it(input1, input1->GetLargestPossibleRegion());
110
111   // typedef input2
112   typename ImageType::Pointer input2 = NULL;
113   IteratorType it2;
114
115   // Special case for normalisation
116   if (mTypeOfOperation == 12) {
117     typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
118     typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
119     ff->SetImage(input1);
120     ff->ComputeMaximum();
121     mScalar = ff->GetMaximum();
122     DD(mScalar);
123     DD("normalisation");
124     mTypeOfOperation = 11; // divide
125   }
126
127   if (mIsOperationUseASecondImage) {
128       // Read input2
129       input2 = this->template GetInput<ImageType>(1);
130       // Set input image iterator
131       it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
132       // Check dimension
133       if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input1, input2)) {
134           std::cerr << "* ERROR * the images (input and input2) must have the same size & spacing";
135           return;
136       }
137   }
138
139   // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
140   if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
141     // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
142     //                     << typeid(PixelType).name()
143     //                     << std::endl;
144     mOverwriteInputImage = false;
145   }
146
147   // ---------------- Overwrite input Image ---------------------
148   if (mOverwriteInputImage) {
149     // Set output iterator (to input1)
150     IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
151     if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
152     else ComputeImage(it, ito);
153     this->template SetNextOutput<ImageType>(input1);
154   }
155
156   // ---------------- Create new output Image ---------------------
157   else {
158     if (mOutputIsFloat) {
159       // Create output image
160       typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
161       typename OutputImageType::Pointer output = OutputImageType::New();
162       output->SetRegions(input1->GetLargestPossibleRegion());
163       output->SetOrigin(input1->GetOrigin());
164       output->SetSpacing(input1->GetSpacing());
165       output->Allocate();
166       // Set output iterator
167       typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
168       IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
169       if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
170       else ComputeImage(it, ito);
171       this->template SetNextOutput<OutputImageType>(output);
172     } else {
173       // Create output image
174       typedef ImageType OutputImageType;
175       typename OutputImageType::Pointer output = OutputImageType::New();
176       output->SetRegions(input1->GetLargestPossibleRegion());
177       output->SetOrigin(input1->GetOrigin());
178       output->SetSpacing(input1->GetSpacing());
179       output->Allocate();
180       // Set output iterator
181       typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
182       IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
183       if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
184       else ComputeImage(it, ito);
185       this->template SetNextOutput<OutputImageType>(output);
186     }
187   }
188 }
189 //--------------------------------------------------------------------
190
191 //--------------------------------------------------------------------
192 template<class args_info_type>
193 template<class Iter1, class Iter2, class Iter3>
194 void  ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
195 {
196   it1.GoToBegin();
197   it2.GoToBegin();
198   ito.GoToBegin();
199   typedef typename Iter3::PixelType PixelType;
200
201   switch (mTypeOfOperation) {
202   case 0: // Addition
203     while (!ito.IsAtEnd()) {
204       ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
205       ++it1;
206       ++it2;
207       ++ito;
208     }
209     break;
210   case 1: // Multiply
211     while (!ito.IsAtEnd()) {
212       ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
213       ++it1;
214       ++it2;
215       ++ito;
216     }
217     break;
218   case 2: // Divide
219     while (!ito.IsAtEnd()) {
220       if (it1.Get() != 0)
221         ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
222       else ito.Set(mDefaultPixelValue);
223       ++it1;
224       ++it2;
225       ++ito;
226     }
227     break;
228   case 3: // Max
229     while (!ito.IsAtEnd()) {
230       if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
231       ++it1;
232       ++it2;
233       ++ito;
234     }
235     break;
236   case 4: // Min
237     while (!ito.IsAtEnd()) {
238       if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
239       ++it1;
240       ++it2;
241       ++ito;
242     }
243     break;
244   case 5: // Absolute difference
245     while (!ito.IsAtEnd()) {
246       ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
247       ++it1;
248       ++it2;
249       ++ito;
250     }
251     break;
252   case 6: // Squared differences
253     while (!ito.IsAtEnd()) {
254       ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
255       ++it1;
256       ++it2;
257       ++ito;
258     }
259     break;
260   case 7: // Difference
261     while (!ito.IsAtEnd()) {
262       ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
263       ++it1;
264       ++it2;
265       ++ito;
266     }
267     break;
268   case 8: // Relative Difference
269     while (!ito.IsAtEnd()) {
270       if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
271       else ito.Set(0.0);
272       ++it1;
273       ++it2;
274       ++ito;
275     }
276     break;
277   default: // error ?
278     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
279     exit(-1);
280   }
281 }
282 //--------------------------------------------------------------------
283
284
285 //--------------------------------------------------------------------
286 template<class args_info_type>
287 template<class Iter1, class Iter2>
288 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
289 {
290   ito.GoToBegin();
291   it.GoToBegin();
292   typedef typename Iter2::PixelType PixelType;
293
294   // Perform operation
295   switch (mTypeOfOperation) {
296   case 0: // Addition
297     while (!it.IsAtEnd()) {
298       ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
299       ++it;
300       ++ito;
301     }
302     break;
303   case 1: // Multiply
304     while (!it.IsAtEnd()) {
305       ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
306       ++it;
307       ++ito;
308     }
309     break;
310   case 2: // Inverse
311     while (!it.IsAtEnd()) {
312       if (it.Get() != 0)
313         ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
314       else ito.Set(mDefaultPixelValue);
315       ++it;
316       ++ito;
317     }
318     break;
319   case 3: // Max
320     while (!it.IsAtEnd()) {
321       if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
322       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
323       ++it;
324       ++ito;
325     }
326     break;
327   case 4: // Min
328     while (!it.IsAtEnd()) {
329       if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
330       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
331       ++it;
332       ++ito;
333     }
334     break;
335   case 5: // Absolute value
336     while (!it.IsAtEnd()) {
337       if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
338       // <= zero to avoid warning for unsigned types
339       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
340       ++it;
341       ++ito;
342     }
343     break;
344   case 6: // Squared value
345     while (!it.IsAtEnd()) {
346       ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
347       ++it;
348       ++ito;
349     }
350     break;
351   case 7: // Log
352     while (!it.IsAtEnd()) {
353       if (it.Get() > 0)
354         ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
355       else ito.Set(mDefaultPixelValue);
356       ++it;
357       ++ito;
358     }
359     break;
360   case 8: // exp
361     while (!it.IsAtEnd()) {
362       ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
363       ++it;
364       ++ito;
365     }
366     break;
367   case 9: // sqrt
368     while (!it.IsAtEnd()) {
369       if (it.Get() > 0)
370         ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
371       else {
372         if (it.Get() ==0) ito.Set(0);
373         else ito.Set(mDefaultPixelValue);
374       }
375       ++it;
376       ++ito;
377     }
378     break;
379   case 10: // exp
380     while (!it.IsAtEnd()) {
381       ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
382       ++it;
383       ++ito;
384     }
385     break;
386   case 11: // divide
387     while (!it.IsAtEnd()) {
388       ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() / mScalar) );
389       ++it;
390       ++ito;
391     }
392     break;
393   default: // error ?
394     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
395     exit(-1);
396   }
397 }
398 //--------------------------------------------------------------------
399
400 } // end namespace
401
402 #endif  //#define CLITKIMAGEARITHMGENERICFILTER_TXX