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();
124 mTypeOfOperation = 11; // divide
127 if (mIsOperationUseASecondImage) {
129 input2 = this->template GetInput<ImageType>(1);
130 // Set input image iterator
131 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
133 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input1, input2)) {
134 std::cerr << "* ERROR * the images (input and input2) must have the same size & spacing";
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()
144 mOverwriteInputImage = false;
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);
156 // ---------------- Create new output Image ---------------------
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());
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);
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());
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);
189 //--------------------------------------------------------------------
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)
199 typedef typename Iter3::PixelType PixelType;
201 switch (mTypeOfOperation) {
203 while (!ito.IsAtEnd()) {
204 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
211 while (!ito.IsAtEnd()) {
212 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
219 while (!ito.IsAtEnd()) {
221 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
222 else ito.Set(mDefaultPixelValue);
229 while (!ito.IsAtEnd()) {
230 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
237 while (!ito.IsAtEnd()) {
238 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
244 case 5: // Absolute difference
245 while (!ito.IsAtEnd()) {
246 ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
252 case 6: // Squared differences
253 while (!ito.IsAtEnd()) {
254 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
260 case 7: // Difference
261 while (!ito.IsAtEnd()) {
262 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
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());
278 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
282 //--------------------------------------------------------------------
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)
292 typedef typename Iter2::PixelType PixelType;
295 switch (mTypeOfOperation) {
297 while (!it.IsAtEnd()) {
298 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
304 while (!it.IsAtEnd()) {
305 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
311 while (!it.IsAtEnd()) {
313 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
314 else ito.Set(mDefaultPixelValue);
320 while (!it.IsAtEnd()) {
321 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
322 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
328 while (!it.IsAtEnd()) {
329 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
330 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
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()));
344 case 6: // Squared value
345 while (!it.IsAtEnd()) {
346 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
352 while (!it.IsAtEnd()) {
354 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
355 else ito.Set(mDefaultPixelValue);
361 while (!it.IsAtEnd()) {
362 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
368 while (!it.IsAtEnd()) {
370 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
372 if (it.Get() ==0) ito.Set(0);
373 else ito.Set(mDefaultPixelValue);
380 while (!it.IsAtEnd()) {
381 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
387 while (!it.IsAtEnd()) {
388 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() / mScalar) );
394 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
398 //--------------------------------------------------------------------
402 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX