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 == 1)
83 mIsOutputScalar = true;
85 else if (mArgsInfo.operation_arg == 5 || mArgsInfo.operation_arg == 6)
86 mIsOutputScalar = true;
88 if (mArgsInfo.output_given) this->SetOutputFilename(mArgsInfo.output_arg);
90 // Check type of operation (with scalar or with other image)
91 if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
92 std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
95 if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
96 if (mArgsInfo.operation_arg < 5) {
97 std::cerr << "Such operation need the --scalar option." << std::endl;
102 //--------------------------------------------------------------------
105 //--------------------------------------------------------------------
106 template<class args_info_type>
107 template<class ImageType>
108 void VectorArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
111 typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
113 // Set input image iterator
114 typedef itk::ImageRegionIterator<ImageType> IteratorType;
115 IteratorType it(input1, input1->GetLargestPossibleRegion());
118 typename ImageType::Pointer input2 = NULL;
121 // Special case for normalisation
123 if (mTypeOfOperation == 12) {
124 typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
125 typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
126 ff->SetImage(input1);
127 ff->ComputeMaximum();
128 mScalar = ff->GetMaximum();
129 mTypeOfOperation = 11; // divide
133 if (mIsOperationUseASecondImage) {
135 input2 = this->template GetInput<ImageType>(1);
136 // Set input image iterator
137 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
139 if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
140 itkExceptionMacro(<< "The images (input and input2) must have the same size");
142 if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
143 itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
144 << "Using first input's information.");
148 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
149 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
150 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
151 // << typeid(PixelType).name()
153 mOverwriteInputImage = false;
156 // ---------------- Overwrite input Image ---------------------
157 if (mOverwriteInputImage && !mIsOutputScalar) {
158 // Set output iterator (to input1)
159 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
160 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
161 else ComputeImage(it, ito);
162 this->template SetNextOutput<ImageType>(input1);
165 // ---------------- Create new output Image ---------------------
167 // Create output image
168 if (!mIsOutputScalar) {
169 typedef ImageType OutputImageType;
170 typename OutputImageType::Pointer output = OutputImageType::New();
171 output->SetRegions(input1->GetLargestPossibleRegion());
172 output->SetOrigin(input1->GetOrigin());
173 output->SetSpacing(input1->GetSpacing());
175 // Set output iterator
176 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
177 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
178 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
179 else ComputeImage(it, ito);
180 this->template SetNextOutput<OutputImageType>(output);
183 // Create scalar output image
184 typedef itk::Image<typename ImageType::PixelType::ValueType, ImageType::ImageDimension> OutputImageType;
185 typename OutputImageType::Pointer output = OutputImageType::New();
186 output->SetRegions(input1->GetLargestPossibleRegion());
187 output->SetOrigin(input1->GetOrigin());
188 output->SetSpacing(input1->GetSpacing());
190 // Set output iterator
191 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
192 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
193 if (mIsOperationUseASecondImage) ComputeScalarImage(it, it2, ito);
194 else ComputeScalarImage(it, ito);
195 this->template SetNextOutput<OutputImageType>(output);
199 //--------------------------------------------------------------------
201 //--------------------------------------------------------------------
202 template<class args_info_type>
203 template<class Iter1, class Iter2, class Iter3>
204 void VectorArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
209 typedef typename Iter3::PixelType PixelType;
211 switch (mTypeOfOperation) {
213 while (!ito.IsAtEnd()) {
214 ito.Set(it1.Get() + it2.Get());
222 while (!ito.IsAtEnd()) {
223 ito.Set(it1.Get() * it2.Get()) );
232 while (!ito.IsAtEnd()) {
234 ito.Set(it1.Get() / it2.Get()));
235 else ito.Set(mDefaultPixelValue);
242 while (!ito.IsAtEnd()) {
243 if (it1.Get() < it2.Get()) ito.Set(it2.Get());
250 while (!ito.IsAtEnd()) {
251 if (it1.Get() > it2.Get()) ito.Set(it2.Get());
259 case 5: // Absolute difference
260 while (!ito.IsAtEnd()) {
261 ito.Set(it2.Get()-it1.Get());
269 case 6: // Squared differences
270 while (!ito.IsAtEnd()) {
271 ito.Set(pow(it1.Get()-it2.Get(),2)));
278 case 7: // Difference
279 while (!ito.IsAtEnd()) {
280 ito.Set(it1.Get()-it2.Get());
287 case 8: // Relative Difference
288 while (!ito.IsAtEnd()) {
289 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
298 case 9: // CrossProduct
299 while (!ito.IsAtEnd()) {
300 ito.Set(it1.Get()^it2.Get());
308 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
312 //--------------------------------------------------------------------
315 //--------------------------------------------------------------------
316 template<class args_info_type>
317 template<class Iter1, class Iter2>
318 void clitk::VectorArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
323 typedef typename Iter1::PixelType PixelType;
325 PixelType scalar_vector;
326 scalar_vector.Fill(mScalar);
329 switch (mTypeOfOperation) {
331 while (!it.IsAtEnd()) {
332 ito.Set(it.Get() + scalar_vector);
338 while (!it.IsAtEnd()) {
339 ito.Set(it.Get() * mScalar);
346 while (!it.IsAtEnd()) {
348 ito.Set(mScalar / it.Get()));
349 else ito.Set(mDefaultPixelValue);
355 while (!it.IsAtEnd()) {
356 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
357 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
363 while (!it.IsAtEnd()) {
364 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
365 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
372 case 5: // Absolute value
373 while (!it.IsAtEnd()) {
374 ito.Set(PixelTypeDownCast<double, PixelType>(it.GetNorm()));
379 case 6: // Squared value
380 while (!it.IsAtEnd()) {
381 ito.Set(PixelTypeDownCast<double, PixelType>(it.GetSquaredNorm());
389 while (!it.IsAtEnd()) {
391 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
392 else ito.Set(mDefaultPixelValue);
398 while (!it.IsAtEnd()) {
399 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
405 while (!it.IsAtEnd()) {
407 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
409 if (it.Get() ==0) ito.Set(0);
410 else ito.Set(mDefaultPixelValue);
417 while (!it.IsAtEnd()) {
418 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
425 while (!it.IsAtEnd()) {
426 ito.Set(it.Get() / mScalar);
431 case 12: // normalize
432 while (!it.IsAtEnd()) {
433 PixelType n = it.Get();
441 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
445 //--------------------------------------------------------------------
447 //--------------------------------------------------------------------
448 template<class args_info_type>
449 template<class Iter1, class Iter2, class Iter3>
450 void VectorArithmGenericFilter<args_info_type>::ComputeScalarImage(Iter1 it1, Iter2 it2, Iter3 ito)
455 typedef typename Iter3::PixelType PixelType;
457 switch (mTypeOfOperation) {
459 while (!ito.IsAtEnd()) {
460 ito.Set(it1.Get() * it2.Get());
467 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
471 //--------------------------------------------------------------------
473 //--------------------------------------------------------------------
474 template<class args_info_type>
475 template<class Iter1, class Iter2>
476 void clitk::VectorArithmGenericFilter<args_info_type>::ComputeScalarImage(Iter1 it, Iter2 ito)
481 typedef typename Iter2::PixelType PixelType;
485 switch (mTypeOfOperation) {
486 case 5: // Absolute value
487 while (!it.IsAtEnd()) {
488 ito.Set(PixelTypeDownCast<double, PixelType>(it.Get().GetNorm()));
493 case 6: // Squared value
494 while (!it.IsAtEnd()) {
495 ito.Set(PixelTypeDownCast<double, PixelType>(it.Get().GetSquaredNorm()));
501 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
505 //--------------------------------------------------------------------
510 #endif //#define CLITKVectorArithmGENERICFILTER_TXX