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_CXX
19 #define CLITKIMAGEARITHMGENERICFILTER_CXX
21 #include "clitkImageArithmGenericFilter.h"
26 // class ImageArithmGenericFilter<args_info_clitkImageArithm>;
28 ////////////////////////////////////////////////////////////////////
29 // specializations for itk::Vector<float, 3u>, 3u
32 void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >()
34 typedef itk::Image< itk::Vector<float, 3u>, 3u > ImageType;
37 ImageType::Pointer input1 = this->GetInput<ImageType>(0);
39 // Set input image iterator
40 typedef itk::ImageRegionIterator<ImageType> IteratorType;
41 IteratorType it(input1, input1->GetLargestPossibleRegion());
44 ImageType::Pointer input2 = NULL;
48 // Special case for normalisation
49 if (mTypeOfOperation == 12) {
50 typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
51 typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
54 mScalar = ff->GetMaximum();
55 mTypeOfOperation = 11; // divide
59 if (mIsOperationUseASecondImage) {
61 input2 = this->GetInput<ImageType>(1);
62 // Set input image iterator
63 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
65 if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
66 itkExceptionMacro(<< "The images (input and input2) must have the same size");
68 if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
69 itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
70 << "Using first input's information.");
75 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
76 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
77 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
78 // << typeid(PixelType).name()
80 mOverwriteInputImage = false;
84 // ---------------- Overwrite input Image ---------------------
85 if (mOverwriteInputImage) {
86 // Set output iterator (to input1)
87 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
88 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
89 else ComputeImage(it, ito);
90 this->SetNextOutput<ImageType>(input1);
92 // ---------------- Create new output Image ---------------------
94 /*if (mOutputIsFloat) {
95 // Create output image
96 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
97 typename OutputImageType::Pointer output = OutputImageType::New();
98 output->SetRegions(input1->GetLargestPossibleRegion());
99 output->SetOrigin(input1->GetOrigin());
100 output->SetSpacing(input1->GetSpacing());
102 // Set output iterator
103 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
104 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
105 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
106 else ComputeImage(it, ito);
107 this->template SetNextOutput<OutputImageType>(output);
109 // Create output image
110 typedef ImageType OutputImageType;
111 OutputImageType::Pointer output = OutputImageType::New();
112 output->SetRegions(input1->GetLargestPossibleRegion());
113 output->SetOrigin(input1->GetOrigin());
114 output->SetSpacing(input1->GetSpacing());
116 // Set output iterator
117 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
118 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
119 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
120 else ComputeImage(it, ito);
121 this->SetNextOutput<OutputImageType>(output);
128 void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
129 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
130 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > >
131 (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it,
132 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
134 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
139 typedef Iter2::PixelType PixelType;
141 PixelType scalar_vector;
142 scalar_vector.Fill(mScalar);
145 switch (mTypeOfOperation) {
147 while (!it.IsAtEnd()) {
148 ito.Set(it.Get() + scalar_vector);
154 while (!it.IsAtEnd()) {
155 ito.Set(it.Get() * mScalar);
162 while (!it.IsAtEnd()) {
164 ito.Set(mScalar / it.Get()));
165 else ito.Set(mDefaultPixelValue);
171 while (!it.IsAtEnd()) {
172 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
173 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
179 while (!it.IsAtEnd()) {
180 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
181 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
186 case 5: // Absolute value
187 while (!it.IsAtEnd()) {
188 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
189 // <= zero to avoid warning for unsigned types
190 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
195 case 6: // Squared value
196 while (!it.IsAtEnd()) {
197 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
203 while (!it.IsAtEnd()) {
205 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
206 else ito.Set(mDefaultPixelValue);
212 while (!it.IsAtEnd()) {
213 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
219 while (!it.IsAtEnd()) {
221 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
223 if (it.Get() ==0) ito.Set(0);
224 else ito.Set(mDefaultPixelValue);
231 while (!it.IsAtEnd()) {
232 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
239 while (!it.IsAtEnd()) {
240 ito.Set(it.Get() / mScalar);
246 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
254 void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
255 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
256 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
257 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >
259 (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it1,
260 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2,
261 itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
263 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter1;
264 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
265 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter3;
270 typedef Iter3::PixelType PixelType;
272 switch (mTypeOfOperation) {
274 while (!ito.IsAtEnd()) {
275 ito.Set(it1.Get() + it2.Get());
283 while (!ito.IsAtEnd()) {
284 ito.Set(it1.Get() * it2.Get()) );
291 while (!ito.IsAtEnd()) {
293 ito.Set(it1.Get() / it2.Get()));
294 else ito.Set(mDefaultPixelValue);
301 while (!ito.IsAtEnd()) {
302 if (it1.Get() < it2.Get()) ito.Set(it2.Get());
309 while (!ito.IsAtEnd()) {
310 if (it1.Get() > it2.Get()) ito.Set(it2.Get());
317 case 5: // Absolute difference
318 while (!ito.IsAtEnd()) {
319 ito.Set(it2.Get()-it1.Get());
326 case 6: // Squared differences
327 while (!ito.IsAtEnd()) {
328 ito.Set(pow(it1.Get()-it2.Get(),2)));
335 case 7: // Difference
336 while (!ito.IsAtEnd()) {
337 ito.Set(it1.Get()-it2.Get());
344 case 8: // Relative Difference
345 while (!ito.IsAtEnd()) {
346 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
355 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
360 ////////////////////////////////////////////////////////////////////
361 // specializations for itk::Vector<double, 3u>, 3u
364 void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<double, 3u>, 3u > >()
366 typedef itk::Image< itk::Vector<double, 3u>, 3u > ImageType;
369 ImageType::Pointer input1 = this->GetInput<ImageType>(0);
371 // Set input image iterator
372 typedef itk::ImageRegionIterator<ImageType> IteratorType;
373 IteratorType it(input1, input1->GetLargestPossibleRegion());
376 ImageType::Pointer input2 = NULL;
380 // Special case for normalisation
381 if (mTypeOfOperation == 12) {
382 typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
383 typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
384 ff->SetImage(input1);
385 ff->ComputeMaximum();
386 mScalar = ff->GetMaximum();
387 mTypeOfOperation = 11; // divide
391 if (mIsOperationUseASecondImage) {
393 input2 = this->GetInput<ImageType>(1);
394 // Set input image iterator
395 it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
397 if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
398 itkExceptionMacro(<< "The images (input and input2) must have the same size");
400 if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
401 itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
402 << "Using first input's information.");
407 // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
408 if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
409 // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
410 // << typeid(PixelType).name()
412 mOverwriteInputImage = false;
416 // ---------------- Overwrite input Image ---------------------
417 if (mOverwriteInputImage) {
418 // Set output iterator (to input1)
419 IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
420 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
421 else ComputeImage(it, ito);
422 this->SetNextOutput<ImageType>(input1);
424 // ---------------- Create new output Image ---------------------
426 /*if (mOutputIsFloat) {
427 // Create output image
428 typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
429 typename OutputImageType::Pointer output = OutputImageType::New();
430 output->SetRegions(input1->GetLargestPossibleRegion());
431 output->SetOrigin(input1->GetOrigin());
432 output->SetSpacing(input1->GetSpacing());
434 // Set output iterator
435 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
436 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
437 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
438 else ComputeImage(it, ito);
439 this->template SetNextOutput<OutputImageType>(output);
441 // Create output image
442 typedef ImageType OutputImageType;
443 OutputImageType::Pointer output = OutputImageType::New();
444 output->SetRegions(input1->GetLargestPossibleRegion());
445 output->SetOrigin(input1->GetOrigin());
446 output->SetSpacing(input1->GetSpacing());
448 // Set output iterator
449 typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
450 IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
451 if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
452 else ComputeImage(it, ito);
453 this->SetNextOutput<OutputImageType>(output);
460 void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
461 itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
462 itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > >
463 (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it,
464 itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito)
466 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter2;
471 typedef Iter2::PixelType PixelType;
473 PixelType scalar_vector;
474 scalar_vector.Fill(mScalar);
477 switch (mTypeOfOperation) {
479 while (!it.IsAtEnd()) {
480 ito.Set(it.Get() + scalar_vector);
486 while (!it.IsAtEnd()) {
487 ito.Set(it.Get() * mScalar);
494 while (!it.IsAtEnd()) {
496 ito.Set(mScalar / it.Get()));
497 else ito.Set(mDefaultPixelValue);
503 while (!it.IsAtEnd()) {
504 if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
505 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
511 while (!it.IsAtEnd()) {
512 if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
513 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
518 case 5: // Absolute value
519 while (!it.IsAtEnd()) {
520 if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
521 // <= zero to avoid warning for unsigned types
522 else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
527 case 6: // Squared value
528 while (!it.IsAtEnd()) {
529 ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
535 while (!it.IsAtEnd()) {
537 ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
538 else ito.Set(mDefaultPixelValue);
544 while (!it.IsAtEnd()) {
545 ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
551 while (!it.IsAtEnd()) {
553 ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
555 if (it.Get() ==0) ito.Set(0);
556 else ito.Set(mDefaultPixelValue);
563 while (!it.IsAtEnd()) {
564 ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
571 while (!it.IsAtEnd()) {
572 ito.Set(it.Get() / mScalar);
578 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
586 void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
587 itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
588 itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
589 itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >
591 (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it1,
592 itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it2,
593 itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito)
595 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter1;
596 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter2;
597 typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter3;
602 typedef Iter3::PixelType PixelType;
604 switch (mTypeOfOperation) {
606 while (!ito.IsAtEnd()) {
607 ito.Set(it1.Get() + it2.Get());
615 while (!ito.IsAtEnd()) {
616 ito.Set(it1.Get() * it2.Get()) );
623 while (!ito.IsAtEnd()) {
625 ito.Set(it1.Get() / it2.Get()));
626 else ito.Set(mDefaultPixelValue);
633 while (!ito.IsAtEnd()) {
634 if (it1.Get() < it2.Get()) ito.Set(it2.Get());
641 while (!ito.IsAtEnd()) {
642 if (it1.Get() > it2.Get()) ito.Set(it2.Get());
649 case 5: // Absolute difference
650 while (!ito.IsAtEnd()) {
651 ito.Set(it2.Get()-it1.Get());
658 case 6: // Squared differences
659 while (!ito.IsAtEnd()) {
660 ito.Set(pow(it1.Get()-it2.Get(),2)));
667 case 7: // Difference
668 while (!ito.IsAtEnd()) {
669 ito.Set(it1.Get()-it2.Get());
676 case 8: // Relative Difference
677 while (!ito.IsAtEnd()) {
678 if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
687 std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
694 #endif //define CLITKIMAGEARITHMGENERICFILTER_CXX