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_IMAGE_TYPE(Dim, char);
41 ADD_IMAGE_TYPE(Dim, uchar);
42 ADD_IMAGE_TYPE(Dim, short);
43 ADD_IMAGE_TYPE(Dim, ushort);
44 ADD_IMAGE_TYPE(Dim, int);
45 ADD_IMAGE_TYPE(Dim, float);
46 ADD_IMAGE_TYPE(Dim, double);
48 //--------------------------------------------------------------------
51 //--------------------------------------------------------------------
52 template<class args_info_type>
53 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) {
65 SetIOVerbose(mArgsInfo.verbose_flag);
66 mTypeOfOperation = mArgsInfo.operation_arg;
67 mDefaultPixelValue = mArgsInfo.pixelValue_arg;
68 mScalar = mArgsInfo.scalar_arg;
69 mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
71 if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
73 if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
74 if (mArgsInfo.input2_given) {
75 mIsOperationUseASecondImage = true;
76 AddInputFilename(mArgsInfo.input2_arg);
79 if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
81 // Check type of operation (with scalar or with other image)
82 if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
83 std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
86 if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
87 if (mArgsInfo.operation_arg < 5) {
88 std::cerr << "Such operation need the --scalar option." << std::endl;
93 //--------------------------------------------------------------------
96 //--------------------------------------------------------------------
97 template<class args_info_type>
98 template<class ImageType>
99 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType() {
101 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
102 typename ImageType::PixelType PixelType;
104 // Set input image iterator
105 typedef itk::ImageRegionIterator<ImageType> IteratorType;
106 IteratorType it(input1, input1->GetLargestPossibleRegion());
109 typename ImageType::Pointer input2 = this->template GetInput<ImageType>(1);
112 if (mIsOperationUseASecondImage) {
114 input2 = this->template GetInput<ImageType>(1);
115 // Set input image iterator
116 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
119 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
120 if (mOverwriteInputImage && mOutputIsFloat && (typeid(PixelType) != typeid(float))) {
121 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
122 // << typeid(PixelType).name()
124 mOverwriteInputImage = false;
127 // ---------------- Overwrite input Image ---------------------
128 if (mOverwriteInputImage) {
129 // Set output iterator (to input1)
130 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
131 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
132 else ComputeImage(it, ito);
133 this->template SetNextOutput<ImageType>(input1);
136 // ---------------- Create new output Image ---------------------
138 if (mOutputIsFloat) {
139 // Create output image
140 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
141 typename OutputImageType::Pointer output = OutputImageType::New();
142 output->SetRegions(input1->GetLargestPossibleRegion());
143 output->SetOrigin(input1->GetOrigin());
144 output->SetSpacing(input1->GetSpacing());
146 // Set output iterator
147 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
148 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
149 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
150 else ComputeImage(it, ito);
151 this->template SetNextOutput<OutputImageType>(output);
154 // Create output image
155 typedef ImageType OutputImageType;
156 typename OutputImageType::Pointer output = OutputImageType::New();
157 output->SetRegions(input1->GetLargestPossibleRegion());
158 output->SetOrigin(input1->GetOrigin());
159 output->SetSpacing(input1->GetSpacing());
161 // Set output iterator
162 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
163 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
164 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
165 else ComputeImage(it, ito);
166 this->template SetNextOutput<OutputImageType>(output);
170 //--------------------------------------------------------------------
172 //--------------------------------------------------------------------
173 template<class args_info_type>
174 template<class Iter1, class Iter2, class Iter3>
175 void ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito) {
179 typedef typename Iter3::PixelType PixelType;
181 switch (mTypeOfOperation) {
183 while (!ito.IsAtEnd()) {
184 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()) );
195 while (!ito.IsAtEnd()) {
197 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
198 else ito.Set(mDefaultPixelValue);
203 while (!ito.IsAtEnd()) {
204 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
209 while (!ito.IsAtEnd()) {
210 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
214 case 5: // Absolute difference
215 while (!ito.IsAtEnd()) {
216 ito.Set(PixelTypeDownCast<double, PixelType>(fabs(it2.Get()-it1.Get())));
220 case 6: // Squared differences
221 while (!ito.IsAtEnd()) {
222 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
226 case 7: // Difference
227 while (!ito.IsAtEnd()) {
228 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
232 case 8: // Relative Difference
233 while (!ito.IsAtEnd()) {
234 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
240 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
244 //--------------------------------------------------------------------
247 //--------------------------------------------------------------------
248 template<class args_info_type>
249 template<class Iter1, class Iter2>
250 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito) {
253 typedef typename Iter2::PixelType PixelType;
256 switch (mTypeOfOperation) {
258 while (!it.IsAtEnd()) {
259 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
264 while (!it.IsAtEnd()) {
265 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
270 while (!it.IsAtEnd()) {
272 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
273 else ito.Set(mDefaultPixelValue);
278 while (!it.IsAtEnd()) {
279 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
280 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
285 while (!it.IsAtEnd()) {
286 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
287 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
291 case 5: // Absolute value
292 while (!it.IsAtEnd()) {
293 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
294 // <= zero to avoid warning for unsigned types
295 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
299 case 6: // Squared value
300 while (!it.IsAtEnd()) {
301 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
306 while (!it.IsAtEnd()) {
308 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
309 else ito.Set(mDefaultPixelValue);
314 while (!it.IsAtEnd()) {
315 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
320 while (!it.IsAtEnd()) {
322 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
324 if (it.Get() ==0) ito.Set(0);
325 else ito.Set(mDefaultPixelValue);
331 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
335 //--------------------------------------------------------------------
339 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX