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);
49 //--------------------------------------------------------------------
52 //--------------------------------------------------------------------
53 template<class args_info_type>
54 void ImageArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b)
56 mOverwriteInputImage = b;
58 //--------------------------------------------------------------------
61 //--------------------------------------------------------------------
62 template<class args_info_type>
63 void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
68 this->SetIOVerbose(mArgsInfo.verbose_flag);
69 mTypeOfOperation = mArgsInfo.operation_arg;
70 mDefaultPixelValue = mArgsInfo.pixelValue_arg;
71 mScalar = mArgsInfo.scalar_arg;
72 mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
74 if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
76 if (mArgsInfo.input1_given) this->AddInputFilename(mArgsInfo.input1_arg);
77 if (mArgsInfo.input2_given) {
78 mIsOperationUseASecondImage = true;
79 this->AddInputFilename(mArgsInfo.input2_arg);
82 if (mArgsInfo.output_given) this->SetOutputFilename(mArgsInfo.output_arg);
84 // Check type of operation (with scalar or with other image)
85 if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
86 std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
89 if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
90 if (mArgsInfo.operation_arg < 5) {
91 std::cerr << "Such operation need the --scalar option." << std::endl;
96 //--------------------------------------------------------------------
99 //--------------------------------------------------------------------
100 template<class args_info_type>
101 template<class ImageType>
102 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
105 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
107 // Set input image iterator
108 typedef itk::ImageRegionIterator<ImageType> IteratorType;
109 IteratorType it(input1, input1->GetLargestPossibleRegion());
112 typename ImageType::Pointer input2 = NULL;
115 // Special case for normalisation
116 if (mTypeOfOperation == 12) {
117 typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
118 typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
119 ff->SetImage(input1);
120 ff->ComputeMaximum();
121 mScalar = ff->GetMaximum();
122 mTypeOfOperation = 11; // divide
125 if (mIsOperationUseASecondImage) {
127 input2 = this->template GetInput<ImageType>(1);
128 // Set input image iterator
129 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
131 if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
132 itkExceptionMacro(<< "The images (input and input2) must have the same size");
134 if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
135 itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
136 << "Using first input's information.");
140 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
141 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
142 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
143 // << typeid(PixelType).name()
145 mOverwriteInputImage = false;
148 // ---------------- Overwrite input Image ---------------------
149 if (mOverwriteInputImage) {
150 // Set output iterator (to input1)
151 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
152 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
153 else ComputeImage(it, ito);
154 this->template SetNextOutput<ImageType>(input1);
157 // ---------------- Create new output Image ---------------------
159 if (mOutputIsFloat) {
160 // Create output image
161 typedef itk::Image<float,ImageType::ImageDimension> 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);
174 // Create output image
175 typedef ImageType OutputImageType;
176 typename OutputImageType::Pointer output = OutputImageType::New();
177 output->SetRegions(input1->GetLargestPossibleRegion());
178 output->SetOrigin(input1->GetOrigin());
179 output->SetSpacing(input1->GetSpacing());
181 // Set output iterator
182 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
183 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
184 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
185 else ComputeImage(it, ito);
186 this->template SetNextOutput<OutputImageType>(output);
190 //--------------------------------------------------------------------
192 //--------------------------------------------------------------------
193 template<class args_info_type>
194 template<class Iter1, class Iter2, class Iter3>
195 void ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
200 typedef typename Iter3::PixelType PixelType;
202 switch (mTypeOfOperation) {
204 while (!ito.IsAtEnd()) {
205 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
212 while (!ito.IsAtEnd()) {
213 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
220 while (!ito.IsAtEnd()) {
222 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
223 else ito.Set(mDefaultPixelValue);
230 while (!ito.IsAtEnd()) {
231 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
238 while (!ito.IsAtEnd()) {
239 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
245 case 5: // Absolute difference
246 while (!ito.IsAtEnd()) {
247 ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
253 case 6: // Squared differences
254 while (!ito.IsAtEnd()) {
255 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
261 case 7: // Difference
262 while (!ito.IsAtEnd()) {
263 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
269 case 8: // Relative Difference
270 while (!ito.IsAtEnd()) {
271 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
279 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
283 //--------------------------------------------------------------------
286 //--------------------------------------------------------------------
287 template<class args_info_type>
288 template<class Iter1, class Iter2>
289 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
293 typedef typename Iter2::PixelType PixelType;
296 switch (mTypeOfOperation) {
298 while (!it.IsAtEnd()) {
299 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
305 while (!it.IsAtEnd()) {
306 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
312 while (!it.IsAtEnd()) {
314 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
315 else ito.Set(mDefaultPixelValue);
321 while (!it.IsAtEnd()) {
322 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
323 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
329 while (!it.IsAtEnd()) {
330 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
331 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
336 case 5: // Absolute value
337 while (!it.IsAtEnd()) {
338 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
339 // <= zero to avoid warning for unsigned types
340 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
345 case 6: // Squared value
346 while (!it.IsAtEnd()) {
347 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
353 while (!it.IsAtEnd()) {
355 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
356 else ito.Set(mDefaultPixelValue);
362 while (!it.IsAtEnd()) {
363 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
369 while (!it.IsAtEnd()) {
371 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
373 if (it.Get() ==0) ito.Set(0);
374 else ito.Set(mDefaultPixelValue);
381 while (!it.IsAtEnd()) {
382 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
388 while (!it.IsAtEnd()) {
389 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() / mScalar) );
395 while (!it.IsAtEnd()) {
396 if (it.Get() == 0) { // special case for fluence image with 0 value in a pixel -> consider 0.5
397 ito.Set(-log(0.5 / mScalar) );
400 ito.Set(PixelTypeDownCast<double, PixelType>(-log((double)it.Get() / mScalar)) );
407 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
411 //--------------------------------------------------------------------
417 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX