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, float);
45 //--------------------------------------------------------------------
48 //--------------------------------------------------------------------
49 template<class args_info_type>
50 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) {
62 SetIOVerbose(mArgsInfo.verbose_flag);
63 mTypeOfOperation = mArgsInfo.operation_arg;
64 mDefaultPixelValue = mArgsInfo.pixelValue_arg;
65 mScalar = mArgsInfo.scalar_arg;
66 mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
68 if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
70 if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
71 if (mArgsInfo.input2_given) {
72 mIsOperationUseASecondImage = true;
73 AddInputFilename(mArgsInfo.input2_arg);
76 if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
78 // Check type of operation (with scalar or with other image)
79 if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
80 std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
83 if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
84 if (mArgsInfo.operation_arg < 5) {
85 std::cerr << "Such operation need the --scalar option." << std::endl;
90 //--------------------------------------------------------------------
93 //--------------------------------------------------------------------
94 template<class args_info_type>
95 template<class ImageType>
96 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType() {
98 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
99 typename ImageType::PixelType PixelType;
101 // Set input image iterator
102 typedef itk::ImageRegionIterator<ImageType> IteratorType;
103 IteratorType it(input1, input1->GetLargestPossibleRegion());
106 typename ImageType::Pointer input2 = this->template GetInput<ImageType>(1);
109 if (mIsOperationUseASecondImage) {
111 input2 = this->template GetInput<ImageType>(1);
112 // Set input image iterator
113 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
116 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
117 if (mOverwriteInputImage && mOutputIsFloat && (typeid(PixelType) != typeid(float))) {
118 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
119 // << typeid(PixelType).name()
121 mOverwriteInputImage = false;
124 // ---------------- Overwrite input Image ---------------------
125 if (mOverwriteInputImage) {
126 // Set output iterator (to input1)
127 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
128 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
129 else ComputeImage(it, ito);
130 this->template SetNextOutput<ImageType>(input1);
133 // ---------------- Create new output Image ---------------------
135 if (mOutputIsFloat) {
136 // Create output image
137 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
138 typename OutputImageType::Pointer output = OutputImageType::New();
139 output->SetRegions(input1->GetLargestPossibleRegion());
140 output->SetOrigin(input1->GetOrigin());
141 output->SetSpacing(input1->GetSpacing());
143 // Set output iterator
144 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
145 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
146 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
147 else ComputeImage(it, ito);
148 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) {
176 typedef typename Iter3::PixelType PixelType;
178 switch (mTypeOfOperation) {
180 while (!ito.IsAtEnd()) {
181 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
186 while (!ito.IsAtEnd()) {
187 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
192 while (!ito.IsAtEnd()) {
194 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
195 else ito.Set(mDefaultPixelValue);
200 while (!ito.IsAtEnd()) {
201 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
206 while (!ito.IsAtEnd()) {
207 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
211 case 5: // Absolute difference
212 while (!ito.IsAtEnd()) {
213 ito.Set(PixelTypeDownCast<double, PixelType>(fabs(it2.Get()-it1.Get())));
217 case 6: // Squared differences
218 while (!ito.IsAtEnd()) {
219 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
223 case 7: // Difference
224 while (!ito.IsAtEnd()) {
225 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
229 case 8: // Relative Difference
230 while (!ito.IsAtEnd()) {
231 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
237 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
241 //--------------------------------------------------------------------
244 //--------------------------------------------------------------------
245 template<class args_info_type>
246 template<class Iter1, class Iter2>
247 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito) {
250 typedef typename Iter2::PixelType PixelType;
253 switch (mTypeOfOperation) {
255 while (!it.IsAtEnd()) {
256 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
261 while (!it.IsAtEnd()) {
262 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
267 while (!it.IsAtEnd()) {
269 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
270 else ito.Set(mDefaultPixelValue);
275 while (!it.IsAtEnd()) {
276 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
277 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
282 while (!it.IsAtEnd()) {
283 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
284 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
288 case 5: // Absolute value
289 while (!it.IsAtEnd()) {
290 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
291 // <= zero to avoid warning for unsigned types
292 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
296 case 6: // Squared value
297 while (!it.IsAtEnd()) {
298 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
303 while (!it.IsAtEnd()) {
305 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
306 else ito.Set(mDefaultPixelValue);
311 while (!it.IsAtEnd()) {
312 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
317 while (!it.IsAtEnd()) {
319 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
321 if (it.Get() ==0) ito.Set(0);
322 else ito.Set(mDefaultPixelValue);
328 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
332 //--------------------------------------------------------------------
336 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX