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) {
27 InitializeImageType<2>();
28 InitializeImageType<3>();
29 InitializeImageType<4>();
30 mIsOperationUseASecondImage = false;
31 mOverwriteInputImage = true;
33 //--------------------------------------------------------------------
36 //--------------------------------------------------------------------
37 template<class args_info_type>
38 template<unsigned int Dim>
39 void ImageArithmGenericFilter<args_info_type>::InitializeImageType() {
40 ADD_DEFAULT_IMAGE_TYPES(Dim);
42 //--------------------------------------------------------------------
45 //--------------------------------------------------------------------
46 template<class args_info_type>
47 void ImageArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b) {
48 mOverwriteInputImage = b;
50 //--------------------------------------------------------------------
53 //--------------------------------------------------------------------
54 template<class args_info_type>
55 void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a) {
59 SetIOVerbose(mArgsInfo.verbose_flag);
60 mTypeOfOperation = mArgsInfo.operation_arg;
61 mDefaultPixelValue = mArgsInfo.pixelValue_arg;
62 mScalar = mArgsInfo.scalar_arg;
63 mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
65 if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
67 if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
68 if (mArgsInfo.input2_given) {
69 mIsOperationUseASecondImage = true;
70 AddInputFilename(mArgsInfo.input2_arg);
73 if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
75 // Check type of operation (with scalar or with other image)
76 if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
77 std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
80 if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
81 if (mArgsInfo.operation_arg < 5) {
82 std::cerr << "Such operation need the --scalar option." << std::endl;
87 //--------------------------------------------------------------------
90 //--------------------------------------------------------------------
91 template<class args_info_type>
92 template<class ImageType>
93 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType() {
95 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
97 // Set input image iterator
98 typedef itk::ImageRegionIterator<ImageType> IteratorType;
99 IteratorType it(input1, input1->GetLargestPossibleRegion());
102 typename ImageType::Pointer input2 = this->template GetInput<ImageType>(1);
105 if (mIsOperationUseASecondImage) {
107 input2 = this->template GetInput<ImageType>(1);
108 // Set input image iterator
109 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
112 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
113 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
114 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
115 // << typeid(PixelType).name()
117 mOverwriteInputImage = false;
120 // ---------------- Overwrite input Image ---------------------
121 if (mOverwriteInputImage) {
122 // Set output iterator (to input1)
123 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
124 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
125 else ComputeImage(it, ito);
126 this->template SetNextOutput<ImageType>(input1);
129 // ---------------- Create new output Image ---------------------
131 if (mOutputIsFloat) {
132 // Create output image
133 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
134 typename OutputImageType::Pointer output = OutputImageType::New();
135 output->SetRegions(input1->GetLargestPossibleRegion());
136 output->SetOrigin(input1->GetOrigin());
137 output->SetSpacing(input1->GetSpacing());
139 // Set output iterator
140 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
141 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
142 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
143 else ComputeImage(it, ito);
144 this->template SetNextOutput<OutputImageType>(output);
147 // Create output image
148 typedef ImageType OutputImageType;
149 typename OutputImageType::Pointer output = OutputImageType::New();
150 output->SetRegions(input1->GetLargestPossibleRegion());
151 output->SetOrigin(input1->GetOrigin());
152 output->SetSpacing(input1->GetSpacing());
154 // Set output iterator
155 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
156 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
157 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
158 else ComputeImage(it, ito);
159 this->template SetNextOutput<OutputImageType>(output);
163 //--------------------------------------------------------------------
165 //--------------------------------------------------------------------
166 template<class args_info_type>
167 template<class Iter1, class Iter2, class Iter3>
168 void ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito) {
172 typedef typename Iter3::PixelType PixelType;
174 switch (mTypeOfOperation) {
176 while (!ito.IsAtEnd()) {
177 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
182 while (!ito.IsAtEnd()) {
183 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
188 while (!ito.IsAtEnd()) {
190 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
191 else ito.Set(mDefaultPixelValue);
196 while (!ito.IsAtEnd()) {
197 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
202 while (!ito.IsAtEnd()) {
203 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
207 case 5: // Absolute difference
208 while (!ito.IsAtEnd()) {
209 ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
213 case 6: // Squared differences
214 while (!ito.IsAtEnd()) {
215 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
219 case 7: // Difference
220 while (!ito.IsAtEnd()) {
221 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
225 case 8: // Relative Difference
226 while (!ito.IsAtEnd()) {
227 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
233 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
237 //--------------------------------------------------------------------
240 //--------------------------------------------------------------------
241 template<class args_info_type>
242 template<class Iter1, class Iter2>
243 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito) {
246 typedef typename Iter2::PixelType PixelType;
249 switch (mTypeOfOperation) {
251 while (!it.IsAtEnd()) {
252 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
257 while (!it.IsAtEnd()) {
258 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
263 while (!it.IsAtEnd()) {
265 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
266 else ito.Set(mDefaultPixelValue);
271 while (!it.IsAtEnd()) {
272 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
273 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
278 while (!it.IsAtEnd()) {
279 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
280 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
284 case 5: // Absolute value
285 while (!it.IsAtEnd()) {
286 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
287 // <= zero to avoid warning for unsigned types
288 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
292 case 6: // Squared value
293 while (!it.IsAtEnd()) {
294 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
299 while (!it.IsAtEnd()) {
301 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
302 else ito.Set(mDefaultPixelValue);
307 while (!it.IsAtEnd()) {
308 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
313 while (!it.IsAtEnd()) {
315 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
317 if (it.Get() ==0) ito.Set(0);
318 else ito.Set(mDefaultPixelValue);
324 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
328 //--------------------------------------------------------------------
332 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX