diff --git a/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp b/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp index 0f2a6dd54f..bbd98ebebb 100644 --- a/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp +++ b/Modules/BasicImageProcessing/src/mitkTransformationOperation.cpp @@ -1,756 +1,758 @@ /*=================================================================== 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) + int signedVImageDimension = static_cast(VImageDimension); + + for (int i = 0; i < signedVImageDimension; ++i) { - for (int j = 0; j < VImageDimension; ++j) + for (int j = 0; j < signedVImageDimension; ++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; + //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]; + //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] - 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/CLUtilities/include/mitkGIFGreyLevelDistanceZone.h b/Modules/Classification/CLUtilities/include/mitkGIFGreyLevelDistanceZone.h index 63b6a7b64c..f3ad178922 100644 --- a/Modules/Classification/CLUtilities/include/mitkGIFGreyLevelDistanceZone.h +++ b/Modules/Classification/CLUtilities/include/mitkGIFGreyLevelDistanceZone.h @@ -1,183 +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; + bool m_SmallNeighbourhoodForDistance=true; }; } #endif //mitkGIFGreyLevelDistanceZone_h diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp index 4fc7fabe14..79fa6c14bd 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp @@ -1,492 +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 1) && (static_cast(direction) < VDimension + 2)) { radius[direction - 2] = 0; } itk::NeighborhoodIterator 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) + if ( useSmallNeighbourhood) { auto cIndex = distanceIter.GetIndex(i); int differentIndexes = 0; - for (int j = 0; j < VDimension; ++j) + for (int j = 0; static_cast(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, bool useSmallNeighbourhood, int direction, int &maxDistance) { AccessByItk_n(input, itkErode2, (output, useSmallNeighbourhood, direction, maxDistance)); } void erodeAndAdd(mitk::Image::Pointer input, mitk::Image::Pointer& finalOutput, bool useSmallNeighbourhood, int direction, int &maxDistance) { maxDistance = 0; 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, 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() : 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()==NULL) { config.distanceMask = mask->Clone(); } else { 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 42b3e49707..eec786aaaa 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyToneDifferenceFeatures.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyToneDifferenceFeatures.cpp @@ -1,228 +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 1) && (static_cast(params.direction) < VImageDimension + 2)) { regionSize[params.direction - 2] = 0; } itk::NeighborhoodIterator 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 b02b1f1417..e7c52db715 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricDensityStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricDensityStatistics.cpp @@ -1,620 +1,620 @@ /*=================================================================== 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; }; static double LegendrePn(unsigned int n, double x) { if (n == 0) { return 1.0; } else if (n == 1) { return x; } else if (n == 2) { return ((3.0 * x*x) - 1.0) * 0.5; } if (x == 1.0) { return 1.0; } if (x == -1.0) { return ((n % 2 == 0) ? 1.0 : -1.0); } if ((x == 0.0) && (n % 2)) { return 0.0; } /* We could simply do this: return (double(((2 * n) - 1)) * x * Pn(n - 1, x) - (double(n - 1)) * Pn(n - 2, x)) / (double)n ; but it could be slow for large n */ double pnm1(LegendrePn(2,x)); double pnm2(LegendrePn(1,x)); double pn(pnm1); for (unsigned int l = 3; l <= n; l++) { pn = (((2.0 * (double)l) - 1.0) * x * pnm1 - (((double)l - 1.0) * pnm2)) / (double)l; pnm2 = pnm1; pnm1 = pn; } return pn; } template void 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; 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 maskA(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; + std::vector< typename ImageType::PointType > vectorPointsInMask; + std::vector< TPixel > vectorValuesInMask; vectorPointsInMask.reserve(Nv); vectorValuesInMask.reserve(Nv); while (!imgA.IsAtEnd()) { if (maskA.Get() > 0) { auto maskAIndex = maskA.GetIndex(); itkImage->TransformIndexToPhysicalPoint(maskAIndex, pointA); 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.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(); if (original_mask.IsNull()) { original_mask = mask->Clone(); } std::string prefix = FeatureDescriptionPrefix(); vtkSmartPointer mesher = vtkSmartPointer::New(); vtkSmartPointer stats = vtkSmartPointer::New(); vtkSmartPointer stats2 = vtkSmartPointer::New(); 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=0; double surf_mvee=0; calculateMEE(mesher->GetOutput(), vol_mvee, surf_mvee); - double vol_mobb=0; - double surf_mobb=0; + //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_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 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, original_mask->GetChannelDescriptor().GetPixelType(), original_mask, original_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 (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); } if (pxMask > 0) { points->InsertNextPoint(x*spacingX, y*spacingY, z*spacingZ); if (pxImage == pxImage) { 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) * 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) { double legendreValue = LegendrePn(i, (alpha*alpha + beta * beta) / (2 * alpha*beta)); ad_aee += 4 * pi*major*minor* legendreValue* (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 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/Wrapping/CMakeLists.txt b/Wrapping/CMakeLists.txt index c757046536..022aa5d433 100644 --- a/Wrapping/CMakeLists.txt +++ b/Wrapping/CMakeLists.txt @@ -1,70 +1,70 @@ find_package(SWIG QUIET) if(SWIG_FOUND) include(mitkLanguageOptions) include(UseSWIG) include(mitkSwigAddLibraryDependencies) include(mitkSwigPrepareFiles) # Path to common files set(MITK_WRAPPING_COMMON_DIR ${MITK_SOURCE_DIR}/Wrapping/Common) # make a manual list of dependencies for the Swig.i files - list( APPEND SWIG_EXTRA_DEPS - "${MITK_WRAPPING_COMMON_DIR}/MITK_Common.i" - ) +# list( APPEND SWIG_EXTRA_DEPS +# "${MITK_WRAPPING_COMMON_DIR}/MITK_Common.i" +# ) # A general packaging target, not built by default, to build packages for each # language. This should depend on all language specific targets. add_custom_target( dist ${CMAKE_COMMAND} -E echo "Finished generating wrapped packages for distribution..." ) # # lua SWIG configuration # #if ( WRAP_LUA ) # add_subdirectory ( Lua ) #endif() # # python SWIG configuration # if ( WRAP_PYTHON ) add_subdirectory ( Python ) endif() # # ruby SWIG configuration # #if ( WRAP_RUBY ) # add_subdirectory ( Ruby ) #endif() # # JAVA SWIG configuration # #if ( WRAP_JAVA ) # add_subdirectory( Java ) #endif() # # C# SWIG configuration # #if ( WRAP_CSHARP ) # add_subdirectory ( CSharp ) #endif() # # TCL SWIG configuration # #if ( WRAP_TCL ) # add_subdirectory ( Tcl ) #endif() # # R SWIG configuration # #if ( WRAP_R ) # add_subdirectory( R ) #endif() endif()