From a0209d730acea6ba9f39a93e04632353f7f1401f Mon Sep 17 00:00:00 2001 From: Romulo Pinho Date: Wed, 22 Feb 2012 10:48:26 +0100 Subject: [PATCH] added support to new vector type arithmetics --- tools/clitkImageArithmGenericFilter.cxx | 334 ++++++++++++++++++++++++ tools/clitkImageArithmGenericFilter.h | 22 ++ tools/clitkImageArithmGenericFilter.txx | 1 + 3 files changed, 357 insertions(+) diff --git a/tools/clitkImageArithmGenericFilter.cxx b/tools/clitkImageArithmGenericFilter.cxx index 71da90e..f0d55c2 100644 --- a/tools/clitkImageArithmGenericFilter.cxx +++ b/tools/clitkImageArithmGenericFilter.cxx @@ -25,6 +25,8 @@ namespace clitk { // template<> // class ImageArithmGenericFilter; +//////////////////////////////////////////////////////////////////// + // specializations for itk::Vector, 3u template<> template<> void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >() @@ -355,6 +357,338 @@ namespace clitk { } } +//////////////////////////////////////////////////////////////////// + // specializations for itk::Vector, 3u + template<> + template<> + void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >() + { + typedef itk::Image< itk::Vector, 3u > ImageType; + + // Read input1 + ImageType::Pointer input1 = this->GetInput(0); + + // Set input image iterator + typedef itk::ImageRegionIterator IteratorType; + IteratorType it(input1, input1->GetLargestPossibleRegion()); + + // typedef input2 + ImageType::Pointer input2 = NULL; + IteratorType it2; + + /* + // Special case for normalisation + if (mTypeOfOperation == 12) { + typedef itk::MinimumMaximumImageCalculator MinMaxFilterType; + typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New(); + ff->SetImage(input1); + ff->ComputeMaximum(); + mScalar = ff->GetMaximum(); + mTypeOfOperation = 11; // divide + } + */ + + if (mIsOperationUseASecondImage) { + // Read input2 + input2 = this->GetInput(1); + // Set input image iterator + it2 = IteratorType(input2, input2->GetLargestPossibleRegion()); + // Check dimension + if (!clitk::HaveSameSize(input1, input2)) { + itkExceptionMacro(<< "The images (input and input2) must have the same size"); + } + if(!clitk::HaveSameSpacing(input1, input2)) { + itkWarningMacro(<< "The images (input and input2) do not have the same spacing. " + << "Using first input's information."); + } + } + + /* + // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite + if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) { + // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is " + // << typeid(PixelType).name() + // << std::endl; + mOverwriteInputImage = false; + } + */ + + // ---------------- Overwrite input Image --------------------- + if (mOverwriteInputImage) { + // Set output iterator (to input1) + IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->SetNextOutput(input1); + } + // ---------------- Create new output Image --------------------- + else { + /*if (mOutputIsFloat) { + // Create output image + typedef itk::Image OutputImageType; + typename OutputImageType::Pointer output = OutputImageType::New(); + output->SetRegions(input1->GetLargestPossibleRegion()); + output->SetOrigin(input1->GetOrigin()); + output->SetSpacing(input1->GetSpacing()); + output->Allocate(); + // Set output iterator + typedef itk::ImageRegionIterator IteratorOutputType; + IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->template SetNextOutput(output); + } else*/ { + // Create output image + typedef ImageType OutputImageType; + OutputImageType::Pointer output = OutputImageType::New(); + output->SetRegions(input1->GetLargestPossibleRegion()); + output->SetOrigin(input1->GetOrigin()); + output->SetSpacing(input1->GetSpacing()); + output->Allocate(); + // Set output iterator + typedef itk::ImageRegionIterator IteratorOutputType; + IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion()); + if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito); + else ComputeImage(it, ito); + this->SetNextOutput(output); + } + } + } + + template<> + template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito) + { + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter2; + + ito.GoToBegin(); + it.GoToBegin(); + + typedef Iter2::PixelType PixelType; + + PixelType scalar_vector; + scalar_vector.Fill(mScalar); + + // Perform operation + switch (mTypeOfOperation) { + case 0: // Addition + while (!it.IsAtEnd()) { + ito.Set(it.Get() + scalar_vector); + ++it; + ++ito; + } + break; + case 1: // Multiply + while (!it.IsAtEnd()) { + ito.Set(it.Get() * mScalar); + ++it; + ++ito; + } + break; + /* + case 2: // Inverse + while (!it.IsAtEnd()) { + if (it.Get() != 0) + ito.Set(mScalar / it.Get())); + else ito.Set(mDefaultPixelValue); + ++it; + ++ito; + } + break; + case 3: // Max + while (!it.IsAtEnd()) { + if (it.Get() < mScalar) ito.Set(PixelTypeDownCast(mScalar)); + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + case 4: // Min + while (!it.IsAtEnd()) { + if (it.Get() > mScalar) ito.Set(PixelTypeDownCast(mScalar)); + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + case 5: // Absolute value + while (!it.IsAtEnd()) { + if (it.Get() <= 0) ito.Set(PixelTypeDownCast(-it.Get())); + // <= zero to avoid warning for unsigned types + else ito.Set(PixelTypeDownCast(it.Get())); + ++it; + ++ito; + } + break; + case 6: // Squared value + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast((double)it.Get()*(double)it.Get())); + ++it; + ++ito; + } + break; + case 7: // Log + while (!it.IsAtEnd()) { + if (it.Get() > 0) + ito.Set(PixelTypeDownCast(log((double)it.Get()))); + else ito.Set(mDefaultPixelValue); + ++it; + ++ito; + } + break; + case 8: // exp + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast(exp((double)it.Get()))); + ++it; + ++ito; + } + break; + case 9: // sqrt + while (!it.IsAtEnd()) { + if (it.Get() > 0) + ito.Set(PixelTypeDownCast(sqrt((double)it.Get()))); + else { + if (it.Get() ==0) ito.Set(0); + else ito.Set(mDefaultPixelValue); + } + ++it; + ++ito; + } + break; + case 10: // exp + while (!it.IsAtEnd()) { + ito.Set(PixelTypeDownCast((0x10000 - (double)it.Get())/mScalar)); + ++it; + ++ito; + } + break; + */ + case 11: // divide + while (!it.IsAtEnd()) { + ito.Set(it.Get() / mScalar); + ++it; + ++ito; + } + break; + default: // error ? + std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; + exit(-1); + } + + } + + template<> + template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > + > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it1, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito) + { + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter1; + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter2; + typedef itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > Iter3; + + it1.GoToBegin(); + it2.GoToBegin(); + ito.GoToBegin(); + typedef Iter3::PixelType PixelType; + + switch (mTypeOfOperation) { + case 0: // Addition + while (!ito.IsAtEnd()) { + ito.Set(it1.Get() + it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 1: // Multiply + while (!ito.IsAtEnd()) { + ito.Set(it1.Get() * it2.Get()) ); + ++it1; + ++it2; + ++ito; + } + break; + case 2: // Divide + while (!ito.IsAtEnd()) { + if (it1.Get() != 0) + ito.Set(it1.Get() / it2.Get())); + else ito.Set(mDefaultPixelValue); + ++it1; + ++it2; + ++ito; + } + break; + case 3: // Max + while (!ito.IsAtEnd()) { + if (it1.Get() < it2.Get()) ito.Set(it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + case 4: // Min + while (!ito.IsAtEnd()) { + if (it1.Get() > it2.Get()) ito.Set(it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + */ + case 5: // Absolute difference + while (!ito.IsAtEnd()) { + ito.Set(it2.Get()-it1.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 6: // Squared differences + while (!ito.IsAtEnd()) { + ito.Set(pow(it1.Get()-it2.Get(),2))); + ++it1; + ++it2; + ++ito; + } + break; + */ + case 7: // Difference + while (!ito.IsAtEnd()) { + ito.Set(it1.Get()-it2.Get()); + ++it1; + ++it2; + ++ito; + } + break; + /* + case 8: // Relative Difference + while (!ito.IsAtEnd()) { + if (it1.Get() != 0) ito.Set(PixelTypeDownCast(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get()); + else ito.Set(0.0); + ++it1; + ++it2; + ++ito; + } + break; + */ + default: // error ? + std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl; + exit(-1); + } + } + } #endif //define CLITKIMAGEARITHMGENERICFILTER_CXX diff --git a/tools/clitkImageArithmGenericFilter.h b/tools/clitkImageArithmGenericFilter.h index 45e0959..39bfeb9 100644 --- a/tools/clitkImageArithmGenericFilter.h +++ b/tools/clitkImageArithmGenericFilter.h @@ -97,6 +97,7 @@ namespace clitk { }; // end class ImageArithmGenericFilter + // specializations for itk::Vector, 3u template<> template<> void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >(); @@ -118,6 +119,27 @@ namespace clitk { itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); + // specializations for itk::Vector, 3u + template<> template<> + void ImageArithmGenericFilter::UpdateWithInputImageType< itk::Image< itk::Vector, 3u > >(); + + template<> template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > + > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); + + template<> template<> + void ImageArithmGenericFilter::ComputeImage< + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > >, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > + > + (itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it1, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > it2, + itk::ImageRegionIterator< itk::Image< itk::Vector, 3u > > ito); } // end namespace //-------------------------------------------------------------------- diff --git a/tools/clitkImageArithmGenericFilter.txx b/tools/clitkImageArithmGenericFilter.txx index 8efd6a7..1afb77f 100644 --- a/tools/clitkImageArithmGenericFilter.txx +++ b/tools/clitkImageArithmGenericFilter.txx @@ -46,6 +46,7 @@ void ImageArithmGenericFilter::InitializeImageType() { ADD_DEFAULT_IMAGE_TYPES(Dim); ADD_VEC_IMAGE_TYPE(3u,3u,float); + ADD_VEC_IMAGE_TYPE(3u,3u,double); } //-------------------------------------------------------------------- -- 2.47.1