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)
28 InitializeImageType<2>();
29 InitializeImageType<3>();
30 InitializeImageType<4>();
31 mIsOperationUseASecondImage = false;
32 mOverwriteInputImage = true;
34 //--------------------------------------------------------------------
37 //--------------------------------------------------------------------
38 template<class args_info_type>
39 template<unsigned int Dim>
40 void ImageArithmGenericFilter<args_info_type>::InitializeImageType()
42 ADD_DEFAULT_IMAGE_TYPES(Dim);
44 //--------------------------------------------------------------------
47 //--------------------------------------------------------------------
48 template<class args_info_type>
49 void ImageArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b)
51 mOverwriteInputImage = b;
53 //--------------------------------------------------------------------
56 //--------------------------------------------------------------------
57 template<class args_info_type>
58 void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
63 SetIOVerbose(mArgsInfo.verbose_flag);
64 mTypeOfOperation = mArgsInfo.operation_arg;
65 mDefaultPixelValue = mArgsInfo.pixelValue_arg;
66 mScalar = mArgsInfo.scalar_arg;
67 mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
69 if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
71 if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
72 if (mArgsInfo.input2_given) {
73 mIsOperationUseASecondImage = true;
74 AddInputFilename(mArgsInfo.input2_arg);
77 if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
79 // Check type of operation (with scalar or with other image)
80 if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
81 std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
84 if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
85 if (mArgsInfo.operation_arg < 5) {
86 std::cerr << "Such operation need the --scalar option." << std::endl;
91 //--------------------------------------------------------------------
94 //--------------------------------------------------------------------
95 template<class args_info_type>
96 template<class ImageType>
97 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
100 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
102 // Set input image iterator
103 typedef itk::ImageRegionIterator<ImageType> IteratorType;
104 IteratorType it(input1, input1->GetLargestPossibleRegion());
107 typename ImageType::Pointer input2 = this->template GetInput<ImageType>(1);
110 if (mIsOperationUseASecondImage) {
112 input2 = this->template GetInput<ImageType>(1);
113 // Set input image iterator
114 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
117 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
118 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
119 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
120 // << typeid(PixelType).name()
122 mOverwriteInputImage = false;
125 // ---------------- Overwrite input Image ---------------------
126 if (mOverwriteInputImage) {
127 // Set output iterator (to input1)
128 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
129 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
130 else ComputeImage(it, ito);
131 this->template SetNextOutput<ImageType>(input1);
134 // ---------------- Create new output Image ---------------------
136 if (mOutputIsFloat) {
137 // Create output image
138 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
139 typename OutputImageType::Pointer output = OutputImageType::New();
140 output->SetRegions(input1->GetLargestPossibleRegion());
141 output->SetOrigin(input1->GetOrigin());
142 output->SetSpacing(input1->GetSpacing());
144 // Set output iterator
145 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
146 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
147 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
148 else ComputeImage(it, ito);
149 this->template SetNextOutput<OutputImageType>(output);
151 // Create output image
152 typedef ImageType OutputImageType;
153 typename OutputImageType::Pointer output = OutputImageType::New();
154 output->SetRegions(input1->GetLargestPossibleRegion());
155 output->SetOrigin(input1->GetOrigin());
156 output->SetSpacing(input1->GetSpacing());
158 // Set output iterator
159 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
160 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
161 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
162 else ComputeImage(it, ito);
163 this->template SetNextOutput<OutputImageType>(output);
167 //--------------------------------------------------------------------
169 //--------------------------------------------------------------------
170 template<class args_info_type>
171 template<class Iter1, class Iter2, class Iter3>
172 void ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
177 typedef typename Iter3::PixelType PixelType;
179 switch (mTypeOfOperation) {
181 while (!ito.IsAtEnd()) {
182 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
189 while (!ito.IsAtEnd()) {
190 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
197 while (!ito.IsAtEnd()) {
199 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
200 else ito.Set(mDefaultPixelValue);
207 while (!ito.IsAtEnd()) {
208 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
215 while (!ito.IsAtEnd()) {
216 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
222 case 5: // Absolute difference
223 while (!ito.IsAtEnd()) {
224 ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
230 case 6: // Squared differences
231 while (!ito.IsAtEnd()) {
232 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
238 case 7: // Difference
239 while (!ito.IsAtEnd()) {
240 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
246 case 8: // Relative Difference
247 while (!ito.IsAtEnd()) {
248 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
256 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
260 //--------------------------------------------------------------------
263 //--------------------------------------------------------------------
264 template<class args_info_type>
265 template<class Iter1, class Iter2>
266 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
270 typedef typename Iter2::PixelType PixelType;
273 switch (mTypeOfOperation) {
275 while (!it.IsAtEnd()) {
276 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
282 while (!it.IsAtEnd()) {
283 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
289 while (!it.IsAtEnd()) {
291 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
292 else ito.Set(mDefaultPixelValue);
298 while (!it.IsAtEnd()) {
299 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
300 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
306 while (!it.IsAtEnd()) {
307 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
308 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
313 case 5: // Absolute value
314 while (!it.IsAtEnd()) {
315 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
316 // <= zero to avoid warning for unsigned types
317 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
322 case 6: // Squared value
323 while (!it.IsAtEnd()) {
324 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
330 while (!it.IsAtEnd()) {
332 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
333 else ito.Set(mDefaultPixelValue);
339 while (!it.IsAtEnd()) {
340 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
346 while (!it.IsAtEnd()) {
348 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
350 if (it.Get() ==0) ito.Set(0);
351 else ito.Set(mDefaultPixelValue);
358 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
362 //--------------------------------------------------------------------
366 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX