-//--------------------------------------------------------------------
-template<class ImageType1, class ImageType2>
-typename ImageType1::Pointer
-ImageArithmGenericFilter::ComputeImage(typename ImageType1::Pointer inputImage1,
- typename ImageType2::Pointer inputImage2) {
-
- typedef typename ImageType1::PixelType PixelType;
- typedef itk::ImageRegionIterator<ImageType1> IteratorType;
- IteratorType it1(inputImage1, inputImage1->GetLargestPossibleRegion());
- it1.GoToBegin();
-
- typedef itk::ImageRegionConstIterator<ImageType2> ConstIteratorType;
- ConstIteratorType it2(inputImage2, inputImage2->GetLargestPossibleRegion());
- it2.GoToBegin();
-
- switch (mTypeOfOperation) {
- case 0: // Addition
- while (!it1.IsAtEnd()) {
- it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
- ++it1; ++it2;
- }
- break;
- case 1: // Multiply
- while (!it1.IsAtEnd()) {
- it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
- ++it1; ++it2;
- }
- break;
- case 2: // Divide
- while (!it1.IsAtEnd()) {
- if (it1.Get() != 0)
- it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
- else it1.Set(mDefaultPixelValue);
- ++it1; ++it2;
- }
- break;
- case 3: // Max
- while (!it1.IsAtEnd()) {
- if (it1.Get() < it2.Get()) it1.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
- ++it1; ++it2;
- }
- break;
- case 4: // Min
- while (!it1.IsAtEnd()) {
- if (it1.Get() > it2.Get()) it1.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
- ++it1; ++it2;
- }
- break;
- case 5: // Absolute difference
- while (!it1.IsAtEnd()) {
- it1.Set(PixelTypeDownCast<double, PixelType>(fabs(it2.Get()-it1.Get())));
- ++it1; ++it2;
- }
- break;
- case 6: // Squared differences
- while (!it1.IsAtEnd()) {
- it1.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
- ++it1; ++it2;
- }
- break;
- case 7: // Difference
- while (!it1.IsAtEnd()) {
- it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
- ++it1; ++it2;
- }
- break;
- case 8: // Relative Difference
- while (!it1.IsAtEnd()) {
- if (it1.Get() != 0) it1.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
- else it1.Set(0.0);
- ++it1; ++it2;
- }
- break;
- default: // error ?
- std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
- exit(-1);
- }
-
- return inputImage1;
-}
-//--------------------------------------------------------------------