diff --git a/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp b/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp index 87269991c4..29dfb3d4d7 100644 --- a/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp +++ b/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp @@ -1,500 +1,500 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkTransformationOperation.h" #include #include #include #include #include #include // Wavelet #include #include #include #include #include #include #include #include -#include +#include #include "itkZeroFluxNeumannBoundaryCondition.h" #include "itkPeriodicBoundaryCondition.h" #include "itkConstantBoundaryCondition.h" //#include #include "itkCastImageFilter.h" #include "itkUnaryFunctorImageFilter.h" #include #include #include namespace mitk { namespace Functor { template< class TInput> class ThresholdValue { public: ThresholdValue() {}; ~ThresholdValue() {}; bool operator!=(const ThresholdValue &) const { return false; } bool operator==(const ThresholdValue & other) const { return !(*this != other); } inline unsigned short operator()(const TInput & A) const { if (A < value) return 0; else return 1; } double value = 0.0; }; template< class TInput, class TOutput> class RoundValue { public: RoundValue() {}; ~RoundValue() {}; bool operator!=(const RoundValue &) const { return false; } bool operator==(const RoundValue & other) const { return !(*this != other); } inline TOutput operator()(const TInput & A) const { return std::round(A); } }; } } template static void ExecuteMultiResolution(itk::Image* image, unsigned int numberOfLevels, bool outputAsDouble, std::vector &resultImages) { typedef itk::Image ImageType; typedef itk::Image DoubleOutputType; typedef itk::RecursiveMultiResolutionPyramidImageFilter ImageTypeFilterType; typedef itk::RecursiveMultiResolutionPyramidImageFilter DoubleTypeFilterType; if (outputAsDouble) { typename DoubleTypeFilterType::Pointer recursiveMultiResolutionPyramidImageFilter = DoubleTypeFilterType::New(); recursiveMultiResolutionPyramidImageFilter->SetInput(image); recursiveMultiResolutionPyramidImageFilter->SetNumberOfLevels(numberOfLevels); recursiveMultiResolutionPyramidImageFilter->Update(); // This outputs the levels (0 is the lowest resolution) for (unsigned int i = 0; i < numberOfLevels; ++i) { mitk::Image::Pointer outputImage = mitk::Image::New(); CastToMitkImage(recursiveMultiResolutionPyramidImageFilter->GetOutput(i), outputImage); resultImages.push_back(outputImage); } } else { typename ImageTypeFilterType::Pointer recursiveMultiResolutionPyramidImageFilter = ImageTypeFilterType::New(); recursiveMultiResolutionPyramidImageFilter->SetInput(image); recursiveMultiResolutionPyramidImageFilter->SetNumberOfLevels(numberOfLevels); recursiveMultiResolutionPyramidImageFilter->Update(); // This outputs the levels (0 is the lowest resolution) for (unsigned int i = 0; i < numberOfLevels; ++i) { mitk::Image::Pointer outputImage = mitk::Image::New(); CastToMitkImage(recursiveMultiResolutionPyramidImageFilter->GetOutput(i), outputImage); resultImages.push_back(outputImage); } } } std::vector mitk::TransformationOperation::MultiResolution(Image::Pointer & image, unsigned int numberOfLevels, bool outputAsDouble) { std::vector resultImages; AccessByItk_n(image, ExecuteMultiResolution, (numberOfLevels, outputAsDouble, resultImages)); return resultImages; } template static void ExecuteLaplacianOfGaussian(itk::Image* image, double sigma, bool outputAsDouble, mitk::Image::Pointer &resultImage) { typedef itk::Image ImageType; typedef itk::Image DoubleOutputType; typedef itk::LaplacianRecursiveGaussianImageFilter ImageTypeFilterType; typedef itk::LaplacianRecursiveGaussianImageFilter DoubleTypeFilterType; if (outputAsDouble) { typename DoubleTypeFilterType::Pointer filter = DoubleTypeFilterType::New(); filter->SetInput(image); filter->SetSigma(sigma); filter->Update(); CastToMitkImage(filter->GetOutput(), resultImage); } else { typename ImageTypeFilterType::Pointer filter = ImageTypeFilterType::New(); filter->SetInput(image); filter->SetSigma(sigma); filter->Update(); CastToMitkImage(filter->GetOutput(), resultImage); } } mitk::Image::Pointer mitk::TransformationOperation::LaplacianOfGaussian(Image::Pointer & image, double sigma, bool outputAsDouble) { Image::Pointer resultImage; AccessByItk_n(image, ExecuteLaplacianOfGaussian, (sigma, outputAsDouble, resultImage)); return resultImage; } template static void ExecuteSpecificWaveletTransformation(itk::Image* image, unsigned int numberOfLevels, unsigned int numberOfBands, mitk::BorderCondition condition, std::vector &resultImages) { const unsigned int Dimension = VImageDimension; typedef TInputPixel PixelType; typedef TOutputPixel OutputPixelType; typedef itk::Image< PixelType, Dimension > ImageType; typedef itk::Image< double, Dimension > DoubleImageType; typedef itk::Image< OutputPixelType, Dimension > OutputImageType; typedef itk::CastImageFilter< ImageType, DoubleImageType > CastFilterType; - typedef itk::FFTPadPositiveIndexImageFilter< DoubleImageType > FFTPadType; + typedef itk::FFTPadImageFilter< DoubleImageType > FFTPadType; typedef itk::ForwardFFTImageFilter< DoubleImageType, itk::Image< std::complex, Dimension> > FFTFilterType; typedef typename FFTFilterType::OutputImageType ComplexImageType; typedef TWaveletFunction WaveletFunctionType; typedef itk::WaveletFrequencyFilterBankGenerator< ComplexImageType, WaveletFunctionType > WaveletFilterBankType; typedef itk::WaveletFrequencyForward< ComplexImageType, ComplexImageType, WaveletFilterBankType > ForwardWaveletType; typedef itk::InverseFFTImageFilter< ComplexImageType, OutputImageType > InverseFFTFilterType; // Convert input parameter unsigned int highSubBands = numberOfBands; //inputBands; unsigned int levels = numberOfLevels; // Perform FFT on input image typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput(image); // Pad Image so it fits the expect typename FFTPadType::Pointer fftpad = FFTPadType::New(); fftpad->SetSizeGreatestPrimeFactor(4); itk::ConstantBoundaryCondition< DoubleImageType > constantBoundaryCondition; itk::PeriodicBoundaryCondition< DoubleImageType > periodicBoundaryCondition; itk::ZeroFluxNeumannBoundaryCondition< DoubleImageType > zeroFluxNeumannBoundaryCondition; switch (condition) { case mitk::BorderCondition::Constant: fftpad->SetBoundaryCondition(&constantBoundaryCondition); break; case mitk::BorderCondition::Periodic: fftpad->SetBoundaryCondition(&periodicBoundaryCondition); break; case mitk::BorderCondition::ZeroFluxNeumann: fftpad->SetBoundaryCondition(&zeroFluxNeumannBoundaryCondition); break; default: break; } fftpad->SetInput(castFilter->GetOutput()); typename FFTFilterType::Pointer fftFilter = FFTFilterType::New(); fftFilter->SetInput(fftpad->GetOutput()); // Calculate forward transformation typename ForwardWaveletType::Pointer forwardWavelet = ForwardWaveletType::New(); forwardWavelet->SetHighPassSubBands(highSubBands); forwardWavelet->SetLevels(levels); forwardWavelet->SetInput(fftFilter->GetOutput()); forwardWavelet->Update(); // Obtain target spacing, size and origin typename ComplexImageType::SpacingType inputSpacing; for (unsigned int i = 0; i < Dimension; ++i) { inputSpacing[i] = image->GetLargestPossibleRegion().GetSize()[i]; } typename ComplexImageType::SpacingType expectedSpacing = inputSpacing; typename ComplexImageType::PointType inputOrigin = image->GetOrigin(); typename ComplexImageType::PointType expectedOrigin = inputOrigin; typename ComplexImageType::SizeType inputSize = fftFilter->GetOutput()->GetLargestPossibleRegion().GetSize(); typename ComplexImageType::SizeType expectedSize = inputSize; // Inverse FFT to obtain filtered images typename InverseFFTFilterType::Pointer inverseFFT = InverseFFTFilterType::New(); for (unsigned int level = 0; level < numberOfLevels + 1; ++level) { double scaleFactorPerLevel = std::pow(static_cast< double >(forwardWavelet->GetScaleFactor()),static_cast< double >(level)); for (unsigned int i = 0; i < Dimension; ++i) { expectedSize[i] = inputSize[i] / scaleFactorPerLevel; expectedOrigin[i] = inputOrigin[i]; expectedSpacing[i] = inputSpacing[i] * scaleFactorPerLevel; } for (unsigned int band = 0; band < highSubBands; ++band) { unsigned int nOutput = level * forwardWavelet->GetHighPassSubBands() + band; // Do not compute bands in low-pass level. if (level == numberOfLevels && band == 0) { nOutput = forwardWavelet->GetTotalOutputs() - 1; } else if (level == numberOfLevels && band != 0) { break; } inverseFFT->SetInput(forwardWavelet->GetOutput(nOutput)); inverseFFT->Update(); auto itkOutputImage = inverseFFT->GetOutput(); itkOutputImage->SetSpacing(expectedSpacing); mitk::Image::Pointer outputImage = mitk::Image::New(); CastToMitkImage(itkOutputImage, outputImage); resultImages.push_back(outputImage); } } } template static void ExecuteWaveletTransformation(itk::Image* image, unsigned int numberOfLevels, unsigned int numberOfBands, mitk::BorderCondition condition, mitk::WaveletType waveletType, std::vector &resultImages) { typedef itk::Point< double, VImageDimension > PointType; typedef itk::HeldIsotropicWavelet< double, VImageDimension, PointType > HeldIsotropicWaveletType; typedef itk::VowIsotropicWavelet< double, VImageDimension, PointType > VowIsotropicWaveletType; typedef itk::SimoncelliIsotropicWavelet< double, VImageDimension, PointType > SimoncelliIsotropicWaveletType; typedef itk::ShannonIsotropicWavelet< double, VImageDimension, PointType > ShannonIsotropicWaveletType; switch (waveletType) { case mitk::WaveletType::Held: ExecuteSpecificWaveletTransformation(image, numberOfLevels, numberOfBands, condition, resultImages); break; case mitk::WaveletType::Shannon: ExecuteSpecificWaveletTransformation(image, numberOfLevels, numberOfBands, condition, resultImages); break; case mitk::WaveletType::Simoncelli: ExecuteSpecificWaveletTransformation(image, numberOfLevels, numberOfBands, condition, resultImages); break; case mitk::WaveletType::Vow: ExecuteSpecificWaveletTransformation(image, numberOfLevels, numberOfBands, condition, resultImages); break; default: ExecuteSpecificWaveletTransformation(image, numberOfLevels, numberOfBands, condition, resultImages); break; } } std::vector mitk::TransformationOperation::WaveletForward(Image::Pointer & image, unsigned int numberOfLevels, unsigned int numberOfBands, mitk::BorderCondition condition, mitk::WaveletType waveletType) { std::vector resultImages; AccessByItk_n(image, ExecuteWaveletTransformation, (numberOfLevels, numberOfBands, condition, waveletType, resultImages)); return resultImages; } template static void ExecuteImageTypeToDouble(itk::Image* image, mitk::Image::Pointer &outputImage) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< double, VImageDimension > DoubleImageType; typedef itk::CastImageFilter< ImageType, DoubleImageType > CastFilterType; typedef itk::ImageDuplicator< DoubleImageType > DuplicatorType; // Perform FFT on input image typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput(image); castFilter->Update(); typename DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage(castFilter->GetOutput()); duplicator->Update(); CastToMitkImage(duplicator->GetOutput(), outputImage); } template static void ExecuteRoundImage(itk::Image* /*image*/, mitk::Image::Pointer resampledImage, mitk::Image::Pointer &outputImage) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< double, VImageDimension > DoubleImageType; typedef itk::UnaryFunctorImageFilter< DoubleImageType, ImageType, mitk::Functor::RoundValue > DefaultFilterType; typename DoubleImageType::Pointer itkImage = DoubleImageType::New(); mitk::CastToItkImage(resampledImage, itkImage); typename DefaultFilterType::Pointer filter = DefaultFilterType::New(); filter->SetInput(itkImage); filter->Update(); CastToMitkImage(filter->GetOutput(), outputImage); } mitk::Image::Pointer mitk::TransformationOperation::ResampleImage(Image::Pointer &image, mitk::Vector3D spacingVector, mitk::ImageMappingInterpolator::Type interpolator, GridInterpolationPositionType position, bool returnAsDouble, bool roundOutput) { // Convert image to double if required mitk::Image::Pointer tmpImage = image; if (returnAsDouble) { AccessByItk_n(image, ExecuteImageTypeToDouble, (tmpImage)); } auto newGeometry = image->GetGeometry()->Clone(); mitk::Vector3D spacing; mitk::BaseGeometry::BoundsArrayType bounds = newGeometry->GetBounds(); for (int i = 0; i < 3; ++i) { spacing[i] = newGeometry->GetSpacing()[i]; //bounds[i*2+1] = newGeometry->GetBounds()[i * 2 + 1]; if (spacingVector[i] > 0) { spacing[i] = spacingVector[i]; if (position == mitk::GridInterpolationPositionType::SameSize) { unsigned int samples = image->GetDimensions()[i]; double currentSpacing = newGeometry->GetSpacing()[i]; double newFactor = std::floor(samples*currentSpacing / spacingVector[i]); spacing[i] = samples * currentSpacing / newFactor; } } bounds[i * 2] = 0; bounds[i*2+1] = std::ceil(bounds[i*2+1] * newGeometry->GetSpacing()[i] *1.0 / spacing[i]); } mitk::Point3D origin = newGeometry->GetOrigin(); if (position == mitk::GridInterpolationPositionType::CenterAligned) { for (int i = 0; i < 3; ++i) { double oldLength = newGeometry->GetSpacing()[i] * newGeometry->GetBounds()[i*2+1]; double newLength = spacing[i] * bounds[i*2+1]; origin[i] = origin[i] - (newLength - oldLength) / 2; } } newGeometry->SetSpacing(spacing); newGeometry->SetOrigin(origin); newGeometry->SetBounds(bounds); mitk::Image::Pointer tmpResult = ImageMappingHelper::map( tmpImage, mitk::GenerateIdentityRegistration3D().GetPointer(), false, 0.0, //Padding Value newGeometry.GetPointer(), false, 0, //Error Value interpolator ); mitk::Image::Pointer result = mitk::Image::New(); if (roundOutput) { AccessByItk_n(tmpImage, ExecuteRoundImage, (tmpResult, result)); } else { result = tmpResult; } return result; } template static void ExecuteImageThresholding(itk::Image* image, mitk::Image::Pointer &resultImage) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::UnaryFunctorImageFilter< ImageType, MaskType, mitk::Functor::ThresholdValue > DefaultFilterType; typename DefaultFilterType::Pointer filter = DefaultFilterType::New(); filter->SetInput(image); filter->GetFunctor().value = 0.5; filter->Update(); CastToMitkImage(filter->GetOutput(), resultImage); } mitk::Image::Pointer mitk::TransformationOperation::ResampleMask(Image::Pointer &image, mitk::Vector3D spacingVector, mitk::ImageMappingInterpolator::Type interpolator, GridInterpolationPositionType position) { mitk::Image::Pointer result; if (interpolator == mitk::ImageMappingInterpolator::NearestNeighbor) { result = TransformationOperation::ResampleImage(image, spacingVector, interpolator, position, false, false); } else { auto tmpResult = TransformationOperation::ResampleImage(image, spacingVector, interpolator, position, true, false); AccessByItk_n(tmpResult, ExecuteImageThresholding, (result)); } return result; } namespace itk { namespace utils { IndexPairType IndexToLevelBandSteerablePyramid(unsigned int linearIndex, unsigned int levels, unsigned int bands) { unsigned int totalOutputs = 1 + levels * bands; if (linearIndex > totalOutputs - 1) { itkGenericExceptionMacro(<< "Failed converting linearIndex " << linearIndex << " with levels: " << levels << " bands: " << bands << " to Level,Band pair : out of bounds"); } // Low pass (band = 0). if (linearIndex == totalOutputs - 1) { return std::make_pair(levels - 1, 0); } unsigned int band = (linearIndex) % bands + 1; // note integer division ahead. unsigned int level = (linearIndex) / bands; itkAssertInDebugAndIgnoreInReleaseMacro(level < levels); return std::make_pair(level, band); } // Instantiation template unsigned int ComputeMaxNumberOfLevels<3>(const Size< 3 >& inputSize, const unsigned int & scaleFactor); template unsigned int ComputeMaxNumberOfLevels<2>(const Size< 2 >& inputSize, const unsigned int & scaleFactor); } // end namespace utils } // end namespace itk diff --git a/Modules/Classification/CLUtilities/CMakeLists.txt b/Modules/Classification/CLUtilities/CMakeLists.txt index 0f05dd8069..577381576a 100644 --- a/Modules/Classification/CLUtilities/CMakeLists.txt +++ b/Modules/Classification/CLUtilities/CMakeLists.txt @@ -1,10 +1,10 @@ mitk_create_module( DEPENDS MitkCore MitkCLCore MitkCommandLine MitkDICOM - PACKAGE_DEPENDS PUBLIC Eigen OpenMP PRIVATE tinyxml2 VTK|FiltersStatistics + PACKAGE_DEPENDS PUBLIC Eigen OpenMP PRIVATE tinyxml2 ITK|MathematicalMorphology+Smoothing VTK|FiltersStatistics ) if(TARGET ${MODULE_TARGET}) if(BUILD_TESTING) add_subdirectory(test) endif() endif() diff --git a/Modules/Classification/CLUtilities/include/itkLocalIntensityFilter.hxx b/Modules/Classification/CLUtilities/include/itkLocalIntensityFilter.hxx index a5f14bd8c3..ac6a61b5b6 100644 --- a/Modules/Classification/CLUtilities/include/itkLocalIntensityFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkLocalIntensityFilter.hxx @@ -1,314 +1,316 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef itkLocalIntensityFilter_cpp #define itkLocalIntensityFilter_cpp #include #include #include #include #include #include "itkImageScanlineIterator.h" #include "itkProgressReporter.h" namespace itk { template< typename TInputImage > LocalIntensityFilter< TInputImage > ::LocalIntensityFilter() :m_ThreadLocalMaximum(1), m_ThreadLocalPeakValue(1), m_ThreadGlobalPeakValue(1) { + this->DynamicMultiThreadingOff(); + // first output is a copy of the image, DataObject created by // superclass // allocate the data objects for the outputs which are // just decorators around real types for (int i = 1; i < 4; ++i) { typename RealObjectType::Pointer output = static_cast< RealObjectType * >(this->MakeOutput(i).GetPointer()); this->ProcessObject::SetNthOutput(i, output.GetPointer()); } } template< typename TInputImage > DataObject::Pointer LocalIntensityFilter< TInputImage > ::MakeOutput(DataObjectPointerArraySizeType output) { switch (output) { case 0: return TInputImage::New().GetPointer(); break; case 1: case 2: case 3: case 4: case 5: case 6: return RealObjectType::New().GetPointer(); break; default: // might as well make an image return TInputImage::New().GetPointer(); break; } } template< typename TInputImage > typename LocalIntensityFilter< TInputImage >::RealObjectType * LocalIntensityFilter< TInputImage > ::GetLocalPeakOutput() { return static_cast< RealObjectType * >(this->ProcessObject::GetOutput(1)); } template< typename TInputImage > const typename LocalIntensityFilter< TInputImage >::RealObjectType * LocalIntensityFilter< TInputImage > ::GetLocalPeakOutput() const { return static_cast< const RealObjectType * >(this->ProcessObject::GetOutput(1)); } template< typename TInputImage > typename LocalIntensityFilter< TInputImage >::RealObjectType * LocalIntensityFilter< TInputImage > ::GetGlobalPeakOutput() { return static_cast< RealObjectType * >(this->ProcessObject::GetOutput(2)); } template< typename TInputImage > const typename LocalIntensityFilter< TInputImage >::RealObjectType * LocalIntensityFilter< TInputImage > ::GetGlobalPeakOutput() const { return static_cast< const RealObjectType * >(this->ProcessObject::GetOutput(2)); } template< typename TInputImage > typename LocalIntensityFilter< TInputImage >::RealObjectType * LocalIntensityFilter< TInputImage > ::GetLocalMaximumOutput() { return static_cast< RealObjectType * >(this->ProcessObject::GetOutput(3)); } template< typename TInputImage > const typename LocalIntensityFilter< TInputImage >::RealObjectType * LocalIntensityFilter< TInputImage > ::GetLocalMaximumOutput() const { return static_cast< const RealObjectType * >(this->ProcessObject::GetOutput(3)); } template< typename TInputImage > void LocalIntensityFilter< TInputImage > ::GenerateInputRequestedRegion() { Superclass::GenerateInputRequestedRegion(); if (this->GetInput()) { InputImagePointer image = const_cast< typename Superclass::InputImageType * >(this->GetInput()); image->SetRequestedRegionToLargestPossibleRegion(); } } template< typename TInputImage > void LocalIntensityFilter< TInputImage > ::EnlargeOutputRequestedRegion(DataObject *data) { Superclass::EnlargeOutputRequestedRegion(data); data->SetRequestedRegionToLargestPossibleRegion(); } template< typename TInputImage > void LocalIntensityFilter< TInputImage > ::AllocateOutputs() { // Pass the input through as the output InputImagePointer image = const_cast< TInputImage * >(this->GetInput()); this->GraftOutput(image); // Nothing that needs to be allocated for the remaining outputs } template< typename TInputImage > void LocalIntensityFilter< TInputImage > ::BeforeThreadedGenerateData() { - ThreadIdType numberOfThreads = this->GetNumberOfThreads(); + ThreadIdType numberOfThreads = this->GetNumberOfWorkUnits(); // Resize the thread temporaries m_ThreadLocalMaximum.SetSize(numberOfThreads); m_ThreadLocalPeakValue.SetSize(numberOfThreads); m_ThreadGlobalPeakValue.SetSize(numberOfThreads); // Initialize the temporaries m_ThreadLocalMaximum.Fill(std::numeric_limits< RealType>::lowest()); m_ThreadLocalPeakValue.Fill(std::numeric_limits< RealType>::lowest()); m_ThreadGlobalPeakValue.Fill(std::numeric_limits< RealType>::lowest()); } template< typename TInputImage > void LocalIntensityFilter< TInputImage > ::AfterThreadedGenerateData() { ThreadIdType i; - ThreadIdType numberOfThreads = this->GetNumberOfThreads(); + ThreadIdType numberOfThreads = this->GetNumberOfWorkUnits(); RealType localMaximum = std::numeric_limits< RealType >::lowest(); RealType localPeakValue = std::numeric_limits< RealType >::lowest(); RealType globalPeakValue = std::numeric_limits< RealType >::lowest(); for (i = 0; i < numberOfThreads; i++) { globalPeakValue = std::max(globalPeakValue, m_ThreadGlobalPeakValue[i]); if (localMaximum == m_ThreadLocalMaximum[i]) { localPeakValue = std::max< RealType >(m_ThreadLocalPeakValue[i], localPeakValue); } else if (localMaximum < m_ThreadLocalMaximum[i]) { localMaximum = m_ThreadLocalMaximum[i]; localPeakValue = m_ThreadLocalPeakValue[i]; } } // Set the outputs this->GetLocalPeakOutput()->Set(localPeakValue); this->GetGlobalPeakOutput()->Set(globalPeakValue); this->GetLocalMaximumOutput()->Set(localMaximum); } template< typename TInputImage > void LocalIntensityFilter< TInputImage > ::ThreadedGenerateData(const RegionType & outputRegionForThread, ThreadIdType threadId) { typename TInputImage::ConstPointer itkImage = this->GetInput(); typename MaskImageType::Pointer itkMask = m_Mask; double range = m_Range; double minimumSpacing = 1; typename TInputImage::SizeType regionSize; int offset = std::ceil(range / minimumSpacing); regionSize.Fill(1); for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i) { minimumSpacing = itkImage->GetSpacing()[i]; offset = std::ceil(range / minimumSpacing); regionSize[i] = offset; } itk::ConstNeighborhoodIterator iter(regionSize, itkImage, outputRegionForThread); itk::ConstNeighborhoodIterator iterMask(regionSize, itkMask, outputRegionForThread); typename TInputImage::PointType origin; typename TInputImage::PointType localPoint; itk::Index index; double tmpPeakValue; double globalPeakValue = std::numeric_limits::lowest(); double localPeakValue = std::numeric_limits::lowest(); PixelType localMaximum = std::numeric_limits::lowest(); std::vector vectorIsInRange; index = iter.GetIndex(); itkImage->TransformIndexToPhysicalPoint(index, origin); for (itk::SizeValueType i = 0; i < iter.Size(); ++i) { itkImage->TransformIndexToPhysicalPoint(iter.GetIndex(i), localPoint); double dist = origin.EuclideanDistanceTo(localPoint); vectorIsInRange.push_back((dist < range)); } int count = 0; iter.NeedToUseBoundaryConditionOff(); iterMask.NeedToUseBoundaryConditionOff(); auto imageSize = itkImage->GetLargestPossibleRegion().GetSize(); unsigned int imageDimension = itkImage->GetImageDimension(); while (!iter.IsAtEnd()) { if (iterMask.GetCenterPixel() > 0) { tmpPeakValue = 0; count = 0; for (itk::SizeValueType i = 0; i < iter.Size(); ++i) { if (vectorIsInRange[i]) { auto localIndex = iter.GetIndex(i); bool calculatePoint = true; for (unsigned int dimension = 0; dimension < imageDimension; ++dimension) { calculatePoint &= (localIndex[dimension] < static_cast(imageSize[dimension])); calculatePoint &= (0 <= localIndex[dimension]); } if (calculatePoint) { tmpPeakValue += iter.GetPixel(i); ++count; } } } tmpPeakValue /= count; globalPeakValue = std::max(tmpPeakValue, globalPeakValue); auto currentCenterPixelValue = iter.GetCenterPixel(); if (localMaximum == currentCenterPixelValue) { localPeakValue = std::max(tmpPeakValue, localPeakValue); } else if (localMaximum < currentCenterPixelValue) { localMaximum = currentCenterPixelValue; localPeakValue = tmpPeakValue; } } ++iterMask; ++iter; } m_ThreadLocalMaximum[threadId] = localMaximum; m_ThreadLocalPeakValue[threadId] = localPeakValue; m_ThreadGlobalPeakValue[threadId] = globalPeakValue; } template< typename TImage > void LocalIntensityFilter< TImage > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "Local Peak: " << this->GetLocalPeak() << std::endl; os << indent << "Global Peak: " << this->GetGlobalPeak() << std::endl; os << indent << "Local Maximum: " << this->GetLocalMaximum() << std::endl; } } // end namespace itk #endif diff --git a/Modules/DICOMTesting/include/mitkTestDICOMLoading.h b/Modules/DICOMTesting/include/mitkTestDICOMLoading.h index 63289abc0b..ad041a45a6 100644 --- a/Modules/DICOMTesting/include/mitkTestDICOMLoading.h +++ b/Modules/DICOMTesting/include/mitkTestDICOMLoading.h @@ -1,111 +1,111 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkTestDICOMLoading_h #define mitkTestDICOMLoading_h #include "mitkClassicDICOMSeriesReader.h" #include "mitkPropertyKeyPath.h" #include "MitkDICOMTestingExports.h" namespace mitk { class MITKDICOMTESTING_EXPORT TestDICOMLoading { public: typedef std::list ImageList; TestDICOMLoading(); ImageList LoadFiles( const StringList & files ); Image::Pointer DecorateVerifyCachedImage( const StringList& files, mitk::Image::Pointer cachedImage ); Image::Pointer DecorateVerifyCachedImage( const StringList& files, DICOMTagCache*, mitk::Image::Pointer cachedImage ); /** \brief Dump relevant image information for later comparison. \sa CompareImageInformationDumps */ std::string DumpImageInformation( const Image* image ); /** \brief Compare two image information dumps. \return true, if dumps are sufficiently equal (see parameters) \sa DumpImageInformation */ bool CompareImageInformationDumps( const std::string& reference, const std::string& test ); private: typedef std::map KeyValueMap; ClassicDICOMSeriesReader::Pointer BuildDICOMReader(); void SetDefaultLocale(); void ResetUserLocale(); - std::string ComponentTypeToString( int type ); + std::string ComponentTypeToString( itk::IOComponentEnum type ); KeyValueMap ParseDump( const std::string& dump ); bool CompareSpacedValueFields( const std::string& reference, const std::string& test, double eps = mitk::eps ); /** Compress whitespace in string \param pString input string \param pFill replacement whitespace (only whitespace in string after reduction) \param pWhitespace characters handled as whitespace */ std::string reduce(const std::string& pString, const std::string& pFill = " ", const std::string& pWhitespace = " \t"); /** Remove leading and trailing whitespace \param pString input string \param pWhitespace characters handled as whitespace */ std::string trim(const std::string& pString, const std::string& pWhitespace = " \t"); template bool StringToNumber(const std::string& s, T& value) { std::stringstream stream(s); stream >> value; return (!stream.fail()) && (std::abs(value) <= std::numeric_limits::max()); } static void AddPropertyToDump(const mitk::PropertyKeyPath& key, const mitk::Image* image, std::stringstream& result); const char* m_PreviousCLocale; std::locale m_PreviousCppLocale; }; } #endif diff --git a/Modules/DICOMTesting/src/mitkTestDICOMLoading.cpp b/Modules/DICOMTesting/src/mitkTestDICOMLoading.cpp index 9d63c65aa9..d6a44162f7 100644 --- a/Modules/DICOMTesting/src/mitkTestDICOMLoading.cpp +++ b/Modules/DICOMTesting/src/mitkTestDICOMLoading.cpp @@ -1,585 +1,585 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "mitkTestDICOMLoading.h" #include "mitkDICOMIOMetaInformationPropertyConstants.h" #include "mitkDICOMProperty.h" #include "mitkArbitraryTimeGeometry.h" #include #include #include #include "itksys/SystemTools.hxx" mitk::TestDICOMLoading::TestDICOMLoading() :m_PreviousCLocale(nullptr) { } void mitk::TestDICOMLoading::SetDefaultLocale() { // remember old locale only once if (m_PreviousCLocale == nullptr) { m_PreviousCLocale = setlocale(LC_NUMERIC, nullptr); // set to "C" setlocale(LC_NUMERIC, "C"); m_PreviousCppLocale = std::cin.getloc(); std::locale l( "C" ); std::cin.imbue(l); std::cout.imbue(l); } } void mitk::TestDICOMLoading::ResetUserLocale() { if (m_PreviousCLocale) { setlocale(LC_NUMERIC, m_PreviousCLocale); std::cin.imbue(m_PreviousCppLocale); std::cout.imbue(m_PreviousCppLocale); m_PreviousCLocale = nullptr; } } mitk::TestDICOMLoading::ImageList mitk::TestDICOMLoading ::LoadFiles( const StringList& files ) { for (auto iter = files.begin(); iter != files.end(); ++iter) { MITK_DEBUG << "File " << *iter; } ImageList result; ClassicDICOMSeriesReader::Pointer reader = this->BuildDICOMReader(); reader->SetTagLookupTableToPropertyFunctor(mitk::GetDICOMPropertyForDICOMValuesFunctor); reader->SetInputFiles( files ); reader->AnalyzeInputFiles(); reader->PrintOutputs(std::cout,true); reader->LoadImages(); unsigned int numberOfImages = reader->GetNumberOfOutputs(); for (unsigned imageIndex = 0; imageIndex < numberOfImages; ++imageIndex) { const DICOMImageBlockDescriptor& block = reader->GetOutput(imageIndex); result.push_back( block.GetMitkImage() ); } return result; } mitk::ClassicDICOMSeriesReader::Pointer mitk::TestDICOMLoading ::BuildDICOMReader() { ClassicDICOMSeriesReader::Pointer reader = ClassicDICOMSeriesReader::New(); reader->SetFixTiltByShearing(true); return reader; } mitk::Image::Pointer mitk::TestDICOMLoading ::DecorateVerifyCachedImage( const StringList& files, mitk::DICOMTagCache* tagCache, mitk::Image::Pointer cachedImage ) { DICOMImageBlockDescriptor block; DICOMImageFrameList framelist; for (auto iter = files.begin(); iter != files.end(); ++iter) { framelist.push_back( DICOMImageFrameInfo::New(*iter) ); } block.SetImageFrameList( framelist ); block.SetTagCache( tagCache ); block.SetMitkImage( cachedImage ); // this should/will create a propertylist describing the image slices return block.GetMitkImage(); } mitk::Image::Pointer mitk::TestDICOMLoading ::DecorateVerifyCachedImage( const StringList& files, mitk::Image::Pointer cachedImage ) { ClassicDICOMSeriesReader::Pointer reader = this->BuildDICOMReader(); reader->SetTagLookupTableToPropertyFunctor(mitk::GetDICOMPropertyForDICOMValuesFunctor); reader->SetInputFiles( files ); reader->AnalyzeInputFiles(); // This just creates a "tag cache and a nice DICOMImageBlockDescriptor. // Both of these could also be produced in a different way. The only // important thing is, that the DICOMImageBlockDescriptor knows a // tag-cache object when PropertyDecorateCachedMitkImageForImageBlockDescriptor // is called. if ( reader->GetNumberOfOutputs() != 1 ) { MITK_ERROR << "Reader produce " << reader->GetNumberOfOutputs() << " images instead of 1 expected.."; return nullptr; } DICOMImageBlockDescriptor block = reader->GetOutput(0); // creates a block copy block.SetMitkImage( cachedImage ); // this should/will create a propertylist describing the image slices return block.GetMitkImage(); } std::string -mitk::TestDICOMLoading::ComponentTypeToString(int type) +mitk::TestDICOMLoading::ComponentTypeToString(itk::IOComponentEnum type) { - if (type == itk::ImageIOBase::UCHAR) + if (type == itk::IOComponentEnum::UCHAR) return "UCHAR"; - else if (type == itk::ImageIOBase::CHAR) + else if (type == itk::IOComponentEnum::CHAR) return "CHAR"; - else if (type == itk::ImageIOBase::USHORT) + else if (type == itk::IOComponentEnum::USHORT) return "USHORT"; - else if (type == itk::ImageIOBase::SHORT) + else if (type == itk::IOComponentEnum::SHORT) return "SHORT"; - else if (type == itk::ImageIOBase::UINT) + else if (type == itk::IOComponentEnum::UINT) return "UINT"; - else if (type == itk::ImageIOBase::INT) + else if (type == itk::IOComponentEnum::INT) return "INT"; - else if (type == itk::ImageIOBase::ULONG) + else if (type == itk::IOComponentEnum::ULONG) return "ULONG"; - else if (type == itk::ImageIOBase::LONG) + else if (type == itk::IOComponentEnum::LONG) return "LONG"; - else if (type == itk::ImageIOBase::FLOAT) + else if (type == itk::IOComponentEnum::FLOAT) return "FLOAT"; - else if (type == itk::ImageIOBase::DOUBLE) + else if (type == itk::IOComponentEnum::DOUBLE) return "DOUBLE"; else return "UNKNOWN"; } // add a line to stringstream result (see DumpImageInformation #define DumpLine(field, data) DumpILine(0, field, data) // add an indented(!) line to stringstream result (see DumpImageInformation #define DumpILine(indent, field, data) \ { \ std::string DumpLine_INDENT; DumpLine_INDENT.resize(indent, ' ' ); \ result << DumpLine_INDENT << field << ": " << data << "\n"; \ } std::string mitk::TestDICOMLoading::DumpImageInformation( const Image* image ) { std::stringstream result; if (image == nullptr) return result.str(); SetDefaultLocale(); // basic image data DumpLine( "Pixeltype", ComponentTypeToString(image->GetPixelType().GetComponentType()) ); DumpLine( "BitsPerPixel", image->GetPixelType().GetBpe() ); DumpLine( "Dimension", image->GetDimension() ); result << "Dimensions: "; for (unsigned int dim = 0; dim < image->GetDimension(); ++dim) result << image->GetDimension(dim) << " "; result << "\n"; // geometry data result << "Geometry: \n"; const TimeGeometry* timeGeometry = image->GetTimeGeometry(); BaseGeometry* geometry = timeGeometry->GetGeometryForTimeStep(0); if (geometry) { AffineTransform3D* transform = geometry->GetIndexToWorldTransform(); if (transform) { result << " " << "Matrix: "; const AffineTransform3D::MatrixType& matrix = transform->GetMatrix(); for (unsigned int i = 0; i < 3; ++i) for (unsigned int j = 0; j < 3; ++j) result << matrix[i][j] << " "; result << "\n"; result << " " << "Offset: "; const AffineTransform3D::OutputVectorType& offset = transform->GetOffset(); for (unsigned int i = 0; i < 3; ++i) result << offset[i] << " "; result << "\n"; result << " " << "Center: "; const AffineTransform3D::InputPointType& center = transform->GetCenter(); for (unsigned int i = 0; i < 3; ++i) result << center[i] << " "; result << "\n"; result << " " << "Translation: "; const AffineTransform3D::OutputVectorType& translation = transform->GetTranslation(); for (unsigned int i = 0; i < 3; ++i) result << translation[i] << " "; result << "\n"; result << " " << "Scale: "; const double* scale = transform->GetScale(); for (unsigned int i = 0; i < 3; ++i) result << scale[i] << " "; result << "\n"; result << " " << "Origin: "; const Point3D& origin = geometry->GetOrigin(); for (unsigned int i = 0; i < 3; ++i) result << origin[i] << " "; result << "\n"; result << " " << "Spacing: "; const Vector3D& spacing = geometry->GetSpacing(); for (unsigned int i = 0; i < 3; ++i) result << spacing[i] << " "; result << "\n"; result << " " << "TimeBounds: "; /////////////////////////////////////// // Workarround T27883. See https://phabricator.mitk.org/T27883#219473 for more details. // This workarround should be removed as soon as T28262 is solved! TimeBounds timeBounds = timeGeometry->GetTimeBounds(); auto atg = dynamic_cast(timeGeometry); if (atg && atg->HasCollapsedFinalTimeStep()) { timeBounds[1] = timeBounds[1] - 1.; } //Original code: //const TimeBounds timeBounds = timeGeometry->GetTimeBounds(); // // End of workarround for T27883 ////////////////////////////////////// for (unsigned int i = 0; i < 2; ++i) result << timeBounds[i] << " "; result << "\n"; } } // io dicom meta information AddPropertyToDump(mitk::DICOMIOMetaInformationPropertyConstants::READER_CONFIGURATION(), image, result); AddPropertyToDump(mitk::DICOMIOMetaInformationPropertyConstants::READER_FILES(), image, result); AddPropertyToDump(mitk::DICOMIOMetaInformationPropertyConstants::READER_GANTRY_TILT_CORRECTED(), image, result); AddPropertyToDump(mitk::DICOMIOMetaInformationPropertyConstants::READER_IMPLEMENTATION_LEVEL(), image, result); AddPropertyToDump(mitk::DICOMIOMetaInformationPropertyConstants::READER_IMPLEMENTATION_LEVEL_STRING(), image, result); AddPropertyToDump(mitk::DICOMIOMetaInformationPropertyConstants::READER_PIXEL_SPACING_INTERPRETATION(), image, result); AddPropertyToDump(mitk::DICOMIOMetaInformationPropertyConstants::READER_PIXEL_SPACING_INTERPRETATION_STRING(), image, result); AddPropertyToDump(mitk::DICOMIOMetaInformationPropertyConstants::READER_3D_plus_t(), image, result); AddPropertyToDump(mitk::DICOMIOMetaInformationPropertyConstants::READER_DCMTK(), image, result); AddPropertyToDump(mitk::DICOMIOMetaInformationPropertyConstants::READER_GDCM(), image, result); ResetUserLocale(); return result.str(); } void mitk::TestDICOMLoading::AddPropertyToDump(const mitk::PropertyKeyPath& key, const mitk::Image* image, std::stringstream& result) { auto propKey = mitk::PropertyKeyPathToPropertyName(key); auto prop = image->GetProperty(propKey.c_str()); if (prop.IsNotNull()) { auto value = prop->GetValueAsString(); auto dicomProp = dynamic_cast< mitk::DICOMProperty*>(prop.GetPointer()); if (dicomProp != nullptr) { auto strippedProp = dicomProp->Clone(); if (key == mitk::DICOMIOMetaInformationPropertyConstants::READER_FILES()) {//strip dicom file information from path to ensure generalized dump files auto timePoints = strippedProp->GetAvailableTimeSteps(); for (auto timePoint : timePoints) { auto slices = strippedProp->GetAvailableSlices(timePoint); for (auto slice : slices) { auto value = strippedProp->GetValue(timePoint, slice); value = itksys::SystemTools::GetFilenameName(value); strippedProp->SetValue(timePoint, slice, value); } } } value = mitk::PropertyPersistenceSerialization::serializeTemporoSpatialStringPropertyToJSON(strippedProp); } result << propKey << ": " << value << "\n"; } } std::string mitk::TestDICOMLoading::trim(const std::string& pString, const std::string& pWhitespace) { const size_t beginStr = pString.find_first_not_of(pWhitespace); if (beginStr == std::string::npos) { // no content return ""; } const size_t endStr = pString.find_last_not_of(pWhitespace); const size_t range = endStr - beginStr + 1; return pString.substr(beginStr, range); } std::string mitk::TestDICOMLoading::reduce(const std::string& pString, const std::string& pFill, const std::string& pWhitespace) { // trim first std::string result(trim(pString, pWhitespace)); // replace sub ranges size_t beginSpace = result.find_first_of(pWhitespace); while (beginSpace != std::string::npos) { const size_t endSpace = result.find_first_not_of(pWhitespace, beginSpace); const size_t range = endSpace - beginSpace; result.replace(beginSpace, range, pFill); const size_t newStart = beginSpace + pFill.length(); beginSpace = result.find_first_of(pWhitespace, newStart); } return result; } bool mitk::TestDICOMLoading::CompareSpacedValueFields( const std::string& reference, const std::string& test, double /*eps*/ ) { bool result(true); // tokenize string, compare each token, if possible by float comparison std::stringstream referenceStream(reduce(reference)); std::stringstream testStream(reduce(test)); std::string refToken; std::string testToken; while ( std::getline( referenceStream, refToken, ' ' ) && std::getline ( testStream, testToken, ' ' ) ) { float refNumber; float testNumber; if ( this->StringToNumber(refToken, refNumber) ) { if ( this->StringToNumber(testToken, testNumber) ) { // print-out compared tokens if DEBUG output allowed MITK_DEBUG << "Reference Token '" << refToken << "'" << " value " << refNumber << ", test Token '" << testToken << "'" << " value " << testNumber; bool old_result = result; result &= ( std::abs(refNumber - testNumber) < 0.0001f /*mitk::eps*/ ); // log the token/number which causes the test to fail if( old_result != result) { MITK_ERROR << std::setprecision(16) << "Reference Token '" << refToken << "'" << " value " << refNumber << ", test Token '" << testToken << "'" << " value " << testNumber; MITK_ERROR << "[FALSE] - difference: " << std::setprecision(16) << std::abs(refNumber - testNumber) << " EPS: " << 0.0001f; //mitk::eps; } } else { MITK_ERROR << refNumber << " cannot be compared to '" << testToken << "'"; } } else { MITK_DEBUG << "Token '" << refToken << "'" << " handled as string"; result &= refToken == testToken; } } if ( std::getline( referenceStream, refToken, ' ' ) ) { MITK_ERROR << "Reference string still had values when test string was already parsed: ref '" << reference << "', test '" << test << "'"; result = false; } else if ( std::getline( testStream, testToken, ' ' ) ) { MITK_ERROR << "Test string still had values when reference string was already parsed: ref '" << reference << "', test '" << test << "'"; result = false; } return result; } bool mitk::TestDICOMLoading::CompareImageInformationDumps( const std::string& referenceDump, const std::string& testDump ) { KeyValueMap reference = ParseDump(referenceDump); KeyValueMap test = ParseDump(testDump); bool testResult(true); // verify all expected values for (KeyValueMap::const_iterator refIter = reference.begin(); refIter != reference.end(); ++refIter) { const std::string& refKey = refIter->first; const std::string& refValue = refIter->second; if ( test.find(refKey) != test.end() ) { const std::string& testValue = test[refKey]; if (refKey == mitk::PropertyKeyPathToPropertyName(mitk::DICOMIOMetaInformationPropertyConstants::READER_DCMTK())) { //check dcmtk version always against the current version of the system bool thisTestResult = testValue == std::string(" ") + PACKAGE_VERSION; testResult &= thisTestResult; MITK_DEBUG << refKey << ": '" << PACKAGE_VERSION << "' == '" << testValue << "' ? " << (thisTestResult ? "YES" : "NO"); } else if (refKey == mitk::PropertyKeyPathToPropertyName(mitk::DICOMIOMetaInformationPropertyConstants::READER_GDCM())) {//check gdcm version always against the current version of the system bool thisTestResult = testValue == std::string(" ") + gdcm::Version::GetVersion(); testResult &= thisTestResult; MITK_DEBUG << refKey << ": '" << gdcm::Version::GetVersion() << "' == '" << testValue << "' ? " << (thisTestResult ? "YES" : "NO"); } else { bool thisTestResult = CompareSpacedValueFields(refValue, testValue); testResult &= thisTestResult; MITK_DEBUG << refKey << ": '" << refValue << "' == '" << testValue << "' ? " << (thisTestResult ? "YES" : "NO"); } } else { MITK_ERROR << "Reference dump contains a key'" << refKey << "' (value '" << refValue << "')." ; MITK_ERROR << "This key is expected to be generated for tests (but was not). Most probably you need to update your test data."; return false; } } // now check test dump does not contain any additional keys for (KeyValueMap::const_iterator testIter = test.begin(); testIter != test.end(); ++testIter) { const std::string& key = testIter->first; const std::string& value = testIter->second; if (key == mitk::PropertyKeyPathToPropertyName(mitk::DICOMIOMetaInformationPropertyConstants::READER_DCMTK())) {//check dcmtk version always against the current version of the system bool thisTestResult = value == std::string(" ")+PACKAGE_VERSION; testResult &= thisTestResult; MITK_DEBUG << key << ": '" << PACKAGE_VERSION << "' == '" << value << "' ? " << (thisTestResult ? "YES" : "NO"); } else if (key == mitk::PropertyKeyPathToPropertyName(mitk::DICOMIOMetaInformationPropertyConstants::READER_GDCM())) {//check gdcm version always against the current version of the system bool thisTestResult = value == std::string(" ") + gdcm::Version::GetVersion(); testResult &= thisTestResult; MITK_DEBUG << key << ": '" << gdcm::Version::GetVersion() << "' == '" << value << "' ? " << (thisTestResult ? "YES" : "NO"); } else if ( reference.find(key) == reference.end() ) { MITK_ERROR << "Test dump contains an unexpected key'" << key << "' (value '" << value << "')." ; MITK_ERROR << "This key is not expected. Most probably you need to update your test data."; return false; } } return testResult; } mitk::TestDICOMLoading::KeyValueMap mitk::TestDICOMLoading::ParseDump( const std::string& dump ) { KeyValueMap parsedResult; std::string shredder(dump); std::stack surroundingKeys; std::stack expectedIndents; expectedIndents.push(0); while (true) { std::string::size_type newLinePos = shredder.find( '\n' ); if (newLinePos == std::string::npos || newLinePos == 0) break; std::string line = shredder.substr( 0, newLinePos ); shredder = shredder.erase( 0, newLinePos+1 ); std::string::size_type keyPosition = line.find_first_not_of( ' ' ); std::string::size_type colonPosition = line.find( ':' ); std::string key = line.substr(keyPosition, colonPosition - keyPosition); std::string::size_type firstSpacePosition = key.find_first_of(" "); if (firstSpacePosition != std::string::npos) { key.erase(firstSpacePosition); } if ( keyPosition > expectedIndents.top() ) { // more indent than before expectedIndents.push(keyPosition); } else { if (!surroundingKeys.empty()) { surroundingKeys.pop(); // last of same length } while (expectedIndents.top() != keyPosition) { expectedIndents.pop(); if (!surroundingKeys.empty()) { surroundingKeys.pop(); } }; // unwind until current indent is found } if (!surroundingKeys.empty()) { key = surroundingKeys.top() + "." + key; // construct current key name } surroundingKeys.push(key); // this is the new embracing key std::string value = line.substr(colonPosition+1); MITK_DEBUG << " Key: '" << key << "' value '" << value << "'" ; parsedResult[key] = value; // store parsing result } return parsedResult; } diff --git a/Modules/ImageDenoising/Testing/itkTotalVariationDenoisingImageFilterTest.cpp b/Modules/ImageDenoising/Testing/itkTotalVariationDenoisingImageFilterTest.cpp index ac6c9ca227..2c8888cc59 100644 --- a/Modules/ImageDenoising/Testing/itkTotalVariationDenoisingImageFilterTest.cpp +++ b/Modules/ImageDenoising/Testing/itkTotalVariationDenoisingImageFilterTest.cpp @@ -1,271 +1,271 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "itkImageRegionIterator.h" #include "itkLocalVariationImageFilter.h" #include "itkTotalVariationDenoisingImageFilter.h" #include "itkTotalVariationSingleIterationImageFilter.h" // image typedefs typedef itk::Image ImageType; typedef itk::ImageRegionIterator IteratorType; // vector image typedefs typedef itk::Vector VectorPixelType; typedef itk::Image VectorImageType; typedef itk::ImageRegionIterator VectorIteratorType; /** * 3x3x3 test image */ ImageType::Pointer GenerateTestImage() { // init ImageType::Pointer image = ImageType::New(); ; // spacing ImageType::SpacingType spacing; spacing[0] = 1; spacing[1] = 1; spacing[2] = 1; image->SetSpacing(spacing); // extent ImageType::RegionType largestPossibleRegion; ImageType::SizeType size = {{3, 3, 1}}; largestPossibleRegion.SetSize(size); ImageType::IndexType index = {{0, 0, 0}}; largestPossibleRegion.SetIndex(index); image->SetLargestPossibleRegion(largestPossibleRegion); image->SetBufferedRegion(largestPossibleRegion); // allocate memory image->Allocate(); int i = 0; IteratorType it(image, largestPossibleRegion); it.GoToBegin(); while (!it.IsAtEnd()) { it.Set((float)i++); ++it; } return image; } VectorImageType::Pointer GenerateVectorTestImage() { // init VectorImageType::Pointer image = VectorImageType::New(); ; // spacing VectorImageType::SpacingType spacing; spacing[0] = 1; spacing[1] = 1; spacing[2] = 1; image->SetSpacing(spacing); // extent VectorImageType::RegionType largestPossibleRegion; VectorImageType::SizeType size = {{3, 3, 1}}; largestPossibleRegion.SetSize(size); VectorImageType::IndexType index = {{0, 0, 0}}; largestPossibleRegion.SetIndex(index); image->SetLargestPossibleRegion(largestPossibleRegion); image->SetBufferedRegion(largestPossibleRegion); // allocate memory image->Allocate(); int i = 0; VectorIteratorType it(image, largestPossibleRegion); it.GoToBegin(); while (!it.IsAtEnd()) { VectorPixelType vec; vec[0] = (float)i; vec[1] = (float)i++; it.Set(vec); ++it; } return image; } void PrintImage(ImageType::Pointer image) { IteratorType it(image, image->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) { std::cout << it.Get() << " "; } std::cout << std::endl; } void PrintVectorImage(VectorImageType::Pointer image) { VectorIteratorType it(image, image->GetLargestPossibleRegion()); for (it.GoToBegin(); !it.IsAtEnd(); ++it) { std::cout << it.Get() << " "; } std::cout << std::endl; } /** * todo */ int itkTotalVariationDenoisingImageFilterTest(int /*argc*/, char * /*argv*/ []) { ImageType::Pointer image = GenerateTestImage(); PrintImage(image); double precision = 0.01; ImageType::IndexType index = {{1, 1, 0}}; VectorImageType::IndexType vecIndex = {{1, 1, 0}}; try { typedef itk::LocalVariationImageFilter LocalFilterType; LocalFilterType::Pointer filter = LocalFilterType::New(); filter->SetInput(image); - filter->SetNumberOfThreads(1); + filter->SetNumberOfWorkUnits(1); filter->Update(); ImageType::Pointer outImage = filter->GetOutput(); PrintImage(outImage); if (fabs(outImage->GetPixel(index) - 4.472) > precision) { return EXIT_FAILURE; } } catch (...) { return EXIT_FAILURE; } try { typedef itk::TotalVariationSingleIterationImageFilter SingleFilterType; SingleFilterType::Pointer sFilter = SingleFilterType::New(); sFilter->SetInput(image); sFilter->SetOriginalImage(GenerateTestImage()); sFilter->SetLambda(0.5); - sFilter->SetNumberOfThreads(1); + sFilter->SetNumberOfWorkUnits(1); sFilter->Update(); ImageType::Pointer outImageS = sFilter->GetOutput(); PrintImage(outImageS); if (fabs(outImageS->GetPixel(index) - 4.0) > precision) { return EXIT_FAILURE; } } catch (...) { return EXIT_FAILURE; } try { typedef itk::TotalVariationDenoisingImageFilter TVFilterType; TVFilterType::Pointer tvFilter = TVFilterType::New(); tvFilter->SetInput(image); tvFilter->SetNumberIterations(30); - tvFilter->SetNumberOfThreads(1); + tvFilter->SetNumberOfWorkUnits(1); tvFilter->SetLambda(0.1); tvFilter->Update(); ImageType::Pointer outImageTV = tvFilter->GetOutput(); PrintImage(outImageTV); if (fabs(outImageTV->GetPixel(index) - 4.0) > precision) { return EXIT_FAILURE; } } catch (...) { return EXIT_FAILURE; } VectorImageType::Pointer vecImage = GenerateVectorTestImage(); PrintVectorImage(vecImage); try { typedef itk::LocalVariationImageFilter LocalVecFilterType; LocalVecFilterType::Pointer vecFilter = LocalVecFilterType::New(); vecFilter->SetInput(vecImage); - vecFilter->SetNumberOfThreads(1); + vecFilter->SetNumberOfWorkUnits(1); vecFilter->Update(); ImageType::Pointer outVecImage = vecFilter->GetOutput(); PrintImage(outVecImage); if (fabs(outVecImage->GetPixel(index) - 6.324) > precision) { return EXIT_FAILURE; } } catch (...) { return EXIT_FAILURE; } try { typedef itk::TotalVariationSingleIterationImageFilter SingleVecFilterType; SingleVecFilterType::Pointer sVecFilter = SingleVecFilterType::New(); sVecFilter->SetInput(vecImage); sVecFilter->SetOriginalImage(vecImage); sVecFilter->SetLambda(0.5); - sVecFilter->SetNumberOfThreads(1); + sVecFilter->SetNumberOfWorkUnits(1); sVecFilter->UpdateLargestPossibleRegion(); VectorImageType::Pointer outVecImageS = sVecFilter->GetOutput(); PrintVectorImage(outVecImageS); if (fabs(outVecImageS->GetPixel(vecIndex)[1] - 4.0) > precision) { return EXIT_FAILURE; } } catch (...) { return EXIT_FAILURE; } try { typedef itk::TotalVariationDenoisingImageFilter TVVectorFilterType; TVVectorFilterType::Pointer tvVecFilter = TVVectorFilterType::New(); tvVecFilter->SetInput(vecImage); tvVecFilter->SetNumberIterations(30); - tvVecFilter->SetNumberOfThreads(1); + tvVecFilter->SetNumberOfWorkUnits(1); tvVecFilter->SetLambda(0.1); tvVecFilter->Update(); VectorImageType::Pointer outVecImageTV = tvVecFilter->GetOutput(); PrintVectorImage(outVecImageTV); if (fabs(outVecImageTV->GetPixel(vecIndex)[1] - 4.0) > precision) { return EXIT_FAILURE; } } catch (...) { return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/ImageDenoising/itkTotalVariationDenoisingImageFilter.txx b/Modules/ImageDenoising/itkTotalVariationDenoisingImageFilter.txx index d5b776752b..25d077d79e 100644 --- a/Modules/ImageDenoising/itkTotalVariationDenoisingImageFilter.txx +++ b/Modules/ImageDenoising/itkTotalVariationDenoisingImageFilter.txx @@ -1,105 +1,105 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef _itkTotalVariationDenoisingImageFilter_txx #define _itkTotalVariationDenoisingImageFilter_txx #include "itkTotalVariationDenoisingImageFilter.h" #include "itkConstShapedNeighborhoodIterator.h" #include "itkImageRegionConstIterator.h" #include "itkImageRegionIterator.h" #include "itkLocalVariationImageFilter.h" #include "itkNeighborhoodAlgorithm.h" #include "itkNeighborhoodInnerProduct.h" #include "itkOffset.h" #include "itkProgressReporter.h" #include "itkZeroFluxNeumannBoundaryCondition.h" #include #include namespace itk { template TotalVariationDenoisingImageFilter::TotalVariationDenoisingImageFilter() : m_Lambda(1.0), m_NumberIterations(0) { } template void TotalVariationDenoisingImageFilter::GenerateData() { // first we cast the input image to match output type typename CastType::Pointer infilter = CastType::New(); infilter->SetInput(this->GetInput()); infilter->Update(); typename TOutputImage::Pointer image = infilter->GetOutput(); // a second copy of the input image is saved as reference infilter = CastType::New(); infilter->SetInput(this->GetInput()); infilter->Update(); typename TOutputImage::Pointer origImage = infilter->GetOutput(); typename SingleIterationFilterType::Pointer filter; for (int i = 0; i < m_NumberIterations; i++) { filter = SingleIterationFilterType::New(); filter->SetInput(image); filter->SetOriginalImage(origImage); filter->SetLambda(m_Lambda); - filter->SetNumberOfThreads(this->GetNumberOfThreads()); + filter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); filter->UpdateLargestPossibleRegion(); image = filter->GetOutput(); std::cout << "Iteration " << i + 1 << "/" << m_NumberIterations << std::endl; } typename OutputImageType::Pointer output = this->GetOutput(); output->SetSpacing(image->GetSpacing()); typename OutputImageType::RegionType largestPossibleRegion; largestPossibleRegion.SetSize(image->GetLargestPossibleRegion().GetSize()); largestPossibleRegion.SetIndex(image->GetLargestPossibleRegion().GetIndex()); output->SetLargestPossibleRegion(image->GetLargestPossibleRegion()); output->SetBufferedRegion(image->GetLargestPossibleRegion()); output->Allocate(); itk::ImageRegionIterator oit(output, output->GetLargestPossibleRegion()); oit.GoToBegin(); itk::ImageRegionConstIterator iit(image, image->GetLargestPossibleRegion()); iit.GoToBegin(); while (!oit.IsAtEnd()) { oit.Set(iit.Get()); ++iit; ++oit; } } /** * Standard "PrintSelf" method */ template void TotalVariationDenoisingImageFilter::PrintSelf(std::ostream &os, Indent indent) const { Superclass::PrintSelf(os, indent); } } // end namespace itk #endif diff --git a/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx b/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx index a8fa6f6063..b409b2500f 100644 --- a/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx +++ b/Modules/ImageDenoising/itkTotalVariationSingleIterationImageFilter.txx @@ -1,251 +1,251 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef _itkTotalVariationSingleIterationImageFilter_txx #define _itkTotalVariationSingleIterationImageFilter_txx #include "itkTotalVariationSingleIterationImageFilter.h" // itk includes #include "itkConstShapedNeighborhoodIterator.h" #include "itkImageRegionIterator.h" #include "itkLocalVariationImageFilter.h" #include "itkNeighborhoodAlgorithm.h" #include "itkNeighborhoodInnerProduct.h" #include "itkOffset.h" #include "itkProgressReporter.h" #include "itkZeroFluxNeumannBoundaryCondition.h" // other includes #include #include namespace itk { /** * constructor */ template TotalVariationSingleIterationImageFilter::TotalVariationSingleIterationImageFilter() { m_Lambda = 1.0; m_LocalVariation = LocalVariationImageType::New(); } /** * generate requested region */ template void TotalVariationSingleIterationImageFilter::GenerateInputRequestedRegion() { // call the superclass' implementation of this method Superclass::GenerateInputRequestedRegion(); // get pointers to the input and output typename Superclass::InputImagePointer inputPtr = const_cast(this->GetInput()); typename Superclass::OutputImagePointer outputPtr = this->GetOutput(); if (!inputPtr || !outputPtr) { return; } // get a copy of the input requested region (should equal the output // requested region) typename TInputImage::RegionType inputRequestedRegion; inputRequestedRegion = inputPtr->GetRequestedRegion(); // pad the input requested region by 1 inputRequestedRegion.PadByRadius(1); // crop the input requested region at the input's largest possible region if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion())) { inputPtr->SetRequestedRegion(inputRequestedRegion); return; } else { // Couldn't crop the region (requested region is outside the largest // possible region). Throw an exception. // store what we tried to request (prior to trying to crop) inputPtr->SetRequestedRegion(inputRequestedRegion); // build an exception InvalidRequestedRegionError e(__FILE__, __LINE__); e.SetLocation(ITK_LOCATION); e.SetDescription("Requested region outside possible region."); e.SetDataObject(inputPtr); throw e; } } /** * generate output */ template void TotalVariationSingleIterationImageFilter::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId) { typename OutputImageType::Pointer output = this->GetOutput(); typename InputImageType::ConstPointer input = this->GetInput(); // Find the data-set boundary "faces" itk::Size size; for (unsigned int i = 0; i < InputImageDimension; i++) size[i] = 1; NeighborhoodAlgorithm::ImageBoundaryFacesCalculator bC; typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType faceList = bC(input, outputRegionForThread, size); NeighborhoodAlgorithm::ImageBoundaryFacesCalculator lv_bC; typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType lv_faceList = lv_bC(m_LocalVariation, outputRegionForThread, size); // support progress methods/callbacks ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); ZeroFluxNeumannBoundaryCondition nbc; ZeroFluxNeumannBoundaryCondition lv_nbc; std::vector ws; std::vector hs; auto lv_fit = lv_faceList.begin(); // Process each of the boundary faces. These are N-d regions which border // the edge of the buffer. for (auto fit = faceList.begin(); fit != faceList.end(); ++fit) { // iterators over output, input, original and local variation image ImageRegionIterator output_image_it = ImageRegionIterator(output, *fit); ImageRegionConstIterator input_image_it = ImageRegionConstIterator(input, *fit); ImageRegionConstIterator orig_image_it = ImageRegionConstIterator(m_OriginalImage, *fit); ImageRegionConstIterator loc_var_image_it = ImageRegionConstIterator(m_LocalVariation, *fit); // neighborhood in input image ConstShapedNeighborhoodIterator input_image_neighbors_it(size, input, *fit); typename ConstShapedNeighborhoodIterator::OffsetType offset; input_image_neighbors_it.OverrideBoundaryCondition(&nbc); input_image_neighbors_it.ClearActiveList(); for (unsigned int i = 0; i < InputImageDimension; i++) { offset.Fill(0); offset[i] = -1; input_image_neighbors_it.ActivateOffset(offset); offset[i] = 1; input_image_neighbors_it.ActivateOffset(offset); } input_image_neighbors_it.GoToBegin(); // neighborhood in local variation image ConstShapedNeighborhoodIterator loc_var_image_neighbors_it( size, m_LocalVariation, *lv_fit); loc_var_image_neighbors_it.OverrideBoundaryCondition(&lv_nbc); loc_var_image_neighbors_it.ClearActiveList(); for (unsigned int i = 0; i < InputImageDimension; i++) { offset.Fill(0); offset[i] = -1; loc_var_image_neighbors_it.ActivateOffset(offset); offset[i] = 1; loc_var_image_neighbors_it.ActivateOffset(offset); } loc_var_image_neighbors_it.GoToBegin(); const unsigned int neighborhoodSize = InputImageDimension * 2; ws.resize(neighborhoodSize); while (!output_image_it.IsAtEnd()) { // 1 / ||nabla_alpha(u)||_a double locvar_alpha_inv = 1.0 / loc_var_image_it.Get(); // compute w_alphabetas int count = 0; double wsum = 0; typename ConstShapedNeighborhoodIterator::ConstIterator loc_var_neighbors_it; for (loc_var_neighbors_it = loc_var_image_neighbors_it.Begin(); !loc_var_neighbors_it.IsAtEnd(); loc_var_neighbors_it++) { // w_alphabeta(u) = // 1 / ||nabla_alpha(u)||_a + 1 / ||nabla_beta(u)||_a ws[count] = locvar_alpha_inv + (1.0 / (double)loc_var_neighbors_it.Get()); wsum += ws[count++]; } // h_alphaalpha * u_alpha^zero typename OutputImageType::PixelType res = static_cast( ((typename OutputImageType::PixelType)orig_image_it.Get()) * (m_Lambda / (m_Lambda + wsum))); // add the different h_alphabeta * u_beta count = 0; typename ConstShapedNeighborhoodIterator::ConstIterator input_neighbors_it; for (input_neighbors_it = input_image_neighbors_it.Begin(); !input_neighbors_it.IsAtEnd(); input_neighbors_it++) { res += input_neighbors_it.Get() * (ws[count++] / (m_Lambda + wsum)); } // set output result output_image_it.Set(res); // increment iterators ++output_image_it; ++input_image_it; ++orig_image_it; ++loc_var_image_it; ++input_image_neighbors_it; ++loc_var_image_neighbors_it; // report progress progress.CompletedPixel(); } ++lv_fit; } } /** * first calculate local variation in the image */ template void TotalVariationSingleIterationImageFilter::BeforeThreadedGenerateData() { typedef typename itk::LocalVariationImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); filter->SetInput(this->GetInput(0)); - filter->SetNumberOfThreads(this->GetNumberOfThreads()); + filter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); filter->Update(); this->m_LocalVariation = filter->GetOutput(); } /** * Standard "PrintSelf" method */ template void TotalVariationSingleIterationImageFilter::PrintSelf(std::ostream &os, Indent indent) const { Superclass::PrintSelf(os, indent); } } // end namespace itk #endif diff --git a/Modules/ModelFit/test/itkMaskedStatisticsImageFilterTest.cpp b/Modules/ModelFit/test/itkMaskedStatisticsImageFilterTest.cpp index 5a00710a8e..590937f13f 100644 --- a/Modules/ModelFit/test/itkMaskedStatisticsImageFilterTest.cpp +++ b/Modules/ModelFit/test/itkMaskedStatisticsImageFilterTest.cpp @@ -1,76 +1,76 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "itkImage.h" #include "itkImageRegionIterator.h" #include "itkMaskedStatisticsImageFilter.h" #include "mitkTestingMacros.h" #include "mitkVector.h" #include "mitkTestDynamicImageGenerator.h" int itkMaskedStatisticsImageFilterTest(int /*argc*/, char*[] /*argv[]*/) { // always start with this! MITK_TEST_BEGIN("itkMaskedStatisticsImageFilterTest") //Prepare test artifacts and helper mitk::TestImageType::Pointer img1 = mitk::GenerateTestImage(); typedef itk::MaskedStatisticsImageFilter FilterType; FilterType::Pointer testFilter = FilterType::New(); testFilter->SetInput(img1); - testFilter->SetNumberOfThreads(2); + testFilter->SetNumberOfWorkUnits(2); testFilter->Update(); FilterType::PixelType max = testFilter->GetMaximum(); FilterType::PixelType min = testFilter->GetMinimum(); FilterType::RealType mean = testFilter->GetMean(); FilterType::RealType sig = testFilter->GetSigma(); FilterType::RealType variance = testFilter->GetVariance(); FilterType::RealType sum = testFilter->GetSum(); CPPUNIT_ASSERT_MESSAGE("Check computed maximum",9 == max); CPPUNIT_ASSERT_MESSAGE("Check computed minimum",1 == min); CPPUNIT_ASSERT_MESSAGE("Check computed mean",5 == mean); CPPUNIT_ASSERT_MESSAGE("Check computed sigma",sqrt(7.5) == sig); CPPUNIT_ASSERT_MESSAGE("Check computed variance",7.5 == variance); CPPUNIT_ASSERT_MESSAGE("Check computed sum",45 == sum); //Test with mask set mitk::TestMaskType::Pointer mask = mitk::GenerateTestMask(); testFilter->SetMask(mask); testFilter->Update(); max = testFilter->GetMaximum(); min = testFilter->GetMinimum(); mean = testFilter->GetMean(); sig = testFilter->GetSigma(); variance = testFilter->GetVariance(); sum = testFilter->GetSum(); CPPUNIT_ASSERT_MESSAGE("Check computed maximum",4 == max); CPPUNIT_ASSERT_MESSAGE("Check computed minimum",2 == min); CPPUNIT_ASSERT_MESSAGE("Check computed mean",3 == mean); CPPUNIT_ASSERT_MESSAGE("Check computed sigma",1 == sig); CPPUNIT_ASSERT_MESSAGE("Check computed variance",1 == variance); CPPUNIT_ASSERT_MESSAGE("Check computed sum",9 == sum); MITK_TEST_END() } diff --git a/Modules/ModelFit/test/itkMultiOutputNaryFunctorImageFilterTest.cpp b/Modules/ModelFit/test/itkMultiOutputNaryFunctorImageFilterTest.cpp index f38866b609..4ca0102a9e 100644 --- a/Modules/ModelFit/test/itkMultiOutputNaryFunctorImageFilterTest.cpp +++ b/Modules/ModelFit/test/itkMultiOutputNaryFunctorImageFilterTest.cpp @@ -1,218 +1,218 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "itkImage.h" #include "itkImageRegionIterator.h" #include "itkMultiOutputNaryFunctorImageFilter.h" #include "mitkTestingMacros.h" #include "mitkVector.h" #include "mitkTestDynamicImageGenerator.h" class TestFunctor { public: typedef std::vector InputPixelArrayType; typedef std::vector OutputPixelArrayType; typedef itk::Index<2> IndexType; TestFunctor() { secondOutputSelection = 0; }; ~TestFunctor() {}; int secondOutputSelection; unsigned int GetNumberOfOutputs() const { return 4; } bool operator!=( const TestFunctor & other) const { return !(*this == other); } bool operator==( const TestFunctor & other ) const { return secondOutputSelection == other.secondOutputSelection; } inline OutputPixelArrayType operator()( const InputPixelArrayType & value, const IndexType& currentIndex ) const { OutputPixelArrayType result; int sum = 0; for (InputPixelArrayType::const_iterator pos = value.begin(); pos != value.end(); ++pos) { sum += *pos; } result.push_back(sum); result.push_back(value[secondOutputSelection]); result.push_back(currentIndex[0]); result.push_back(currentIndex[1]); return result; } }; int itkMultiOutputNaryFunctorImageFilterTest(int /*argc*/, char*[] /*argv[]*/) { // always start with this! MITK_TEST_BEGIN("itkMultiOutputNaryFunctorImageFilter") //Prepare test artifacts and helper mitk::TestImageType::Pointer img1 = mitk::GenerateTestImage(); mitk::TestImageType::Pointer img2 = mitk::GenerateTestImage(10); mitk::TestImageType::Pointer img3 = mitk::GenerateTestImage(100); mitk::TestImageType::IndexType testIndex1; testIndex1[0] = 0; testIndex1[1] = 0; mitk::TestImageType::IndexType testIndex2; testIndex2[0] = 2; testIndex2[1] = 0; mitk::TestImageType::IndexType testIndex3; testIndex3[0] = 0; testIndex3[1] = 1; mitk::TestImageType::IndexType testIndex4; testIndex4[0] = 1; testIndex4[1] = 1; mitk::TestImageType::IndexType testIndex5; testIndex5[0] = 2; testIndex5[1] = 2; //Test default usage of filter typedef itk::MultiOutputNaryFunctorImageFilter FilterType; FilterType::Pointer testFilter = FilterType::New(); testFilter->SetInput(0,img1); testFilter->SetInput(1,img2); testFilter->SetInput(2,img3); - testFilter->SetNumberOfThreads(2); + testFilter->SetNumberOfWorkUnits(2); testFilter->Update(); mitk::TestImageType::Pointer out1 = testFilter->GetOutput(0); mitk::TestImageType::Pointer out2 = testFilter->GetOutput(1); mitk::TestImageType::Pointer out3 = testFilter->GetOutput(2); mitk::TestImageType::Pointer out4 = testFilter->GetOutput(3); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #1 index #1 (functor #1)",111 == out1->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #1 index #2 (functor #1)",333 == out1->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #1 index #3 (functor #1)",444 == out1->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #1 index #4 (functor #1)",555 == out1->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #1 index #5 (functor #1)",999 == out1->GetPixel(testIndex5)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #2 index #1 (functor #1)",1 == out2->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #2 index #2 (functor #1)",3 == out2->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #2 index #3 (functor #1)",4 == out2->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #2 index #4 (functor #1)",5 == out2->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #2 index #5 (functor #1)",9 == out2->GetPixel(testIndex5)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #3 index #1 (functor #1)",0 == out3->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #3 index #2 (functor #1)",2 == out3->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #3 index #3 (functor #1)",0 == out3->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #3 index #4 (functor #1)",1 == out3->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #3 index #5 (functor #1)",2 == out3->GetPixel(testIndex5)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #4 index #1 (functor #1)",0 == out4->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #4 index #2 (functor #1)",0 == out4->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #4 index #3 (functor #1)",1 == out4->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #4 index #4 (functor #1)",1 == out4->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #4 index #5 (functor #1)",2 == out4->GetPixel(testIndex5)); //Test with functor set by user TestFunctor funct2; funct2.secondOutputSelection = 1; testFilter->SetFunctor(funct2); testFilter->Update(); out1 = testFilter->GetOutput(0); out2 = testFilter->GetOutput(1); out3 = testFilter->GetOutput(2); out4 = testFilter->GetOutput(3); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #1 index #1 (functor #2)",111 == out1->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #1 index #2 (functor #2)",333 == out1->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #1 index #3 (functor #2)",444 == out1->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #1 index #4 (functor #2)",555 == out1->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #1 index #5 (functor #2)",999 == out1->GetPixel(testIndex5)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #2 index #1 (functor #2)",10 == out2->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #2 index #2 (functor #2)",30 == out2->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #2 index #3 (functor #2)",40 == out2->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #2 index #4 (functor #2)",50 == out2->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #2 index #5 (functor #2)",90 == out2->GetPixel(testIndex5)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #3 index #1 (functor #2)",0 == out3->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #3 index #2 (functor #2)",2 == out3->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #3 index #3 (functor #2)",0 == out3->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #3 index #4 (functor #2)",1 == out3->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #3 index #5 (functor #2)",2 == out3->GetPixel(testIndex5)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #4 index #1 (functor #2)",0 == out4->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #4 index #2 (functor #2)",0 == out4->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #4 index #3 (functor #2)",1 == out4->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #4 index #4 (functor #2)",1 == out4->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of output #4 index #5 (functor #2)",2 == out4->GetPixel(testIndex5)); //Test with mask set mitk::TestMaskType::Pointer mask = mitk::GenerateTestMask(); testFilter->SetMask(mask); testFilter->Update(); out1 = testFilter->GetOutput(0); out2 = testFilter->GetOutput(1); out3 = testFilter->GetOutput(2); out4 = testFilter->GetOutput(3); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #1 index #1 (functor #2)",0 == out1->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #1 index #2 (functor #2)",333 == out1->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #1 index #3 (functor #2)",444 == out1->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #1 index #4 (functor #2)",0 == out1->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #1 index #5 (functor #2)",0 == out1->GetPixel(testIndex5)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #2 index #1 (functor #2)",0 == out2->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #2 index #2 (functor #2)",30 == out2->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #2 index #3 (functor #2)",40 == out2->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #2 index #4 (functor #2)",0 == out2->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #2 index #5 (functor #2)",0 == out2->GetPixel(testIndex5)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #3 index #1 (functor #2)",0 == out3->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #3 index #2 (functor #2)",2 == out3->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #3 index #3 (functor #2)",0 == out3->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #3 index #4 (functor #2)",0 == out3->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #3 index #5 (functor #2)",0 == out3->GetPixel(testIndex5)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #4 index #1 (functor #2)",0 == out4->GetPixel(testIndex1)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #4 index #2 (functor #2)",0 == out4->GetPixel(testIndex2)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #4 index #3 (functor #2)",1 == out4->GetPixel(testIndex3)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #4 index #4 (functor #2)",0 == out4->GetPixel(testIndex4)); CPPUNIT_ASSERT_MESSAGE("Check pixel of masked output #4 index #5 (functor #2)",0 == out4->GetPixel(testIndex5)); MITK_TEST_END() } diff --git a/Modules/TubeGraph/include/mitkTubeGraphDataInteractor.h b/Modules/TubeGraph/include/mitkTubeGraphDataInteractor.h index 8d07e4f3d8..622ed78781 100644 --- a/Modules/TubeGraph/include/mitkTubeGraphDataInteractor.h +++ b/Modules/TubeGraph/include/mitkTubeGraphDataInteractor.h @@ -1,122 +1,122 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkTubeGraphDataInteractor3D_h_ #define mitkTubeGraphDataInteractor3D_h_ #include #include #include #include "mitkTubeGraph.h" #include "mitkTubeGraphProperty.h" namespace mitk { // Define events for TubeGraph interaction notifications - itkEventMacro(SelectionChangedTubeGraphEvent, itk::AnyEvent); + itkEventMacroDeclaration(SelectionChangedTubeGraphEvent, itk::AnyEvent); /** * \brief * * \ingroup Interaction */ // Inherit from DataInteratcor, this provides functionality of a state machine and configurable inputs. class MITKTUBEGRAPH_EXPORT TubeGraphDataInteractor : public DataInteractor { public: mitkClassMacro(TubeGraphDataInteractor, DataInteractor); itkNewMacro(Self); /** * Describes, which activation modes are available based on the * currently picked tube: * * \li None means "no tube is active" * \li Single means "only the picked tube is active" * \li ToRoot means "all tubes from the picked on down to the root of the tube graph are active" * \li ToPeriphery means "all tubes included in the subgraph of the currently picked vessel are active" * \li Points means "shortes path between two picked tubes are active" * \li Multiple means "all picked tubes are active" */ enum ActivationMode { None = 0, Single, ToRoot, ToPeriphery, Points, Multiple }; enum ActionMode { AttributationMode = 0, AnnotationMode, EditMode, RootMode, InformationMode }; void SetActivationMode(const ActivationMode &activationMode); ActivationMode GetActivationMode(); void SetActionMode(const ActionMode &actionMode); ActionMode GetActionMode(); void ResetPickedTubes(); mitk::Point3D GetLastPickedPosition(); protected: TubeGraphDataInteractor(); ~TubeGraphDataInteractor() override; /** * Here actions strings from the loaded state machine pattern are mapped to functions of * the DataInteractor. These functions are called when an action from the state machine pattern is executed. */ void ConnectActionsAndFunctions() override; /** * This function is called when a DataNode has been set/changed. */ void DataNodeChanged() override; /** * Initializes the movement, stores starting position. */ virtual bool CheckOverTube(const InteractionEvent *); virtual void SelectTube(StateMachineAction *, InteractionEvent *); virtual void DeselectTube(StateMachineAction *, InteractionEvent *); void SelectTubesByActivationModus(); void UpdateActivation(); private: std::vector GetTubesToRoot(); std::vector GetTubesBetweenPoints(); std::vector GetPathToPeriphery(); std::vector GetPathBetweenTubes(const TubeGraph::TubeDescriptorType &start, const TubeGraph::TubeDescriptorType &end); TubeGraph::Pointer m_TubeGraph; TubeGraphProperty::Pointer m_TubeGraphProperty; TubeGraph::TubeDescriptorType m_LastPickedTube; TubeGraph::TubeDescriptorType m_SecondLastPickedTube; ActivationMode m_ActivationMode; ActionMode m_ActionMode; mitk::TubeElement *m_LastPickedElement = nullptr; }; } #endif diff --git a/Modules/TubeGraph/src/Interactions/mitkTubeGraphDataInteractor.cpp b/Modules/TubeGraph/src/Interactions/mitkTubeGraphDataInteractor.cpp index 55cab1dbcb..bd0af564dd 100644 --- a/Modules/TubeGraph/src/Interactions/mitkTubeGraphDataInteractor.cpp +++ b/Modules/TubeGraph/src/Interactions/mitkTubeGraphDataInteractor.cpp @@ -1,285 +1,290 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkTubeGraphDataInteractor.h" #include #include #include #include "mitkTubeGraphPicker.h" #include #include #include #include #include +namespace mitk +{ + itkEventMacroDefinition(SelectionChangedTubeGraphEvent, itk::AnyEvent); +} + mitk::TubeGraphDataInteractor::TubeGraphDataInteractor() : m_LastPickedTube(TubeGraph::ErrorId), m_SecondLastPickedTube(TubeGraph::ErrorId), m_ActivationMode(None), m_ActionMode(AttributationMode) { } mitk::TubeGraphDataInteractor::~TubeGraphDataInteractor() {} void mitk::TubeGraphDataInteractor::ConnectActionsAndFunctions() { // **Conditions** that can be used in the state machine, to ensure that certain conditions are met, before actually // executing an action CONNECT_CONDITION("isOverTube", CheckOverTube); // **Function** in the statmachine patterns also referred to as **Actions** CONNECT_FUNCTION("selectTube", SelectTube); CONNECT_FUNCTION("deselectTube", DeselectTube); } void mitk::TubeGraphDataInteractor::DataNodeChanged() { if (GetDataNode() != nullptr) { if (GetDataNode()->GetData() != nullptr) { m_TubeGraph = dynamic_cast(GetDataNode()->GetData()); m_TubeGraphProperty = dynamic_cast( m_TubeGraph->GetProperty("Tube Graph.Visualization Information").GetPointer()); if (m_TubeGraphProperty.IsNull()) MITK_ERROR << "Something went wrong! No tube graph property!"; } else m_TubeGraph = nullptr; } else m_TubeGraph = nullptr; } bool mitk::TubeGraphDataInteractor::CheckOverTube(const InteractionEvent *interactionEvent) { const auto *positionEvent = dynamic_cast(interactionEvent); if (positionEvent == nullptr) return false; auto *picker = new mitk::TubeGraphPicker(); picker->SetTubeGraph(m_TubeGraph); auto pickedTube = picker->GetPickedTube(positionEvent->GetPositionInWorld()); TubeGraph::TubeDescriptorType tubeDescriptor = pickedTube.first; if (tubeDescriptor != TubeGraph::ErrorId) { m_LastPickedElement = pickedTube.second; m_SecondLastPickedTube = m_LastPickedTube; m_LastPickedTube = tubeDescriptor; return true; } else // nothing picked return false; } void mitk::TubeGraphDataInteractor::SelectTube(StateMachineAction *, InteractionEvent *) { if (m_TubeGraph.IsNull()) return; this->SelectTubesByActivationModus(); RenderingManager::GetInstance()->RequestUpdateAll(); if (m_ActivationMode != None) { // show tube id on status bar std::stringstream displayText; displayText << "Picked tube: ID [" << m_LastPickedTube.first << "," << m_LastPickedTube.second << "]"; StatusBar::GetInstance()->DisplayText(displayText.str().c_str()); // TODO!!! this->InvokeEvent(SelectionChangedTubeGraphEvent()); } } void mitk::TubeGraphDataInteractor::DeselectTube(StateMachineAction *, InteractionEvent *) { if (m_TubeGraph.IsNull()) return; if ((m_ActivationMode != Multiple) && (m_ActivationMode != Points)) { m_TubeGraphProperty->DeactivateAllTubes(); RenderingManager::GetInstance()->RequestUpdateAll(); // TODO!!!this->InvokeEvent(SelectionChangedTubeGraphEvent()); } // show info on status bar StatusBar::GetInstance()->DisplayText("No tube hit!"); } void mitk::TubeGraphDataInteractor::SetActivationMode(const ActivationMode &activationMode) { m_ActivationMode = activationMode; if (m_TubeGraph.IsNotNull()) if (m_LastPickedTube != mitk::TubeGraph::ErrorId) this->UpdateActivation(); } mitk::TubeGraphDataInteractor::ActivationMode mitk::TubeGraphDataInteractor::GetActivationMode() { return m_ActivationMode; } void mitk::TubeGraphDataInteractor::SetActionMode(const ActionMode &actionMode) { m_ActionMode = actionMode; } mitk::TubeGraphDataInteractor::ActionMode mitk::TubeGraphDataInteractor::GetActionMode() { return m_ActionMode; } void mitk::TubeGraphDataInteractor::SelectTubesByActivationModus() { if (m_LastPickedTube != mitk::TubeGraph::ErrorId) { this->UpdateActivation(); } } void mitk::TubeGraphDataInteractor::UpdateActivation() { if (m_ActionMode == RootMode) { m_TubeGraphProperty->DeactivateAllTubes(); m_TubeGraphProperty->SetTubeActive(m_LastPickedTube, true); // QmitkTubeGraphSelectRootDialog* dialog = new QmitkTubeGraphSelectRootDialog(m_Parent); // int dialogReturnValue = dialog->exec(); // delete dialog; // if ( dialogReturnValue != QDialog::Rejected ) // user doesn't clicked cancel or pressed Esc or something similar //{ m_TubeGraph->SetRootTube(m_LastPickedTube); //} m_TubeGraphProperty->DeactivateAllTubes(); RenderingManager::GetInstance()->RequestUpdateAll(); } else { switch (m_ActivationMode) { case None: { m_TubeGraphProperty->DeactivateAllTubes(); } break; case Single: { m_TubeGraphProperty->DeactivateAllTubes(); m_TubeGraphProperty->SetTubeActive(m_LastPickedTube, true); } break; case Multiple: { // special deactivation for multiple modus // if activated--> deactivate; if not activated--> activate if (m_TubeGraphProperty->IsTubeActive(m_LastPickedTube)) m_TubeGraphProperty->SetTubeActive(m_LastPickedTube, false); else m_TubeGraphProperty->SetTubeActive(m_LastPickedTube, true); } break; case ToRoot: { m_TubeGraphProperty->DeactivateAllTubes(); std::vector activeTubes = this->GetTubesToRoot(); m_TubeGraphProperty->SetTubesActive(activeTubes); } break; case ToPeriphery: { m_TubeGraphProperty->DeactivateAllTubes(); std::vector activeTubes = this->GetPathToPeriphery(); m_TubeGraphProperty->SetTubesActive(activeTubes); } break; case Points: { m_TubeGraphProperty->DeactivateAllTubes(); std::vector activeTubes = this->GetTubesBetweenPoints(); m_TubeGraphProperty->SetTubesActive(activeTubes); } break; default: MITK_WARN << "Unknown tube graph interaction mode!"; break; } } } std::vector mitk::TubeGraphDataInteractor::GetTubesToRoot() { TubeGraph::TubeDescriptorType root = m_TubeGraph->GetRootTube(); if (root == TubeGraph::ErrorId) { root = m_TubeGraph->GetThickestTube(); m_TubeGraph->SetRootTube(root); } return this->GetPathBetweenTubes(m_LastPickedTube, root); } std::vector mitk::TubeGraphDataInteractor::GetTubesBetweenPoints() { return this->GetPathBetweenTubes(m_LastPickedTube, m_SecondLastPickedTube); } std::vector mitk::TubeGraphDataInteractor::GetPathBetweenTubes( const mitk::TubeGraph::TubeDescriptorType &start, const mitk::TubeGraph::TubeDescriptorType &end) { std::vector solutionPath; if ((start != TubeGraph::ErrorId) && (end != TubeGraph::ErrorId)) { if (start != end) solutionPath = m_TubeGraph->SearchAllPathBetweenVertices(start, end); else solutionPath.push_back(start); } return solutionPath; } std::vector mitk::TubeGraphDataInteractor::GetPathToPeriphery() { std::vector solutionPath; if (m_LastPickedTube != TubeGraph::ErrorId) solutionPath = m_TubeGraph->SearchPathToPeriphery(m_LastPickedTube); return solutionPath; } void mitk::TubeGraphDataInteractor::ResetPickedTubes() { m_LastPickedTube = TubeGraph::ErrorId; m_SecondLastPickedTube = TubeGraph::ErrorId; } mitk::Point3D mitk::TubeGraphDataInteractor::GetLastPickedPosition() { if (m_LastPickedElement) return m_LastPickedElement->GetCoordinates(); else return mitk::Point3D(); }