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);
48 ADD_VEC_IMAGE_TYPE(3u,3u,float);
49 ADD_VEC_IMAGE_TYPE(3u,3u,double);
51 //--------------------------------------------------------------------
54 //--------------------------------------------------------------------
55 template<class args_info_type>
56 void ImageArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b)
58 mOverwriteInputImage = b;
60 //--------------------------------------------------------------------
63 //--------------------------------------------------------------------
64 template<class args_info_type>
65 void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
70 SetIOVerbose(mArgsInfo.verbose_flag);
71 mTypeOfOperation = mArgsInfo.operation_arg;
72 mDefaultPixelValue = mArgsInfo.pixelValue_arg;
73 mScalar = mArgsInfo.scalar_arg;
74 mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
76 if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
78 if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
79 if (mArgsInfo.input2_given) {
80 mIsOperationUseASecondImage = true;
81 AddInputFilename(mArgsInfo.input2_arg);
84 if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
86 // Check type of operation (with scalar or with other image)
87 if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
88 std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
91 if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
92 if (mArgsInfo.operation_arg < 5) {
93 std::cerr << "Such operation need the --scalar option." << std::endl;
98 //--------------------------------------------------------------------
101 //--------------------------------------------------------------------
102 template<class args_info_type>
103 template<class ImageType>
104 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
107 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
109 // Set input image iterator
110 typedef itk::ImageRegionIterator<ImageType> IteratorType;
111 IteratorType it(input1, input1->GetLargestPossibleRegion());
114 typename ImageType::Pointer input2 = NULL;
117 // Special case for normalisation
118 if (mTypeOfOperation == 12) {
119 typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
120 typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
121 ff->SetImage(input1);
122 ff->ComputeMaximum();
123 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::HaveSameSize<ImageType, ImageType>(input1, input2)) {
134 itkExceptionMacro(<< "The images (input and input2) must have the same size");
136 if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
137 itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
138 << "Using first input's information.");
142 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
143 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
144 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
145 // << typeid(PixelType).name()
147 mOverwriteInputImage = false;
150 // ---------------- Overwrite input Image ---------------------
151 if (mOverwriteInputImage) {
152 // Set output iterator (to input1)
153 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
154 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
155 else ComputeImage(it, ito);
156 this->template SetNextOutput<ImageType>(input1);
159 // ---------------- Create new output Image ---------------------
161 if (mOutputIsFloat) {
162 // Create output image
163 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
164 typename OutputImageType::Pointer output = OutputImageType::New();
165 output->SetRegions(input1->GetLargestPossibleRegion());
166 output->SetOrigin(input1->GetOrigin());
167 output->SetSpacing(input1->GetSpacing());
169 // Set output iterator
170 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
171 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
172 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
173 else ComputeImage(it, ito);
174 this->template SetNextOutput<OutputImageType>(output);
176 // Create output image
177 typedef ImageType OutputImageType;
178 typename OutputImageType::Pointer output = OutputImageType::New();
179 output->SetRegions(input1->GetLargestPossibleRegion());
180 output->SetOrigin(input1->GetOrigin());
181 output->SetSpacing(input1->GetSpacing());
183 // Set output iterator
184 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
185 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
186 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
187 else ComputeImage(it, ito);
188 this->template SetNextOutput<OutputImageType>(output);
192 //--------------------------------------------------------------------
194 //--------------------------------------------------------------------
195 template<class args_info_type>
196 template<class Iter1, class Iter2, class Iter3>
197 void ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
202 typedef typename Iter3::PixelType PixelType;
204 switch (mTypeOfOperation) {
206 while (!ito.IsAtEnd()) {
207 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
214 while (!ito.IsAtEnd()) {
215 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
222 while (!ito.IsAtEnd()) {
224 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
225 else ito.Set(mDefaultPixelValue);
232 while (!ito.IsAtEnd()) {
233 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
240 while (!ito.IsAtEnd()) {
241 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
247 case 5: // Absolute difference
248 while (!ito.IsAtEnd()) {
249 ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
255 case 6: // Squared differences
256 while (!ito.IsAtEnd()) {
257 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
263 case 7: // Difference
264 while (!ito.IsAtEnd()) {
265 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
271 case 8: // Relative Difference
272 while (!ito.IsAtEnd()) {
273 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
281 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
285 //--------------------------------------------------------------------
288 //--------------------------------------------------------------------
289 template<class args_info_type>
290 template<class Iter1, class Iter2>
291 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
295 typedef typename Iter2::PixelType PixelType;
298 switch (mTypeOfOperation) {
300 while (!it.IsAtEnd()) {
301 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
307 while (!it.IsAtEnd()) {
308 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
314 while (!it.IsAtEnd()) {
316 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
317 else ito.Set(mDefaultPixelValue);
323 while (!it.IsAtEnd()) {
324 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
325 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
331 while (!it.IsAtEnd()) {
332 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
333 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
338 case 5: // Absolute value
339 while (!it.IsAtEnd()) {
340 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
341 // <= zero to avoid warning for unsigned types
342 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
347 case 6: // Squared value
348 while (!it.IsAtEnd()) {
349 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
355 while (!it.IsAtEnd()) {
357 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
358 else ito.Set(mDefaultPixelValue);
364 while (!it.IsAtEnd()) {
365 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
371 while (!it.IsAtEnd()) {
373 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
375 if (it.Get() ==0) ito.Set(0);
376 else ito.Set(mDefaultPixelValue);
383 while (!it.IsAtEnd()) {
384 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
390 while (!it.IsAtEnd()) {
391 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() / mScalar) );
397 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
401 //--------------------------------------------------------------------
407 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX