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 CLITKIMAGEARITHMGENERICFILTER_TXX
19 #define CLITKIMAGEARITHMGENERICFILTER_TXX
21 #include "clitkImageCommon.h"
23 #include "itkMinimumMaximumImageCalculator.h"
28 //--------------------------------------------------------------------
29 template<class args_info_type>
30 ImageArithmGenericFilter<args_info_type>::ImageArithmGenericFilter()
31 :ImageToImageGenericFilter<Self>("ImageArithmGenericFilter"),mTypeOfOperation(0)
33 InitializeImageType<2>();
34 InitializeImageType<3>();
35 InitializeImageType<4>();
36 mIsOperationUseASecondImage = false;
37 mOverwriteInputImage = true;
39 //--------------------------------------------------------------------
42 //--------------------------------------------------------------------
43 template<class args_info_type>
44 template<unsigned int Dim>
45 void ImageArithmGenericFilter<args_info_type>::InitializeImageType()
47 ADD_DEFAULT_IMAGE_TYPES(Dim);
49 //--------------------------------------------------------------------
52 //--------------------------------------------------------------------
53 template<class args_info_type>
54 void ImageArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b)
56 mOverwriteInputImage = b;
58 //--------------------------------------------------------------------
61 //--------------------------------------------------------------------
62 template<class args_info_type>
63 void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
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;
74 if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
76 if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
77 if (mArgsInfo.input2_given) {
78 mIsOperationUseASecondImage = true;
79 AddInputFilename(mArgsInfo.input2_arg);
82 if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
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;
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;
96 //--------------------------------------------------------------------
99 //--------------------------------------------------------------------
100 template<class args_info_type>
101 template<class ImageType>
102 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
105 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
107 // Set input image iterator
108 typedef itk::ImageRegionIterator<ImageType> IteratorType;
109 IteratorType it(input1, input1->GetLargestPossibleRegion());
112 typename ImageType::Pointer input2 = NULL;
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 mTypeOfOperation = 11; // divide
125 if (mIsOperationUseASecondImage) {
127 input2 = this->template GetInput<ImageType>(1);
128 // Set input image iterator
129 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
131 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input1, input2)) {
132 std::cerr << "* ERROR * the images (input and input2) must have the same size & spacing";
137 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
138 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
139 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
140 // << typeid(PixelType).name()
142 mOverwriteInputImage = false;
145 // ---------------- Overwrite input Image ---------------------
146 if (mOverwriteInputImage) {
147 // Set output iterator (to input1)
148 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
149 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
150 else ComputeImage(it, ito);
151 this->template SetNextOutput<ImageType>(input1);
154 // ---------------- Create new output Image ---------------------
156 if (mOutputIsFloat) {
157 // Create output image
158 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
159 typename OutputImageType::Pointer output = OutputImageType::New();
160 output->SetRegions(input1->GetLargestPossibleRegion());
161 output->SetOrigin(input1->GetOrigin());
162 output->SetSpacing(input1->GetSpacing());
164 // Set output iterator
165 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
166 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
167 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
168 else ComputeImage(it, ito);
169 this->template SetNextOutput<OutputImageType>(output);
171 // Create output image
172 typedef ImageType OutputImageType;
173 typename OutputImageType::Pointer output = OutputImageType::New();
174 output->SetRegions(input1->GetLargestPossibleRegion());
175 output->SetOrigin(input1->GetOrigin());
176 output->SetSpacing(input1->GetSpacing());
178 // Set output iterator
179 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
180 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
181 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
182 else ComputeImage(it, ito);
183 this->template SetNextOutput<OutputImageType>(output);
187 //--------------------------------------------------------------------
189 //--------------------------------------------------------------------
190 template<class args_info_type>
191 template<class Iter1, class Iter2, class Iter3>
192 void ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
197 typedef typename Iter3::PixelType PixelType;
199 switch (mTypeOfOperation) {
201 while (!ito.IsAtEnd()) {
202 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
209 while (!ito.IsAtEnd()) {
210 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
217 while (!ito.IsAtEnd()) {
219 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
220 else ito.Set(mDefaultPixelValue);
227 while (!ito.IsAtEnd()) {
228 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
235 while (!ito.IsAtEnd()) {
236 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
242 case 5: // Absolute difference
243 while (!ito.IsAtEnd()) {
244 ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
250 case 6: // Squared differences
251 while (!ito.IsAtEnd()) {
252 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
258 case 7: // Difference
259 while (!ito.IsAtEnd()) {
260 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
266 case 8: // Relative Difference
267 while (!ito.IsAtEnd()) {
268 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
276 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
280 //--------------------------------------------------------------------
283 //--------------------------------------------------------------------
284 template<class args_info_type>
285 template<class Iter1, class Iter2>
286 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
290 typedef typename Iter2::PixelType PixelType;
293 switch (mTypeOfOperation) {
295 while (!it.IsAtEnd()) {
296 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
302 while (!it.IsAtEnd()) {
303 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
309 while (!it.IsAtEnd()) {
311 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
312 else ito.Set(mDefaultPixelValue);
318 while (!it.IsAtEnd()) {
319 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
320 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
326 while (!it.IsAtEnd()) {
327 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
328 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
333 case 5: // Absolute value
334 while (!it.IsAtEnd()) {
335 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
336 // <= zero to avoid warning for unsigned types
337 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
342 case 6: // Squared value
343 while (!it.IsAtEnd()) {
344 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
350 while (!it.IsAtEnd()) {
352 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
353 else ito.Set(mDefaultPixelValue);
359 while (!it.IsAtEnd()) {
360 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
366 while (!it.IsAtEnd()) {
368 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
370 if (it.Get() ==0) ito.Set(0);
371 else ito.Set(mDefaultPixelValue);
378 while (!it.IsAtEnd()) {
379 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
385 while (!it.IsAtEnd()) {
386 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() / mScalar) );
392 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
396 //--------------------------------------------------------------------
400 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX