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 CLITKVectorArithmGENERICFILTER_TXX
19 #define CLITKVectorArithmGENERICFILTER_TXX
21 #include "clitkImageCommon.h"
23 #include "itkMinimumMaximumImageCalculator.h"
28 //--------------------------------------------------------------------
29 template<class args_info_type>
30 VectorArithmGenericFilter<args_info_type>::VectorArithmGenericFilter()
31 :ImageToImageGenericFilter<Self>("VectorArithmGenericFilter"),mTypeOfOperation(0)
33 InitializeImageType<3>();
34 mIsOperationUseASecondImage = false;
35 mIsOutputScalar = false;
36 mOverwriteInputImage = true;
38 //--------------------------------------------------------------------
41 //--------------------------------------------------------------------
42 template<class args_info_type>
43 template<unsigned int Dim>
44 void VectorArithmGenericFilter<args_info_type>::InitializeImageType()
46 ADD_VEC_IMAGE_TYPE(Dim,3u,float);
47 ADD_VEC_IMAGE_TYPE(Dim,3u,double);
48 ADD_VEC_IMAGE_TYPE(Dim,2u,float);
49 ADD_VEC_IMAGE_TYPE(Dim,2u,double);
51 //--------------------------------------------------------------------
54 //--------------------------------------------------------------------
55 template<class args_info_type>
56 void VectorArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b)
58 mOverwriteInputImage = b;
60 //--------------------------------------------------------------------
63 //--------------------------------------------------------------------
64 template<class args_info_type>
65 void VectorArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
70 this->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) this->AddInputFilename(mArgsInfo.input1_arg);
79 if (mArgsInfo.input2_given) {
80 mIsOperationUseASecondImage = true;
81 this->AddInputFilename(mArgsInfo.input2_arg);
82 if (mArgsInfo.operation_arg == 2)
83 mIsOutputScalar = true;
85 else if (mArgsInfo.operation_arg == 5 || mArgsInfo.operation_arg == 6)
86 mIsOutputScalar = true;
88 if (mArgsInfo.output_given) {
89 this->SetOutputFilename(mArgsInfo.output_arg);
90 mOverwriteInputImage = false;
93 // Check type of operation (with scalar or with other image)
94 if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
95 std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
98 if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
99 if (mArgsInfo.operation_arg < 5) {
100 std::cerr << "Such operation need the --scalar option." << std::endl;
105 //--------------------------------------------------------------------
108 //--------------------------------------------------------------------
109 template<class args_info_type>
110 template<class ImageType>
111 void VectorArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
114 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
116 // Set input image iterator
117 typedef itk::ImageRegionIterator<ImageType> IteratorType;
118 IteratorType it(input1, input1->GetLargestPossibleRegion());
121 typename ImageType::Pointer input2 = ITK_NULLPTR;
124 // Special case for normalisation
126 if (mTypeOfOperation == 12) {
127 typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
128 typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
129 ff->SetImage(input1);
130 ff->ComputeMaximum();
131 mScalar = ff->GetMaximum();
132 mTypeOfOperation = 11; // divide
136 if (mIsOperationUseASecondImage) {
138 input2 = this->template GetInput<ImageType>(1);
139 // Set input image iterator
140 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
142 if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
143 itkExceptionMacro(<< "The images (input and input2) must have the same size");
145 if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
146 itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
147 << "Using first input's information.");
151 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
152 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
153 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
154 // << typeid(PixelType).name()
156 mOverwriteInputImage = false;
159 // ---------------- Overwrite input Image ---------------------
160 if (mOverwriteInputImage && !mIsOutputScalar) {
161 // Set output iterator (to input1)
162 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
163 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
164 else ComputeImage(it, ito);
165 this->template SetNextOutput<ImageType>(input1);
168 // ---------------- Create new output Image ---------------------
170 // Create output image
171 if (!mIsOutputScalar) {
172 typedef ImageType OutputImageType;
173 typename OutputImageType::Pointer output = OutputImageType::New();
174 output->SetRegions(input1->GetLargestPossibleRegion());
175 output->SetOrigin(input1->GetOrigin());
176 output->SetSpacing(input1->GetSpacing());
178 // Set output iterator
179 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
180 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
181 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
182 else ComputeImage(it, ito);
183 this->template SetNextOutput<OutputImageType>(output);
186 // Create scalar output image
187 typedef itk::Image<typename ImageType::PixelType::ValueType, ImageType::ImageDimension> OutputImageType;
188 typename OutputImageType::Pointer output = OutputImageType::New();
189 output->SetRegions(input1->GetLargestPossibleRegion());
190 output->SetOrigin(input1->GetOrigin());
191 output->SetSpacing(input1->GetSpacing());
193 // Set output iterator
194 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
195 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
196 if (mIsOperationUseASecondImage) ComputeScalarImage(it, it2, ito);
197 else ComputeScalarImage(it, ito);
198 this->template SetNextOutput<OutputImageType>(output);
202 //--------------------------------------------------------------------
204 //--------------------------------------------------------------------
205 template<class args_info_type>
206 template<class Iter1, class Iter2, class Iter3>
207 void VectorArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
212 typedef typename Iter3::PixelType PixelType;
214 switch (mTypeOfOperation) {
216 while (!ito.IsAtEnd()) {
217 ito.Set(it1.Get() + it2.Get());
224 case 1: // term to term Multiply
225 while (!ito.IsAtEnd()) {
226 typename Iter1::PixelType outputPixel(ito.Get());
227 outputPixel.SetVnlVector(element_product(it1.Get().GetVnlVector(),it2.Get().GetVnlVector()));
228 ito.Set(outputPixel);
237 while (!ito.IsAtEnd()) {
239 ito.Set(it1.Get() / it2.Get()));
240 else ito.Set(mDefaultPixelValue);
247 while (!ito.IsAtEnd()) {
248 if (it1.Get() < it2.Get()) ito.Set(it2.Get());
255 while (!ito.IsAtEnd()) {
256 if (it1.Get() > it2.Get()) ito.Set(it2.Get());
264 case 5: // Absolute difference
265 while (!ito.IsAtEnd()) {
266 ito.Set(it2.Get()-it1.Get());
274 case 6: // Squared differences
275 while (!ito.IsAtEnd()) {
276 ito.Set(pow(it1.Get()-it2.Get(),2)));
283 case 7: // Difference
284 while (!ito.IsAtEnd()) {
285 ito.Set(it1.Get()-it2.Get());
292 case 8: // Relative Difference
293 while (!ito.IsAtEnd()) {
294 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
303 case 9: // CrossProduct
304 while (!ito.IsAtEnd()) {
305 ito.Set(it1.Get()^it2.Get());
313 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
317 //--------------------------------------------------------------------
320 //--------------------------------------------------------------------
321 template<class args_info_type>
322 template<class Iter1, class Iter2>
323 void clitk::VectorArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
328 typedef typename Iter1::PixelType PixelType;
330 PixelType scalar_vector;
331 scalar_vector.Fill(mScalar);
334 switch (mTypeOfOperation) {
336 while (!it.IsAtEnd()) {
337 ito.Set(it.Get() + scalar_vector);
343 while (!it.IsAtEnd()) {
344 ito.Set(it.Get() * mScalar);
351 while (!it.IsAtEnd()) {
353 ito.Set(mScalar / it.Get()));
354 else ito.Set(mDefaultPixelValue);
360 while (!it.IsAtEnd()) {
361 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
362 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
368 while (!it.IsAtEnd()) {
369 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
370 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
377 case 5: // Absolute value
378 while (!it.IsAtEnd()) {
379 ito.Set(PixelTypeDownCast<double, PixelType>(it.GetNorm()));
384 case 6: // Squared value
385 while (!it.IsAtEnd()) {
386 ito.Set(PixelTypeDownCast<double, PixelType>(it.GetSquaredNorm());
394 while (!it.IsAtEnd()) {
396 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
397 else ito.Set(mDefaultPixelValue);
403 while (!it.IsAtEnd()) {
404 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
410 while (!it.IsAtEnd()) {
412 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
414 if (it.Get() ==0) ito.Set(0);
415 else ito.Set(mDefaultPixelValue);
422 while (!it.IsAtEnd()) {
423 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
430 while (!it.IsAtEnd()) {
431 ito.Set(it.Get() / mScalar);
436 case 12: // normalize
437 while (!it.IsAtEnd()) {
438 PixelType n = it.Get();
439 if (n.GetNorm() != 0)
448 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
452 //--------------------------------------------------------------------
454 //--------------------------------------------------------------------
455 template<class args_info_type>
456 template<class Iter1, class Iter2, class Iter3>
457 void VectorArithmGenericFilter<args_info_type>::ComputeScalarImage(Iter1 it1, Iter2 it2, Iter3 ito)
462 typedef typename Iter3::PixelType PixelType;
464 switch (mTypeOfOperation) {
466 while (!ito.IsAtEnd()) {
467 ito.Set(it1.Get() * it2.Get());
474 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
478 //--------------------------------------------------------------------
480 //--------------------------------------------------------------------
481 template<class args_info_type>
482 template<class Iter1, class Iter2>
483 void clitk::VectorArithmGenericFilter<args_info_type>::ComputeScalarImage(Iter1 it, Iter2 ito)
488 typedef typename Iter2::PixelType PixelType;
492 switch (mTypeOfOperation) {
493 case 5: // Absolute value
494 while (!it.IsAtEnd()) {
495 ito.Set(PixelTypeDownCast<double, PixelType>(it.Get().GetNorm()));
500 case 6: // Squared value
501 while (!it.IsAtEnd()) {
502 ito.Set(PixelTypeDownCast<double, PixelType>(it.Get().GetSquaredNorm()));
508 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
512 //--------------------------------------------------------------------
517 #endif //#define CLITKVectorArithmGENERICFILTER_TXX