diff --git a/Modules/BasicImageProcessing/include/mitkTransformationOperation.h b/Modules/BasicImageProcessing/include/mitkTransformationOperation.h index da75c42d7c..c36e4478c8 100644 --- a/Modules/BasicImageProcessing/include/mitkTransformationOperation.h +++ b/Modules/BasicImageProcessing/include/mitkTransformationOperation.h @@ -1,79 +1,79 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ -#ifndef mitkArithmeticOperation_h -#define mitkArithmeticOperation_h +#ifndef mitkTransformationOperation_h +#define mitkTransformationOperation_h #include #include #include namespace mitk { enum BorderCondition { Constant, Periodic, ZeroFluxNeumann }; enum WaveletType { Held, Vow, Simoncelli, Shannon }; enum GridInterpolationPositionType { SameSize, OriginAligned, CenterAligned }; /** \brief Executes a transformation operations on one or two images * * All parameters of the arithmetic operations must be specified during construction. * The actual operation is executed when calling GetResult(). */ class MITKBASICIMAGEPROCESSING_EXPORT TransformationOperation { public: static std::vector MultiResolution(Image::Pointer & image, unsigned int numberOfLevels, bool outputAsDouble = false); static Image::Pointer LaplacianOfGaussian(Image::Pointer & image, double sigma, bool outputAsDouble = false); static std::vector WaveletForward(Image::Pointer image, unsigned int numberOfLevels, unsigned int numberOfBands, BorderCondition condition, WaveletType waveletType); static Image::Pointer ResampleImage(Image::Pointer &image, mitk::Vector3D spacing, mitk::ImageMappingInterpolator::Type interpolator, GridInterpolationPositionType position, bool returnAsDouble, bool roundOutput); static Image::Pointer ResampleMask(Image::Pointer &image, mitk::Vector3D spacing, mitk::ImageMappingInterpolator::Type interpolator, GridInterpolationPositionType position); static Image::Pointer CastToUnsignedChar(Image::Pointer image, int &error); static Image::Pointer CastToSignedChar(Image::Pointer image, int &error); static Image::Pointer CastToChar(Image::Pointer image, int &error); static Image::Pointer CastToUnsignedShort(Image::Pointer image, int &error); static Image::Pointer CastToShort(Image::Pointer image, int &error); static Image::Pointer CastToUnsignedInt(Image::Pointer image, int &error); static Image::Pointer CastToInt(Image::Pointer image, int &error); static Image::Pointer CastToUnsignedLong(Image::Pointer image, int &error); static Image::Pointer CastToLong(Image::Pointer image, int &error); static Image::Pointer CastToFloat(Image::Pointer image, int &error); static Image::Pointer CastToDouble(Image::Pointer image, int &error); }; } -#endif // mitkArithmeticOperation_h \ No newline at end of file +#endif // mitkTransformationOperation_h \ No newline at end of file diff --git a/Modules/BasicImageProcessing/src/mitkMaskCleaningOperation.cpp b/Modules/BasicImageProcessing/src/mitkMaskCleaningOperation.cpp index ce6d2bb4b2..75d5320b8c 100644 --- a/Modules/BasicImageProcessing/src/mitkMaskCleaningOperation.cpp +++ b/Modules/BasicImageProcessing/src/mitkMaskCleaningOperation.cpp @@ -1,141 +1,142 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkMaskCleaningOperation.h" #include #include #include #include #include #include #include "itkCastImageFilter.h" namespace mitk { namespace Functor { template< class TInput, class TOutput> class RangeBasedMasking { public: RangeBasedMasking() {}; ~RangeBasedMasking() {}; bool operator!=(const RangeBasedMasking &) const { return false; } bool operator==(const RangeBasedMasking & other) const { return !(*this != other); } inline TOutput operator()(const TInput & A, const TOutput & B) const { unsigned short erg = B; if (lowerLimitOn && (A < lowerLimit)) { erg = 0; } if (upperLimitOn && (upperLimit < A)) { erg = 0; } return erg; } bool lowerLimitOn = false; bool upperLimitOn = false; double lowerLimit = 0; double upperLimit = 1; }; } } template static void ExecuteRangeBasedMasking(itk::Image* image, mitk::Image::Pointer & mask, bool lowerLimitOn, double lowerLimit, bool upperLimitOn, double upperLimit, mitk::Image::Pointer & resultImage) { typedef itk::Image ImageType; typedef itk::Image MaskImageType; typedef itk::BinaryFunctorImageFilter< ImageType, MaskImageType, MaskImageType, mitk::Functor::RangeBasedMasking > DefaultFilterType; typename MaskImageType::Pointer itkMask = MaskImageType::New(); mitk::CastToItkImage(mask.GetPointer(), itkMask); typename DefaultFilterType::Pointer filter = DefaultFilterType::New(); filter->SetInput1(image); filter->SetInput2(itkMask); filter->GetFunctor().lowerLimitOn = lowerLimitOn; filter->GetFunctor().upperLimitOn = upperLimitOn; filter->GetFunctor().lowerLimit = lowerLimit; filter->GetFunctor().upperLimit = upperLimit; filter->Update(); CastToMitkImage(filter->GetOutput(), resultImage); } mitk::Image::Pointer mitk::MaskCleaningOperation::RangeBasedMasking(Image::Pointer & image, Image::Pointer & mask, bool lowerLimitOn, double lowerLimit, bool upperLimitOn, double upperLimit) { Image::Pointer resultImage; AccessByItk_n(image, ExecuteRangeBasedMasking, (mask, lowerLimitOn, lowerLimit, upperLimitOn, upperLimit, resultImage)); return resultImage; } template static void CalculateMeanAndStd(itk::Image* image, mitk::Image::Pointer & mask, double &mean, double &std) { typedef itk::Image ImageType; typedef itk::Image MaskImageType; typename MaskImageType::Pointer itkMask = MaskImageType::New(); mitk::CastToItkImage(mask.GetPointer(), itkMask); itk::ImageRegionConstIterator image_iter(image, image->GetLargestPossibleRegion()); itk::ImageRegionConstIterator mask_iter(itkMask, itkMask->GetLargestPossibleRegion()); unsigned int number_of_voxels = 0; while (!image_iter.IsAtEnd()) { if (mask_iter.Get() > 0) { mean += image_iter.Get(); std += image_iter.Get() * image_iter.Get(); ++number_of_voxels; } ++image_iter; ++mask_iter; } mean /= number_of_voxels; std /= number_of_voxels; - std -= mean; + std -= mean*mean; std = sqrt(std); + MITK_INFO << "Mean: " << mean << " Sigma: " << std; } mitk::Image::Pointer mitk::MaskCleaningOperation::MaskOutlierFiltering(Image::Pointer & image, Image::Pointer & mask) { Image::Pointer resultImage; double mean = 0; double std = 0; AccessByItk_n(image, CalculateMeanAndStd, (mask, mean, std)); return MaskCleaningOperation::RangeBasedMasking(image, mask, true, mean - 3 * std, true, mean + 3 * std); } \ No newline at end of file diff --git a/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp b/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp index 00aeebf9c5..0f2a6dd54f 100644 --- a/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp +++ b/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp @@ -1,756 +1,756 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTransformationOperation.h" #include #include #include #include #include #include #include #include // Wavelet #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); } }; } } // Error Codes: 0: Successfull // 1: General Error (Not Specific) // 2: Error when duplicating images // 3: Error when casting image template static void ExecuteConvertToDataType(itk::Image* image, mitk::Image::Pointer &resultImage, int & isSuccessful) { typedef itk::Image InputImageType; typedef itk::Image OutputImageType; isSuccessful = 1; if (std::is_same::value) // If input has already the "right" data type, only clone it. { try { typedef itk::ImageDuplicator DuplicateFilterType; typename DuplicateFilterType::Pointer duplicateFilter = DuplicateFilterType::New(); duplicateFilter->SetInputImage(image); duplicateFilter->Update(); CastToMitkImage(duplicateFilter->GetOutput(), resultImage); } catch (...) { resultImage = nullptr; isSuccessful = 2; return; } } else // Convert the image data type if it has the "wrong" datatype. { try { typedef itk::CastImageFilter CastFilterType; typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput(image); castFilter->Update(); CastToMitkImage(castFilter->GetOutput(), resultImage); } catch (...) { resultImage = nullptr; isSuccessful = 3; return; } } isSuccessful = 0; return; } // Error Codes: 0: Successfull // 1: General Error (Not Specific) // 2: Error when duplicating images // 3: Error when casting image // 4: General Error (really early) // 5: Unsupported image type // // Commented types are not template static void ExecuteConvertToGivenType(itk::Image* image, std::type_index &type, mitk::Image::Pointer &resultImage, int & isSuccessful) { isSuccessful = 4; if (type == typeid(unsigned char)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } else if (type == typeid(signed char)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } else if (type == typeid(char)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } //else if (type == typeid(wchar_t)) //{ // ExecuteConvertToDataType(image, resultImage, isSuccessful); //} //else if (type == typeid(char16_t)) //{ // ExecuteConvertToDataType(image, resultImage, isSuccessful); //} //else if (type == typeid(char32_t)) //{ // ExecuteConvertToDataType(image, resultImage, isSuccessful); //} //else if (type == typeid(bool)) //{ // ExecuteConvertToDataType(image, resultImage, isSuccessful); //} else if (type == typeid(unsigned short int)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } else if (type == typeid(short int)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } else if (type == typeid(unsigned int)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } else if (type == typeid(int)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } else if (type == typeid(unsigned long int)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } else if (type == typeid(long int)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } //else if (type == typeid(unsigned long long int)) //{ // ExecuteConvertToDataType(image, resultImage, isSuccessful); //} //else if (type == typeid(long long int)) //{ // ExecuteConvertToDataType(image, resultImage, isSuccessful); //} else if (type == typeid(float)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } else if (type == typeid(double)) { ExecuteConvertToDataType(image, resultImage, isSuccessful); } //else if (type == typeid(long double)) //{ // ExecuteConvertToDataType(image, resultImage, isSuccessful); //} else { resultImage = NULL; isSuccessful = 5; } } 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::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; auto oldDirection = image->GetDirection(); typename ImageType::DirectionType oldDir = image->GetDirection(); typename ImageType::DirectionType newDir; for (int i = 0; i < VImageDimension; ++i) { for (int j = 0; j < VImageDimension; ++j) { newDir[i][j] = ((i == j) ? itk::SpacePrecisionType(1.0) : itk::SpacePrecisionType(0)); } } image->SetDirection(newDir); // 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->GetSpacing()[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) { 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->SetOrigin(inputOrigin); itkOutputImage->SetSpacing(expectedSpacing); itkOutputImage->SetDirection(oldDir); 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]; + double oldLength = newGeometry->GetSpacing()[i] * (newGeometry->GetBounds()[i*2+1] - 1); + double newLength = spacing[i] * (bounds[i*2+1] - 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; } mitk::Image::Pointer mitk::TransformationOperation::CastToUnsignedChar(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(unsigned char); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); return result; } mitk::Image::Pointer mitk::TransformationOperation::CastToSignedChar(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(signed char); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); return result; } mitk::Image::Pointer mitk::TransformationOperation::CastToChar(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(char); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); return result; } mitk::Image::Pointer mitk::TransformationOperation::CastToUnsignedShort(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(unsigned short int); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); return result; } mitk::Image::Pointer mitk::TransformationOperation::CastToShort(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(short int); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); return result; } mitk::Image::Pointer mitk::TransformationOperation::CastToUnsignedInt(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(unsigned int); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); return result; } mitk::Image::Pointer mitk::TransformationOperation::CastToInt(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(int); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); return result; } mitk::Image::Pointer mitk::TransformationOperation::CastToUnsignedLong(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(unsigned long int); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); return result; } mitk::Image::Pointer mitk::TransformationOperation::CastToLong(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(long int); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); return result; } mitk::Image::Pointer mitk::TransformationOperation::CastToFloat(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(float); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); return result; } mitk::Image::Pointer mitk::TransformationOperation::CastToDouble(Image::Pointer image, int &error) { mitk::Image::Pointer result; std::type_index type = typeid(double); AccessByItk_n(image, ExecuteConvertToGivenType, (type, result, error)); 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/CLCore/include/mitkAbstractGlobalImageFeature.h b/Modules/Classification/CLCore/include/mitkAbstractGlobalImageFeature.h index 362dcd73bf..fba5996563 100644 --- a/Modules/Classification/CLCore/include/mitkAbstractGlobalImageFeature.h +++ b/Modules/Classification/CLCore/include/mitkAbstractGlobalImageFeature.h @@ -1,295 +1,295 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkAbstractGlobalImageFeature_h #define mitkAbstractGlobalImageFeature_h #include #include #include #include #include // STD Includes // Eigen #include // MITK includes #include namespace mitk { /** * * * ## Histogram Configuration ## * Most Feature Generation Classes that use histograms use the same parameters and * initialization logic. In general, all information can be passed either by the corresponding * Setter (which does not differenciate between global setting and feature specific setting) and * a parameter object which can be obtained from the command line arguments, for example. * * If the image values are used for the initializiation of the histogram, it can be defined * whether the whole image is used or only the masked areas to find minima and maxima. This is * done by the option SetIgnoreMask or the corrsponding options * -NAME::ignore-mask-for-histogram and -ignore-mask-for-histogram. If these are * true, the whole image is used for the calculation. * * Depending on the passed arguments, different initialization methods are used. The initialization * is in the following order: * - If Minimum Intensity, Maximum Intensity, and Binsize: The histogram is * initialized between the minimum and maximum intensity. the number of bins is determined by the * binsize. If the distance between minimum and maximum is not a multiple of the binsize, the maximum * is increase so that it is. * - Minimum Intensity, Bins, and Binsize: The histogram is initialized with the * given binsize, and the intensity range from the minimum to \f$maximum = minimum + binsize*bins\f$. * - Minimum Intensity, Maximum Intensity, and Bins: The histogram is initialized * between the given minimum and maximum intensity. The binsize is calculated so that the number * of bins is equal to the given number of bins. * - Binsize, and Minimum Intensity: The maximum is set to the maximum that * occur in the given image. Depending if the mask is considered or not, either only masked voxels or * the whole image is used for the calculation. The initialization is then equal as if the minimum * and maximum would have been given right from the beginning. * - Binsize, and Maximum Intensity: The minimum intensity is set to the minimum that * occur in the given image. Depending if the mask is considered or not, either only masked voxels or * the whole image is used for the calculation. The initialization is then equal as if the minimum * and maximum would have been given right from the beginning. * - Binsize: The maximum and the minimum intensity is set to the minimum and maximum that * occur in the given image. Depending if the mask is considered or not, either only masked voxels or * the whole image is used for the calculation. The initialization is then equal as if the minimum * and maximum would have been given right from the beginning. * - Bins, and Minimum Intensity: The maximum is calculated from the image. Depending * if the mask is considered or not, either only masked voxels or the whole image is used for the calculation. The histogram is * then initialized as if these values would have been given as minimum and maximum intensity. * - Bins, and Maximum Intensity: The minimum is calculated from the image. Depending * if the mask is considered or not, either only masked voxels or the whole image is used for the calculation. The histogram is * then initialized as if these values would have been given as minimum and maximum intensity. * - Bins: The minimum and the maximum is calculated from the image. Depending * if the mask is considered or not, either only masked voxels or * the whole image is used for the calculation. The histogram is * then initialized as if these values would have been given as minimum and maximum intensity. * - No Parameter given:The minimum and maximum intensity from the whole image or masked image is calculated and * the histogram then initialized to this with a standard number of bins (Is set by each filter on its own.) * * ### Remark about command line parameter#### * There are generally two options to set a parameter via the command line. A global one that works for * all filters that use histograms and a local one that set this parameter specific for this filter. The * local parameters start with the filter name (Indiciated by NAME) followed by two colons, for example * vol::min to set the minimum intensity for the volume filter. The global parameter is overwritten * by the local parameter, if it is specified. Otherwise, it is still valid. If this prevents the specification * of an histogram initialization method (for example, because the binsize is globally specified but the histogram * should be initialized using a fixed numbe of bins), the parameter NAME::ignore-global-histogram can be passed. * Then, all global histogram parameters are ignored and only local ones are used. * * The maximum intensity can be set by different command line parameters: global for all filters that use histograms * by -minimum-intensity and -minimum. Alternative it can be set only for this filter by * -NAME::minimum and -NAME::min. * * The minimum intensity can be set by different command line parameters: global for all filters that use histograms * by -maximum-intensity and -maximum. Alternative it can be set only for this filter by * \NAME::maximum and NAME::max. * * The binsize can be set by different command line parameters: global for all filters that use histograms * by -binsize. Alternative it can be set only for this filter by * \NAME::binsize. * * The number of bins can be set by different command line parameters: global for all filters that use histograms * by -bins. Alternative it can be set only for this filter by * \NAME::bins. * ### Note to the developers ### * All features are supposed to work the same way if a histogram is used somewhere in * the code. For this, each derived class that makes use of a histogram should use * the Quantifier object. In order to use this object correctly, the AddArguments-Function should * contain the line AddQuantifierArguments(parser);, the CalculateFeaturesUsingParameters function * should contain the line InitializeQuantifierFromParameters(feature, mask); and the CalculateFeatures function * sould contain the line InitializeQuantifier(image, mask);. These function * calls ensure that the necessary options are given to the configuration file, and that the initialization * of the quantifier is done correctly. This ensures an consistend behavior over all FeatureGeneration Classes. * */ class MITKCLCORE_EXPORT AbstractGlobalImageFeature : public BaseData { public: mitkClassMacro(AbstractGlobalImageFeature, BaseData) typedef std::vector< std::pair > FeatureListType; typedef std::vector< std::string> FeatureNameListType; typedef std::map ParameterTypes; /** * \brief Calculates the feature of this abstact interface. Does not necessarily considers the parameter settings. */ virtual FeatureListType CalculateFeatures(const Image::Pointer & feature, const Image::Pointer &mask) = 0; /** * \brief Calculates the given feature Slice-wise. Might not be availble for an individual filter! */ FeatureListType CalculateFeaturesSlicewise(const Image::Pointer & feature, const Image::Pointer &mask, int sliceID); /** * \brief Calculates the feature of this abstact interface. Does not necessarily considers the parameter settings. */ virtual void CalculateFeaturesSliceWiseUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, int sliceID, FeatureListType &featureList); /** * \brief Calculates the feature of this abstact interface. Does not necessarily considers the parameter settings. */ virtual void CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList) = 0; /** * \brief Returns a list of the names of all features that are calculated from this class */ virtual FeatureNameListType GetFeatureNames() = 0; /** * \brief Adds an additional Separator to the name of the feature, which encodes the used parameters */ virtual std::string GetCurrentFeatureEncoding(); /** * \brief Returns a string that encodes the feature class name. * * This Feature returns a string that should be put before every other feature * value. It has the format [::]::[::]. * The options and are only added, if the corresponding * options are set. */ std::string FeatureDescriptionPrefix(); itkSetMacro(Prefix, std::string); itkSetMacro(ShortName, std::string); itkSetMacro(LongName, std::string); itkSetMacro(FeatureClassName, std::string); itkSetMacro(Direction, int); void SetParameter(ParameterTypes param) { m_Parameter=param; }; itkGetConstMacro(Prefix, std::string); itkGetConstMacro(ShortName, std::string); itkGetConstMacro(LongName, std::string); itkGetConstMacro(FeatureClassName, std::string); itkGetConstMacro(Parameter, ParameterTypes); itkSetMacro(UseQuantifier, bool); itkGetConstMacro(UseQuantifier, bool); itkSetMacro(Quantifier, IntensityQuantifier::Pointer); itkGetMacro(Quantifier, IntensityQuantifier::Pointer); itkGetConstMacro(Direction, int); itkSetMacro(MinimumIntensity, double); itkSetMacro(UseMinimumIntensity, bool); itkSetMacro(MaximumIntensity, double); itkSetMacro(UseMaximumIntensity, bool); itkGetConstMacro(MinimumIntensity, double); itkGetConstMacro(UseMinimumIntensity, bool); itkGetConstMacro(MaximumIntensity, double); itkGetConstMacro(UseMaximumIntensity, bool); itkSetMacro(Binsize, double); itkSetMacro(UseBinsize, bool); itkGetConstMacro(Binsize, double); itkGetConstMacro(UseBinsize, bool); - itkSetMacro(MorphMask, mitk::Image::Pointer); - itkGetConstMacro(MorphMask, mitk::Image::Pointer); + itkSetObjectMacro(MorphMask, mitk::Image); + itkGetObjectMacro(MorphMask, mitk::Image); itkSetMacro(Bins, int); itkSetMacro(UseBins, bool); itkGetConstMacro(UseBins, bool); itkGetConstMacro(Bins, int); itkSetMacro(IgnoreMask, bool); itkGetConstMacro(IgnoreMask, bool); itkSetMacro(EncodeParameters, bool); itkGetConstMacro(EncodeParameters, bool); std::string GetOptionPrefix() const { if (m_Prefix.length() > 0) return m_Prefix + "::" + m_ShortName; return m_ShortName; } virtual void AddArguments(mitkCommandLineParser &parser) = 0; std::vector SplitDouble(std::string str, char delimiter); void AddQuantifierArguments(mitkCommandLineParser &parser); void InitializeQuantifierFromParameters(const Image::Pointer & feature, const Image::Pointer &mask,unsigned int defaultBins = 256); void InitializeQuantifier(const Image::Pointer & feature, const Image::Pointer &mask, unsigned int defaultBins = 256); std::string QuantifierParameterString(); public: //#ifndef DOXYGEN_SKIP void SetRequestedRegionToLargestPossibleRegion() override {}; bool RequestedRegionIsOutsideOfTheBufferedRegion() override { return true; }; bool VerifyRequestedRegion() override { return false; }; void SetRequestedRegion (const itk::DataObject * /*data*/) override {}; // Override bool IsEmpty() const override { if(IsInitialized() == false) return true; const TimeGeometry* timeGeometry = const_cast(this)->GetUpdatedTimeGeometry(); if(timeGeometry == nullptr) return true; return false; } private: std::string m_Prefix; // Prefix before all input parameters std::string m_ShortName; // Name of all variables std::string m_LongName; // Long version of the name (For turning on) std::string m_FeatureClassName; ParameterTypes m_Parameter; // Parameter setting bool m_UseQuantifier = false; IntensityQuantifier::Pointer m_Quantifier; double m_MinimumIntensity = 0; bool m_UseMinimumIntensity = false; double m_MaximumIntensity = 100; bool m_UseMaximumIntensity = false; bool m_EncodeParameters = false; double m_Binsize = 1; bool m_UseBinsize = false; int m_Bins = 256; bool m_UseBins = true; int m_Direction = 0; bool m_IgnoreMask = false; bool m_CalculateWithParameter = false; mitk::Image::Pointer m_MorphMask = nullptr; //#endif // Skip Doxygen }; } #endif //mitkAbstractGlobalImageFeature_h diff --git a/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp b/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp index 5920b99ce2..64fc4cc9ce 100644 --- a/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp +++ b/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp @@ -1,406 +1,425 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include static void ExtractSlicesFromImages(mitk::Image::Pointer image, mitk::Image::Pointer mask, int direction, std::vector &imageVector, std::vector &maskVector) { typedef itk::Image< double, 2 > FloatImage2DType; typedef itk::Image< unsigned short, 2 > MaskImage2DType; typedef itk::Image< double, 3 > FloatImageType; typedef itk::Image< unsigned short, 3 > MaskImageType; FloatImageType::Pointer itkFloat = FloatImageType::New(); MaskImageType::Pointer itkMask = MaskImageType::New(); mitk::CastToItkImage(mask, itkMask); mitk::CastToItkImage(image, itkFloat); int idxA, idxB, idxC; switch (direction) { case 0: idxA = 1; idxB = 2; idxC = 0; break; case 1: idxA = 0; idxB = 2; idxC = 1; break; case 2: idxA = 0; idxB = 1; idxC = 2; break; default: idxA = 1; idxB = 2; idxC = 0; break; } auto imageSize = image->GetLargestPossibleRegion().GetSize(); FloatImageType::IndexType index3D; FloatImage2DType::IndexType index2D; FloatImage2DType::SpacingType spacing2D; spacing2D[0] = itkFloat->GetSpacing()[idxA]; spacing2D[1] = itkFloat->GetSpacing()[idxB]; for (unsigned int i = 0; i < imageSize[idxC]; ++i) { FloatImage2DType::RegionType region; FloatImage2DType::IndexType start; FloatImage2DType::SizeType size; start[0] = 0; start[1] = 0; size[0] = imageSize[idxA]; size[1] = imageSize[idxB]; region.SetIndex(start); region.SetSize(size); FloatImage2DType::Pointer image2D = FloatImage2DType::New(); image2D->SetRegions(region); image2D->Allocate(); MaskImage2DType::Pointer mask2D = MaskImage2DType::New(); mask2D->SetRegions(region); mask2D->Allocate(); unsigned long voxelsInMask = 0; for (unsigned int a = 0; a < imageSize[idxA]; ++a) { for (unsigned int b = 0; b < imageSize[idxB]; ++b) { index3D[idxA] = a; index3D[idxB] = b; index3D[idxC] = i; index2D[0] = a; index2D[1] = b; image2D->SetPixel(index2D, itkFloat->GetPixel(index3D)); mask2D->SetPixel(index2D, itkMask->GetPixel(index3D)); voxelsInMask += (itkMask->GetPixel(index3D) > 0) ? 1 : 0; } } image2D->SetSpacing(spacing2D); mask2D->SetSpacing(spacing2D); mitk::Image::Pointer tmpFloatImage = mitk::Image::New(); tmpFloatImage->InitializeByItk(image2D.GetPointer()); mitk::GrabItkImageMemory(image2D, tmpFloatImage); mitk::Image::Pointer tmpMaskImage = mitk::Image::New(); tmpMaskImage->InitializeByItk(mask2D.GetPointer()); mitk::GrabItkImageMemory(mask2D, tmpMaskImage); if (voxelsInMask > 0) { imageVector.push_back(tmpFloatImage); maskVector.push_back(tmpMaskImage); } } } std::vector mitk::AbstractGlobalImageFeature::SplitDouble(std::string str, char delimiter) { std::vector internal; std::stringstream ss(str); // Turn the string into a stream. std::string tok; double val; while (std::getline(ss, tok, delimiter)) { std::stringstream s2(tok); s2 >> val; internal.push_back(val); } return internal; } void mitk::AbstractGlobalImageFeature::AddQuantifierArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(name + "::minimum", name + "::min", mitkCommandLineParser::Float, "Minium Intensity for Quantification", "Defines the minimum Intensity used for Quantification", us::Any()); parser.addArgument(name + "::maximum", name + "::max", mitkCommandLineParser::Float, "Maximum Intensity for Quantification", "Defines the maximum Intensity used for Quantification", us::Any()); parser.addArgument(name + "::bins", name + "::bins", mitkCommandLineParser::Int, "Number of Bins", "Define the number of bins that is used ", us::Any()); parser.addArgument(name + "::binsize", name + "::binsize", mitkCommandLineParser::Float, "Binsize", "Define the size of the used bins", us::Any()); parser.addArgument(name + "::ignore-global-histogram", name + "::ignore-global-histogram", mitkCommandLineParser::Bool, "Ignore the global histogram Parameters", "Ignores the global histogram parameters", us::Any()); parser.addArgument(name + "::ignore-mask-for-histogram", name + "::ignore-mask", mitkCommandLineParser::Bool, "Ignore the global histogram Parameters", "Ignores the global histogram parameters", us::Any()); } void mitk::AbstractGlobalImageFeature::InitializeQuantifierFromParameters(const Image::Pointer & feature, const Image::Pointer &mask, unsigned int defaultBins) { unsigned int bins = 0; double binsize = 0; double minimum = 0; double maximum = 0; auto parsedArgs = GetParameter(); std::string name = GetOptionPrefix(); bool useGlobal = true; if (parsedArgs.count(name + "::ignore-global-histogram")) { useGlobal = false; SetUseMinimumIntensity(false); SetUseMaximumIntensity(false); SetUseBinsize(false); SetUseBins(false); } if (useGlobal) { if (parsedArgs.count("ignore-mask-for-histogram")) { bool tmp = us::any_cast(parsedArgs["ignore-mask-for-histogram"]); SetIgnoreMask(tmp); } if (parsedArgs.count("minimum-intensity")) { minimum = us::any_cast(parsedArgs["minimum-intensity"]); SetMinimumIntensity(minimum); SetUseMinimumIntensity(true); } if (parsedArgs.count("maximum-intensity")) { maximum = us::any_cast(parsedArgs["maximum-intensity"]); SetMaximumIntensity(maximum); SetUseMaximumIntensity(true); } if (parsedArgs.count("bins")) { bins = us::any_cast(parsedArgs["bins"]); SetBins(bins); SetUseBins(true); } if (parsedArgs.count("binsize")) { binsize = us::any_cast(parsedArgs["binsize"]); SetBinsize(binsize); SetUseBinsize(true); } } if (parsedArgs.count(name+"::ignore-mask-for-histogram")) { bool tmp = us::any_cast(parsedArgs[name+"::ignore-mask-for-histogram"]); SetIgnoreMask(tmp); } if (parsedArgs.count(name + "::minimum")) { minimum = us::any_cast(parsedArgs[name + "::minimum"]); SetMinimumIntensity(minimum); SetUseMinimumIntensity(true); } if (parsedArgs.count(name + "::maximum")) { maximum = us::any_cast(parsedArgs[name + "::maximum"]); SetMaximumIntensity(maximum); SetUseMaximumIntensity(true); } if (parsedArgs.count(name + "::bins")) { bins = us::any_cast(parsedArgs[name + "::bins"]); SetBins(bins); } if (parsedArgs.count(name + "::binsize")) { binsize = us::any_cast(parsedArgs[name + "::binsize"]); SetBinsize(binsize); SetUseBinsize(true); } InitializeQuantifier(feature, mask, defaultBins); } void mitk::AbstractGlobalImageFeature::InitializeQuantifier(const Image::Pointer & feature, const Image::Pointer &mask, unsigned int defaultBins) { m_Quantifier = IntensityQuantifier::New(); if (GetUseMinimumIntensity() && GetUseMaximumIntensity() && GetUseBinsize()) m_Quantifier->InitializeByBinsizeAndMaximum(GetMinimumIntensity(), GetMaximumIntensity(), GetBinsize()); else if (GetUseMinimumIntensity() && GetUseBins() && GetUseBinsize()) m_Quantifier->InitializeByBinsizeAndBins(GetMinimumIntensity(), GetBins(), GetBinsize()); else if (GetUseMinimumIntensity() && GetUseMaximumIntensity() && GetUseBins()) m_Quantifier->InitializeByMinimumMaximum(GetMinimumIntensity(), GetMaximumIntensity(), GetBins()); // Intialize from Image and Binsize else if (GetUseBinsize() && GetIgnoreMask() && GetUseMinimumIntensity()) m_Quantifier->InitializeByImageAndBinsizeAndMinimum(feature, GetMinimumIntensity(), GetBinsize()); else if (GetUseBinsize() && GetIgnoreMask() && GetUseMaximumIntensity()) m_Quantifier->InitializeByImageAndBinsizeAndMaximum(feature, GetMaximumIntensity(), GetBinsize()); else if (GetUseBinsize() && GetIgnoreMask()) m_Quantifier->InitializeByImageAndBinsize(feature, GetBinsize()); // Initialize form Image, Mask and Binsize else if (GetUseBinsize() && GetUseMinimumIntensity()) m_Quantifier->InitializeByImageRegionAndBinsizeAndMinimum(feature, mask, GetMinimumIntensity(), GetBinsize()); else if (GetUseBinsize() && GetUseMaximumIntensity()) m_Quantifier->InitializeByImageRegionAndBinsizeAndMaximum(feature, mask, GetMaximumIntensity(), GetBinsize()); else if (GetUseBinsize()) m_Quantifier->InitializeByImageRegionAndBinsize(feature, mask, GetBinsize()); // Intialize from Image and Bins else if (GetUseBins() && GetIgnoreMask() && GetUseMinimumIntensity()) m_Quantifier->InitializeByImageAndMinimum(feature, GetMinimumIntensity(), GetBins()); else if (GetUseBins() && GetIgnoreMask() && GetUseMaximumIntensity()) m_Quantifier->InitializeByImageAndMaximum(feature, GetMaximumIntensity(), GetBins()); - else if (GetUseBins()) + else if (GetUseBins() && GetIgnoreMask()) m_Quantifier->InitializeByImage(feature, GetBins()); // Intialize from Image, Mask and Bins else if (GetUseBins() && GetUseMinimumIntensity()) m_Quantifier->InitializeByImageRegionAndMinimum(feature, mask, GetMinimumIntensity(), GetBins()); else if (GetUseBins() && GetUseMaximumIntensity()) m_Quantifier->InitializeByImageRegionAndMaximum(feature, mask, GetMaximumIntensity(), GetBins()); - else if (GetUseBins()) + else if (GetUseBins() ) m_Quantifier->InitializeByImageRegion(feature, mask, GetBins()); // Default else if (GetIgnoreMask()) - m_Quantifier->InitializeByImage(feature, GetBins()); + m_Quantifier->InitializeByImage(feature, defaultBins); else m_Quantifier->InitializeByImageRegion(feature, mask, defaultBins); } std::string mitk::AbstractGlobalImageFeature::GetCurrentFeatureEncoding() { return ""; } std::string mitk::AbstractGlobalImageFeature::FeatureDescriptionPrefix() { std::string output; output = m_FeatureClassName + "::"; if (m_EncodeParameters) { output += GetCurrentFeatureEncoding() + "::"; } return output; } std::string mitk::AbstractGlobalImageFeature::QuantifierParameterString() { std::stringstream ss; if (GetUseMinimumIntensity() && GetUseMaximumIntensity() && GetUseBinsize()) ss << "Min-" << GetMinimumIntensity() << "_Max-" << GetMaximumIntensity() << "_BS-" << GetBinsize(); else if (GetUseMinimumIntensity() && GetUseBins() && GetUseBinsize()) ss << "Min-" << GetMinimumIntensity() << "_Bins-" << GetBins() << "_BS-" << GetBinsize(); else if (GetUseMinimumIntensity() && GetUseMaximumIntensity() && GetUseBins()) ss << "Min-" << GetMinimumIntensity() << "_Max-" << GetMaximumIntensity() << "_Bins-" << GetBins(); // Intialize from Image and Binsize else if (GetUseBinsize() && GetIgnoreMask() && GetUseMinimumIntensity()) ss << "Min-" << GetMinimumIntensity() << "_BS-" << GetBinsize() << "_FullImage"; else if (GetUseBinsize() && GetIgnoreMask() && GetUseMaximumIntensity()) ss << "Max-" << GetMaximumIntensity() << "_BS-" << GetBinsize() << "_FullImage"; else if (GetUseBinsize() && GetIgnoreMask()) ss << "BS-" << GetBinsize() << "_FullImage"; // Initialize form Image, Mask and Binsize else if (GetUseBinsize() && GetUseMinimumIntensity()) ss << "Min-" << GetMinimumIntensity() << "_BS-" << GetBinsize(); else if (GetUseBinsize() && GetUseMaximumIntensity()) ss << "Max-" << GetMaximumIntensity() << "_BS-" << GetBinsize(); else if (GetUseBinsize()) ss << "BS-" << GetBinsize(); // Intialize from Image and Bins else if (GetUseBins() && GetIgnoreMask() && GetUseMinimumIntensity()) ss << "Min-" << GetMinimumIntensity() << "_Bins-" << GetBins() << "_FullImage"; else if (GetUseBins() && GetIgnoreMask() && GetUseMaximumIntensity()) ss << "Max-" << GetMaximumIntensity() << "_Bins-" << GetBins() << "_FullImage"; else if (GetUseBins()) ss << "Bins-" << GetBins() << "_FullImage"; // Intialize from Image, Mask and Bins else if (GetUseBins() && GetUseMinimumIntensity()) ss << "Min-" << GetMinimumIntensity() << "_Bins-" << GetBins(); else if (GetUseBins() && GetUseMaximumIntensity()) ss << "Max-" << GetMaximumIntensity() << "_Bins-" << GetBins(); else if (GetUseBins()) ss << "Bins-" << GetBins(); // Default else if (GetIgnoreMask()) ss << "Bins-" << GetBins() << "_FullImage"; else ss << "Bins-" << GetBins(); return ss.str(); } void mitk::AbstractGlobalImageFeature::CalculateFeaturesSliceWiseUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, int sliceID, FeatureListType &featureList) { m_CalculateWithParameter = true; auto result = CalculateFeaturesSlicewise(feature, mask, sliceID); featureList.insert(featureList.end(), result.begin(), result.end()); } mitk::AbstractGlobalImageFeature::FeatureListType mitk::AbstractGlobalImageFeature::CalculateFeaturesSlicewise(const Image::Pointer & feature, const Image::Pointer &mask, int sliceID) { std::vector imageVector; std::vector maskVector; + std::vector morphMaskVector; ExtractSlicesFromImages(feature, mask,sliceID, imageVector, maskVector); std::vector statVector; + bool morphMaskIsNotNull = (GetMorphMask() != NULL); + mitk::Image::Pointer morphMask; + if (morphMaskIsNotNull) + { + morphMask = this->GetMorphMask(); + std::vector tmpVector; + ExtractSlicesFromImages(feature, morphMask, sliceID, tmpVector, morphMaskVector); + } + for (std::size_t index = 0; index < imageVector.size(); ++index) { + if (morphMaskIsNotNull) + { + this->SetMorphMask(morphMaskVector[index].GetPointer()); + } if (m_CalculateWithParameter) { FeatureListType stat; this->CalculateFeaturesUsingParameters(imageVector[index], maskVector[index], maskVector[index], stat); statVector.push_back(stat); } else { auto stat = this->CalculateFeatures(imageVector[index], maskVector[index]); statVector.push_back(stat); } } + if (morphMaskIsNotNull) + { + this->SetMorphMask(morphMask.GetPointer()); + } + if (statVector.size() < 1) return FeatureListType(); FeatureListType statMean, statStd, result; for (std::size_t i = 0; i < statVector[0].size(); ++i) { auto cElement1 = statVector[0][i]; cElement1.first = "SliceWise Mean " + cElement1.first; cElement1.second = 0.0; auto cElement2 = statVector[0][i]; cElement2.first = "SliceWise Var. " + cElement2.first; cElement2.second = 0.0; statMean.push_back(cElement1); statStd.push_back(cElement2); } for (auto cStat : statVector) { for (std::size_t i = 0; i < cStat.size(); ++i) { statMean[i].second += cStat[i].second / (1.0*statVector.size()); } } for (auto cStat : statVector) { for (std::size_t i = 0; i < cStat.size(); ++i) { statStd[i].second += (cStat[i].second - statMean[i].second)*(cStat[i].second - statMean[i].second) / (1.0*statVector.size()); } } for (auto cStat : statVector) { std::copy(cStat.begin(), cStat.end(), std::back_inserter(result)); } std::copy(statMean.begin(), statMean.end(), std::back_inserter(result)); std::copy(statStd.begin(), statStd.end(), std::back_inserter(result)); return result; } \ No newline at end of file diff --git a/Modules/Classification/CLUtilities/include/mitkGIFGreyLevelDistanceZone.h b/Modules/Classification/CLUtilities/include/mitkGIFGreyLevelDistanceZone.h index 870e556978..63b6a7b64c 100644 --- a/Modules/Classification/CLUtilities/include/mitkGIFGreyLevelDistanceZone.h +++ b/Modules/Classification/CLUtilities/include/mitkGIFGreyLevelDistanceZone.h @@ -1,178 +1,183 @@ #ifndef mitkGIFGreyLevelDistanceZone_h #define mitkGIFGreyLevelDistanceZone_h #include #include #include #include namespace mitk { struct GreyLevelDistanceZoneFeatures { GreyLevelDistanceZoneFeatures() : SmallDistanceEmphasis(0), LargeDistanceEmphasis(0), LowGreyLevelEmphasis(0), HighGreyLevelEmphasis(0), SmallDistanceLowGreyLevelEmphasis(0), SmallDistanceHighGreyLevelEmphasis(0), LargeDistanceLowGreyLevelEmphasis(0), LargeDistanceHighGreyLevelEmphasis(0), GreyLevelNonUniformity(0), GreyLevelNonUniformityNormalized(0), ZoneDistanceNonUniformity(0), ZoneDistanceNoneUniformityNormalized(0), ZonePercentage(0), GreyLevelMean(0), GreyLevelVariance(0), ZoneDistanceMean(0), ZoneDistanceVariance(0), ZoneDistanceEntropy(0) { } public: double SmallDistanceEmphasis; double LargeDistanceEmphasis; double LowGreyLevelEmphasis; double HighGreyLevelEmphasis; double SmallDistanceLowGreyLevelEmphasis; double SmallDistanceHighGreyLevelEmphasis; double LargeDistanceLowGreyLevelEmphasis; double LargeDistanceHighGreyLevelEmphasis; double GreyLevelNonUniformity; double GreyLevelNonUniformityNormalized; double ZoneDistanceNonUniformity; double ZoneDistanceNoneUniformityNormalized; double ZonePercentage; double GreyLevelMean; double GreyLevelVariance; double ZoneDistanceMean; double ZoneDistanceVariance; double ZoneDistanceEntropy; }; class MITKCLUTILITIES_EXPORT GIFGreyLevelDistanceZone : public AbstractGlobalImageFeature { /** * \brief Calculates the Grey Level Distance Zone * * This class can be used to calculate Grey Level Distance Zone as presented in Thibault et al. 2014. * * The basic idea behind the Grey Level Distance Zone based features is to count the connected areas * with a given intensity value \f$x_i\f$ and a given distance to the border of each segmentation \f$d_i\f$. * Several features are then calculated based on a matrix, that gives the number of occurence for each * combination of \f$x_i\f$ and \f$ d_i \f$ as \f$m_{x,d}\f$. * * This feature calculator is activated by the option -grey-level-distance-zone or -gldz. * * The connected areas are based on the binned image, the binning parameters can be set via the default * parameters as described in AbstractGlobalImageFeature. It is also possible to determine the * dimensionality of the neighbourhood using direction-related commands as described in AbstractGlobalImageFeature. * No other options are possible beside these two options. * * The features are calculated based on a mask. It is assumed that the mask is * of the type of an unsigned short image. It is expected that the image contains only voxels with value 0 and 1, * of which all voxels with an value equal to one are treated as masked. * * The features depend on the distance to the border of the segmentation ROI. In some cases, the border * definition might be different from the definition of the masked area, for example, if small openings * in the mask should not influence the distance. Due to that, it is possible to submit a second mask, * named Morphological Mask to the features that is then used to calculate the distance of each voxel to * border of the segmented area. The morpological mask can be either set by the function call SetMorphMask() * or by the corresponding global option. (Not parsed by the filter itself, but by the command line tool). * * Beside the complete matrix, which is represented by its individual elements \f$m_{x,d}\f$, som eadditional * values are used for the definition. \f$N_g\f$ is the number of discrete grey levels, \f$ N_d\f$ the number * (or maximum value) of possible distances, and \f$N_s\f$ the total number of zones. * \f$m_{x,d}\f$ gives the number of connected areas with the discrete * grey level x and distance to the boarder of d. Corresponding, \f$p_{x,d} = \frac{m_{x,d}}{N_s} gives the relativ * probability of this matrix cell. \f$ \f$N_v\f$ is the number of voxels. In addition, the marginal * sums \f$m_{x,\cdot} = m_x = \sum_d m_{x,d} \f$ , representing the sum of all zones with a given intensity, and * sums \f$m_{\cdot, d} = m_d = \sum_x m_{x,d} \f$ , representing the sum of all zones with a given distance, are used. * The distance are given as the number of voxels to the border and the voxel intensity is given as the * bin number of the binned image, starting with 1. * * The following features are then defined: * - Grey Level Distance Zone::Small Distance Emphasis: * \f[ \textup{Small Distance Emphasis}= \frac{1}{N_s} \sum_d \frac{m_d}{d^2} \f] * - Grey Level Distance Zone::Large Distance Emphasis: * \f[ \textup{Large Distance Emphasis}= \frac{1}{N_s} \sum_d d^2 m_d \f] * - Grey Level Distance Zone::Low Grey Level Emphasis: * \f[ \textup{Low Grey Level Emphasis}= \frac{1}{N_s} \sum_x \frac{m_x}{x^2} \f] * - Grey Level Distance Zone::High Grey Level Emphasis: * \f[ \textup{High Grey Level Emphasis}= \frac{1}{N_s} \sum_x x^2 m_x \f] * - Grey Level Distance Zone::Small Distance Low Grey Level Emphasis: * \f[ \textup{Small Distance Low Grey Level Emphasis}= \frac{1}{N_s} \sum_x \sum_d \frac{ m_{x,d}}{x^2 d^2}\f] * - Grey Level Distance Zone::Small Distance High Grey Level Emphasis: * \f[ \textup{Small Distance High Grey Level Emphasis}= \frac{1}{N_s} \sum_x \sum_d \frac{x^2 m_{x,d}}{d^2}\f] * - Grey Level Distance Zone::Large Distance Low Grey Level Emphasis: * \f[ \textup{Large Distance Low Grey Level Emphasis}= \frac{1}{N_s} \sum_x \sum_d \frac{d^2 m_{x,d}}{x^2}\f] * - Grey Level Distance Zone::Large Distance High Grey Level Emphasis: * \f[ \textup{Large Distance High Grey Level Emphasis}= \frac{1}{N_s} \sum_x \sum_d \x^2 d^2 m_{x,d} \f] * - Grey Level Distance Zone::Grey Level Non-Uniformity: * \f[ \textup{Grey Level Non-Uniformity}= \frac{1}{N_s} \sum_x m_x^2 \f] * - Grey Level Distance Zone::Grey Level Non-Uniformity Normalized: * \f[ \textup{Grey Level Non-Uniformity Normalized}= \frac{1}{N_s^2} \sum_x m_x^2 \f] * - Grey Level Distance Zone::Zone Distance Non-Uniformity: * \f[ \textup{Grey Level Non-Uniformity}= \frac{1}{N_s} \sum_d m_d^2 \f] * - Grey Level Distance Zone::Zone Distance Non-Uniformity Normalized: * \f[ \textup{Grey Level Non-Uniformity Normalized}= \frac{1}{N_s^2} \sum_d m_d^2 \f] * - Grey Level Distance Zone::Zone Percentage: The ratio of zones to the possible zones: * \f[ \textup{Zone Percentage}= \frac{N_s}{N_v} \f] * - Grey Level Distance Zone::Grey Level Mean: * \f[ \textup{Grey Level Mean}= \mu_x = \sum_x \sum_d x p_{x,d} \f] * - Grey Level Distance Zone::Grey Level Variance: * \f[ \textup{Grey Level Variance} = \sum_x \sum_d \left(x - \mu_x \right)^2 p_{x,d} \f] * - Grey Level Distance Zone::Zone Distance Mean: * \f[ \textup{Zone Distance Mean}= \mu_d = \sum_x \sum_d d p_{x,d} \f] * - Grey Level Distance Zone::Zone Distance Variance: * \f[ \textup{Zone Distance Variance} = \sum_x \sum_d \left(d - \mu_d \right)^2 p_{x,d} \f] * - Grey Level Distance Zone::Zone Distance Entropy : * \f[ \textup{Zone Distance Entropy} = - \sum_x \sum_d p_{x,d} \textup{log}_2 ( p_{x,d} ) \f] * - Grey Level Distance Zone::Grey Level Entropy : * \f[ \textup{Grey Level Entropy} = - \sum_x \sum_d p_{x,d} \textup{log}_2 ( p_{x,d} ) \f] */ public: mitkClassMacro(GIFGreyLevelDistanceZone, AbstractGlobalImageFeature); itkFactorylessNewMacro(Self) itkCloneMacro(Self) GIFGreyLevelDistanceZone(); /** * \brief Calculates the Cooccurence-Matrix based features for this class. */ virtual FeatureListType CalculateFeatures(const Image::Pointer & image, const Image::Pointer &feature) override; /** * \brief Returns a list of the names of all features that are calculated from this class */ virtual FeatureNameListType GetFeatureNames() override; virtual void CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList); virtual void AddArguments(mitkCommandLineParser &parser); virtual std::string GetCurrentFeatureEncoding() override; + itkSetMacro(SmallNeighbourhoodForDistance, bool); + itkGetConstMacro(SmallNeighbourhoodForDistance, bool); + struct GIFGreyLevelDistanceZoneConfiguration { mitk::Image::Pointer distanceMask; mitk::IntensityQuantifier::Pointer Quantifier; unsigned int direction; double MinimumIntensity; double MaximumIntensity; int Bins; + bool SmallNeighbourhoodForDistance; std::string prefix; }; private: + bool m_SmallNeighbourhoodForDistance; }; } #endif //mitkGIFGreyLevelDistanceZone_h diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp index 74e875ab26..4fc7fabe14 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp @@ -1,470 +1,492 @@ #include // MITK #include #include #include #include #include // ITK #include #include #include #include #include namespace mitk{ struct GreyLevelDistanceZoneMatrixHolder { public: GreyLevelDistanceZoneMatrixHolder(mitk::IntensityQuantifier::Pointer quantifier, int number, int maxSize); int IntensityToIndex(double intensity); int m_NumberOfBins; int m_MaximumSize; int m_NumerOfVoxels; Eigen::MatrixXd m_Matrix; mitk::IntensityQuantifier::Pointer m_Quantifier; }; } static void MatrixFeaturesTo(mitk::GreyLevelDistanceZoneFeatures features, std::string prefix, mitk::GIFGreyLevelDistanceZone::FeatureListType &featureList); mitk::GreyLevelDistanceZoneMatrixHolder::GreyLevelDistanceZoneMatrixHolder(mitk::IntensityQuantifier::Pointer quantifier, int number, int maxSize) : m_NumberOfBins(number), m_MaximumSize(maxSize), m_NumerOfVoxels(0), m_Quantifier(quantifier) { m_Matrix.resize(number, maxSize); m_Matrix.fill(0); } int mitk::GreyLevelDistanceZoneMatrixHolder::IntensityToIndex(double intensity) { return m_Quantifier->IntensityToIndex(intensity); } template int CalculateGlSZMatrix(itk::Image* itkImage, itk::Image* mask, itk::Image* distanceImage, std::vector > offsets, bool estimateLargestRegion, mitk::GreyLevelDistanceZoneMatrixHolder &holder) { typedef itk::Image ImageType; typedef itk::Image MaskImageType; typedef typename ImageType::IndexType IndexType; typedef itk::ImageRegionIteratorWithIndex ConstIterType; typedef itk::ImageRegionIteratorWithIndex ConstMaskIterType; auto region = mask->GetLargestPossibleRegion(); typename MaskImageType::RegionType newRegion; newRegion.SetSize(region.GetSize()); newRegion.SetIndex(region.GetIndex()); ConstIterType imageIter(itkImage, itkImage->GetLargestPossibleRegion()); ConstMaskIterType maskIter(mask, mask->GetLargestPossibleRegion()); typename MaskImageType::Pointer visitedImage = MaskImageType::New(); visitedImage->SetRegions(newRegion); visitedImage->Allocate(); visitedImage->FillBuffer(0); int largestRegion = 0; holder.m_NumberOfBins = 0; while (!maskIter.IsAtEnd()) { if (maskIter.Value() > 0 ) { auto startIntensityIndex = holder.IntensityToIndex(imageIter.Value()); std::vector indices; indices.push_back(maskIter.GetIndex()); unsigned int steps = 0; int smallestDistance = 500; while (indices.size() > 0) { auto currentIndex = indices.back(); indices.pop_back(); if (!region.IsInside(currentIndex)) { continue; } auto wasVisited = visitedImage->GetPixel(currentIndex); auto newIntensityIndex = holder.IntensityToIndex(itkImage->GetPixel(currentIndex)); auto isInMask = mask->GetPixel(currentIndex); if ((isInMask > 0) && (newIntensityIndex == startIntensityIndex) && (wasVisited < 1)) { ++(holder.m_NumerOfVoxels); smallestDistance = (smallestDistance > distanceImage->GetPixel(currentIndex)) ? distanceImage->GetPixel(currentIndex) : smallestDistance; ++steps; visitedImage->SetPixel(currentIndex, 1); for (auto offset : offsets) { auto newIndex = currentIndex + offset; indices.push_back(newIndex); newIndex = currentIndex - offset; indices.push_back(newIndex); } } } if (steps > 0) { largestRegion = std::max(steps, largestRegion); steps = std::min(steps, holder.m_MaximumSize); if (!estimateLargestRegion) { holder.m_Matrix(startIntensityIndex, smallestDistance-1) += 1; } } } ++imageIter; ++maskIter; } return largestRegion; } template void itkErode2( itk::Image *sourceImage, mitk::Image::Pointer &resultImage, + bool useSmallNeighbourhood, + int direction, int &maxDistance) { typedef itk::Image ImageType; typedef unsigned short MaskType; typedef itk::Image MaskImageType; typename MaskImageType::Pointer distanceImage = MaskImageType::New(); distanceImage->SetRegions(sourceImage->GetLargestPossibleRegion()); distanceImage->SetOrigin(sourceImage->GetOrigin()); distanceImage->SetSpacing(sourceImage->GetSpacing()); distanceImage->SetDirection(sourceImage->GetDirection()); distanceImage->Allocate(); distanceImage->FillBuffer(std::numeric_limits::max()-1); typename ImageType::SizeType radius; radius.Fill(1); + if ((direction > 1) && (direction - 2 neighbourIter(radius, sourceImage, sourceImage->GetLargestPossibleRegion()); itk::NeighborhoodIterator distanceIter(radius, distanceImage, distanceImage->GetLargestPossibleRegion()); bool imageChanged = true; while (imageChanged) { imageChanged = false; maxDistance = 0; neighbourIter.GoToBegin(); distanceIter.GoToBegin(); while (!neighbourIter.IsAtEnd()) { MaskType oldDistance = distanceIter.GetCenterPixel(); maxDistance = std::max(maxDistance, oldDistance); if (neighbourIter.GetCenterPixel() < 1) { if (oldDistance > 0) { distanceIter.SetCenterPixel(0); imageChanged = true; } } else if (oldDistance>0) { MaskType minimumDistance = oldDistance; + auto index = distanceIter.GetIndex(); for (unsigned int i = 0; i < distanceIter.Size(); ++i) { + if (true) + { + auto cIndex = distanceIter.GetIndex(i); + int differentIndexes = 0; + for (int j = 0; j < VDimension; ++j) + { + differentIndexes += (index[j] == cIndex[j]) ? 0 : 1; + } + if (differentIndexes > 1) + { + continue; + } + } minimumDistance = std::min(minimumDistance, 1+distanceIter.GetPixel(i)); } if (minimumDistance != oldDistance) { distanceIter.SetCenterPixel(minimumDistance); imageChanged = true; } } ++neighbourIter; ++distanceIter; } } mitk::CastToMitkImage(distanceImage, resultImage); } -void erode(mitk::Image::Pointer input, mitk::Image::Pointer &output, int &maxDistance) +void erode(mitk::Image::Pointer input, mitk::Image::Pointer &output, bool useSmallNeighbourhood, int direction, int &maxDistance) { - AccessByItk_2(input, itkErode2, output, maxDistance); + AccessByItk_n(input, itkErode2, (output, useSmallNeighbourhood, direction, maxDistance)); } -void erodeAndAdd(mitk::Image::Pointer input, mitk::Image::Pointer& finalOutput, int &maxDistance) +void erodeAndAdd(mitk::Image::Pointer input, mitk::Image::Pointer& finalOutput, bool useSmallNeighbourhood, int direction, int &maxDistance) { maxDistance = 0; - erode(input, finalOutput, maxDistance); + erode(input, finalOutput, useSmallNeighbourhood, direction, maxDistance); } void static CalculateFeatures( mitk::GreyLevelDistanceZoneMatrixHolder &holder, mitk::GreyLevelDistanceZoneFeatures & results ) { auto SgzMatrix = holder.m_Matrix; auto pgzMatrix = holder.m_Matrix; auto pgMatrix = holder.m_Matrix; auto pzMatrix = holder.m_Matrix; double Ns = pgzMatrix.sum(); pgzMatrix /= Ns; pgMatrix.rowwise().normalize(); pzMatrix.colwise().normalize(); for (int i = 0; i < pgzMatrix.rows(); ++i) for (int j = 0; j < pgzMatrix.cols(); ++j) { if (pgzMatrix(i, j) != pgzMatrix(i, j)) pgzMatrix(i, j) = 0; if (pgMatrix(i, j) != pgMatrix(i, j)) pgMatrix(i, j) = 0; if (pzMatrix(i, j) != pzMatrix(i, j)) pzMatrix(i, j) = 0; } Eigen::VectorXd SgVector = SgzMatrix.rowwise().sum(); Eigen::VectorXd SzVector = SgzMatrix.colwise().sum(); for (int j = 0; j < SzVector.size(); ++j) { results.SmallDistanceEmphasis += SzVector(j) / (j+1) / (j+1); results.LargeDistanceEmphasis += SzVector(j) * (j + 1.0) * (j + 1.0); results.ZoneDistanceNonUniformity += SzVector(j) * SzVector(j); results.ZoneDistanceNoneUniformityNormalized += SzVector(j) * SzVector(j); } for (int i = 0; i < SgVector.size(); ++i) { results.LowGreyLevelEmphasis += SgVector(i) / (i + 1) / (i + 1); results.HighGreyLevelEmphasis += SgVector(i) * (i + 1) * (i + 1); results.GreyLevelNonUniformity += SgVector(i)*SgVector(i); results.GreyLevelNonUniformityNormalized += SgVector(i)*SgVector(i); } for (int i = 0; i < SgzMatrix.rows(); ++i) { for (int j = 0; j < SgzMatrix.cols(); ++j) { results.SmallDistanceLowGreyLevelEmphasis += SgzMatrix(i, j) / (i + 1) / (i + 1) / (j + 1) / (j + 1); results.SmallDistanceHighGreyLevelEmphasis += SgzMatrix(i, j) * (i + 1) * (i + 1) / (j + 1) / (j + 1); results.LargeDistanceLowGreyLevelEmphasis += SgzMatrix(i, j) / (i + 1) / (i + 1) * (j + 1.0) * (j + 1.0); results.LargeDistanceHighGreyLevelEmphasis += SgzMatrix(i, j) * (i + 1) * (i + 1) * (j + 1.0) * (j + 1.0); results.ZonePercentage += SgzMatrix(i, j); results.GreyLevelMean += (i + 1)*pgzMatrix(i, j); results.ZoneDistanceMean += (j + 1)*pgzMatrix(i, j); if (pgzMatrix(i, j) > 0) results.ZoneDistanceEntropy -= pgzMatrix(i, j) * std::log(pgzMatrix(i, j)) / std::log(2); } } for (int i = 0; i < SgzMatrix.rows(); ++i) { for (int j = 0; j < SgzMatrix.cols(); ++j) { results.GreyLevelVariance += (i + 1 - results.GreyLevelMean)*(i + 1 - results.GreyLevelMean)*pgzMatrix(i, j); results.ZoneDistanceVariance += (j + 1 - results.ZoneDistanceMean)*(j + 1 - results.ZoneDistanceMean)*pgzMatrix(i, j); } } results.SmallDistanceEmphasis /= Ns; results.LargeDistanceEmphasis /= Ns; results.LowGreyLevelEmphasis /= Ns; results.HighGreyLevelEmphasis /= Ns; results.SmallDistanceLowGreyLevelEmphasis /= Ns; results.SmallDistanceHighGreyLevelEmphasis /= Ns; results.LargeDistanceLowGreyLevelEmphasis /= Ns; results.LargeDistanceHighGreyLevelEmphasis /= Ns; results.GreyLevelNonUniformity /= Ns; results.GreyLevelNonUniformityNormalized /= Ns*Ns; results.ZoneDistanceNonUniformity /= Ns; results.ZoneDistanceNoneUniformityNormalized /= Ns*Ns; results.ZonePercentage = Ns / holder.m_NumerOfVoxels;// results.ZonePercentage; } template static void CalculateGreyLevelDistanceZoneFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFGreyLevelDistanceZone::FeatureListType & featureList, mitk::GIFGreyLevelDistanceZone::GIFGreyLevelDistanceZoneConfiguration config) { typedef itk::Image MaskType; typedef itk::Neighborhood NeighborhoodType; typedef itk::Offset OffsetType; /////////////////////////////////////////////////////////////////////////////////////////////// int maximumDistance = 0; mitk::Image::Pointer mitkDistanceImage = mitk::Image::New(); - erodeAndAdd(config.distanceMask, mitkDistanceImage, maximumDistance); + erodeAndAdd(config.distanceMask, mitkDistanceImage, config.SmallNeighbourhoodForDistance, config.direction, maximumDistance); typename MaskType::Pointer distanceImage = MaskType::New(); mitk::CastToItkImage(mitkDistanceImage, distanceImage); typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); //Find possible directions std::vector < itk::Offset > offsetVector; NeighborhoodType hood; hood.SetRadius(1); unsigned int centerIndex = hood.GetCenterNeighborhoodIndex(); OffsetType offset; for (unsigned int d = 0; d < centerIndex; d++) { offset = hood.GetOffset(d); bool useOffset = true; for (unsigned int i = 0; i < VImageDimension; ++i) { if ((config.direction == i + 2) && offset[i] != 0) { useOffset = false; } } if (useOffset) { offsetVector.push_back(offset); } } if (config.direction == 1) { offsetVector.clear(); offset[0] = 0; offset[1] = 0; offset[2] = 1; offsetVector.push_back(offset); } MITK_INFO << "Maximum Distance: " << maximumDistance; std::vector resultVector; mitk::GreyLevelDistanceZoneMatrixHolder holderOverall(config.Quantifier, config.Bins, maximumDistance + 1); mitk::GreyLevelDistanceZoneFeatures overallFeature; CalculateGlSZMatrix(itkImage, maskImage, distanceImage, offsetVector, false, holderOverall); CalculateFeatures(holderOverall, overallFeature); MatrixFeaturesTo(overallFeature, config.prefix, featureList); } static void MatrixFeaturesTo(mitk::GreyLevelDistanceZoneFeatures features, std::string prefix, mitk::GIFGreyLevelDistanceZone::FeatureListType &featureList) { featureList.push_back(std::make_pair(prefix + "Small Distance Emphasis", features.SmallDistanceEmphasis)); featureList.push_back(std::make_pair(prefix + "Large Distance Emphasis", features.LargeDistanceEmphasis)); featureList.push_back(std::make_pair(prefix + "Low Grey Level Emphasis", features.LowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "High Grey Level Emphasis", features.HighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Small Distance Low Grey Level Emphasis", features.SmallDistanceLowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Small Distance High Grey Level Emphasis", features.SmallDistanceHighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Large Distance Low Grey Level Emphasis", features.LargeDistanceLowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Large Distance High Grey Level Emphasis", features.LargeDistanceHighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity", features.GreyLevelNonUniformity)); featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity Normalized", features.GreyLevelNonUniformityNormalized)); featureList.push_back(std::make_pair(prefix + "Distance Size Non-Uniformity", features.ZoneDistanceNonUniformity)); featureList.push_back(std::make_pair(prefix + "Distance Size Non-Uniformity Normalized", features.ZoneDistanceNoneUniformityNormalized)); featureList.push_back(std::make_pair(prefix + "Zone Percentage", features.ZonePercentage)); featureList.push_back(std::make_pair(prefix + "Grey Level Mean", features.GreyLevelMean)); featureList.push_back(std::make_pair(prefix + "Grey Level Variance", features.GreyLevelVariance)); featureList.push_back(std::make_pair(prefix + "Zone Distance Mean", features.ZoneDistanceMean)); featureList.push_back(std::make_pair(prefix + "Zone Distance Variance", features.ZoneDistanceVariance)); featureList.push_back(std::make_pair(prefix + "Zone Distance Entropy", features.ZoneDistanceEntropy)); featureList.push_back(std::make_pair(prefix + "Grey Level Entropy", features.ZoneDistanceEntropy)); } - mitk::GIFGreyLevelDistanceZone::GIFGreyLevelDistanceZone() + mitk::GIFGreyLevelDistanceZone::GIFGreyLevelDistanceZone() : m_SmallNeighbourhoodForDistance(true) { SetShortName("gldz"); SetLongName("distance-zone"); SetFeatureClassName("Grey Level Distance Zone"); } mitk::GIFGreyLevelDistanceZone::FeatureListType mitk::GIFGreyLevelDistanceZone::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { InitializeQuantifier(image, mask); FeatureListType featureList; GIFGreyLevelDistanceZoneConfiguration config; config.direction = GetDirection(); - if (GetMorphMask().IsNull()) + if (GetMorphMask()==NULL) { config.distanceMask = mask->Clone(); } else { - config.distanceMask = GetMorphMask(); + config.distanceMask = GetMorphMask()->Clone(); } config.MinimumIntensity = GetQuantifier()->GetMinimum(); config.MaximumIntensity = GetQuantifier()->GetMaximum(); config.Bins = GetQuantifier()->GetBins(); config.prefix = FeatureDescriptionPrefix(); config.Quantifier = GetQuantifier(); + config.SmallNeighbourhoodForDistance = GetSmallNeighbourhoodForDistance(); AccessByItk_3(image, CalculateGreyLevelDistanceZoneFeatures, mask, featureList, config); return featureList; } mitk::GIFGreyLevelDistanceZone::FeatureNameListType mitk::GIFGreyLevelDistanceZone::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFGreyLevelDistanceZone::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Grey Level Distance Zone", "Calculates the size zone based features.", us::Any()); + AddQuantifierArguments(parser); } void mitk::GIFGreyLevelDistanceZone::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); std::string name = GetOptionPrefix(); if (parsedArgs.count(GetLongName())) { InitializeQuantifierFromParameters(feature, mask); MITK_INFO << "Start calculating Grey Level Distance Zone ...."; auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating Grey Level Distance Zone."; } } std::string mitk::GIFGreyLevelDistanceZone::GetCurrentFeatureEncoding() { return QuantifierParameterString(); } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyToneDifferenceFeatures.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyToneDifferenceFeatures.cpp index 70a6837f27..42b3e49707 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyToneDifferenceFeatures.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyToneDifferenceFeatures.cpp @@ -1,221 +1,228 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include #include #include // ITK #include #include // STL #include struct GIFNeighbourhoodGreyToneDifferenceParameter { int Range = 1; + int direction = 0; mitk::IntensityQuantifier::Pointer quantifier; std::string prefix; }; template static void CalculateIntensityPeak(itk::Image* itkImage, mitk::Image::Pointer mask, GIFNeighbourhoodGreyToneDifferenceParameter params, mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::FeatureListType & featureList) { typedef itk::Image ImageType; typedef itk::Image MaskType; typename MaskType::Pointer itkMask = MaskType::New(); mitk::CastToItkImage(mask, itkMask); typename ImageType::SizeType regionSize; regionSize.Fill(params.Range); + if ((params.direction > 1) && (params.direction - 2 iter(regionSize, itkImage, itkImage->GetLargestPossibleRegion()); itk::NeighborhoodIterator iterMask(regionSize, itkMask, itkMask->GetLargestPossibleRegion()); std::vector pVector; std::vector sVector; pVector.resize(params.quantifier->GetBins(), 0); sVector.resize(params.quantifier->GetBins(), 0); int count = 0; while (!iter.IsAtEnd()) { if (iterMask.GetCenterPixel() > 0) { int localCount = 0; double localMean = 0; unsigned int localIndex = params.quantifier->IntensityToIndex(iter.GetCenterPixel()); for (itk::SizeValueType i = 0; i < iter.Size(); ++i) { if (i == (iter.Size() / 2)) continue; if (iterMask.GetPixel(i) > 0) { ++localCount; localMean += params.quantifier->IntensityToIndex(iter.GetPixel(i)) + 1; } } if (localCount > 0) { localMean /= localCount; } localMean = std::abs(localIndex + 1 - localMean); pVector[localIndex] += 1; sVector[localIndex] += localMean; ++count; } ++iterMask; ++iter; } unsigned int Ngp = 0; for (unsigned int i = 0; i < params.quantifier->GetBins(); ++i) { if (pVector[i] > 0.1) { ++Ngp; } pVector[i] /= count; } double sumS = 0; double sumStimesP = 0; double contrastA = 0; double busynessA = 0; double complexity = 0; double strengthA = 0; for (unsigned int i = 0; i < params.quantifier->GetBins(); ++i) { sumS += sVector[i]; sumStimesP += pVector[i] * sVector[i]; for (unsigned int j = 0; j < params.quantifier->GetBins(); ++j) { double iMinusj = 1.0*i - 1.0*j; contrastA += pVector[i] * pVector[j] * iMinusj*iMinusj; if ((pVector[i] > 0) && (pVector[j] > 0)) { busynessA += std::abs((i + 1.0)*pVector[i] - (j + 1.0)*pVector[j]); complexity += std::abs(iMinusj)*(pVector[i] * sVector[i] + pVector[j] * sVector[j]) / (pVector[i] + pVector[j]); strengthA += (pVector[i] + pVector[j])*iMinusj*iMinusj; } } } double coarsness = 1.0 / std::min(sumStimesP, 1000000); double contrast = 0; double busyness = 0; if (Ngp > 1) { contrast = contrastA / Ngp / (Ngp - 1) / count * sumS; busyness = sumStimesP / busynessA; } complexity /= count; double strength = 0; if (sumS > 0) { strength = strengthA / sumS; } std::string prefix = params.prefix; featureList.push_back(std::make_pair(prefix + "Coarsness", coarsness)); featureList.push_back(std::make_pair(prefix + "Contrast", contrast)); featureList.push_back(std::make_pair(prefix + "Busyness", busyness)); featureList.push_back(std::make_pair(prefix + "Complexity", complexity)); featureList.push_back(std::make_pair(prefix + "Strength", strength)); } mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::GIFNeighbourhoodGreyToneDifferenceFeatures() : m_Range(1) { SetLongName("neighbourhood-grey-tone-difference"); SetShortName("ngtd"); SetFeatureClassName("Neighbourhood Grey Tone Difference"); } mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::FeatureListType mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; InitializeQuantifierFromParameters(image, mask); GIFNeighbourhoodGreyToneDifferenceParameter params; params.Range = GetRange(); params.quantifier = GetQuantifier(); params.prefix = FeatureDescriptionPrefix(); + params.direction = GetDirection(); AccessByItk_3(image, CalculateIntensityPeak, mask, params, featureList); return featureList; } mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::FeatureNameListType mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::AddArguments(mitkCommandLineParser &parser) { AddQuantifierArguments(parser); std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Neighbourhood Grey Tone Difference", "calculates Neighborhood Grey Tone based features", us::Any()); parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::Int, "Range for the local intensity", "Give the range that should be used for the local intensity in mm", us::Any()); } void mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &, FeatureListType &featureList) { InitializeQuantifierFromParameters(feature, mask); std::string name = GetOptionPrefix(); auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { MITK_INFO << "Start calculating Neighbourhood Grey Tone Difference features ...."; if (parsedArgs.count(name + "::range")) { int range = us::any_cast(parsedArgs[name + "::range"]); this->SetRange(range); } auto localResults = this->CalculateFeatures(feature, mask); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating Neighbourhood Grey Tone Difference features...."; } } std::string mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::GetCurrentFeatureEncoding() { std::ostringstream ss; ss << m_Range; std::string strRange = ss.str(); return QuantifierParameterString() + "_Range-" + ss.str(); } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricDensityStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricDensityStatistics.cpp index 9e184ae7db..37a2b91744 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricDensityStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricDensityStatistics.cpp @@ -1,531 +1,565 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include #include #include // ITK #include #include #include #include // VTK #include #include #include #include #include #include #include #include // STL #include #include // Eigen #include struct GIFVolumetricDensityStatisticsParameters { double volume; std::string prefix; }; template void -CalculateVolumeDensityStatistic(itk::Image* itkImage, mitk::Image::Pointer mask, GIFVolumetricDensityStatisticsParameters params, mitk::GIFVolumetricDensityStatistics::FeatureListType & featureList) +CalculateVolumeDensityStatistic(itk::Image* itkImage, mitk::Image::Pointer originalMask, mitk::Image::Pointer mask, GIFVolumetricDensityStatisticsParameters params, mitk::GIFVolumetricDensityStatistics::FeatureListType & featureList) { typedef itk::Image ImageType; typedef itk::Image MaskType; - double volume = params.volume; std::string prefix = params.prefix; typename MaskType::Pointer maskImage = MaskType::New(); + typename MaskType::Pointer originalMaskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); + mitk::CastToItkImage(originalMask, originalMaskImage); + itk::ImageRegionConstIteratorWithIndex imgA(itkImage, itkImage->GetLargestPossibleRegion()); - itk::ImageRegionConstIteratorWithIndex imgB(itkImage, itkImage->GetLargestPossibleRegion()); itk::ImageRegionConstIteratorWithIndex maskA(maskImage, maskImage->GetLargestPossibleRegion()); - itk::ImageRegionConstIteratorWithIndex maskB(maskImage, maskImage->GetLargestPossibleRegion()); + itk::ImageRegionConstIteratorWithIndex orignalMaskIter(originalMaskImage, originalMaskImage->GetLargestPossibleRegion()); double moranA = 0; double moranB = 0; double geary = 0; double Nv = 0; double w_ij = 0; double mean = 0; + unsigned int originalMaskVoxelCount = 0; typename ImageType::PointType pointA; typename ImageType::PointType pointB; while (!imgA.IsAtEnd()) { if (maskA.Get() > 0) { Nv += 1; mean += imgA.Get(); } + if (orignalMaskIter.Get() > 0) + { + ++originalMaskVoxelCount; + } ++imgA; ++maskA; + ++orignalMaskIter; } mean /= Nv; + double volume = originalMaskVoxelCount * originalMaskImage->GetSpacing()[0] * originalMaskImage->GetSpacing()[1] * originalMaskImage->GetSpacing()[2]; + imgA.GoToBegin(); maskA.GoToBegin(); + std::vector vectorPointsInMask; + std::vector vectorValuesInMask; + + vectorPointsInMask.reserve(Nv); + vectorValuesInMask.reserve(Nv); + while (!imgA.IsAtEnd()) { if (maskA.Get() > 0) { - imgB.GoToBegin(); - maskB.GoToBegin(); - while (!imgB.IsAtEnd()) - { - if ((imgA.GetIndex() == imgB.GetIndex()) || - (maskB.Get() < 1)) - { - ++imgB; - ++maskB; - continue; - } - itkImage->TransformIndexToPhysicalPoint(maskA.GetIndex(), pointA); - itkImage->TransformIndexToPhysicalPoint(maskB.GetIndex(), pointB); + auto maskAIndex = maskA.GetIndex(); + itkImage->TransformIndexToPhysicalPoint(maskAIndex, pointA); - double w = 1 / pointA.EuclideanDistanceTo(pointB); - moranA += w*(imgA.Get() - mean)* (imgB.Get() - mean); - geary += w * (imgA.Get() - imgB.Get()) * (imgA.Get() - imgB.Get()); - - w_ij += w; - - ++imgB; - ++maskB; - } - moranB += (imgA.Get() - mean)* (imgA.Get() - mean); + vectorPointsInMask.push_back(pointA); + vectorValuesInMask.push_back(imgA.Get()); } ++imgA; ++maskA; } + std::size_t elementsInVectors = vectorPointsInMask.size(); + + for (std::size_t idxA = 0; idxA < elementsInVectors; ++idxA) + { + pointA = vectorPointsInMask[idxA]; + double valueA = vectorValuesInMask[idxA]; + double valueAMinusMean = valueA - mean; + + for (std::size_t idxB = idxA+1; idxB < elementsInVectors; ++idxB) + { + pointB = vectorPointsInMask[idxB]; + double valueB = vectorValuesInMask[idxB]; + + double w = 2 * 1 / pointA.EuclideanDistanceTo(pointB); // Using symetry to avoid double calculation + moranA += w * valueAMinusMean * (valueB - mean); + geary += w * (valueA - valueB) * (valueA - valueB); + + w_ij += w; + } + moranB += valueAMinusMean * valueAMinusMean; + } MITK_INFO << "Volume: " << volume; MITK_INFO << " Mean: " << mean; featureList.push_back(std::make_pair(prefix + "Volume integrated intensity", volume* mean)); featureList.push_back(std::make_pair(prefix + "Volume Moran's I index", Nv / w_ij * moranA / moranB)); featureList.push_back(std::make_pair(prefix + "Volume Geary's C measure", ( Nv -1 ) / 2 / w_ij * geary/ moranB)); } void calculateMOBB(vtkPointSet *pointset, double &volume, double &surface) { volume = std::numeric_limits::max(); for (int cellID = 0; cellID < pointset->GetNumberOfCells(); ++cellID) { auto cell = pointset->GetCell(cellID); for (int edgeID = 0; edgeID < 3; ++edgeID) { auto edge = cell->GetEdge(edgeID); double pA[3], pB[3]; double pAA[3], pBB[3]; vtkSmartPointer transform = vtkSmartPointer::New(); transform->PostMultiply(); pointset->GetPoint(edge->GetPointId(0), pA); pointset->GetPoint(edge->GetPointId(1), pB); double angleZ = std::atan2((- pA[2] + pB[2]) ,(pA[1] - pB[1])); angleZ *= 180 / vnl_math::pi; if (pA[2] == pB[2]) angleZ = 0; transform->RotateX(angleZ); transform->TransformPoint(pA, pAA); transform->TransformPoint(pB, pBB); double angleY = std::atan2((pAA[1] -pBB[1]) ,-(pAA[0] - pBB[0])); angleY *= 180 / vnl_math::pi; if (pAA[1] == pBB[1]) angleY = 0; transform->RotateZ(angleY); double p0[3]; pointset->GetPoint(edge->GetPointId(0), p0); double curMinX = std::numeric_limits::max(); double curMaxX = std::numeric_limits::lowest(); double curMinY = std::numeric_limits::max(); double curMaxY = std::numeric_limits::lowest(); double curMinZ = std::numeric_limits::max(); double curMaxZ = std::numeric_limits::lowest(); for (int pointID = 0; pointID < pointset->GetNumberOfPoints(); ++pointID) { double p[3]; double p2[3]; pointset->GetPoint(pointID, p); p[0] -= p0[0]; p[1] -= p0[1]; p[2] -= p0[2]; transform->TransformPoint(p, p2); curMinX = std::min(p2[0], curMinX); curMaxX = std::max(p2[0], curMaxX); curMinY = std::min(p2[1], curMinY); curMaxY = std::max(p2[1], curMaxY); curMinZ = std::min(p2[2], curMinZ); curMaxZ = std::max(p2[2], curMaxZ); } if ((curMaxX - curMinX)*(curMaxY - curMinY)*(curMaxZ - curMinZ) < volume) { volume = (curMaxX - curMinX)*(curMaxY - curMinY)*(curMaxZ - curMinZ); surface = (curMaxX - curMinX)*(curMaxX - curMinX) + (curMaxY - curMinY)*(curMaxY - curMinY) + (curMaxZ - curMinZ)*(curMaxZ - curMinZ); surface *= 2; } } } } -void calculateMEE(vtkPointSet *pointset, double &vol, double &surf, double tolerance=0.0001) +void calculateMEE(vtkPointSet *pointset, double &vol, double &surf, double tolerance=0.05) { // Inspired by https://github.com/smdabdoub/ProkaryMetrics/blob/master/calc/fitting.py int numberOfPoints = pointset->GetNumberOfPoints(); int dimension = 3; Eigen::MatrixXd points(3, numberOfPoints); Eigen::MatrixXd Q(3+1, numberOfPoints); double p[3]; std::cout << "Initialize Q " << std::endl; for (int i = 0; i < numberOfPoints; ++i) { pointset->GetPoint(i, p); points(0, i) = p[0]; points(1, i) = p[1]; points(2, i) = p[2]; Q(0, i) = p[0]; Q(1, i) = p[1]; Q(2, i) = p[2]; Q(3, i) = 1.0; } int count = 1; double error = 1; Eigen::VectorXd u_vector(numberOfPoints); u_vector.fill(1.0 / numberOfPoints); Eigen::DiagonalMatrix u = u_vector.asDiagonal(); Eigen::VectorXd ones(dimension + 1); ones.fill(1); Eigen::MatrixXd Ones = ones.asDiagonal(); // Khachiyan Algorithm while (error > tolerance) { auto Qt = Q.transpose(); Eigen::MatrixXd X = Q*u*Qt; Eigen::FullPivHouseholderQR qr(X); Eigen::MatrixXd Xi = qr.solve(Ones); Eigen::MatrixXd M = Qt * Xi * Q; double maximumValue = M(0, 0); int maximumPosition = 0; for (int i = 0; i < numberOfPoints; ++i) { if (maximumValue < M(i, i)) { maximumValue = M(i, i); maximumPosition = i; } } double stepsize = (maximumValue - dimension - 1) / ((dimension + 1) * (maximumValue - 1)); Eigen::DiagonalMatrix new_u = (1.0 - stepsize) * u; new_u.diagonal()[maximumPosition] = (new_u.diagonal())(maximumPosition) + stepsize; ++count; error = (new_u.diagonal() - u.diagonal()).norm(); u.diagonal() = new_u.diagonal(); } // U = u Eigen::MatrixXd Ai = points * u * points.transpose() - points * u *(points * u).transpose(); Eigen::FullPivHouseholderQR qr(Ai); Eigen::VectorXd ones2(dimension); ones2.fill(1); Eigen::MatrixXd Ones2 = ones2.asDiagonal(); Eigen::MatrixXd A = qr.solve(Ones2)*1.0/dimension; Eigen::JacobiSVD svd(A); double c = 1 / sqrt(svd.singularValues()[0]); double b = 1 / sqrt(svd.singularValues()[1]); double a = 1 / sqrt(svd.singularValues()[2]); double V = 4 * vnl_math::pi*a*b*c / 3; double ad_mvee= 0; double alpha = std::sqrt(1 - b*b / a / a); double beta = std::sqrt(1 - c*c / a / a); for (int i = 0; i < 20; ++i) { ad_mvee += 4 * vnl_math::pi*a*b*(alpha*alpha + beta*beta) / (2 * alpha*beta) * (std::pow(alpha*beta, i)) / (1 - 4 * i*i); } vol = V; surf = ad_mvee; } mitk::GIFVolumetricDensityStatistics::FeatureListType mitk::GIFVolumetricDensityStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { + // Check if Image is at least of dimension 3, otherwise return without calculating + // as features only defined in 3D FeatureListType featureList; if (image->GetDimension() < 3) { return featureList; } + mitk::Image::Pointer original_mask = GetMorphMask(); + + std::string prefix = FeatureDescriptionPrefix(); vtkSmartPointer mesher = vtkSmartPointer::New(); vtkSmartPointer stats = vtkSmartPointer::New(); vtkSmartPointer stats2 = vtkSmartPointer::New(); - mesher->SetInputData(mask->GetVtkImageData()); + mesher->SetInputData(original_mask->GetVtkImageData()); mesher->SetValue(0, 0.5); stats->SetInputConnection(mesher->GetOutputPort()); stats->Update(); vtkSmartPointer delaunay = vtkSmartPointer< vtkDelaunay3D >::New(); delaunay->SetInputConnection(mesher->GetOutputPort()); delaunay->SetAlpha(0); delaunay->Update(); vtkSmartPointer geometryFilter = vtkSmartPointer::New(); geometryFilter->SetInputConnection(delaunay->GetOutputPort()); geometryFilter->Update(); stats2->SetInputConnection(geometryFilter->GetOutputPort()); stats2->Update(); - double vol_mvee; - double surf_mvee; + double vol_mvee=0; + double surf_mvee=0; calculateMEE(mesher->GetOutput(), vol_mvee, surf_mvee); - double vol_mobb; - double surf_mobb; - calculateMOBB(geometryFilter->GetOutput(), vol_mobb, surf_mobb); + double vol_mobb=0; + double surf_mobb=0; + //calculateMOBB(geometryFilter->GetOutput(), vol_mobb, surf_mobb); double pi = vnl_math::pi; double meshVolume = stats->GetVolume(); double meshSurf = stats->GetSurfaceArea(); GIFVolumetricDensityStatisticsParameters params; params.volume = meshVolume; params.prefix = prefix; - AccessByItk_3(image, CalculateVolumeDensityStatistic, mask, params, featureList); + AccessByItk_n(image, CalculateVolumeDensityStatistic, (original_mask, mask, params, featureList)); //Calculate center of mass shift int xx = mask->GetDimensions()[0]; int yy = mask->GetDimensions()[1]; int zz = mask->GetDimensions()[2]; - double xd = mask->GetGeometry()->GetSpacing()[0]; - double yd = mask->GetGeometry()->GetSpacing()[1]; - double zd = mask->GetGeometry()->GetSpacing()[2]; + double spacingX = mask->GetGeometry()->GetSpacing()[0]; + double spacingY = mask->GetGeometry()->GetSpacing()[1]; + double spacingZ = mask->GetGeometry()->GetSpacing()[2]; int minimumX=xx; int maximumX=0; int minimumY=yy; int maximumY=0; int minimumZ=zz; int maximumZ=0; vtkSmartPointer dataset1Arr = vtkSmartPointer::New(); vtkSmartPointer dataset2Arr = vtkSmartPointer::New(); vtkSmartPointer dataset3Arr = vtkSmartPointer::New(); dataset1Arr->SetNumberOfComponents(1); dataset2Arr->SetNumberOfComponents(1); dataset3Arr->SetNumberOfComponents(1); dataset1Arr->SetName("M1"); dataset2Arr->SetName("M2"); dataset3Arr->SetName("M3"); vtkSmartPointer dataset1ArrU = vtkSmartPointer::New(); vtkSmartPointer dataset2ArrU = vtkSmartPointer::New(); vtkSmartPointer dataset3ArrU = vtkSmartPointer::New(); dataset1ArrU->SetNumberOfComponents(1); dataset2ArrU->SetNumberOfComponents(1); dataset3ArrU->SetNumberOfComponents(1); dataset1ArrU->SetName("M1"); dataset2ArrU->SetName("M2"); dataset3ArrU->SetName("M3"); vtkSmartPointer points = vtkSmartPointer< vtkPoints >::New(); for (int x = 0; x < xx; x++) { for (int y = 0; y < yy; y++) { for (int z = 0; z < zz; z++) { itk::Image::IndexType index; index[0] = x; index[1] = y; index[2] = z; mitk::ScalarType pxImage; mitk::ScalarType pxMask; + mitk::ScalarType valueOriginalMask; mitkPixelTypeMultiplex5( mitk::FastSinglePixelAccess, image->GetChannelDescriptor().GetPixelType(), image, image->GetVolumeData(), index, pxImage, 0); mitkPixelTypeMultiplex5( mitk::FastSinglePixelAccess, mask->GetChannelDescriptor().GetPixelType(), mask, mask->GetVolumeData(), index, pxMask, 0); + mitkPixelTypeMultiplex5( + mitk::FastSinglePixelAccess, + mask->GetChannelDescriptor().GetPixelType(), + mask, + mask->GetVolumeData(), + index, + valueOriginalMask, + 0); + //Check if voxel is contained in segmentation - if (pxMask > 0) + if (valueOriginalMask > 0) { minimumX = std::min(x, minimumX); minimumY = std::min(y, minimumY); minimumZ = std::min(z, minimumZ); maximumX = std::max(x, maximumX); maximumY = std::max(y, maximumY); maximumZ = std::max(z, maximumZ); - points->InsertNextPoint(x*xd, y*yd, z*zd); - + } + if (pxMask > 0) + { + points->InsertNextPoint(x*spacingX, y*spacingY, z*spacingZ); if (pxImage == pxImage) { - dataset1Arr->InsertNextValue(x*xd); - dataset2Arr->InsertNextValue(y*yd); - dataset3Arr->InsertNextValue(z*zd); + dataset1Arr->InsertNextValue(x*spacingX); + dataset2Arr->InsertNextValue(y*spacingY); + dataset3Arr->InsertNextValue(z*spacingZ); } } } } } vtkSmartPointer datasetTable = vtkSmartPointer::New(); datasetTable->AddColumn(dataset1Arr); datasetTable->AddColumn(dataset2Arr); datasetTable->AddColumn(dataset3Arr); vtkSmartPointer pcaStatistics = vtkSmartPointer::New(); pcaStatistics->SetInputData(vtkStatisticsAlgorithm::INPUT_DATA, datasetTable); pcaStatistics->SetColumnStatus("M1", 1); pcaStatistics->SetColumnStatus("M2", 1); pcaStatistics->SetColumnStatus("M3", 1); pcaStatistics->RequestSelectedColumns(); pcaStatistics->SetDeriveOption(true); pcaStatistics->Update(); vtkSmartPointer eigenvalues = vtkSmartPointer::New(); pcaStatistics->GetEigenvalues(eigenvalues); std::vector eigen_val(3); eigen_val[2] = eigenvalues->GetValue(0); eigen_val[1] = eigenvalues->GetValue(1); eigen_val[0] = eigenvalues->GetValue(2); double major = 2*sqrt(eigen_val[2]); double minor = 2*sqrt(eigen_val[1]); double least = 2*sqrt(eigen_val[0]); double alpha = std::sqrt(1 - minor*minor / major / major); double beta = std::sqrt(1 - least*least / major / major); - double a = (maximumX - minimumX+1) * xd; - double b = (maximumY - minimumY+1) * yd; - double c = (maximumZ - minimumZ+1) * zd; + double a = (maximumX - minimumX+1) * spacingX; + double b = (maximumY - minimumY+1) * spacingY; + double c = (maximumZ - minimumZ+1) * spacingZ; double vd_aabb = meshVolume / (a*b*c); double ad_aabb = meshSurf / (2 * a*b + 2 * a*c + 2 * b*c); double vd_aee = 3 * meshVolume / (4.0*pi*major*minor*least); double ad_aee = 0; for (int i = 0; i < 20; ++i) { ad_aee += 4 * pi*major*minor*(alpha*alpha + beta*beta) / (2 * alpha*beta) * (std::pow(alpha*beta, i)) / (1 - 4 * i*i); } ad_aee = meshSurf / ad_aee; double vd_ch = meshVolume / stats2->GetVolume(); double ad_ch = meshSurf / stats2->GetSurfaceArea(); featureList.push_back(std::make_pair(prefix + "Volume density axis-aligned bounding box", vd_aabb)); featureList.push_back(std::make_pair(prefix + "Surface density axis-aligned bounding box", ad_aabb)); - featureList.push_back(std::make_pair(prefix + "Volume density oriented minimum bounding box", meshVolume / vol_mobb)); - featureList.push_back(std::make_pair(prefix + "Surface density oriented minimum bounding box", meshSurf / surf_mobb)); +// featureList.push_back(std::make_pair(prefix + "Volume density oriented minimum bounding box", meshVolume / vol_mobb)); +// featureList.push_back(std::make_pair(prefix + "Surface density oriented minimum bounding box", meshSurf / surf_mobb)); featureList.push_back(std::make_pair(prefix + "Volume density approx. enclosing ellipsoid", vd_aee)); featureList.push_back(std::make_pair(prefix + "Surface density approx. enclosing ellipsoid", ad_aee)); featureList.push_back(std::make_pair(prefix + "Volume density approx. minimum volume enclosing ellipsoid", meshVolume / vol_mvee)); featureList.push_back(std::make_pair(prefix + "Surface density approx. minimum volume enclosing ellipsoid", meshSurf / surf_mvee)); featureList.push_back(std::make_pair(prefix + "Volume density convex hull", vd_ch)); featureList.push_back(std::make_pair(prefix + "Surface density convex hull", ad_ch)); return featureList; } mitk::GIFVolumetricDensityStatistics::GIFVolumetricDensityStatistics() { SetLongName("volume-density"); SetShortName("volden"); SetFeatureClassName("Morphological Density"); } mitk::GIFVolumetricDensityStatistics::FeatureNameListType mitk::GIFVolumetricDensityStatistics::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFVolumetricDensityStatistics::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Volume-Density Statistic", "calculates volume density based features", us::Any()); } void mitk::GIFVolumetricDensityStatistics::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &, FeatureListType &featureList) { auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { MITK_INFO << "Start calculating volumetric density features ...."; auto localResults = this->CalculateFeatures(feature, mask); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating volumetric density features...."; } } diff --git a/Modules/MatchPointRegistration/Helper/mitkImageMappingHelper.cpp b/Modules/MatchPointRegistration/Helper/mitkImageMappingHelper.cpp index eb7b353ef3..b755e551e3 100644 --- a/Modules/MatchPointRegistration/Helper/mitkImageMappingHelper.cpp +++ b/Modules/MatchPointRegistration/Helper/mitkImageMappingHelper.cpp @@ -1,393 +1,393 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include "mapRegistration.h" #include "mitkImageMappingHelper.h" #include "mitkRegistrationHelper.h" template typename ::itk::InterpolateImageFunction< TImage >::Pointer generateInterpolator(mitk::ImageMappingInterpolator::Type interpolatorType) { typedef ::itk::InterpolateImageFunction< TImage > BaseInterpolatorType; typename BaseInterpolatorType::Pointer result; switch (interpolatorType) { case mitk::ImageMappingInterpolator::NearestNeighbor: { result = ::itk::NearestNeighborInterpolateImageFunction::New(); break; } case mitk::ImageMappingInterpolator::BSpline_3: { typename ::itk::BSplineInterpolateImageFunction::Pointer spInterpolator = ::itk::BSplineInterpolateImageFunction::New(); spInterpolator->SetSplineOrder(3); result = spInterpolator; break; } case mitk::ImageMappingInterpolator::WSinc_Hamming: { result = ::itk::WindowedSincInterpolateImageFunction::New(); break; } case mitk::ImageMappingInterpolator::WSinc_Welch: { result = ::itk::WindowedSincInterpolateImageFunction >::New(); break; } default: { result = ::itk::LinearInterpolateImageFunction::New(); break; } } return result; }; template void doMITKMap(const ::itk::Image* input, mitk::ImageMappingHelper::ResultImageType::Pointer& result, const mitk::ImageMappingHelper::RegistrationType*& registration, bool throwOnOutOfInputAreaError, const double& paddingValue, const mitk::ImageMappingHelper::ResultImageGeometryType*& resultGeometry, bool throwOnMappingError, const double& errorValue, mitk::ImageMappingInterpolator::Type interpolatorType) { typedef ::map::core::Registration ConcreteRegistrationType; typedef ::map::core::ImageMappingTask, ::itk::Image > MappingTaskType; typename MappingTaskType::Pointer spTask = MappingTaskType::New(); typedef typename MappingTaskType::ResultImageDescriptorType ResultImageDescriptorType; typename ResultImageDescriptorType::Pointer resultDescriptor; //check if image and result geometry fits the passed registration ///////////////////////////////////////////////////////////////// if (registration->getMovingDimensions()!=VImageDimension) { map::core::OStringStream str; str << "Dimension of MITK image ("<getMovingDimensions()<<")."; throw mitk::AccessByItkException(str.str()); } if (registration->getTargetDimensions()!=VImageDimension) { map::core::OStringStream str; str << "Dimension of MITK image ("<getTargetDimensions()<<")."; throw mitk::AccessByItkException(str.str()); } const ConcreteRegistrationType* castedReg = dynamic_cast(registration); if (registration->getTargetDimensions()==2 && resultGeometry) { mitk::ImageMappingHelper::ResultImageGeometryType::BoundsArrayType bounds = resultGeometry->GetBounds(); if (bounds[4]!=0 || bounds[5]!=0) { //array "bounds" is constructed as [min Dim1, max Dim1, min Dim2, max Dim2, min Dim3, max Dim3] //therfore [4] and [5] must be 0 map::core::OStringStream str; str << "Dimension of defined result geometry does not equal the target dimension of the registration object ("<getTargetDimensions()<<")."; throw mitk::AccessByItkException(str.str()); } } //check/create resultDescriptor ///////////////////////// if (resultGeometry) { resultDescriptor = ResultImageDescriptorType::New(); typename ResultImageDescriptorType::PointType origin; typename ResultImageDescriptorType::SizeType size; typename ResultImageDescriptorType::SpacingType fieldSpacing; typename ResultImageDescriptorType::DirectionType matrix; mitk::ImageMappingHelper::ResultImageGeometryType::BoundsArrayType geoBounds = resultGeometry->GetBounds(); mitk::Vector3D geoSpacing = resultGeometry->GetSpacing(); mitk::Point3D geoOrigin = resultGeometry->GetOrigin(); mitk::AffineTransform3D::MatrixType geoMatrix = resultGeometry->GetIndexToWorldTransform()->GetMatrix(); for (unsigned int i = 0; i(geoOrigin[i]); fieldSpacing[i] = static_cast(geoSpacing[i]); size[i] = static_cast(geoBounds[(2*i)+1]-geoBounds[2*i])*fieldSpacing[i]; } //Matrix extraction matrix.SetIdentity(); unsigned int i; unsigned int j; /// \warning 2D MITK images could have a 3D rotation, since they have a 3x3 geometry matrix. /// If it is only a rotation around the transversal plane normal, it can be express with a 2x2 matrix. /// In this case, the ITK image conservs this information and is identical to the MITK image! /// If the MITK image contains any other rotation, the ITK image will have no rotation at all. /// Spacing is of course conserved in both cases. // the following loop devides by spacing now to normalize columns. // counterpart of InitializeByItk in mitkImage.h line 372 of revision 15092. // Check if information is lost if ( VImageDimension == 2) { if ( ( geoMatrix[0][2] != 0) || ( geoMatrix[1][2] != 0) || ( geoMatrix[2][0] != 0) || ( geoMatrix[2][1] != 0) || (( geoMatrix[2][2] != 1) && ( geoMatrix[2][2] != -1) )) { // The 2D MITK image contains 3D rotation information. // This cannot be expressed in a 2D ITK image, so the ITK image will have no rotation } else { // The 2D MITK image can be converted to an 2D ITK image without information loss! for ( i=0; i < 2; ++i) { for( j=0; j < 2; ++j ) { matrix[i][j] = geoMatrix[i][j]/fieldSpacing[j]; } } } } else if (VImageDimension == 3) { // Normal 3D image. Conversion possible without problem! for ( i=0; i < 3; ++i) { for( j=0; j < 3; ++j ) { matrix[i][j] = geoMatrix[i][j]/fieldSpacing[j]; } } } else { assert(0); throw mitk::AccessByItkException("Usage of resultGeometry for 2D images is not yet implemented."); /**@TODO Implement extraction of 2D-Rotation-Matrix out of 3D-Rotation-Matrix * to cover this case as well. * matrix = extract2DRotationMatrix(resultGeometry)*/ } resultDescriptor->setOrigin(origin); resultDescriptor->setSize(size); resultDescriptor->setSpacing(fieldSpacing); resultDescriptor->setDirection(matrix); } //do the mapping ///////////////////////// typedef ::itk::InterpolateImageFunction< ::itk::Image > BaseInterpolatorType; typename BaseInterpolatorType::Pointer interpolator = generateInterpolator< ::itk::Image >(interpolatorType); assert(interpolator.IsNotNull()); spTask->setImageInterpolator(interpolator); spTask->setInputImage(input); spTask->setRegistration(castedReg); spTask->setResultImageDescriptor(resultDescriptor); spTask->setThrowOnMappingError(throwOnMappingError); spTask->setErrorValue(errorValue); spTask->setThrowOnPaddingError(throwOnOutOfInputAreaError); spTask->setPaddingValue(paddingValue); spTask->execute(); mitk::CastToMitkImage<>(spTask->getResultImage(),result); } mitk::ImageMappingHelper::ResultImageType::Pointer mitk::ImageMappingHelper::map(const InputImageType* input, const RegistrationType* registration, bool throwOnOutOfInputAreaError, const double& paddingValue, const ResultImageGeometryType* resultGeometry, bool throwOnMappingError, const double& errorValue, mitk::ImageMappingInterpolator::Type interpolatorType) { if (!registration) { mitkThrow() << "Cannot map image. Passed registration wrapper pointer is nullptr."; } if (!input) { mitkThrow() << "Cannot map image. Passed image pointer is nullptr."; } ResultImageType::Pointer result; if(input->GetTimeSteps()==1) { //map the image and done AccessByItk_n(input, doMITKMap, (result, registration, throwOnOutOfInputAreaError, paddingValue, resultGeometry, throwOnMappingError, errorValue, interpolatorType)); } else { //map every time step and compose mitk::TimeGeometry::ConstPointer timeGeometry = input->GetTimeGeometry(); mitk::TimeGeometry::Pointer mappedTimeGeometry = timeGeometry->Clone(); for (unsigned int i = 0; iGetTimeSteps(); ++i) { ResultImageGeometryType::Pointer mappedGeometry = resultGeometry->Clone(); mappedTimeGeometry->SetTimeStepGeometry(mappedGeometry,i); } result = mitk::Image::New(); result->Initialize(input->GetPixelType(),*mappedTimeGeometry, 1, input->GetTimeSteps()); for (unsigned int i = 0; iGetTimeSteps(); ++i) { mitk::ImageTimeSelector::Pointer imageTimeSelector = mitk::ImageTimeSelector::New(); imageTimeSelector->SetInput(input); imageTimeSelector->SetTimeNr(i); imageTimeSelector->UpdateLargestPossibleRegion(); InputImageType::Pointer timeStepInput = imageTimeSelector->GetOutput(); ResultImageType::Pointer timeStepResult; AccessByItk_n(timeStepInput, doMITKMap, (timeStepResult, registration, throwOnOutOfInputAreaError, paddingValue, resultGeometry, throwOnMappingError, errorValue, interpolatorType)); mitk::ImageReadAccessor readAccess(timeStepResult); result->SetVolume(readAccess.GetData(),i); } } return result; } mitk::ImageMappingHelper::ResultImageType::Pointer mitk::ImageMappingHelper::map(const InputImageType* input, const MITKRegistrationType* registration, bool throwOnOutOfInputAreaError, const double& paddingValue, const ResultImageGeometryType* resultGeometry, - bool throwOnMappingError, const double& errorValue, mitk::ImageMappingInterpolator::Type) + bool throwOnMappingError, const double& errorValue, mitk::ImageMappingInterpolator::Type interpolator) { if (!registration) { mitkThrow() << "Cannot map image. Passed registration wrapper pointer is nullptr."; } if (!registration->GetRegistration()) { mitkThrow() << "Cannot map image. Passed registration wrapper containes no registration."; } if (!input) { mitkThrow() << "Cannot map image. Passed image pointer is nullptr."; } - ResultImageType::Pointer result = map(input, registration->GetRegistration(), throwOnOutOfInputAreaError, paddingValue, resultGeometry, throwOnMappingError, errorValue); + ResultImageType::Pointer result = map(input, registration->GetRegistration(), throwOnOutOfInputAreaError, paddingValue, resultGeometry, throwOnMappingError, errorValue, interpolator); return result; } mitk::ImageMappingHelper::ResultImageType::Pointer mitk::ImageMappingHelper:: refineGeometry(const InputImageType* input, const RegistrationType* registration, bool throwOnError) { mitk::ImageMappingHelper::ResultImageType::Pointer result = nullptr; if (!registration) { mitkThrow() << "Cannot refine image geometry. Passed registration pointer is nullptr."; } if (!input) { mitkThrow() << "Cannot refine image geometry. Passed image pointer is nullptr."; } mitk::MITKRegistrationHelper::Affine3DTransformType::Pointer spTransform = mitk::MITKRegistrationHelper::getAffineMatrix(registration,false); if(spTransform.IsNull() && throwOnError) { mitkThrow() << "Cannot refine image geometry. Registration does not contain a suitable direct mapping kernel (3D affine transformation or compatible required)."; } if(spTransform.IsNotNull()) { //copy input image result = input->Clone(); //refine geometries for(unsigned int i = 0; i < result->GetTimeSteps(); ++i) { //refine every time step result->GetGeometry(i)->Compose(spTransform); } result->GetTimeGeometry()->Update(); } return result; } mitk::ImageMappingHelper::ResultImageType::Pointer mitk::ImageMappingHelper:: refineGeometry(const InputImageType* input, const MITKRegistrationType* registration, bool throwOnError) { if (!registration) { mitkThrow() << "Cannot refine image geometry. Passed registration wrapper pointer is nullptr."; } if (!registration->GetRegistration()) { mitkThrow() << "Cannot refine image geometry. Passed registration wrapper containes no registration."; } if (!input) { mitkThrow() << "Cannot refine image geometry. Passed image pointer is nullptr."; } ResultImageType::Pointer result = refineGeometry(input, registration->GetRegistration(), throwOnError); return result; } bool mitk::ImageMappingHelper:: canRefineGeometry(const RegistrationType* registration) { bool result = true; if (!registration) { mitkThrow() << "Cannot check refine capability of registration. Passed registration pointer is nullptr."; } //if the helper does not return null, we can refine the geometry. result = mitk::MITKRegistrationHelper::getAffineMatrix(registration,false).IsNotNull(); return result; } bool mitk::ImageMappingHelper:: canRefineGeometry(const MITKRegistrationType* registration) { if (!registration) { mitkThrow() << "Cannot check refine capability of registration. Passed registration wrapper pointer is nullptr."; } if (!registration->GetRegistration()) { mitkThrow() << "Cannot check refine capability of registration. Passed registration wrapper containes no registration."; } return canRefineGeometry(registration->GetRegistration()); } diff --git a/Wrapping/Common/mitk_swig_classes.i b/Wrapping/Common/mitk_swig_classes.i index 088820d94e..784530af17 100644 --- a/Wrapping/Common/mitk_swig_classes.i +++ b/Wrapping/Common/mitk_swig_classes.i @@ -1,158 +1,202 @@ // // Defining some Macros that make problems with SWIG as the // corresponding definitions are not included by default. // Luckely, these includes are not necessary for SWIG. // #define ITK_NOEXCEPT #define ITKCommon_EXPORT #define ITK_OVERRIDE #define MITKCORE_EXPORT #define MITKCLCORE_EXPORT #define MITKCLUTILITIES_EXPORT #define ITKCommon_EXPORT #define MITKMATCHPOINTREGISTRATION_EXPORT #define MAPDeployment_EXPORT #define MAPAlgorithms_EXPORT #define MITKSEGMENTATION_EXPORT #define MITKMULTILABEL_EXPORT +#define MITKBASICIMAGEPROCESSING_EXPORT + +#define VCL_TEMPLATE_EXPORT #define ITKCommon_EXPORT #define ITK_FORWARD_EXPORT #define ITK_OVERRIDE #define ITK_NOEXCEPT %include %include %include %include %include %include #define DEPRECATED(func) func #undef ITK_DISALLOW_COPY_AND_ASSIGN #define ITK_DISALLOW_COPY_AND_ASSIGN(TypeName) +typedef double ScalarType; +typedef Vector3D mitk::Vector3D; + + %pythoncode %{ convertion_list = {} %} SWIG_ADD_NONOBJECT_CLASS(Indent, itkIndent.h, itk) SWIG_ADD_MITK_CLASS(LightObject, itkLightObject.h, itk) SWIG_ADD_MITK_CLASS(Object, itkObject.h, itk) SWIG_ADD_MITK_CLASS(DataObject, itkDataObject.h, itk) SWIG_ADD_MITK_CLASS(TimeGeometry, mitkTimeGeometry.h, mitk) SWIG_ADD_MITK_CLASS(ArbitraryTimeGeometry, mitkArbitraryTimeGeometry.h, mitk) SWIG_ADD_MITK_CLASS(ProportionalTimeGeometry, mitkProportionalTimeGeometry.h, mitk) SWIG_ADD_MITK_CLASS(BaseGeometry, mitkBaseGeometry.h, mitk) SWIG_ADD_MITK_CLASS(Geometry3D, mitkGeometry3D.h, mitk) SWIG_ADD_MITK_CLASS(SlicedGeometry3D, mitkSlicedGeometry3D.h, mitk) SWIG_ADD_MITK_CLASS(PlaneGeometry , mitkPlaneGeometry.h, mitk) SWIG_ADD_NONOBJECT_NOVECTOR_CLASS(BoundingBox, mitkBaseGeometry.h, mitk) SWIG_ADD_NONOBJECT_CLASS(TimeBounds, mitkBaseGeometry.h, mitk) SWIG_ADD_NONOBJECT_CLASS(FixedArrayType, mitkBaseGeometry.h, mitk) SWIG_ADD_NONOBJECT_CLASS(Point2D, mitkPoint.h, mitk) SWIG_ADD_NONOBJECT_CLASS(Point3D, mitkPoint.h, mitk) SWIG_ADD_NONOBJECT_CLASS(Point4D, mitkPoint.h, mitk) SWIG_ADD_NONOBJECT_CLASS(Point2I, mitkPoint.h, mitk) SWIG_ADD_NONOBJECT_CLASS(Point3I, mitkPoint.h, mitk) SWIG_ADD_NONOBJECT_CLASS(Point4I, mitkPoint.h, mitk) SWIG_ADD_NONOBJECT_CLASS(VnlVector, mitkVector.h, mitk) -SWIG_ADD_NONOBJECT_CLASS(Vector2D, mitkVector.h, mitk) -SWIG_ADD_NONOBJECT_CLASS(Vector3D, mitkVector.h, mitk) -SWIG_ADD_NONOBJECT_CLASS(Vector4D, mitkVector.h, mitk) +%ignore rBegin() const; +%ignore rEnd() const; +%ignore rBegin(); +%ignore rEnd(); +%ignore itk::Vector::GetVnlVector; +%ignore itk::Vector::Get_vnl_vector; +%ignore itk::Vector::Get_vnl_vector const; + +%{ + #include + #include + #include + typedef vnl_vector_ref VnlVectorRef; +%} +//%include +//MITKSWIG_ADD_CLASS(VnlVectorRef, mitkVector.h,) +//typedef ::VnlVectorRef VnlVectorRef; +class VnlVectorRef; +%template(VectorVnlVectorRefPointer) std::vector<::VnlVectorRef * >; +//SWIG_ADD_NONOBJECT_NOVECTOR_CLASS(VnlVectorRef, mitkVector.h, ) + +%include +%include +class itk::FixedArray::ReverseIterator; +//%template(VnlVectorRef) vnl_vector_ref; +%template(itkFixedArray3D) itk::FixedArray; +%template(itkVector2D) itk::Vector; +%template(itkVector3D) itk::Vector; +%template(itkVector4D) itk::Vector; + + +SWIG_ADD_NONOBJECT_TEMPLATE_CLASS(Vector2D, mitkVector.h, mitk) +%template(Vector2D) mitk::Vector; +SWIG_ADD_NONOBJECT_TEMPLATE_CLASS(Vector3D, mitkVector.h, mitk) +%template(Vector3D) mitk::Vector; +SWIG_ADD_NONOBJECT_TEMPLATE_CLASS(Vector4D, mitkVector.h, mitk) +%template(Vector4D) mitk::Vector; SWIG_ADD_MITK_CLASS(BaseData, mitkBaseData.h, mitk) SWIG_ADD_MITK_CLASS(SlicedData, mitkSlicedData.h, mitk) SWIG_ADD_MITK_CLASS(Image, mitkImage.h, mitk) SWIG_ADD_MITK_CLASS(LabelSetImage, mitkLabelSetImage.h, mitk) SWIG_ADD_MITK_CLASS(PointSet, mitkPointSet.h, mitk) %{ using mitk::Message; %} // // Phenotyping Related Classes // SWIG_ADD_MITK_CLASS(AbstractGlobalImageFeature, mitkAbstractGlobalImageFeature.h, mitk) SWIG_ADD_MITK_CLASS(GIFImageDescriptionFeatures, mitkGIFImageDescriptionFeatures.h, mitk) SWIG_ADD_MITK_CLASS(GIFFirstOrderStatistics, mitkGIFFirstOrderStatistics.h, mitk) SWIG_ADD_MITK_CLASS(GIFFirstOrderHistogramStatistics, mitkGIFFirstOrderHistogramStatistics.h, mitk) SWIG_ADD_MITK_CLASS(GIFFirstOrderNumericStatistics, mitkGIFFirstOrderNumericStatistics.h, mitk) SWIG_ADD_MITK_CLASS(GIFVolumetricStatistics, mitkGIFVolumetricStatistics.h, mitk) SWIG_ADD_MITK_CLASS(GIFVolumetricDensityStatistics, mitkGIFVolumetricDensityStatistics.h, mitk) SWIG_ADD_MITK_CLASS(GIFCooccurenceMatrix2, mitkGIFCooccurenceMatrix2.h, mitk) SWIG_ADD_MITK_CLASS(GIFNeighbouringGreyLevelDependenceFeature, mitkGIFNeighbouringGreyLevelDependenceFeatures.h, mitk) SWIG_ADD_MITK_CLASS(GIFGreyLevelRunLength, mitkGIFGreyLevelRunLength.h, mitk) SWIG_ADD_MITK_CLASS(GIFGreyLevelSizeZone, mitkGIFGreyLevelSizeZone.h, mitk) SWIG_ADD_MITK_CLASS(GIFGreyLevelDistanceZone, mitkGIFGreyLevelDistanceZone.h, mitk) SWIG_ADD_MITK_CLASS(GIFLocalIntensity, mitkGIFLocalIntensity.h, mitk) SWIG_ADD_MITK_CLASS(GIFIntensityVolumeHistogramFeatures, mitkGIFIntensityVolumeHistogramFeatures.h, mitk) SWIG_ADD_MITK_CLASS(GIFNeighbourhoodGreyToneDifferenceFeatures, mitkGIFNeighbourhoodGreyToneDifferenceFeatures.h, mitk) SWIG_ADD_MITK_CLASS(GIFCurvatureStatistic, mitkGIFCurvatureStatistic.h, mitk) // // Conversion and Segmentation based Classes // SWIG_ADD_MITK_CLASS(ContourModelSetToImageFilter, mitkContourModelSetToImageFilter.h, mitk) SWIG_ADD_NONOBJECT_NOVECTOR_CLASS(BooleanOperation, mitkBooleanOperation.h, mitk) SWIG_ADD_NONOBJECT_NOVECTOR_CLASS(MorphologicalOperations, mitkMorphologicalOperations.h, mitk) +SWIG_ADD_NONOBJECT_NOVECTOR_CLASS(MaskCleaningOperation, mitkMaskCleaningOperation.h, mitk) +SWIG_ADD_NONOBJECT_NOVECTOR_CLASS(TransformationOperation, mitkTransformationOperation.h, mitk) +SWIG_ADD_NONOBJECT_NOVECTOR_CLASS(ArithmeticOperation, mitkArithmeticOperation.h, mitk) %{ #include typedef itk::DataObject::DataObjectIdentifierType DataObjectIdentifierType; typedef itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType; %} // // MatchPoint Related Classes // %ignore DLLHandle::LibraryHandleType; %{ #include namespace core { typedef std::string String; } typedef map::deployment::DLLHandle::LibraryHandleType LibraryHandleType; typedef map::deployment::DLLHandle DLLHandle; %} namespace core { typedef std::string String; } %nodefaultctor LibraryHandleType; struct LibraryHandleType {}; %include SWIG_ADD_MITK_CLASS(UID,mapUID.h, map::algorithm) SWIG_ADD_MITK_CLASS(DLLInfo, mapDeploymentDLLInfo.h, map::deployment) //SWIG_ADD_MITK_CLASS_VECTORFREE(DLLHandle, mapDeploymentDLLHandle.h, map::deployment) SWIG_ADD_MITK_CLASS_VECTORFREE(DLLDirectoryBrowser, mapDeploymentDLLDirectoryBrowser.h, ::map::deployment) MITKSWIG_ADD_HEADERFILE(mapDeploymentDLLAccess.h) %{ DLLHandle const * ConvertDLLHandleSmartPointerToPointer(itk::SmartPointer p) { return p.GetPointer(); } %} DLLHandle const * ConvertDLLHandleSmartPointerToPointer(DLLHandle::Pointer p); MITKSWIG_ADD_CLASS(MITKAlgorithmHelper, mitkAlgorithmHelper.h, mitk) MITKSWIG_ADD_CLASS(RegistrationType, mitkImageMappingHelper.h, mitk::ImageMappingHelper) MITKSWIG_ADD_CLASS(MITKRegistrationType, mitkImageMappingHelper.h, mitk::ImageMappingHelper) // SWIG_ADD_MITK_CLASS(FastSymmetricForcesDemonsMultiResDefaultRegistrationAlgorithm, mitkFastSymmetricForcesDemonsMultiResDefaultRegistrationAlgorithm.h, mitk) // SWIG_ADD_MITK_CLASS(LevelSetMotionMultiResDefaultRegistrationAlgorithm, mitkLevelSetMotionMultiResDefaultRegistrationAlgorithm.h, mitk) // SWIG_ADD_MITK_CLASS(MultiModalAffineDefaultRegistrationAlgorithm, mitkMultiModalAffineDefaultRegistrationAlgorithm.h, mitk) // SWIG_ADD_MITK_CLASS(MultiModalRigidDefaultRegistrationAlgorithm, mitkMultiModalRigidDefaultRegistrationAlgorithm.h, mitk) // SWIG_ADD_MITK_CLASS(MultiModalTransDefaultRegistrationAlgorithm, mitkMultiModalTransDefaultRegistrationAlgorithm.h, mitk) // SWIG_ADD_MITK_CLASS(RigidClosedFormPointsDefaultRegistrationAlgorithm, mitkRigidClosedFormPointsDefaultRegistrationAlgorithm.h, mitk) // SWIG_ADD_MITK_CLASS(RigidICPDefaultRegistrationAlgorithm, mitkRigidICPDefaultRegistrationAlgorithm.h, mitk) \ No newline at end of file diff --git a/Wrapping/Common/mitk_swig_macros.i b/Wrapping/Common/mitk_swig_macros.i index 76d3f761e6..8a9938ac96 100644 --- a/Wrapping/Common/mitk_swig_macros.i +++ b/Wrapping/Common/mitk_swig_macros.i @@ -1,239 +1,257 @@ // // This file contains macros for swig. // // // MITKSWIG_ADD_HEADERFILE includes a header-file into SWIG // %define MITKSWIG_ADD_HEADERFILE( classinclude ) // Include the given header, where the class definition is found %include < ## classinclude ## > // Include the include file in the generated cpp file %{ #include < ## classinclude ## > %} %enddef // // MITKSWIG_ADD_CLASS is a helper macro in order to do // all important stuff in order to wrap an existing // class // %define MITKSWIG_ADD_CLASS(classname, classinclude, nspace) MITKSWIG_ADD_HEADERFILE( classinclude ) // Using class name in order to remove ambigiouties %{ typedef nspace ## :: ## classname classname ## ; using nspace ## :: ## classname ; %} using nspace ##:: ## classname ; // Typedef is necessary to overcome ambigiouties resulting in the fact that SWIG // ignores namespaces. This can lead to some problems with templates. typedef nspace ## :: ## classname classname ## ; %enddef // // MITKSWIG_AUTOMATED_CASTING is a helper macro in order to // provide a convinience interface for up/downcasting of // classes // %define MITKSWIG_AUTOMATED_CASTING(classname, classinclude, nspace) // Define a conversion method to convert from and to types %template(ConvertTo ## classname) ConvertTo< nspace ##:: ## classname ## >; // This extend is necessary to have the automatic cast methods available %extend itk::SmartPointer< nspace ## :: ## classname ## ::Self> { %pythoncode %{ def _GetListOfValidItems(self): return [str(k).replace("class itk::","") for k in self.GetClassHierarchy() if str(k).replace("class itk::","") in convertion_list.keys() ] %} %pythoncode %{ def __getattr__(self, item): if type(item)==str: if (len(item) > 9) and ('ConvertTo' in item): searchString=item[9:] if searchString in self._GetListOfValidItems(): def func_t(): return convertion_list[searchString](self) return func_t %} %pythoncode %{ def __dir__(self): return super().__dir__() + ['ConvertTo'+k for k in self._GetListOfValidItems()] %} } %extend std::vector< nspace ## :: ## classname *>::value_type { %pythoncode %{ def _GetListOfValidItems(self): return [str(k) for k in self.GetClassHierarchy() if str(k) in convertion_list.keys() ] %} %pythoncode %{ def __getattr__(self, item): if type(item)==str: if (len(item) > 9) and ('ConvertTo' in item): searchString=item[9:] if searchString in self._GetListOfValidItems(): def func_t(): return convertion_list[searchString](self) return func_t %} %pythoncode %{ def __dir__(self): return super().__dir__() + ['ConvertTo'+k for k in self._GetListOfValidItems()] %} } %pythoncode %{ convertion_list['classname'] = ConvertTo ## classname %} %enddef // // MITKSWIG_SMARTPOINTERVECTOR : Wrapper for Vectors of Smartpointer-Classes // %define MITKSWIG_SMARTPOINTERVECTOR(classname, classinclude, nspace) // Initianziation of std. vectors containing pointers to these classes. This allows to use // vectors of these types as target language arrays. %template(Vector ## classname ## Pointer) std::vector< nspace ## :: ## classname ## ::Pointer >; %template(Vector ## classname) std::vector< nspace ## :: ## classname ## ::Self *>; %enddef // // MITKSWIG_POINTERVECTOR : Wrapper for Vectors of Classes // %define MITKSWIG_POINTERVECTOR(classname, classinclude, nspace) // Initianziation of std. vectors containing pointers to these classes. This allows to use // vectors of these types as target language arrays. %template(Vector ## classname ## Pointer) std::vector< nspace ## :: ## classname * >; %template(Vector ## classname) std::vector< nspace ## :: ## classname >; %enddef // // MITKSWIG_MITKSMARTPOINTER_INITIALIZATION : Wrapper for Vectors of Smartpointer-Classes // %define MITKSWIG_MITKSMARTPOINTER_INITIALIZATION(classname, classinclude, nspace) // Declaring that this class is a smart-pointer class, in order to handle // online upcasting where necessary (for example python) %feature("smartptr", noblock=1) classname { classname ## ::Pointer } %feature("smartptr", noblock=1) nspace ##:: ## classname { nspace ## :: ## classname ## ::Pointer } %feature("smartptr", noblock=1) nspace ##:: ## classname { classname ## ::Pointer } %feature("smartptr", noblock=1) nspace ##:: ## classname { itk::SmartPointer } %enddef // // MITKSWIG_MITKSMARTPOINTER_TEMPLATE : Wrapper for Vectors of Smartpointer-Classes // %define MITKSWIG_MITKSMARTPOINTER_TEMPLATE(classname, classinclude, nspace) // Defining the Smartpointer, allows easy access in target language %template(classname ## Pointer) itk::SmartPointer; %enddef // // SWIG_ADD_MITK_CLASS is a helper macro in order to do // all important stuff before an mitk::Class is included. // Requires the name of the class as it is in c++ as classname // and the include file, in which the class is defined. // It is assumed that the class is somehow inherited from // mitk::BaseData, and supports smartpointers. // %define SWIG_ADD_MITK_CLASS_VECTORFREE(classname, classinclude, nspace) %include < ## classinclude ## > MITKSWIG_MITKSMARTPOINTER_INITIALIZATION(classname, classinclude, nspace) MITKSWIG_ADD_CLASS( classname, classinclude, nspace ) //class classname ## ; //class classname ## ::Pointer; class nspace ## :: ## classname ## ; class nspace ## :: ## classname ## ::Pointer; MITKSWIG_MITKSMARTPOINTER_TEMPLATE(classname, classinclude, nspace) MITKSWIG_AUTOMATED_CASTING(classname, classinclude, nspace) %enddef // // SWIG_ADD_MITK_CLASS is a helper macro in order to do // all important stuff before an mitk::Class is included. // Requires the name of the class as it is in c++ as classname // and the include file, in which the class is defined. // It is assumed that the class is somehow inherited from // mitk::BaseData, and supports smartpointers. // %define SWIG_ADD_MITK_CLASS(classname, classinclude, nspace) %include < ## classinclude ## > MITKSWIG_MITKSMARTPOINTER_INITIALIZATION(classname, classinclude, nspace) MITKSWIG_ADD_CLASS( classname, classinclude, nspace ) // Defining only the original classname // If the ::Pointer is also defined, smartpointer won't work class nspace ## :: ## classname ## ; //class nspace ## :: ## classname ## ::Pointer; // It is important to first define the Vectors and // then define the Smartpointer. Otherwise a SWIG-bug ... MITKSWIG_SMARTPOINTERVECTOR(classname, classinclude, nspace) MITKSWIG_MITKSMARTPOINTER_TEMPLATE(classname, classinclude, nspace) MITKSWIG_AUTOMATED_CASTING(classname, classinclude, nspace) %enddef // // SWIG_ADD_NONOBJECT_CLASS is a helper macro in order to do // all important stuff before an mitk::Class is included. // Requires the name of the class as it is in c++ as classname // and the include file, in which the class is defined. // It is assumed that the class is somehow inherited from // mitk::BaseData, and supports smartpointers. // %define SWIG_ADD_NONOBJECT_CLASS(classname, classinclude, nspace) %include < ## classinclude ## > MITKSWIG_ADD_CLASS( classname, classinclude, nspace ) // Typedef is necessary to overcome ambigiouties resulting in the fact that SWIG // ignores namespaces. This can lead to some problems with templates. typedef nspace ## :: ## classname classname ## ; class nspace ## :: ## classname ## ; MITKSWIG_POINTERVECTOR(classname, classinclude, nspace) %enddef +// +// SWIG_ADD_NONOBJECT_CLASS is a helper macro in order to do +// all important stuff before an mitk::Class is included. +// Requires the name of the class as it is in c++ as classname +// and the include file, in which the class is defined. +// It is assumed that the class is somehow inherited from +// mitk::BaseData, and supports smartpointers. +// +%define SWIG_ADD_NONOBJECT_TEMPLATE_CLASS(classname, classinclude, nspace) + %include < ## classinclude ## > + MITKSWIG_ADD_CLASS( classname, classinclude, nspace ) + + //%template(classname) templateArgs; + //class nspace ## :: ## classname ## ; + + MITKSWIG_POINTERVECTOR(classname, classinclude, nspace) +%enddef + // // SWIG_ADD_NONOBJECT_CLASS is a helper macro in order to do // all important stuff before an mitk::Class is included. // Requires the name of the class as it is in c++ as classname // and the include file, in which the class is defined. // It is assumed that the class is somehow inherited from // mitk::BaseData, and supports smartpointers. // %define SWIG_ADD_NONOBJECT_NOVECTOR_CLASS(classname, classinclude, nspace) %include < ## classinclude ## > MITKSWIG_ADD_CLASS( classname, classinclude, nspace ) // Typedef is necessary to overcome ambigiouties resulting in the fact that SWIG // ignores namespaces. This can lead to some problems with templates. typedef nspace ## :: ## classname classname ## ; // Initianziation of std. vectors containing pointers to these classes. This allows to use // vectors of these types as target language arrays. %template(Vector ## classname ## Pointer) std::vector< nspace ## :: ## classname * >; %enddef diff --git a/Wrapping/Python/CMakeLists.txt b/Wrapping/Python/CMakeLists.txt index cb9327dc76..724d6cdca1 100644 --- a/Wrapping/Python/CMakeLists.txt +++ b/Wrapping/Python/CMakeLists.txt @@ -1,77 +1,77 @@ # Version 2.8.1 is the minium requirement for this script. # this is lower than the general minimum requirement. #cmake_minimum_required ( VERSION 2.8.1 FATAL_ERROR ) include(mitkTargetLinkLibrariesWithDynamicLookup) project( MITK_Python ) set(CMAKE_SHARED_LINKER_FLAGS "" CACHE INTERNAL "" FORCE) set(CMAKE_MODULE_LINKER_FLAGS "" CACHE INTERNAL "" FORCE) mitk_check_dynamic_lookup(MODULE SHARED MITK_UNDEFINED_SYMBOLS_ALLOWED ) # # Find the necessary libraries etc.. # if ( MITK_UNDEFINED_SYMBOLS_ALLOWED ) set( _QUIET_LIBRARY "QUIET" ) else() set( _QUIET_LIBRARY "REQUIRED" ) endif() find_package ( PythonInterp REQUIRED ) find_package ( PythonLibs ${_QUIET_LIBRARY} ) include_directories ( ${CMAKE_CURRENT_SOURCE_DIR} ) # # Options # option ( MITK_PYTHON_THREADS "Enable threaded python usage by unlocking the GIL." ON ) mark_as_advanced( MITK_PYTHON_THREADS ) option ( MITK_PYTHON_EGG "Add building of python eggs to the dist target." OFF ) mark_as_advanced( MITK_PYTHON_EGG ) option ( MITK_PYTHON_WHEEL "Add building of python wheels to the dist target." ON ) mark_as_advanced( MITK_PYTHON_WHEEL ) # Prepare the SWIG-File, i.e. especially add necessary include folders -mitkSwigPrepareFiles(pyMITK MITK.i "MitkCore;MitkCLCore;MitkCLUtilities;ITKCommon;MitkMatchPointRegistration;MitkSegmentation;MitkMultilabel;MitkDICOMReader;MitkDICOMReaderServices;MitkDicomRT") +mitkSwigPrepareFiles(pyMITK MITK.i "MitkCore;MitkCLCore;MitkCLUtilities;MitkBasicImageProcessing;ITKCommon;MitkMatchPointRegistration;MitkSegmentation;MitkMultilabel;MitkDICOMReader;MitkDICOMReaderServices;MitkDicomRT") # Add additional SWIG Parameters # These parameters depend on the target language set(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_GLOBAL_FLAGS} -features autodoc=1 -keyword ) if( MITK_PYTHON_THREADS ) set(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_FLAGS} -threads) endif() set(CMAKE_SWIG_OUTDIR ${CMAKE_CURRENT_BINARY_DIR}) # Create the actual SWIG project swig_add_module(pyMITK python MITK.i ) -mitkSwigAddLibraryDependencies(pyMITK "MitkCore;MitkCLCore;MitkCLUtilities;ITKCommon;MitkMatchPointRegistration;MitkSegmentation;MitkMultilabel;MitkDICOMReader;MitkDICOMReaderServices;MitkDicomRT") +mitkSwigAddLibraryDependencies(pyMITK "MitkCore;MitkCLCore;MitkCLUtilities;MitkBasicImageProcessing;ITKCommon;MitkMatchPointRegistration;MitkSegmentation;MitkMultilabel;MitkDICOMReader;MitkDICOMReaderServices;MitkDicomRT") mitk_target_link_libraries_with_dynamic_lookup(${SWIG_MODULE_pyMITK_REAL_NAME} ${PYTHON_LIBRARIES}) if(DEFINED SKBUILD) message(WARNING "SKBuild exists") # Currently this installation install(FILES ${CMAKE_CURRENT_BINARY_DIR}/pyMITK.py ${CMAKE_CURRENT_SOURCE_DIR}/Packaging/__init__.py # ${MITK_DOC_FILES} DESTINATION pyMITK COMPONENT Runtime ) install(TARGETS ${SWIG_MODULE_pyMITK_REAL_NAME} RUNTIME DESTINATION pyMITK LIBRARY DESTINATION pyMITK COMPONENT Runtime ) else() message(WARNING "SKBuild missing") include(LegacyPackaging.cmake) endif()