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
21 #include "clitkImageCommon.h"
26 //--------------------------------------------------------------------
27 template<class args_info_type>
28 ImageArithmGenericFilter<args_info_type>::ImageArithmGenericFilter()
29 :ImageToImageGenericFilter<Self>("ImageArithmGenericFilter"),mTypeOfOperation(0)
31 InitializeImageType<2>();
32 InitializeImageType<3>();
33 InitializeImageType<4>();
34 mIsOperationUseASecondImage = false;
35 mOverwriteInputImage = true;
37 //--------------------------------------------------------------------
40 //--------------------------------------------------------------------
41 template<class args_info_type>
42 template<unsigned int Dim>
43 void ImageArithmGenericFilter<args_info_type>::InitializeImageType()
45 ADD_DEFAULT_IMAGE_TYPES(Dim);
47 //--------------------------------------------------------------------
50 //--------------------------------------------------------------------
51 template<class args_info_type>
52 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)
66 SetIOVerbose(mArgsInfo.verbose_flag);
67 mTypeOfOperation = mArgsInfo.operation_arg;
68 mDefaultPixelValue = mArgsInfo.pixelValue_arg;
69 mScalar = mArgsInfo.scalar_arg;
70 mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
72 if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
74 if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
75 if (mArgsInfo.input2_given) {
76 mIsOperationUseASecondImage = true;
77 AddInputFilename(mArgsInfo.input2_arg);
80 if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
82 // Check type of operation (with scalar or with other image)
83 if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
84 std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
87 if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
88 if (mArgsInfo.operation_arg < 5) {
89 std::cerr << "Such operation need the --scalar option." << std::endl;
94 //--------------------------------------------------------------------
97 //--------------------------------------------------------------------
98 template<class args_info_type>
99 template<class ImageType>
100 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
103 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
105 // Set input image iterator
106 typedef itk::ImageRegionIterator<ImageType> IteratorType;
107 IteratorType it(input1, input1->GetLargestPossibleRegion());
110 typename ImageType::Pointer input2 = NULL;
114 if (mIsOperationUseASecondImage) {
116 input2 = this->template GetInput<ImageType>(1);
117 // Set input image iterator
118 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
120 if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input1, input2)) {
121 std::cerr << "* ERROR * the images (input and input2) must have the same size & spacing";
126 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
127 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
128 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
129 // << typeid(PixelType).name()
131 mOverwriteInputImage = false;
134 // ---------------- Overwrite input Image ---------------------
135 if (mOverwriteInputImage) {
136 // Set output iterator (to input1)
137 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
138 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
139 else ComputeImage(it, ito);
140 this->template SetNextOutput<ImageType>(input1);
143 // ---------------- Create new output Image ---------------------
145 if (mOutputIsFloat) {
146 // Create output image
147 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
148 typename OutputImageType::Pointer output = OutputImageType::New();
149 output->SetRegions(input1->GetLargestPossibleRegion());
150 output->SetOrigin(input1->GetOrigin());
151 output->SetSpacing(input1->GetSpacing());
153 // Set output iterator
154 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
155 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
156 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
157 else ComputeImage(it, ito);
158 this->template SetNextOutput<OutputImageType>(output);
160 // Create output image
161 typedef ImageType OutputImageType;
162 typename OutputImageType::Pointer output = OutputImageType::New();
163 output->SetRegions(input1->GetLargestPossibleRegion());
164 output->SetOrigin(input1->GetOrigin());
165 output->SetSpacing(input1->GetSpacing());
167 // Set output iterator
168 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
169 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
170 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
171 else ComputeImage(it, ito);
172 this->template SetNextOutput<OutputImageType>(output);
176 //--------------------------------------------------------------------
178 //--------------------------------------------------------------------
179 template<class args_info_type>
180 template<class Iter1, class Iter2, class Iter3>
181 void ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
186 typedef typename Iter3::PixelType PixelType;
188 switch (mTypeOfOperation) {
190 while (!ito.IsAtEnd()) {
191 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
198 while (!ito.IsAtEnd()) {
199 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
206 while (!ito.IsAtEnd()) {
208 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
209 else ito.Set(mDefaultPixelValue);
216 while (!ito.IsAtEnd()) {
217 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
224 while (!ito.IsAtEnd()) {
225 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
231 case 5: // Absolute difference
232 while (!ito.IsAtEnd()) {
233 ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
239 case 6: // Squared differences
240 while (!ito.IsAtEnd()) {
241 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
247 case 7: // Difference
248 while (!ito.IsAtEnd()) {
249 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
255 case 8: // Relative Difference
256 while (!ito.IsAtEnd()) {
257 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
265 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
269 //--------------------------------------------------------------------
272 //--------------------------------------------------------------------
273 template<class args_info_type>
274 template<class Iter1, class Iter2>
275 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
279 typedef typename Iter2::PixelType PixelType;
282 switch (mTypeOfOperation) {
284 while (!it.IsAtEnd()) {
285 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
291 while (!it.IsAtEnd()) {
292 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
298 while (!it.IsAtEnd()) {
300 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
301 else ito.Set(mDefaultPixelValue);
307 while (!it.IsAtEnd()) {
308 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
309 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
315 while (!it.IsAtEnd()) {
316 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
317 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
322 case 5: // Absolute value
323 while (!it.IsAtEnd()) {
324 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
325 // <= zero to avoid warning for unsigned types
326 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
331 case 6: // Squared value
332 while (!it.IsAtEnd()) {
333 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
339 while (!it.IsAtEnd()) {
341 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
342 else ito.Set(mDefaultPixelValue);
348 while (!it.IsAtEnd()) {
349 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
355 while (!it.IsAtEnd()) {
357 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
359 if (it.Get() ==0) ito.Set(0);
360 else ito.Set(mDefaultPixelValue);
367 while (!it.IsAtEnd()) {
368 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
374 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
378 //--------------------------------------------------------------------
382 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX