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 = ITK_NULLPTR;
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());
166 output->SetDirection(input1->GetDirection());
168 // Set output iterator
169 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
170 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
171 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
172 else ComputeImage(it, ito);
173 this->template SetNextOutput<OutputImageType>(output);
175 // Create output image
176 typedef ImageType OutputImageType;
177 typename OutputImageType::Pointer output = OutputImageType::New();
178 output->SetRegions(input1->GetLargestPossibleRegion());
179 output->SetOrigin(input1->GetOrigin());
180 output->SetSpacing(input1->GetSpacing());
182 // Set output iterator
183 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
184 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
185 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
186 else ComputeImage(it, ito);
187 this->template SetNextOutput<OutputImageType>(output);
191 //--------------------------------------------------------------------
193 //--------------------------------------------------------------------
194 template<class args_info_type>
195 template<class Iter1, class Iter2, class Iter3>
196 void ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
201 typedef typename Iter3::PixelType PixelType;
203 switch (mTypeOfOperation) {
205 while (!ito.IsAtEnd()) {
206 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
213 while (!ito.IsAtEnd()) {
214 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
221 while (!ito.IsAtEnd()) {
223 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
224 else ito.Set(mDefaultPixelValue);
231 while (!ito.IsAtEnd()) {
232 if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
239 while (!ito.IsAtEnd()) {
240 if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
246 case 5: // Absolute difference
247 while (!ito.IsAtEnd()) {
248 ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
254 case 6: // Squared differences
255 while (!ito.IsAtEnd()) {
256 ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
262 case 7: // Difference
263 while (!ito.IsAtEnd()) {
264 ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
270 case 8: // Relative Difference
271 while (!ito.IsAtEnd()) {
272 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
273 else ito.Set(mDefaultPixelValue);
280 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
284 //--------------------------------------------------------------------
287 //--------------------------------------------------------------------
288 template<class args_info_type>
289 template<class Iter1, class Iter2>
290 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
294 typedef typename Iter2::PixelType PixelType;
297 switch (mTypeOfOperation) {
299 while (!it.IsAtEnd()) {
300 ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
306 while (!it.IsAtEnd()) {
307 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
313 while (!it.IsAtEnd()) {
315 ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
316 else ito.Set(mDefaultPixelValue);
322 while (!it.IsAtEnd()) {
323 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
324 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
330 while (!it.IsAtEnd()) {
331 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
332 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
337 case 5: // Absolute value
338 while (!it.IsAtEnd()) {
339 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
340 // <= zero to avoid warning for unsigned types
341 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
346 case 6: // Squared value
347 while (!it.IsAtEnd()) {
348 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
354 while (!it.IsAtEnd()) {
356 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
357 else ito.Set(mDefaultPixelValue);
363 while (!it.IsAtEnd()) {
364 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
370 while (!it.IsAtEnd()) {
372 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
374 if (it.Get() ==0) ito.Set(0);
375 else ito.Set(mDefaultPixelValue);
382 while (!it.IsAtEnd()) {
383 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
389 while (!it.IsAtEnd()) {
390 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() / mScalar) );
396 while (!it.IsAtEnd()) {
397 if (it.Get() == 0) { // special case for fluence image with 0 value in a pixel -> consider 0.5
398 ito.Set(-log(0.5 / mScalar) );
401 ito.Set(PixelTypeDownCast<double, PixelType>(-log((double)it.Get() / mScalar)) );
408 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
412 //--------------------------------------------------------------------
418 #endif //#define CLITKIMAGEARITHMGENERICFILTER_TXX