diff --git a/Modules/Classification/CLCore/CMakeLists.txt b/Modules/Classification/CLCore/CMakeLists.txt index 045293844a..7726503bbe 100644 --- a/Modules/Classification/CLCore/CMakeLists.txt +++ b/Modules/Classification/CLCore/CMakeLists.txt @@ -1,8 +1,8 @@ MITK_CREATE_MODULE( - DEPENDS MitkCore + DEPENDS MitkCore MitkCommandLine PACKAGE_DEPENDS PUBLIC Eigen WARNINGS_AS_ERRORS ) #add_subdirectory(test) diff --git a/Modules/Classification/CLCore/include/mitkAbstractGlobalImageFeature.h b/Modules/Classification/CLCore/include/mitkAbstractGlobalImageFeature.h index 0320ba5abf..99a9ff726e 100644 --- a/Modules/Classification/CLCore/include/mitkAbstractGlobalImageFeature.h +++ b/Modules/Classification/CLCore/include/mitkAbstractGlobalImageFeature.h @@ -1,78 +1,115 @@ /*=================================================================== 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 + + // Eigen #include // STD Includes // MITK includes #include namespace mitk { class MITKCLCORE_EXPORT AbstractGlobalImageFeature : public BaseData { public: typedef std::vector< std::pair > FeatureListType; - typedef std::vector< std::string> FeatureNameListType; + typedef std::vector< std::string> FeatureNameListType; + typedef std::map ParameterTypes; /** - * \brief Calculates the feature of this abstact interface. + * \brief Calculates the feature of this abstact interface. Does not necessarily considers the parameter settings. */ - virtual FeatureListType CalculateFeatures(const Image::Pointer & image, const Image::Pointer &feature) = 0; + virtual FeatureListType CalculateFeatures(const Image::Pointer & feature, const Image::Pointer &mask) = 0; + + /** + * \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; + itkSetMacro(Prefix, std::string); + itkSetMacro(ShortName, std::string); + itkSetMacro(LongName, std::string); + ParameterTypes GetParameter() { return m_Parameter; }; + + itkGetConstMacro(Prefix, std::string); + itkGetConstMacro(ShortName, std::string); + itkGetConstMacro(LongName, std::string); + itkGetConstMacro(Parameter, ParameterTypes); + + 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); + public: //#ifndef DOXYGEN_SKIP virtual void SetRequestedRegionToLargestPossibleRegion() override {}; virtual bool RequestedRegionIsOutsideOfTheBufferedRegion() override { return true; }; virtual bool VerifyRequestedRegion() override { return false; }; virtual void SetRequestedRegion (const itk::DataObject * /*data*/) override {}; // Override virtual bool IsEmpty() const override { if(IsInitialized() == false) return true; const TimeGeometry* timeGeometry = const_cast(this)->GetUpdatedTimeGeometry(); if(timeGeometry == NULL) 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) + ParameterTypes m_Parameter; // Parameter setting + //#endif // Skip Doxygen }; } #endif //mitkAbstractGlobalImageFeature_h diff --git a/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp b/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp index cf41725808..5de783d940 100644 --- a/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp +++ b/Modules/Classification/CLCore/src/mitkAbstractGlobalImageFeature.cpp @@ -1,17 +1,31 @@ /*=================================================================== 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 + +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; +} \ No newline at end of file diff --git a/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp b/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp index fb50e373f9..b18301afca 100644 --- a/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp +++ b/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp @@ -1,375 +1,396 @@ /*=================================================================== 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 mitkCLPolyToNrrd_cpp #define mitkCLPolyToNrrd_cpp #include "time.h" #include #include #include #include "mitkCommandLineParser.h" #include #include #include #include #include +#include #include #include #include #include #include #include "itkNearestNeighborInterpolateImageFunction.h" #include "itkResampleImageFilter.h" typedef itk::Image< double, 3 > FloatImageType; typedef itk::Image< unsigned char, 3 > MaskImageType; static std::vector 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; } template void ResampleImage(itk::Image* itkImage, float resolution, mitk::Image::Pointer& newImage) { typedef itk::Image ImageType; typedef itk::ResampleImageFilter ResampleFilterType; typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); auto spacing = itkImage->GetSpacing(); auto size = itkImage->GetLargestPossibleRegion().GetSize(); for (int i = 0; i < VImageDimension; ++i) { size[i] = size[i] / (1.0*resolution)*(1.0*spacing[i])+1.0; } spacing.Fill(resolution); resampler->SetInput(itkImage); resampler->SetSize(size); resampler->SetOutputSpacing(spacing); resampler->SetOutputOrigin(itkImage->GetOrigin()); resampler->SetOutputDirection(itkImage->GetDirection()); resampler->Update(); newImage->InitializeByItk(resampler->GetOutput()); mitk::GrabItkImageMemory(resampler->GetOutput(), newImage); } static void CreateNoNaNMask(mitk::Image::Pointer mask, mitk::Image::Pointer ref, mitk::Image::Pointer& newMask) { MaskImageType::Pointer itkMask = MaskImageType::New(); FloatImageType::Pointer itkValue = FloatImageType::New(); mitk::CastToItkImage(mask, itkMask); mitk::CastToItkImage(ref, itkValue); typedef itk::ImageDuplicator< MaskImageType > DuplicatorType; DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage(itkMask); duplicator->Update(); auto tmpMask = duplicator->GetOutput(); itk::ImageRegionIterator mask1Iter(itkMask, itkMask->GetLargestPossibleRegion()); itk::ImageRegionIterator mask2Iter(tmpMask, tmpMask->GetLargestPossibleRegion()); itk::ImageRegionIterator imageIter(itkValue, itkValue->GetLargestPossibleRegion()); while (!mask1Iter.IsAtEnd()) { mask2Iter.Set(0); if (mask1Iter.Value() > 0) { // Is not NaN if (imageIter.Value() == imageIter.Value()) { mask2Iter.Set(1); } } ++mask1Iter; ++mask2Iter; ++imageIter; } newMask->InitializeByItk(tmpMask); mitk::GrabItkImageMemory(tmpMask, newMask); } static void ResampleMask(mitk::Image::Pointer mask, mitk::Image::Pointer ref, mitk::Image::Pointer& newMask) { typedef itk::NearestNeighborInterpolateImageFunction< MaskImageType> NearestNeighborInterpolateImageFunctionType; typedef itk::ResampleImageFilter ResampleFilterType; NearestNeighborInterpolateImageFunctionType::Pointer nn_interpolator = NearestNeighborInterpolateImageFunctionType::New(); MaskImageType::Pointer itkMoving = MaskImageType::New(); MaskImageType::Pointer itkRef = MaskImageType::New(); mitk::CastToItkImage(mask, itkMoving); mitk::CastToItkImage(ref, itkRef); ResampleFilterType::Pointer resampler = ResampleFilterType::New(); resampler->SetInput(itkMoving); resampler->SetReferenceImage(itkRef); resampler->UseReferenceImageOn(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); newMask->InitializeByItk(resampler->GetOutput()); mitk::GrabItkImageMemory(resampler->GetOutput(), newMask); } int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); // required params parser.addArgument("image", "i", mitkCommandLineParser::InputImage, "Input Image", "Path to the input VTK polydata", us::Any(), false); parser.addArgument("mask", "m", mitkCommandLineParser::InputImage, "Input Mask", "Mask Image that specifies the area over for the statistic, (Values = 1)", us::Any(), false); parser.addArgument("output", "o", mitkCommandLineParser::OutputFile, "Output text file", "Target file. The output statistic is appended to this file.", us::Any(), false); parser.addArgument("cooccurence", "cooc", mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features", us::Any()); parser.addArgument("cooccurence2", "cooc2", mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features (new implementation)", us::Any()); + parser.addArgument("ngldm", "ngldm", mitkCommandLineParser::String, "Calculate Neighbouring grey level dependence based features", "Calculate Neighbouring grey level dependence based features", us::Any()); parser.addArgument("volume","vol",mitkCommandLineParser::String, "Use Volume-Statistic", "calculates volume based features",us::Any()); parser.addArgument("run-length","rl",mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features",us::Any()); parser.addArgument("first-order","fo",mitkCommandLineParser::String, "Use First Order Features", "calculates First order based features",us::Any()); parser.addArgument("header","head",mitkCommandLineParser::String,"Add Header (Labels) to output","",us::Any()); parser.addArgument("description","d",mitkCommandLineParser::String,"Text","Description that is added to the output",us::Any()); parser.addArgument("same-space", "sp", mitkCommandLineParser::String, "Bool", "Set the spacing of all images to equal. Otherwise an error will be thrown. ", us::Any()); parser.addArgument("resample-mask", "rm", mitkCommandLineParser::Bool, "Bool", "Resamples the mask to the resolution of the input image ", us::Any()); parser.addArgument("save-resample-mask", "srm", mitkCommandLineParser::String, "String", "If specified the resampled mask is saved to this path (if -rm is 1)", us::Any()); parser.addArgument("fixed-isotropic", "fi", mitkCommandLineParser::Float, "Float", "Input image resampled to fixed isotropic resolution given in mm. Should be used with resample-mask ", us::Any()); parser.addArgument("direction", "dir", mitkCommandLineParser::String, "Int", "Allows to specify the direction for Cooc and RL. 0: All directions, 1: Only single direction (Test purpose), 2,3,4... Without dimension 0,1,2... ", us::Any()); // Miniapp Infos parser.setCategory("Classification Tools"); parser.setTitle("Global Image Feature calculator"); parser.setDescription("Calculates different global statistics for a given segmentation / image combination"); parser.setContributor("MBI"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) { return EXIT_FAILURE; } if ( parsedArgs.count("help") || parsedArgs.count("h")) { return EXIT_SUCCESS; } MITK_INFO << "Version: "<< 1.7; //bool useCooc = parsedArgs.count("cooccurence"); bool resampleMask = false; if (parsedArgs.count("resample-mask")) { resampleMask = us::any_cast(parsedArgs["resample-mask"]); } mitk::Image::Pointer image = mitk::IOUtil::LoadImage(parsedArgs["image"].ToString()); mitk::Image::Pointer mask = mitk::IOUtil::LoadImage(parsedArgs["mask"].ToString()); if (parsedArgs.count("fixed-isotropic")) { mitk::Image::Pointer newImage = mitk::Image::New(); float resolution = us::any_cast(parsedArgs["fixed-isotropic"]); AccessByItk_2(image, ResampleImage, resolution, newImage); image = newImage; } if (resampleMask) { mitk::Image::Pointer newMaskImage = mitk::Image::New(); ResampleMask(mask, image, newMaskImage); mask = newMaskImage; if (parsedArgs.count("save-resample-mask")) { mitk::IOUtil::SaveImage(mask, parsedArgs["save-resample-mask"].ToString()); } } bool fixDifferentSpaces = parsedArgs.count("same-space"); if ( ! mitk::Equal(mask->GetGeometry(0)->GetOrigin(), image->GetGeometry(0)->GetOrigin())) { MITK_INFO << "Not equal Origins"; if (fixDifferentSpaces) { image->GetGeometry(0)->SetOrigin(mask->GetGeometry(0)->GetOrigin()); } else { return -1; } } if ( ! mitk::Equal(mask->GetGeometry(0)->GetSpacing(), image->GetGeometry(0)->GetSpacing())) { MITK_INFO << "Not equal Sapcings"; if (fixDifferentSpaces) { image->GetGeometry(0)->SetSpacing(mask->GetGeometry(0)->GetSpacing()); } else { return -1; } } int direction = 0; if (parsedArgs.count("direction")) { direction = splitDouble(parsedArgs["direction"].ToString(), ';')[0]; } mitk::Image::Pointer maskNoNaN = mitk::Image::New(); CreateNoNaNMask(mask, image, maskNoNaN); mitk::AbstractGlobalImageFeature::FeatureListType stats; //////////////////////////////////////////////////////////////// // Calculate First Order Features //////////////////////////////////////////////////////////////// if (parsedArgs.count("first-order")) { MITK_INFO << "Start calculating first order statistics...."; mitk::GIFFirstOrderStatistics::Pointer firstOrderCalculator = mitk::GIFFirstOrderStatistics::New(); auto localResults = firstOrderCalculator->CalculateFeatures(image, maskNoNaN); stats.insert(stats.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating first order statistics...."; } //////////////////////////////////////////////////////////////// // Calculate Volume based Features //////////////////////////////////////////////////////////////// if (parsedArgs.count("volume")) { MITK_INFO << "Start calculating volumetric ...."; mitk::GIFVolumetricStatistics::Pointer volCalculator = mitk::GIFVolumetricStatistics::New(); auto localResults = volCalculator->CalculateFeatures(image, mask); stats.insert(stats.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating volumetric...."; } //////////////////////////////////////////////////////////////// // Calculate Co-occurence Features //////////////////////////////////////////////////////////////// if (parsedArgs.count("cooccurence")) { auto ranges = splitDouble(parsedArgs["cooccurence"].ToString(),';'); for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; mitk::GIFCooccurenceMatrix::Pointer coocCalculator = mitk::GIFCooccurenceMatrix::New(); coocCalculator->SetRange(ranges[i]); coocCalculator->SetDirection(direction); auto localResults = coocCalculator->CalculateFeatures(image, maskNoNaN); stats.insert(stats.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; } } + //////////////////////////////////////////////////////////////// + // Calculate ngldm Features + //////////////////////////////////////////////////////////////// + if (parsedArgs.count("ngldm")) + { + auto ranges = splitDouble(parsedArgs["ngldm"].ToString(), ';'); + + for (std::size_t i = 0; i < ranges.size(); ++i) + { + MITK_INFO << "Start calculating NGLDM with range " << ranges[i] << "...."; + mitk::GIFNeighbouringGreyLevelDependenceFeature::Pointer coocCalculator = mitk::GIFNeighbouringGreyLevelDependenceFeature::New(); + coocCalculator->SetRange(ranges[i]); + coocCalculator->SetDirection(direction); + auto localResults = coocCalculator->CalculateFeatures(image, maskNoNaN); + stats.insert(stats.end(), localResults.begin(), localResults.end()); + MITK_INFO << "Finished calculating NGLDM with range " << ranges[i] << "...."; + } + } + //////////////////////////////////////////////////////////////// // Calculate Co-occurence Features //////////////////////////////////////////////////////////////// if (parsedArgs.count("cooccurence2")) { auto ranges = splitDouble(parsedArgs["cooccurence2"].ToString(), ';'); for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; mitk::GIFCooccurenceMatrix2::Pointer coocCalculator = mitk::GIFCooccurenceMatrix2::New(); coocCalculator->SetRange(ranges[i]); coocCalculator->SetDirection(direction); auto localResults = coocCalculator->CalculateFeatures(image, maskNoNaN); stats.insert(stats.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; } } //////////////////////////////////////////////////////////////// // Calculate Run-Length Features //////////////////////////////////////////////////////////////// if (parsedArgs.count("run-length")) { auto ranges = splitDouble(parsedArgs["run-length"].ToString(),';'); for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating run-length with number of bins " << ranges[i] << "...."; mitk::GIFGrayLevelRunLength::Pointer calculator = mitk::GIFGrayLevelRunLength::New(); calculator->SetRange(ranges[i]); calculator->SetDirection(direction); auto localResults = calculator->CalculateFeatures(image, mask); stats.insert(stats.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating run-length with number of bins " << ranges[i] << "...."; } } for (std::size_t i = 0; i < stats.size(); ++i) { std::cout << stats[i].first << " - " << stats[i].second < #include #include namespace mitk { class MITKCLUTILITIES_EXPORT GIFCooccurenceMatrix : public AbstractGlobalImageFeature { public: mitkClassMacro(GIFCooccurenceMatrix,AbstractGlobalImageFeature) itkFactorylessNewMacro(Self) itkCloneMacro(Self) GIFCooccurenceMatrix(); /** * \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; itkGetConstMacro(Range,double); itkSetMacro(Range, double); itkGetConstMacro(Direction, unsigned int); itkSetMacro(Direction, unsigned int); + virtual void CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList); + virtual void AddArguments(mitkCommandLineParser &parser); + + struct GIFCooccurenceMatrixConfiguration { double range; unsigned int direction; }; private: double m_Range; unsigned int m_Direction; }; } #endif //mitkGIFCooccurenceMatrix_h diff --git a/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h b/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h index 5a228faa6a..62b6a491bb 100644 --- a/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h +++ b/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h @@ -1,138 +1,141 @@ #ifndef mitkGIFCooccurenceMatrix2_h #define mitkGIFCooccurenceMatrix2_h #include #include #include #include namespace mitk { struct CoocurenceMatrixHolder { public: CoocurenceMatrixHolder(double min, double max, int number); int IntensityToIndex(double intensity); double IndexToMinIntensity(int index); double IndexToMeanIntensity(int index); double IndexToMaxIntensity(int index); double m_MinimumRange; double m_MaximumRange; double m_Stepsize; int m_NumberOfBins; Eigen::MatrixXd m_Matrix; }; struct CoocurenceMatrixFeatures { CoocurenceMatrixFeatures(): JointMaximum(0), JointAverage(0), JointVariance(0), JointEntropy(0), RowMaximum(0), RowAverage(0), RowVariance(0), RowEntropy(0), FirstRowColumnEntropy(0), SecondRowColumnEntropy(0), DifferenceAverage(0), DifferenceVariance(0), DifferenceEntropy(0), SumAverage(0), SumVariance(0), SumEntropy(0), AngularSecondMoment(0), Contrast(0), Dissimilarity(0), InverseDifference(0), InverseDifferenceNormalised(0), InverseDifferenceMoment(0), InverseDifferenceMomentNormalised(0), InverseVariance(0), Correlation(0), Autocorrelation(0), ClusterTendency(0), ClusterShade(0), ClusterProminence(0), FirstMeasureOfInformationCorrelation(0), SecondMeasureOfInformationCorrelation(0) { } public: double JointMaximum; double JointAverage; double JointVariance; double JointEntropy; double RowMaximum; double RowAverage; double RowVariance; double RowEntropy; double FirstRowColumnEntropy; double SecondRowColumnEntropy; double DifferenceAverage; double DifferenceVariance; double DifferenceEntropy; double SumAverage; double SumVariance; double SumEntropy; double AngularSecondMoment; double Contrast; double Dissimilarity; double InverseDifference; double InverseDifferenceNormalised; double InverseDifferenceMoment; double InverseDifferenceMomentNormalised; double InverseVariance; double Correlation; double Autocorrelation; double ClusterTendency; double ClusterShade; double ClusterProminence; double FirstMeasureOfInformationCorrelation; double SecondMeasureOfInformationCorrelation; }; class MITKCLUTILITIES_EXPORT GIFCooccurenceMatrix2 : public AbstractGlobalImageFeature { public: - mitkClassMacro(GIFCooccurenceMatrix2,AbstractGlobalImageFeature) + mitkClassMacro(GIFCooccurenceMatrix2, AbstractGlobalImageFeature); itkFactorylessNewMacro(Self) itkCloneMacro(Self) GIFCooccurenceMatrix2(); /** * \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); + itkGetConstMacro(Range,double); itkSetMacro(Range, double); itkGetConstMacro(Direction, unsigned int); itkSetMacro(Direction, unsigned int); struct GIFCooccurenceMatrix2Configuration { double range; unsigned int direction; }; private: double m_Range; unsigned int m_Direction; }; } #endif //mitkGIFCooccurenceMatrix2_h diff --git a/Modules/Classification/CLUtilities/include/mitkGIFFirstOrderStatistics.h b/Modules/Classification/CLUtilities/include/mitkGIFFirstOrderStatistics.h index 1d648411a4..fd3ba52ab4 100644 --- a/Modules/Classification/CLUtilities/include/mitkGIFFirstOrderStatistics.h +++ b/Modules/Classification/CLUtilities/include/mitkGIFFirstOrderStatistics.h @@ -1,63 +1,67 @@ /*=================================================================== 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 mitkGIFFirstOrderStatistics_h #define mitkGIFFirstOrderStatistics_h #include #include #include namespace mitk { class MITKCLUTILITIES_EXPORT GIFFirstOrderStatistics : public AbstractGlobalImageFeature { public: mitkClassMacro(GIFFirstOrderStatistics,AbstractGlobalImageFeature) itkFactorylessNewMacro(Self) itkCloneMacro(Self) GIFFirstOrderStatistics(); /** * \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; itkGetConstMacro(Range,double); itkSetMacro(Range, double); itkGetConstMacro(HistogramSize,int); itkSetMacro(HistogramSize, int); itkGetConstMacro(UseCtRange,bool); itkSetMacro(UseCtRange, bool); + virtual void CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList); + virtual void AddArguments(mitkCommandLineParser &parser); + + struct ParameterStruct { int m_HistogramSize; bool m_UseCtRange; }; private: double m_Range; int m_HistogramSize; bool m_UseCtRange; }; } #endif //mitkGIFFirstOrderStatistics_h diff --git a/Modules/Classification/CLUtilities/include/mitkGIFGrayLevelRunLength.h b/Modules/Classification/CLUtilities/include/mitkGIFGrayLevelRunLength.h index bb0d665bf8..96140b0537 100644 --- a/Modules/Classification/CLUtilities/include/mitkGIFGrayLevelRunLength.h +++ b/Modules/Classification/CLUtilities/include/mitkGIFGrayLevelRunLength.h @@ -1,51 +1,55 @@ #ifndef mitkGIFGrayLevelRunLength_h #define mitkGIFGrayLevelRunLength_h #include #include #include namespace mitk { class MITKCLUTILITIES_EXPORT GIFGrayLevelRunLength : public AbstractGlobalImageFeature { public: mitkClassMacro(GIFGrayLevelRunLength,AbstractGlobalImageFeature) itkFactorylessNewMacro(Self) itkCloneMacro(Self) GIFGrayLevelRunLength(); /** * \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; itkGetConstMacro(Range,double); itkSetMacro(Range, double); itkGetConstMacro(UseCtRange, bool); itkSetMacro(UseCtRange, bool); itkGetConstMacro(Direction, unsigned int); itkSetMacro(Direction, unsigned int); + virtual void CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList); + virtual void AddArguments(mitkCommandLineParser &parser); + + struct ParameterStruct { bool m_UseCtRange; double m_Range; unsigned int m_Direction; }; private: double m_Range; bool m_UseCtRange; unsigned int m_Direction; }; } #endif //mitkGIFGrayLevelRunLength_h diff --git a/Modules/Classification/CLUtilities/include/mitkGIFGrayLevelSizeZone.h b/Modules/Classification/CLUtilities/include/mitkGIFGrayLevelSizeZone.h index 0bf67bf627..c2a3ae9e66 100644 --- a/Modules/Classification/CLUtilities/include/mitkGIFGrayLevelSizeZone.h +++ b/Modules/Classification/CLUtilities/include/mitkGIFGrayLevelSizeZone.h @@ -1,51 +1,55 @@ #ifndef mitkGIFGrayLevelSizeZone_h #define mitkGIFGrayLevelSizeZone_h #include #include #include namespace mitk { class MITKCLUTILITIES_EXPORT GIFGrayLevelSizeZone : public AbstractGlobalImageFeature { public: mitkClassMacro(GIFGrayLevelSizeZone,AbstractGlobalImageFeature) itkFactorylessNewMacro(Self) itkCloneMacro(Self) GIFGrayLevelSizeZone(); /** * \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; itkGetConstMacro(Range,double); itkSetMacro(Range, double); itkGetConstMacro(UseCtRange, bool); itkSetMacro(UseCtRange, bool); itkGetConstMacro(Direction, unsigned int); itkSetMacro(Direction, unsigned int); + virtual void CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList); + virtual void AddArguments(mitkCommandLineParser &parser); + + struct ParameterStruct { bool m_UseCtRange; double m_Range; unsigned int m_Direction; }; private: double m_Range; bool m_UseCtRange; unsigned int m_Direction; }; } #endif //mitkGIFGrayLevelSizeZone_h diff --git a/Modules/Classification/CLUtilities/include/mitkGIFNeighbourhoodGreyLevelDifference.h b/Modules/Classification/CLUtilities/include/mitkGIFNeighbourhoodGreyLevelDifference.h index d730ee5c30..f8118ca102 100644 --- a/Modules/Classification/CLUtilities/include/mitkGIFNeighbourhoodGreyLevelDifference.h +++ b/Modules/Classification/CLUtilities/include/mitkGIFNeighbourhoodGreyLevelDifference.h @@ -1,51 +1,55 @@ #ifndef mitkGIFNeighbourhoodGreyLevelDifference_h #define mitkGIFNeighbourhoodGreyLevelDifference_h #include #include #include namespace mitk { class MITKCLUTILITIES_EXPORT GIFNeighbourhoodGreyLevelDifference : public AbstractGlobalImageFeature { public: mitkClassMacro(GIFNeighbourhoodGreyLevelDifference,AbstractGlobalImageFeature) itkFactorylessNewMacro(Self) itkCloneMacro(Self) GIFNeighbourhoodGreyLevelDifference(); /** * \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; itkGetConstMacro(Range,double); itkSetMacro(Range, double); itkGetConstMacro(UseCtRange, bool); itkSetMacro(UseCtRange, bool); itkGetConstMacro(Direction, unsigned int); itkSetMacro(Direction, unsigned int); + virtual void CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList); + virtual void AddArguments(mitkCommandLineParser &parser); + + struct ParameterStruct { bool m_UseCtRange; double m_Range; unsigned int m_Direction; }; private: double m_Range; bool m_UseCtRange; unsigned int m_Direction; }; } #endif //mitkGIFNeighbourhoodGreyLevelDifference_h diff --git a/Modules/Classification/CLUtilities/include/mitkGIFNeighbouringGreyLevelDependenceFeatures.h b/Modules/Classification/CLUtilities/include/mitkGIFNeighbouringGreyLevelDependenceFeatures.h new file mode 100644 index 0000000000..63a00b12c6 --- /dev/null +++ b/Modules/Classification/CLUtilities/include/mitkGIFNeighbouringGreyLevelDependenceFeatures.h @@ -0,0 +1,144 @@ +#ifndef mitkGIFNeighbouringGreyLevelDependenceFeatures_h +#define mitkGIFNeighbouringGreyLevelDependenceFeatures_h + +#include +#include +#include + +#include + +namespace mitk +{ + struct NGLDMMatrixHolder + { + public: + NGLDMMatrixHolder(double min, double max, int number, int depenence); + + int IntensityToIndex(double intensity); + double IndexToMinIntensity(int index); + double IndexToMeanIntensity(int index); + double IndexToMaxIntensity(int index); + + double m_MinimumRange; + double m_MaximumRange; + double m_Stepsize; + int m_NumberOfDependences; + int m_NumberOfBins; + Eigen::MatrixXd m_Matrix; + + int m_NeighbourhoodSize; + unsigned long m_NumberOfNeighbourVoxels; + unsigned long m_NumberOfDependenceNeighbourVoxels; + unsigned long m_NumberOfNeighbourhoods; + unsigned long m_NumberOfCompleteNeighbourhoods; + }; + + struct NGLDMMatrixFeatures + { + NGLDMMatrixFeatures() : + LowDependenceEmphasis(0), + HighDependenceEmphasis(0), + LowGreyLevelCountEmphasis(0), + HighGreyLevelCountEmphasis(0), + LowDependenceLowGreyLevelEmphasis(0), + LowDependenceHighGreyLevelEmphasis(0), + HighDependenceLowGreyLevelEmphasis(0), + HighDependenceHighGreyLevelEmphasis(0), + GreyLevelNonUniformity(0), + GreyLevelNonUniformityNormalised(0), + DependenceCountNonUniformity(0), + DependenceCountNonUniformityNormalised(0), + DependenceCountPercentage(0), + GreyLevelVariance(0), + DependenceCountVariance(0), + DependenceCountEntropy(0), + MeanGreyLevelCount(0), + MeanDependenceCount(0), + ExpectedNeighbourhoodSize(0), + AverageNeighbourhoodSize(0), + AverageIncompleteNeighbourhoodSize(0), + PercentageOfCompleteNeighbourhoods(0), + PercentageOfDependenceNeighbours(0) + { + } + + public: + double LowDependenceEmphasis; + double HighDependenceEmphasis; + double LowGreyLevelCountEmphasis; + double HighGreyLevelCountEmphasis; + double LowDependenceLowGreyLevelEmphasis; + double LowDependenceHighGreyLevelEmphasis; + double HighDependenceLowGreyLevelEmphasis; + double HighDependenceHighGreyLevelEmphasis; + + double GreyLevelNonUniformity; + double GreyLevelNonUniformityNormalised; + double DependenceCountNonUniformity; + double DependenceCountNonUniformityNormalised; + + double DependenceCountPercentage; + double GreyLevelVariance; + double DependenceCountVariance; + double DependenceCountEntropy; + double MeanGreyLevelCount; + double MeanDependenceCount; + + double ExpectedNeighbourhoodSize; + double AverageNeighbourhoodSize; + double AverageIncompleteNeighbourhoodSize; + double PercentageOfCompleteNeighbourhoods; + double PercentageOfDependenceNeighbours; + + }; + + + class MITKCLUTILITIES_EXPORT GIFNeighbouringGreyLevelDependenceFeature : public AbstractGlobalImageFeature + { + public: + mitkClassMacro(GIFNeighbouringGreyLevelDependenceFeature, AbstractGlobalImageFeature) + itkFactorylessNewMacro(Self) + itkCloneMacro(Self) + + GIFNeighbouringGreyLevelDependenceFeature(); + + /** + * \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; + + itkGetConstMacro(Range,double); + itkSetMacro(Range, double); + itkGetConstMacro(Direction, unsigned int); + itkSetMacro(Direction, unsigned int); + itkGetConstMacro(Alpha, int); + itkSetMacro(Alpha, int); + itkGetConstMacro(Bins, int); + itkSetMacro(Bins, int); + + virtual void CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList); + virtual void AddArguments(mitkCommandLineParser &parser); + + + struct GIFNeighbouringGreyLevelDependenceFeatureConfiguration + { + double range; + unsigned int direction; + int alpha; + int bins; + }; + + private: + double m_Range; + unsigned int m_Direction; + int m_Alpha; + int m_Bins; + }; + +} +#endif //mitkGIFNeighbouringGreyLevelDependenceFeature_h diff --git a/Modules/Classification/CLUtilities/include/mitkGIFVolumetricStatistics.h b/Modules/Classification/CLUtilities/include/mitkGIFVolumetricStatistics.h index ff6ae4bd9f..2310cacf87 100644 --- a/Modules/Classification/CLUtilities/include/mitkGIFVolumetricStatistics.h +++ b/Modules/Classification/CLUtilities/include/mitkGIFVolumetricStatistics.h @@ -1,52 +1,51 @@ /*=================================================================== 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 mitkGIFVolumetricStatistics_h #define mitkGIFVolumetricStatistics_h #include #include #include namespace mitk { class MITKCLUTILITIES_EXPORT GIFVolumetricStatistics : public AbstractGlobalImageFeature { public: mitkClassMacro(GIFVolumetricStatistics,AbstractGlobalImageFeature) itkFactorylessNewMacro(Self) itkCloneMacro(Self) GIFVolumetricStatistics(); /** * \brief Calculates the Cooccurence-Matrix based features for this class. */ virtual FeatureListType CalculateFeatures(const Image::Pointer & image, const Image::Pointer &feature); /** * \brief Returns a list of the names of all features that are calculated from this class */ virtual FeatureNameListType GetFeatureNames(); - itkGetConstMacro(Range,double); - itkSetMacro(Range, double); + virtual void CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList); + virtual void AddArguments(mitkCommandLineParser &parser); private: - double m_Range; }; } #endif //mitkGIFVolumetricStatistics_h diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp index f086dbc940..960cd29e5e 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp @@ -1,312 +1,362 @@ #include // MITK #include #include #include // ITK #include #include // STL #include template void CalculateCoocurenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFCooccurenceMatrix::FeatureListType & featureList, mitk::GIFCooccurenceMatrix::GIFCooccurenceMatrixConfiguration config) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::Statistics::EnhancedScalarImageToTextureFeaturesFilter FilterType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef typename FilterType::TextureFeaturesFilterType TextureFilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename FilterType::Pointer filter = FilterType::New(); typename FilterType::OffsetVector::Pointer newOffset = FilterType::OffsetVector::New(); auto oldOffsets = filter->GetOffsets(); auto oldOffsetsIterator = oldOffsets->Begin(); while(oldOffsetsIterator != oldOffsets->End()) { bool continueOuterLoop = false; typename FilterType::OffsetType offset = oldOffsetsIterator->Value(); for (unsigned int i = 0; i < VImageDimension; ++i) { offset[i] *= config.range; if (config.direction == i + 2 && offset[i] != 0) { continueOuterLoop = true; } } if (config.direction == 1) { offset[0] = 0; offset[1] = 0; offset[2] = 1; newOffset->push_back(offset); break; } oldOffsetsIterator++; if (continueOuterLoop) continue; newOffset->push_back(offset); } filter->SetOffsets(newOffset); // All features are required typename FilterType::FeatureNameVectorPointer requestedFeatures = FilterType::FeatureNameVector::New(); requestedFeatures->push_back(TextureFilterType::Energy); requestedFeatures->push_back(TextureFilterType::Entropy); requestedFeatures->push_back(TextureFilterType::Correlation); requestedFeatures->push_back(TextureFilterType::InverseDifferenceMoment); requestedFeatures->push_back(TextureFilterType::Inertia); requestedFeatures->push_back(TextureFilterType::ClusterShade); requestedFeatures->push_back(TextureFilterType::ClusterProminence); requestedFeatures->push_back(TextureFilterType::HaralickCorrelation); requestedFeatures->push_back(TextureFilterType::Autocorrelation); requestedFeatures->push_back(TextureFilterType::Contrast); requestedFeatures->push_back(TextureFilterType::Dissimilarity); requestedFeatures->push_back(TextureFilterType::MaximumProbability); requestedFeatures->push_back(TextureFilterType::InverseVariance); requestedFeatures->push_back(TextureFilterType::Homogeneity1); requestedFeatures->push_back(TextureFilterType::ClusterTendency); requestedFeatures->push_back(TextureFilterType::Variance); requestedFeatures->push_back(TextureFilterType::SumAverage); requestedFeatures->push_back(TextureFilterType::SumEntropy); requestedFeatures->push_back(TextureFilterType::SumVariance); requestedFeatures->push_back(TextureFilterType::DifferenceAverage); requestedFeatures->push_back(TextureFilterType::DifferenceEntropy); requestedFeatures->push_back(TextureFilterType::DifferenceVariance); requestedFeatures->push_back(TextureFilterType::InverseDifferenceMomentNormalized); requestedFeatures->push_back(TextureFilterType::InverseDifferenceNormalized); requestedFeatures->push_back(TextureFilterType::InverseDifference); requestedFeatures->push_back(TextureFilterType::JointAverage); requestedFeatures->push_back(TextureFilterType::FirstMeasureOfInformationCorrelation); requestedFeatures->push_back(TextureFilterType::SecondMeasureOfInformationCorrelation); typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); filter->SetInput(itkImage); filter->SetMaskImage(maskImage); filter->SetRequestedFeatures(requestedFeatures); filter->SetPixelValueMinMax(minMaxComputer->GetMinimum()-0.5,minMaxComputer->GetMaximum()+0.5); //filter->SetPixelValueMinMax(-1024,3096); //filter->SetNumberOfBinsPerAxis(5); filter->Update(); auto featureMeans = filter->GetFeatureMeans (); auto featureStd = filter->GetFeatureStandardDeviations(); std::ostringstream ss; ss << config.range; std::string strRange = ss.str(); for (std::size_t i = 0; i < featureMeans->size(); ++i) { switch (i) { case TextureFilterType::Energy : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Energy Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Energy Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Entropy : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Entropy Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Entropy Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Correlation : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Correlation Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Correlation Std.",featureStd->ElementAt(i))); break; case TextureFilterType::InverseDifferenceMoment : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMoment Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMoment Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Inertia : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Inertia Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Inertia Std.",featureStd->ElementAt(i))); break; case TextureFilterType::ClusterShade : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterShade Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterShade Std.",featureStd->ElementAt(i))); break; case TextureFilterType::ClusterProminence : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterProminence Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterProminence Std.",featureStd->ElementAt(i))); break; case TextureFilterType::HaralickCorrelation : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") HaralickCorrelation Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") HaralickCorrelation Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Autocorrelation : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Autocorrelation Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Autocorrelation Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Contrast : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Contrast Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Contrast Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Dissimilarity : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Dissimilarity Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Dissimilarity Std.",featureStd->ElementAt(i))); break; case TextureFilterType::MaximumProbability : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") MaximumProbability Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") MaximumProbability Std.",featureStd->ElementAt(i))); break; case TextureFilterType::InverseVariance : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseVariance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseVariance Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Homogeneity1 : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Homogeneity1 Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Homogeneity1 Std.",featureStd->ElementAt(i))); break; case TextureFilterType::ClusterTendency : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterTendency Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterTendency Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Variance : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Variance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Variance Std.",featureStd->ElementAt(i))); break; case TextureFilterType::SumAverage : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumAverage Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumAverage Std.",featureStd->ElementAt(i))); break; case TextureFilterType::SumEntropy : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumEntropy Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumEntropy Std.",featureStd->ElementAt(i))); break; case TextureFilterType::SumVariance : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumVariance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumVariance Std.",featureStd->ElementAt(i))); break; case TextureFilterType::DifferenceAverage : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceAverage Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceAverage Std.",featureStd->ElementAt(i))); break; case TextureFilterType::DifferenceEntropy : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceEntropy Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceEntropy Std.",featureStd->ElementAt(i))); break; case TextureFilterType::DifferenceVariance : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceVariance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceVariance Std.",featureStd->ElementAt(i))); break; case TextureFilterType::InverseDifferenceMomentNormalized : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMomentNormalized Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMomentNormalized Std.",featureStd->ElementAt(i))); break; case TextureFilterType::InverseDifferenceNormalized : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceNormalized Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceNormalized Std.",featureStd->ElementAt(i))); break; case TextureFilterType::InverseDifference : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifference Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifference Std.",featureStd->ElementAt(i))); break; case TextureFilterType::JointAverage : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") JointAverage Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") JointAverage Std.",featureStd->ElementAt(i))); break; case TextureFilterType::FirstMeasureOfInformationCorrelation : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") FirstMeasureOfInformationCorrelation Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") FirstMeasureOfInformationCorrelation Std.",featureStd->ElementAt(i))); break; case TextureFilterType::SecondMeasureOfInformationCorrelation : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SecondMeasureOfInformationCorrelation Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SecondMeasureOfInformationCorrelation Std.",featureStd->ElementAt(i))); break; default: break; } } } mitk::GIFCooccurenceMatrix::GIFCooccurenceMatrix(): m_Range(1.0), m_Direction(0) { + SetShortName("cooc"); + SetLongName("cooccurence"); } mitk::GIFCooccurenceMatrix::FeatureListType mitk::GIFCooccurenceMatrix::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; GIFCooccurenceMatrixConfiguration config; config.direction = m_Direction; config.range = m_Range; AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config); return featureList; } mitk::GIFCooccurenceMatrix::FeatureNameListType mitk::GIFCooccurenceMatrix::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("co-occ. Energy Means"); featureList.push_back("co-occ. Energy Std."); featureList.push_back("co-occ. Entropy Means"); featureList.push_back("co-occ. Entropy Std."); featureList.push_back("co-occ. Correlation Means"); featureList.push_back("co-occ. Correlation Std."); featureList.push_back("co-occ. InverseDifferenceMoment Means"); featureList.push_back("co-occ. InverseDifferenceMoment Std."); featureList.push_back("co-occ. Inertia Means"); featureList.push_back("co-occ. Inertia Std."); featureList.push_back("co-occ. ClusterShade Means"); featureList.push_back("co-occ. ClusterShade Std."); featureList.push_back("co-occ. ClusterProminence Means"); featureList.push_back("co-occ. ClusterProminence Std."); featureList.push_back("co-occ. HaralickCorrelation Means"); featureList.push_back("co-occ. HaralickCorrelation Std."); featureList.push_back("co-occ. Autocorrelation Means"); featureList.push_back("co-occ. Autocorrelation Std."); featureList.push_back("co-occ. Contrast Means"); featureList.push_back("co-occ. Contrast Std."); featureList.push_back("co-occ. Dissimilarity Means"); featureList.push_back("co-occ. Dissimilarity Std."); featureList.push_back("co-occ. MaximumProbability Means"); featureList.push_back("co-occ. MaximumProbability Std."); featureList.push_back("co-occ. InverseVariance Means"); featureList.push_back("co-occ. InverseVariance Std."); featureList.push_back("co-occ. Homogeneity1 Means"); featureList.push_back("co-occ. Homogeneity1 Std."); featureList.push_back("co-occ. ClusterTendency Means"); featureList.push_back("co-occ. ClusterTendency Std."); featureList.push_back("co-occ. Variance Means"); featureList.push_back("co-occ. Variance Std."); featureList.push_back("co-occ. SumAverage Means"); featureList.push_back("co-occ. SumAverage Std."); featureList.push_back("co-occ. SumEntropy Means"); featureList.push_back("co-occ. SumEntropy Std."); featureList.push_back("co-occ. SumVariance Means"); featureList.push_back("co-occ. SumVariance Std."); featureList.push_back("co-occ. DifferenceAverage Means"); featureList.push_back("co-occ. DifferenceAverage Std."); featureList.push_back("co-occ. DifferenceEntropy Means"); featureList.push_back("co-occ. DifferenceEntropy Std."); featureList.push_back("co-occ. DifferenceVariance Means"); featureList.push_back("co-occ. DifferenceVariance Std."); featureList.push_back("co-occ. InverseDifferenceMomentNormalized Means"); featureList.push_back("co-occ. InverseDifferenceMomentNormalized Std."); featureList.push_back("co-occ. InverseDifferenceNormalized Means"); featureList.push_back("co-occ. InverseDifferenceNormalized Std."); featureList.push_back("co-occ. InverseDifference Means"); featureList.push_back("co-occ. InverseDifference Std."); featureList.push_back("co-occ. JointAverage Means"); featureList.push_back("co-occ. JointAverage Std."); featureList.push_back("co-occ. FirstMeasurementOfInformationCorrelation Means"); featureList.push_back("co-occ. FirstMeasurementOfInformationCorrelation Std."); featureList.push_back("co-occ. SecondMeasurementOfInformationCorrelation Means"); featureList.push_back("co-occ. SecondMeasurementOfInformationCorrelation Std."); return featureList; } + + + +void mitk::GIFCooccurenceMatrix::AddArguments(mitkCommandLineParser &parser) +{ + std::string name = GetOptionPrefix(); + + parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features", us::Any()); + parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::String, "Cooc 2 Range", "Define the range that is used (Semicolon-separated)", us::Any()); + parser.addArgument(name + "::direction", name + "::dir", mitkCommandLineParser::String, "Int", "Allows to specify the direction for Cooc and RL. 0: All directions, 1: Only single direction (Test purpose), 2,3,4... without dimension 0,1,2... ", us::Any()); +} + +void +mitk::GIFCooccurenceMatrix::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())) + { + int direction = 0; + if (parsedArgs.count(name + "::direction")) + { + direction = SplitDouble(parsedArgs[name + "::direction"].ToString(), ';')[0]; + } + std::vector ranges; + if (parsedArgs.count(name + "::range")) + { + ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); + } + else + { + ranges.push_back(1); + } + + for (std::size_t i = 0; i < ranges.size(); ++i) + { + MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; + mitk::GIFCooccurenceMatrix::Pointer coocCalculator = mitk::GIFCooccurenceMatrix::New(); + coocCalculator->SetRange(ranges[i]); + coocCalculator->SetDirection(direction); + auto localResults = coocCalculator->CalculateFeatures(feature, maskNoNAN); + featureList.insert(featureList.end(), localResults.begin(), localResults.end()); + MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; + } + } +} + diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp index 9c2a7a559a..32569b9c54 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp @@ -1,591 +1,642 @@ #include // MITK #include #include #include // ITK #include #include #include #include // STL #include static void MatrixFeaturesTo(mitk::CoocurenceMatrixFeatures features, std::string prefix, mitk::GIFCooccurenceMatrix2::FeatureListType &featureList); static void CalculateMeanAndStdDevFeatures(std::vector featureList, mitk::CoocurenceMatrixFeatures &meanFeature, mitk::CoocurenceMatrixFeatures &stdFeature); static void NormalizeMatrixFeature(mitk::CoocurenceMatrixFeatures &features, std::size_t number); mitk::CoocurenceMatrixHolder::CoocurenceMatrixHolder(double min, double max, int number) : m_MinimumRange(min), m_MaximumRange(max), m_NumberOfBins(number) { m_Matrix.resize(number, number); m_Matrix.fill(0); m_Stepsize = (max - min) / (number); } int mitk::CoocurenceMatrixHolder::IntensityToIndex(double intensity) { return std::floor((intensity - m_MinimumRange) / m_Stepsize); } double mitk::CoocurenceMatrixHolder::IndexToMinIntensity(int index) { return m_MinimumRange + index * m_Stepsize; } double mitk::CoocurenceMatrixHolder::IndexToMeanIntensity(int index) { return m_MinimumRange + (index+0.5) * m_Stepsize; } double mitk::CoocurenceMatrixHolder::IndexToMaxIntensity(int index) { return m_MinimumRange + (index + 1) * m_Stepsize; } template void CalculateCoOcMatrix(itk::Image* itkImage, itk::Image* mask, itk::Offset offset, int range, mitk::CoocurenceMatrixHolder &holder) { typedef itk::Image ImageType; typedef itk::Image MaskImageType; typedef itk::ShapedNeighborhoodIterator ShapeIterType; typedef itk::ShapedNeighborhoodIterator ShapeMaskIterType; typedef itk::ImageRegionConstIterator ConstIterType; typedef itk::ImageRegionConstIterator ConstMaskIterType; itk::Size radius; radius.Fill(range+1); ShapeIterType imageOffsetIter(radius, itkImage, itkImage->GetLargestPossibleRegion()); ShapeMaskIterType maskOffsetIter(radius, mask, mask->GetLargestPossibleRegion()); imageOffsetIter.ActivateOffset(offset); maskOffsetIter.ActivateOffset(offset); ConstIterType imageIter(itkImage, itkImage->GetLargestPossibleRegion()); ConstMaskIterType maskIter(mask, mask->GetLargestPossibleRegion()); // iterator.GetIndex() + ci.GetNeighborhoodOffset() auto region = mask->GetLargestPossibleRegion(); while (!maskIter.IsAtEnd()) { auto ciMask = maskOffsetIter.Begin(); auto ciValue = imageOffsetIter.Begin(); if (maskIter.Value() > 0 && ciMask.Get() > 0 && imageIter.Get() == imageIter.Get() && ciValue.Get() == ciValue.Get() && region.IsInside(maskOffsetIter.GetIndex() + ciMask.GetNeighborhoodOffset())) { int i = holder.IntensityToIndex(imageIter.Get()); int j = holder.IntensityToIndex(ciValue.Get()); holder.m_Matrix(i, j) += 1; holder.m_Matrix(j, i) += 1; } ++imageOffsetIter; ++maskOffsetIter; ++imageIter; ++maskIter; } } void CalculateFeatures( mitk::CoocurenceMatrixHolder &holder, mitk::CoocurenceMatrixFeatures & results ) { auto pijMatrix = holder.m_Matrix; auto piMatrix = holder.m_Matrix; auto pjMatrix = holder.m_Matrix; double Ng = holder.m_NumberOfBins; int NgSize = holder.m_NumberOfBins; pijMatrix /= pijMatrix.sum(); piMatrix.rowwise().normalize(); pjMatrix.colwise().normalize(); for (int i = 0; i < holder.m_NumberOfBins; ++i) for (int j = 0; j < holder.m_NumberOfBins; ++j) { if (pijMatrix(i, j) != pijMatrix(i, j)) pijMatrix(i, j) = 0; if (piMatrix(i, j) != piMatrix(i, j)) piMatrix(i, j) = 0; if (pjMatrix(i, j) != pjMatrix(i, j)) pjMatrix(i, j) = 0; } Eigen::VectorXd piVector = pijMatrix.colwise().sum(); Eigen::VectorXd pjVector = pijMatrix.rowwise().sum(); double sigmai = 0;; for (int i = 0; i < holder.m_NumberOfBins; ++i) { double iInt = holder.IndexToMeanIntensity(i); results.RowAverage += iInt * piVector(i); if (piVector(i) > 0) { results.RowEntropy -= piVector(i) * std::log(piVector(i)) / std::log(2); } } for (int i = 0; i < holder.m_NumberOfBins; ++i) { double iInt = holder.IndexToMeanIntensity(i); results.RowVariance += (iInt - results.RowAverage)*(iInt - results.RowAverage) * piVector(i); } results.RowMaximum = piVector.maxCoeff(); sigmai = std::sqrt(results.RowVariance); Eigen::VectorXd pimj(NgSize); pimj.fill(0); Eigen::VectorXd pipj(2*NgSize); pipj.fill(0); results.JointMaximum += pijMatrix.maxCoeff(); for (int i = 0; i < holder.m_NumberOfBins; ++i) { for (int j = 0; j < holder.m_NumberOfBins; ++j) { double iInt = holder.IndexToMeanIntensity(i); double jInt = holder.IndexToMeanIntensity(j); double pij = pijMatrix(i, j); int deltaK = (i - j)>0?(i-j) : (j-i); pimj(deltaK) += pij; pipj(i + j) += pij; results.JointAverage += iInt * pij; if (pij > 0) { results.JointEntropy -= pij * std::log(pij) / std::log(2); results.FirstRowColumnEntropy -= pij * std::log(piVector(i)*pjVector(j)) / std::log(2); } if (piVector(i) > 0 && pjVector(j) > 0 ) { results.SecondRowColumnEntropy -= piVector(i)*pjVector(j) * std::log(piVector(i)*pjVector(j)) / std::log(2); } results.AngularSecondMoment += pij*pij; results.Contrast += (iInt - jInt)* (iInt - jInt) * pij; results.Dissimilarity += std::abs(iInt - jInt) * pij; results.InverseDifference += pij / (1 + (std::abs(iInt - jInt))); results.InverseDifferenceNormalised += pij / (1 + (std::abs(iInt - jInt) / Ng)); results.InverseDifferenceMoment += pij / (1 + (iInt - jInt)*(iInt - jInt)); results.InverseDifferenceMomentNormalised += pij / (1 + (iInt - jInt)*(iInt - jInt)/Ng/Ng); results.Autocorrelation += iInt*jInt * pij; double cluster = (iInt + jInt - 2 * results.RowAverage); results.ClusterTendency += cluster*cluster * pij; results.ClusterShade += cluster*cluster*cluster * pij; results.ClusterProminence += cluster*cluster*cluster*cluster * pij; if (iInt != jInt) { results.InverseVariance += pij / (iInt - jInt) / (iInt - jInt); } } } results.Correlation = 1 / sigmai / sigmai * (-results.RowAverage*results.RowAverage+ results.Autocorrelation); results.FirstMeasureOfInformationCorrelation = (results.JointEntropy - results.FirstRowColumnEntropy) / results.RowEntropy; if (results.JointEntropy < results.SecondRowColumnEntropy) { results.SecondMeasureOfInformationCorrelation = sqrt(1 - exp(-2 * (results.SecondRowColumnEntropy - results.JointEntropy))); } else { results.SecondMeasureOfInformationCorrelation = 0; } for (int i = 0; i < holder.m_NumberOfBins; ++i) { for (int j = 0; j < holder.m_NumberOfBins; ++j) { double iInt = holder.IndexToMeanIntensity(i); double jInt = holder.IndexToMeanIntensity(j); double pij = pijMatrix(i, j); results.JointVariance += (iInt - results.JointAverage)* (iInt - results.JointAverage)*pij; } } for (int k = 0; k < NgSize; ++k) { results.DifferenceAverage += k* pimj(k); if (pimj(k) > 0) { results.DifferenceEntropy -= pimj(k) * log(pimj(k)) / std::log(2); } } for (int k = 0; k < NgSize; ++k) { results.DifferenceVariance += (results.DifferenceAverage-k)* (results.DifferenceAverage-k)*pimj(k); } for (int k = 0; k <2* NgSize ; ++k) { results.SumAverage += (2+k)* pipj(k); if (pipj(k) > 0) { results.SumEntropy -= pipj(k) * log(pipj(k)) / std::log(2); } } for (int k = 0; k < 2*NgSize; ++k) { results.SumVariance += (2+k - results.SumAverage)* (2+k - results.SumAverage)*pipj(k); } //MITK_INFO << std::endl << holder.m_Matrix; //MITK_INFO << std::endl << pijMatrix; //MITK_INFO << std::endl << piMatrix; //MITK_INFO << std::endl << pjMatrix; //for (int i = 0; i < holder.m_NumberOfBins; ++i) //{ // MITK_INFO << "Bin " << i << " Min: " << holder.IndexToMinIntensity(i) << " Max: " << holder.IndexToMaxIntensity(i); //} //MITK_INFO << pimj; //MITK_INFO << pipj; } template void CalculateCoocurenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFCooccurenceMatrix2::FeatureListType & featureList, mitk::GIFCooccurenceMatrix2::GIFCooccurenceMatrix2Configuration config) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef itk::Neighborhood NeighborhoodType; typedef itk::Offset OffsetType; typedef itk::Statistics::EnhancedScalarImageToTextureFeaturesFilter FilterType; typedef typename FilterType::TextureFeaturesFilterType TextureFilterType; /////////////////////////////////////////////////////////////////////////////////////////////// typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); double rangeMin = -0.5 + minMaxComputer->GetMinimum(); double rangeMax = 0.5 + minMaxComputer->GetMaximum(); // Define Range int numberOfBins = 6; 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 (int i = 0; i < VImageDimension; ++i) { offset[i] *= config.range; if (config.direction == i + 2 && offset[i] != 0) { useOffset = false; } } if (useOffset) { offsetVector.push_back(offset); } } std::vector resultVector; mitk::CoocurenceMatrixHolder holderOverall(rangeMin, rangeMax, numberOfBins); mitk::CoocurenceMatrixFeatures overallFeature; for (int i = 0; i < offsetVector.size(); ++i) { offset = offsetVector[i]; MITK_INFO << offset; mitk::CoocurenceMatrixHolder holder(rangeMin, rangeMax, numberOfBins); mitk::CoocurenceMatrixFeatures coocResults; CalculateCoOcMatrix(itkImage, maskImage, offset, config.range, holder); holderOverall.m_Matrix += holder.m_Matrix; CalculateFeatures(holder, coocResults); resultVector.push_back(coocResults); } CalculateFeatures(holderOverall, overallFeature); //NormalizeMatrixFeature(overallFeature, offsetVector.size()); mitk::CoocurenceMatrixFeatures featureMean; mitk::CoocurenceMatrixFeatures featureStd; CalculateMeanAndStdDevFeatures(resultVector, featureMean, featureStd); std::ostringstream ss; ss << config.range; std::string strRange = ss.str(); MatrixFeaturesTo(overallFeature, "co-occ. 2 (" + strRange + ") overall", featureList); MatrixFeaturesTo(featureMean, "co-occ. 2 (" + strRange + ") mean", featureList); MatrixFeaturesTo(featureStd, "co-occ. 2 (" + strRange + ") std.dev.", featureList); } static void MatrixFeaturesTo(mitk::CoocurenceMatrixFeatures features, std::string prefix, mitk::GIFCooccurenceMatrix2::FeatureListType &featureList) { featureList.push_back(std::make_pair(prefix + " Joint Maximum", features.JointMaximum)); featureList.push_back(std::make_pair(prefix + " Joint Average", features.JointAverage)); featureList.push_back(std::make_pair(prefix + " Joint Variance", features.JointVariance)); featureList.push_back(std::make_pair(prefix + " Joint Entropy", features.JointEntropy)); featureList.push_back(std::make_pair(prefix + " Row Maximum", features.RowMaximum)); featureList.push_back(std::make_pair(prefix + " Row Average", features.RowAverage)); featureList.push_back(std::make_pair(prefix + " Row Variance", features.RowVariance)); featureList.push_back(std::make_pair(prefix + " Row Entropy", features.RowEntropy)); featureList.push_back(std::make_pair(prefix + " First Row-Column Entropy", features.FirstRowColumnEntropy)); featureList.push_back(std::make_pair(prefix + " Second Row-Column Entropy", features.SecondRowColumnEntropy)); featureList.push_back(std::make_pair(prefix + " Difference Average", features.DifferenceAverage)); featureList.push_back(std::make_pair(prefix + " Difference Variance", features.DifferenceVariance)); featureList.push_back(std::make_pair(prefix + " Difference Entropy", features.DifferenceEntropy)); featureList.push_back(std::make_pair(prefix + " Sum Average", features.SumAverage)); featureList.push_back(std::make_pair(prefix + " Sum Variance", features.SumVariance)); featureList.push_back(std::make_pair(prefix + " Sum Entropy", features.SumEntropy)); featureList.push_back(std::make_pair(prefix + " Angular Second Moment", features.AngularSecondMoment)); featureList.push_back(std::make_pair(prefix + " Contrast", features.Contrast)); featureList.push_back(std::make_pair(prefix + " Dissimilarity", features.Dissimilarity)); featureList.push_back(std::make_pair(prefix + " Inverse Difference", features.InverseDifference)); featureList.push_back(std::make_pair(prefix + " Inverse Difference Normalized", features.InverseDifferenceNormalised)); featureList.push_back(std::make_pair(prefix + " Inverse Difference Moment", features.InverseDifferenceMoment)); featureList.push_back(std::make_pair(prefix + " Inverse Difference Moment Normalized", features.InverseDifferenceMomentNormalised)); featureList.push_back(std::make_pair(prefix + " Inverse Variance", features.InverseVariance)); featureList.push_back(std::make_pair(prefix + " Correlation", features.Correlation)); featureList.push_back(std::make_pair(prefix + " Autocorrleation", features.Autocorrelation)); featureList.push_back(std::make_pair(prefix + " Cluster Tendency", features.ClusterTendency)); featureList.push_back(std::make_pair(prefix + " Cluster Shade", features.ClusterShade)); featureList.push_back(std::make_pair(prefix + " Cluster Prominence", features.ClusterProminence)); featureList.push_back(std::make_pair(prefix + " First Measure of Information Correlation", features.FirstMeasureOfInformationCorrelation)); featureList.push_back(std::make_pair(prefix + " Second Measure of Information Correlation", features.SecondMeasureOfInformationCorrelation)); } static void CalculateMeanAndStdDevFeatures(std::vector featureList, mitk::CoocurenceMatrixFeatures &meanFeature, mitk::CoocurenceMatrixFeatures &stdFeature) { #define ADDFEATURE(a) meanFeature.a += featureList[i].a;stdFeature.a += featureList[i].a*featureList[i].a #define CALCVARIANCE(a) stdFeature.a =sqrt(stdFeature.a - meanFeature.a*meanFeature.a) for (int i = 0; i < featureList.size(); ++i) { ADDFEATURE(JointMaximum); ADDFEATURE(JointAverage); ADDFEATURE(JointVariance); ADDFEATURE(JointEntropy); ADDFEATURE(RowMaximum); ADDFEATURE(RowAverage); ADDFEATURE(RowVariance); ADDFEATURE(RowEntropy); ADDFEATURE(FirstRowColumnEntropy); ADDFEATURE(SecondRowColumnEntropy); ADDFEATURE(DifferenceAverage); ADDFEATURE(DifferenceVariance); ADDFEATURE(DifferenceEntropy); ADDFEATURE(SumAverage); ADDFEATURE(SumVariance); ADDFEATURE(SumEntropy); ADDFEATURE(AngularSecondMoment); ADDFEATURE(Contrast); ADDFEATURE(Dissimilarity); ADDFEATURE(InverseDifference); ADDFEATURE(InverseDifferenceNormalised); ADDFEATURE(InverseDifferenceMoment); ADDFEATURE(InverseDifferenceMomentNormalised); ADDFEATURE(InverseVariance); ADDFEATURE(Correlation); ADDFEATURE(Autocorrelation); ADDFEATURE(ClusterShade); ADDFEATURE(ClusterTendency); ADDFEATURE(ClusterProminence); ADDFEATURE(FirstMeasureOfInformationCorrelation); ADDFEATURE(SecondMeasureOfInformationCorrelation); } NormalizeMatrixFeature(meanFeature, featureList.size()); NormalizeMatrixFeature(stdFeature, featureList.size()); CALCVARIANCE(JointMaximum); CALCVARIANCE(JointAverage); CALCVARIANCE(JointVariance); CALCVARIANCE(JointEntropy); CALCVARIANCE(RowMaximum); CALCVARIANCE(RowAverage); CALCVARIANCE(RowVariance); CALCVARIANCE(RowEntropy); CALCVARIANCE(FirstRowColumnEntropy); CALCVARIANCE(SecondRowColumnEntropy); CALCVARIANCE(DifferenceAverage); CALCVARIANCE(DifferenceVariance); CALCVARIANCE(DifferenceEntropy); CALCVARIANCE(SumAverage); CALCVARIANCE(SumVariance); CALCVARIANCE(SumEntropy); CALCVARIANCE(AngularSecondMoment); CALCVARIANCE(Contrast); CALCVARIANCE(Dissimilarity); CALCVARIANCE(InverseDifference); CALCVARIANCE(InverseDifferenceNormalised); CALCVARIANCE(InverseDifferenceMoment); CALCVARIANCE(InverseDifferenceMomentNormalised); CALCVARIANCE(InverseVariance); CALCVARIANCE(Correlation); CALCVARIANCE(Autocorrelation); CALCVARIANCE(ClusterShade); CALCVARIANCE(ClusterTendency); CALCVARIANCE(ClusterProminence); CALCVARIANCE(FirstMeasureOfInformationCorrelation); CALCVARIANCE(SecondMeasureOfInformationCorrelation); #undef ADDFEATURE #undef CALCVARIANCE } static void NormalizeMatrixFeature(mitk::CoocurenceMatrixFeatures &features, std::size_t number) { features.JointMaximum = features.JointMaximum / number; features.JointAverage = features.JointAverage / number; features.JointVariance = features.JointVariance / number; features.JointEntropy = features.JointEntropy / number; features.RowMaximum = features.RowMaximum / number; features.RowAverage = features.RowAverage / number; features.RowVariance = features.RowVariance / number; features.RowEntropy = features.RowEntropy / number; features.FirstRowColumnEntropy = features.FirstRowColumnEntropy / number; features.SecondRowColumnEntropy = features.SecondRowColumnEntropy / number; features.DifferenceAverage = features.DifferenceAverage / number; features.DifferenceVariance = features.DifferenceVariance / number; features.DifferenceEntropy = features.DifferenceEntropy / number; features.SumAverage = features.SumAverage / number; features.SumVariance = features.SumVariance / number; features.SumEntropy = features.SumEntropy / number; features.AngularSecondMoment = features.AngularSecondMoment / number; features.Contrast = features.Contrast / number; features.Dissimilarity = features.Dissimilarity / number; features.InverseDifference = features.InverseDifference / number; features.InverseDifferenceNormalised = features.InverseDifferenceNormalised / number; features.InverseDifferenceMoment = features.InverseDifferenceMoment / number; features.InverseDifferenceMomentNormalised = features.InverseDifferenceMomentNormalised / number; features.InverseVariance = features.InverseVariance / number; features.Correlation = features.Correlation / number; features.Autocorrelation = features.Autocorrelation / number; features.ClusterShade = features.ClusterShade / number; features.ClusterTendency = features.ClusterTendency / number; features.ClusterProminence = features.ClusterProminence / number; features.FirstMeasureOfInformationCorrelation = features.FirstMeasureOfInformationCorrelation / number; features.SecondMeasureOfInformationCorrelation = features.SecondMeasureOfInformationCorrelation / number; } mitk::GIFCooccurenceMatrix2::GIFCooccurenceMatrix2(): m_Range(1.0), m_Direction(0) { + SetShortName("cooc2"); + SetLongName("cooccurence2"); } mitk::GIFCooccurenceMatrix2::FeatureListType mitk::GIFCooccurenceMatrix2::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; GIFCooccurenceMatrix2Configuration config; config.direction = m_Direction; config.range = m_Range; AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config); return featureList; } mitk::GIFCooccurenceMatrix2::FeatureNameListType mitk::GIFCooccurenceMatrix2::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("co-occ. Energy Means"); featureList.push_back("co-occ. Energy Std."); featureList.push_back("co-occ. Entropy Means"); featureList.push_back("co-occ. Entropy Std."); featureList.push_back("co-occ. Correlation Means"); featureList.push_back("co-occ. Correlation Std."); featureList.push_back("co-occ. InverseDifferenceMoment Means"); featureList.push_back("co-occ. InverseDifferenceMoment Std."); featureList.push_back("co-occ. Inertia Means"); featureList.push_back("co-occ. Inertia Std."); featureList.push_back("co-occ. ClusterShade Means"); featureList.push_back("co-occ. ClusterShade Std."); featureList.push_back("co-occ. ClusterProminence Means"); featureList.push_back("co-occ. ClusterProminence Std."); featureList.push_back("co-occ. HaralickCorrelation Means"); featureList.push_back("co-occ. HaralickCorrelation Std."); featureList.push_back("co-occ. Autocorrelation Means"); featureList.push_back("co-occ. Autocorrelation Std."); featureList.push_back("co-occ. Contrast Means"); featureList.push_back("co-occ. Contrast Std."); featureList.push_back("co-occ. Dissimilarity Means"); featureList.push_back("co-occ. Dissimilarity Std."); featureList.push_back("co-occ. MaximumProbability Means"); featureList.push_back("co-occ. MaximumProbability Std."); featureList.push_back("co-occ. InverseVariance Means"); featureList.push_back("co-occ. InverseVariance Std."); featureList.push_back("co-occ. Homogeneity1 Means"); featureList.push_back("co-occ. Homogeneity1 Std."); featureList.push_back("co-occ. ClusterTendency Means"); featureList.push_back("co-occ. ClusterTendency Std."); featureList.push_back("co-occ. Variance Means"); featureList.push_back("co-occ. Variance Std."); featureList.push_back("co-occ. SumAverage Means"); featureList.push_back("co-occ. SumAverage Std."); featureList.push_back("co-occ. SumEntropy Means"); featureList.push_back("co-occ. SumEntropy Std."); featureList.push_back("co-occ. SumVariance Means"); featureList.push_back("co-occ. SumVariance Std."); featureList.push_back("co-occ. DifferenceAverage Means"); featureList.push_back("co-occ. DifferenceAverage Std."); featureList.push_back("co-occ. DifferenceEntropy Means"); featureList.push_back("co-occ. DifferenceEntropy Std."); featureList.push_back("co-occ. DifferenceVariance Means"); featureList.push_back("co-occ. DifferenceVariance Std."); featureList.push_back("co-occ. InverseDifferenceMomentNormalized Means"); featureList.push_back("co-occ. InverseDifferenceMomentNormalized Std."); featureList.push_back("co-occ. InverseDifferenceNormalized Means"); featureList.push_back("co-occ. InverseDifferenceNormalized Std."); featureList.push_back("co-occ. InverseDifference Means"); featureList.push_back("co-occ. InverseDifference Std."); featureList.push_back("co-occ. JointAverage Means"); featureList.push_back("co-occ. JointAverage Std."); featureList.push_back("co-occ. FirstMeasurementOfInformationCorrelation Means"); featureList.push_back("co-occ. FirstMeasurementOfInformationCorrelation Std."); featureList.push_back("co-occ. SecondMeasurementOfInformationCorrelation Means"); featureList.push_back("co-occ. SecondMeasurementOfInformationCorrelation Std."); return featureList; } + + + + +void mitk::GIFCooccurenceMatrix2::AddArguments(mitkCommandLineParser &parser) +{ + std::string name = GetOptionPrefix(); + + parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features (new implementation)", us::Any()); + parser.addArgument(name+"::range", name+"::range", mitkCommandLineParser::String, "Cooc 2 Range", "Define the range that is used (Semicolon-separated)", us::Any()); + parser.addArgument(name+"::direction", name+"::dir", mitkCommandLineParser::String, "Int", "Allows to specify the direction for Cooc and RL. 0: All directions, 1: Only single direction (Test purpose), 2,3,4... without dimension 0,1,2... ", us::Any()); +} + +void +mitk::GIFCooccurenceMatrix2::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())) + { + int direction = 0; + if (parsedArgs.count(name + "::direction")) + { + direction = SplitDouble(parsedArgs[name + "::direction"].ToString(), ';')[0]; + } + std::vector ranges; + if (parsedArgs.count(name + "::range")) + { + ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); + } + else + { + ranges.push_back(1); + } + + for (std::size_t i = 0; i < ranges.size(); ++i) + { + MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; + this->SetRange(ranges[i]); + this->SetDirection(direction); + auto localResults = this->CalculateFeatures(feature, maskNoNAN); + featureList.insert(featureList.end(), localResults.begin(), localResults.end()); + MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; + } + } + +} + diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp index 9a7b7f0999..07febc61fd 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp @@ -1,214 +1,238 @@ /*=================================================================== 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 // ITK #include #include // STL #include template void CalculateFirstOrderStatistics(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFFirstOrderStatistics::FeatureListType & featureList, mitk::GIFFirstOrderStatistics::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::LabelStatisticsImageFilter FilterType; typedef typename FilterType::HistogramType HistogramType; typedef typename HistogramType::IndexType HIndexType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); double imageRange = minMaxComputer->GetMaximum() - minMaxComputer->GetMinimum(); typename FilterType::Pointer labelStatisticsImageFilter = FilterType::New(); labelStatisticsImageFilter->SetInput( itkImage ); labelStatisticsImageFilter->SetLabelInput(maskImage); labelStatisticsImageFilter->SetUseHistograms(true); if (params.m_UseCtRange) { labelStatisticsImageFilter->SetHistogramParameters(1024.5+3096.5, -1024.5,3096.5); } else { labelStatisticsImageFilter->SetHistogramParameters(params.m_HistogramSize, minMaxComputer->GetMinimum(),minMaxComputer->GetMaximum()); } labelStatisticsImageFilter->Update(); // --------------- Range -------------------- double range = labelStatisticsImageFilter->GetMaximum(1) - labelStatisticsImageFilter->GetMinimum(1); // --------------- Uniformity, Entropy -------------------- double count = labelStatisticsImageFilter->GetCount(1); //double std_dev = labelStatisticsImageFilter->GetSigma(1); double uncorrected_std_dev = std::sqrt((count - 1) / count * labelStatisticsImageFilter->GetVariance(1)); double mean = labelStatisticsImageFilter->GetMean(1); auto histogram = labelStatisticsImageFilter->GetHistogram(1); HIndexType index; index.SetSize(1); double binWidth = histogram->GetBinMax(0, 0) - histogram->GetBinMin(0, 0); double uniformity = 0; double entropy = 0; double squared_sum = 0; double kurtosis = 0; double mean_absolut_deviation = 0; double skewness = 0; double sum_prob = 0; double p10th = histogram->Quantile(0,0.10); double p25th = histogram->Quantile(0,0.25); double p75th = histogram->Quantile(0,0.75); double p90th = histogram->Quantile(0,0.90); double Log2=log(2); for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; double prob = histogram->GetFrequency(index); if (prob < 0.1) continue; double voxelValue = histogram->GetBinMin(0, i) +binWidth * 0.5; sum_prob += prob; squared_sum += prob * voxelValue*voxelValue; prob /= count; mean_absolut_deviation += prob* std::abs(voxelValue - mean); kurtosis +=prob* (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean); skewness += prob* (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean); uniformity += prob*prob; if (prob > 0) { entropy += prob * std::log(prob) / Log2; } } entropy = -entropy; double rms = std::sqrt(squared_sum / count); kurtosis = kurtosis / (uncorrected_std_dev*uncorrected_std_dev * uncorrected_std_dev*uncorrected_std_dev); skewness = skewness / (uncorrected_std_dev*uncorrected_std_dev * uncorrected_std_dev); //mean_absolut_deviation = mean_absolut_deviation; double coveredGrayValueRange = range / imageRange; featureList.push_back(std::make_pair("FirstOrder Range",range)); featureList.push_back(std::make_pair("FirstOrder Uniformity",uniformity)); featureList.push_back(std::make_pair("FirstOrder Entropy",entropy)); featureList.push_back(std::make_pair("FirstOrder Energy",squared_sum)); featureList.push_back(std::make_pair("FirstOrder RMS",rms)); featureList.push_back(std::make_pair("FirstOrder Kurtosis", kurtosis)); featureList.push_back(std::make_pair("FirstOrder Excess Kurtosis", kurtosis-3)); featureList.push_back(std::make_pair("FirstOrder Skewness",skewness)); featureList.push_back(std::make_pair("FirstOrder Mean absolute deviation",mean_absolut_deviation)); featureList.push_back(std::make_pair("FirstOrder Covered Image Intensity Range",coveredGrayValueRange)); featureList.push_back(std::make_pair("FirstOrder Minimum",labelStatisticsImageFilter->GetMinimum(1))); featureList.push_back(std::make_pair("FirstOrder Maximum",labelStatisticsImageFilter->GetMaximum(1))); featureList.push_back(std::make_pair("FirstOrder Mean",labelStatisticsImageFilter->GetMean(1))); featureList.push_back(std::make_pair("FirstOrder Variance",labelStatisticsImageFilter->GetVariance(1))); featureList.push_back(std::make_pair("FirstOrder Sum",labelStatisticsImageFilter->GetSum(1))); featureList.push_back(std::make_pair("FirstOrder Median",labelStatisticsImageFilter->GetMedian(1))); featureList.push_back(std::make_pair("FirstOrder Standard deviation",labelStatisticsImageFilter->GetSigma(1))); featureList.push_back(std::make_pair("FirstOrder No. of Voxel",labelStatisticsImageFilter->GetCount(1))); featureList.push_back(std::make_pair("FirstOrder 10th Percentile",p10th)); featureList.push_back(std::make_pair("FirstOrder 90th Percentile",p90th)); featureList.push_back(std::make_pair("FirstOrder Interquartile Range",(p75th - p25th))); //Calculate the robus mean absolute deviation //First, set all frequencies to 0 that are <10th or >90th percentile for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; if (histogram->GetBinMax(0, i) < p10th) { histogram->SetFrequencyOfIndex(index, 0); } else if (histogram->GetBinMin(0, i) > p90th) { histogram->SetFrequencyOfIndex(index, 0); } } //Calculate the mean double meanRobust = 0.0; for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; meanRobust += histogram->GetFrequency(index) * 0.5 * (histogram->GetBinMin(0,i) + histogram->GetBinMax(0,i)); } meanRobust = meanRobust / histogram->GetTotalFrequency(); double robustMeanAbsoluteDeviation = 0.0; for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; robustMeanAbsoluteDeviation += std::abs(histogram->GetFrequency(index) * ( (0.5 * (histogram->GetBinMin(0,i) + histogram->GetBinMax(0,i))) - meanRobust )); } robustMeanAbsoluteDeviation = robustMeanAbsoluteDeviation / histogram->GetTotalFrequency(); featureList.push_back(std::make_pair("FirstOrder Robust Mean", meanRobust)); featureList.push_back(std::make_pair("FirstOrder Robust Mean Absolute Deviation",robustMeanAbsoluteDeviation)); } mitk::GIFFirstOrderStatistics::GIFFirstOrderStatistics() : m_HistogramSize(256), m_UseCtRange(false) { + SetShortName("fo"); + SetLongName("first-order"); } mitk::GIFFirstOrderStatistics::FeatureListType mitk::GIFFirstOrderStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; ParameterStruct params; params.m_HistogramSize = this->m_HistogramSize; params.m_UseCtRange = this->m_UseCtRange; AccessByItk_3(image, CalculateFirstOrderStatistics, mask, featureList, params); return featureList; } mitk::GIFFirstOrderStatistics::FeatureNameListType mitk::GIFFirstOrderStatistics::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("FirstOrder Minimum"); featureList.push_back("FirstOrder Maximum"); featureList.push_back("FirstOrder Mean"); featureList.push_back("FirstOrder Variance"); featureList.push_back("FirstOrder Sum"); featureList.push_back("FirstOrder Median"); featureList.push_back("FirstOrder Standard deviation"); featureList.push_back("FirstOrder No. of Voxel"); return featureList; } + + +void mitk::GIFFirstOrderStatistics::AddArguments(mitkCommandLineParser &parser) +{ + std::string name = GetOptionPrefix(); + + parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Volume-Statistic", "calculates volume based features", us::Any()); +} + +void +mitk::GIFFirstOrderStatistics::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList) +{ + auto parsedArgs = GetParameter(); + if (parsedArgs.count(GetLongName())) + { + MITK_INFO << "Start calculating first order features ...."; + auto localResults = this->CalculateFeatures(feature, mask); + featureList.insert(featureList.end(), localResults.begin(), localResults.end()); + MITK_INFO << "Finished calculating first order features...."; + } +} + diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp index cc5b33b65f..35feb2dd63 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp @@ -1,248 +1,296 @@ /*=================================================================== 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 // ITK #include #include // STL #include template void CalculateGrayLevelRunLengthFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFGrayLevelRunLength::FeatureListType & featureList, mitk::GIFGrayLevelRunLength::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::Statistics::EnhancedScalarImageToRunLengthFeaturesFilter FilterType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef typename FilterType::RunLengthFeaturesFilterType TextureFilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename FilterType::Pointer filter = FilterType::New(); typename FilterType::OffsetVector::Pointer newOffset = FilterType::OffsetVector::New(); auto oldOffsets = filter->GetOffsets(); auto oldOffsetsIterator = oldOffsets->Begin(); while (oldOffsetsIterator != oldOffsets->End()) { bool continueOuterLoop = false; typename FilterType::OffsetType offset = oldOffsetsIterator->Value(); for (unsigned int i = 0; i < VImageDimension; ++i) { if (params.m_Direction == i + 2 && offset[i] != 0) { continueOuterLoop = true; } } if (params.m_Direction == 1) { offset[0] = 0; offset[1] = 0; offset[2] = 1; newOffset->push_back(offset); break; } oldOffsetsIterator++; if (continueOuterLoop) continue; newOffset->push_back(offset); } filter->SetOffsets(newOffset); // All features are required typename FilterType::FeatureNameVectorPointer requestedFeatures = FilterType::FeatureNameVector::New(); requestedFeatures->push_back(TextureFilterType::ShortRunEmphasis); requestedFeatures->push_back(TextureFilterType::LongRunEmphasis); requestedFeatures->push_back(TextureFilterType::GreyLevelNonuniformity); requestedFeatures->push_back(TextureFilterType::GreyLevelNonuniformityNormalized); requestedFeatures->push_back(TextureFilterType::RunLengthNonuniformity); requestedFeatures->push_back(TextureFilterType::RunLengthNonuniformityNormalized); requestedFeatures->push_back(TextureFilterType::LowGreyLevelRunEmphasis); requestedFeatures->push_back(TextureFilterType::HighGreyLevelRunEmphasis); requestedFeatures->push_back(TextureFilterType::ShortRunLowGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::ShortRunHighGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::LongRunLowGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::LongRunHighGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::RunPercentage); requestedFeatures->push_back(TextureFilterType::NumberOfRuns); requestedFeatures->push_back(TextureFilterType::GreyLevelVariance); requestedFeatures->push_back(TextureFilterType::RunLengthVariance); requestedFeatures->push_back(TextureFilterType::RunEntropy); typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); filter->SetInput(itkImage); filter->SetMaskImage(maskImage); filter->SetRequestedFeatures(requestedFeatures); int rangeOfPixels = params.m_Range; if (rangeOfPixels < 2) rangeOfPixels = 256; if (params.m_UseCtRange) { filter->SetPixelValueMinMax((TPixel)(-1024.5),(TPixel)(3096.5)); filter->SetNumberOfBinsPerAxis(3096.5+1024.5); } else { MITK_INFO << "Min: " << minMaxComputer->GetMinimum() << " Max: " << minMaxComputer->GetMaximum() << " Range: " << rangeOfPixels; filter->SetPixelValueMinMax(minMaxComputer->GetMinimum(),minMaxComputer->GetMaximum()); filter->SetNumberOfBinsPerAxis(rangeOfPixels); } filter->SetDistanceValueMinMax(0,rangeOfPixels); filter->Update(); auto featureMeans = filter->GetFeatureMeans (); auto featureStd = filter->GetFeatureStandardDeviations(); std::ostringstream ss; ss << rangeOfPixels; std::string strRange = ss.str(); for (std::size_t i = 0; i < featureMeans->size(); ++i) { switch (i) { case TextureFilterType::ShortRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::LongRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::GreyLevelNonuniformity : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelNonuniformity Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelNonuniformity Std.",featureStd->ElementAt(i))); break; case TextureFilterType::GreyLevelNonuniformityNormalized : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelNonuniformityNormalized Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelNonuniformityNormalized Std.",featureStd->ElementAt(i))); break; case TextureFilterType::RunLengthNonuniformity : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthNonuniformity Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthNonuniformity Std.",featureStd->ElementAt(i))); break; case TextureFilterType::RunLengthNonuniformityNormalized : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthNonuniformityNormalized Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthNonuniformityNormalized Std.",featureStd->ElementAt(i))); break; case TextureFilterType::LowGreyLevelRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LowGreyLevelRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LowGreyLevelRunEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::HighGreyLevelRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") HighGreyLevelRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") HighGreyLevelRunEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::ShortRunLowGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunLowGreyLevelEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::ShortRunHighGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunHighGreyLevelEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::LongRunLowGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunLowGreyLevelEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::LongRunHighGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunHighGreyLevelEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::RunPercentage : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunPercentage Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunPercentage Std.",featureStd->ElementAt(i))); break; case TextureFilterType::NumberOfRuns : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") NumberOfRuns Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") NumberOfRuns Std.",featureStd->ElementAt(i))); break; case TextureFilterType::GreyLevelVariance : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelVariance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelVariance Std.",featureStd->ElementAt(i))); break; case TextureFilterType::RunLengthVariance : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthVariance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthVariance Std.",featureStd->ElementAt(i))); break; case TextureFilterType::RunEntropy : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunEntropy Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunEntropy Std.",featureStd->ElementAt(i))); break; default: break; } } } mitk::GIFGrayLevelRunLength::GIFGrayLevelRunLength(): m_Range(1.0), m_UseCtRange(false), m_Direction(0) { + SetShortName("rl"); + SetLongName("run-length"); } mitk::GIFGrayLevelRunLength::FeatureListType mitk::GIFGrayLevelRunLength::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; ParameterStruct params; params.m_UseCtRange=m_UseCtRange; params.m_Range = m_Range; params.m_Direction = m_Direction; AccessByItk_3(image, CalculateGrayLevelRunLengthFeatures, mask, featureList,params); return featureList; } mitk::GIFGrayLevelRunLength::FeatureNameListType mitk::GIFGrayLevelRunLength::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("RunLength. ShortRunEmphasis Means"); featureList.push_back("RunLength. ShortRunEmphasis Std."); featureList.push_back("RunLength. LongRunEmphasis Means"); featureList.push_back("RunLength. LongRunEmphasis Std."); featureList.push_back("RunLength. GreyLevelNonuniformity Means"); featureList.push_back("RunLength. GreyLevelNonuniformity Std."); featureList.push_back("RunLength. RunLengthNonuniformity Means"); featureList.push_back("RunLength. RunLengthNonuniformity Std."); featureList.push_back("RunLength. LowGreyLevelRunEmphasis Means"); featureList.push_back("RunLength. LowGreyLevelRunEmphasis Std."); featureList.push_back("RunLength. HighGreyLevelRunEmphasis Means"); featureList.push_back("RunLength. HighGreyLevelRunEmphasis Std."); featureList.push_back("RunLength. ShortRunLowGreyLevelEmphasis Means"); featureList.push_back("RunLength. ShortRunLowGreyLevelEmphasis Std."); featureList.push_back("RunLength. ShortRunHighGreyLevelEmphasis Means"); featureList.push_back("RunLength. ShortRunHighGreyLevelEmphasis Std."); featureList.push_back("RunLength. LongRunHighGreyLevelEmphasis Means"); featureList.push_back("RunLength. LongRunHighGreyLevelEmphasis Std."); return featureList; } + + +void mitk::GIFGrayLevelRunLength::AddArguments(mitkCommandLineParser &parser) +{ + std::string name = GetOptionPrefix(); + + parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features (new implementation)", us::Any()); + parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::String, "Cooc 2 Range", "Define the range that is used (Semicolon-separated)", us::Any()); + parser.addArgument(name + "::direction", name + "::dir", mitkCommandLineParser::String, "Int", "Allows to specify the direction for Cooc and RL. 0: All directions, 1: Only single direction (Test purpose), 2,3,4... without dimension 0,1,2... ", us::Any()); +} + +void +mitk::GIFGrayLevelRunLength::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())) + { + int direction = 0; + if (parsedArgs.count(name + "::direction")) + { + direction = SplitDouble(parsedArgs[name + "::direction"].ToString(), ';')[0]; + } + std::vector ranges; + if (parsedArgs.count(name + "::range")) + { + ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); + } + else + { + ranges.push_back(1); + } + + for (std::size_t i = 0; i < ranges.size(); ++i) + { + MITK_INFO << "Start calculating Run-length with range " << ranges[i] << "...."; + this->SetRange(ranges[i]); + this->SetDirection(direction); + auto localResults = this->CalculateFeatures(feature, maskNoNAN); + featureList.insert(featureList.end(), localResults.begin(), localResults.end()); + MITK_INFO << "Finished calculating Run-length with range " << ranges[i] << "...."; + } + } + +} diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp index d185e134cf..af3df166b1 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp @@ -1,226 +1,284 @@ /*=================================================================== 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 // ITK #include #include // STL #include template void CalculateGrayLevelSizeZoneFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFGrayLevelSizeZone::FeatureListType & featureList, mitk::GIFGrayLevelSizeZone::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::Statistics::EnhancedScalarImageToSizeZoneFeaturesFilter FilterType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef typename FilterType::SizeZoneFeaturesFilterType TextureFilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename FilterType::Pointer filter = FilterType::New(); typename FilterType::OffsetVector::Pointer newOffset = FilterType::OffsetVector::New(); auto oldOffsets = filter->GetOffsets(); auto oldOffsetsIterator = oldOffsets->Begin(); while (oldOffsetsIterator != oldOffsets->End()) { bool continueOuterLoop = false; typename FilterType::OffsetType offset = oldOffsetsIterator->Value(); for (unsigned int i = 0; i < VImageDimension; ++i) { if (params.m_Direction == i + 2 && offset[i] != 0) { continueOuterLoop = true; } } if (params.m_Direction == 1) { offset[0] = 0; offset[1] = 0; offset[2] = 1; newOffset->push_back(offset); break; } oldOffsetsIterator++; if (continueOuterLoop) continue; newOffset->push_back(offset); } filter->SetOffsets(newOffset); // All features are required typename FilterType::FeatureNameVectorPointer requestedFeatures = FilterType::FeatureNameVector::New(); requestedFeatures->push_back(TextureFilterType::SmallZoneEmphasis); requestedFeatures->push_back(TextureFilterType::LargeZoneEmphasis); requestedFeatures->push_back(TextureFilterType::GreyLevelNonuniformity); requestedFeatures->push_back(TextureFilterType::GreyLevelNonuniformityNormalized); requestedFeatures->push_back(TextureFilterType::SizeZoneNonuniformity); requestedFeatures->push_back(TextureFilterType::SizeZoneNonuniformityNormalized); requestedFeatures->push_back(TextureFilterType::LowGreyLevelZoneEmphasis); requestedFeatures->push_back(TextureFilterType::HighGreyLevelZoneEmphasis); requestedFeatures->push_back(TextureFilterType::SmallZoneLowGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::SmallZoneHighGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::LargeZoneLowGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::LargeZoneHighGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::ZonePercentage); requestedFeatures->push_back(TextureFilterType::GreyLevelVariance); requestedFeatures->push_back(TextureFilterType::SizeZoneVariance); requestedFeatures->push_back(TextureFilterType::ZoneEntropy); typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); filter->SetInput(itkImage); filter->SetMaskImage(maskImage); filter->SetRequestedFeatures(requestedFeatures); int rangeOfPixels = params.m_Range; if (rangeOfPixels < 2) rangeOfPixels = 256; if (params.m_UseCtRange) { filter->SetPixelValueMinMax((TPixel)(-1024.5),(TPixel)(3096.5)); filter->SetNumberOfBinsPerAxis(3096.5+1024.5); } else { filter->SetPixelValueMinMax(minMaxComputer->GetMinimum(),minMaxComputer->GetMaximum()); filter->SetNumberOfBinsPerAxis(rangeOfPixels); } filter->SetDistanceValueMinMax(0,rangeOfPixels); filter->Update(); auto featureMeans = filter->GetFeatureMeans (); auto featureStd = filter->GetFeatureStandardDeviations(); std::ostringstream ss; ss << rangeOfPixels; std::string strRange = ss.str(); for (std::size_t i = 0; i < featureMeans->size(); ++i) { switch (i) { case TextureFilterType::SmallZoneEmphasis : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SmallZoneEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::LargeZoneEmphasis : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") LargeZoneEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::GreyLevelNonuniformity : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") GreyLevelNonuniformity Means",featureMeans->ElementAt(i))); break; case TextureFilterType::GreyLevelNonuniformityNormalized : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") GreyLevelNonuniformityNormalized Means",featureMeans->ElementAt(i))); break; case TextureFilterType::SizeZoneNonuniformity : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SizeZoneNonuniformity Means",featureMeans->ElementAt(i))); break; case TextureFilterType::SizeZoneNonuniformityNormalized : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SizeZoneNonuniformityNormalized Means",featureMeans->ElementAt(i))); break; case TextureFilterType::LowGreyLevelZoneEmphasis : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") LowGreyLevelZoneEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::HighGreyLevelZoneEmphasis : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") HighGreyLevelZoneEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::SmallZoneLowGreyLevelEmphasis : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SmallZoneLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::SmallZoneHighGreyLevelEmphasis : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SmallZoneHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::LargeZoneLowGreyLevelEmphasis : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") LargeZoneLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::LargeZoneHighGreyLevelEmphasis : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") LargeZoneHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::ZonePercentage : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") ZonePercentage Means",featureMeans->ElementAt(i))); break; case TextureFilterType::GreyLevelVariance : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") GreyLevelVariance Means",featureMeans->ElementAt(i))); break; case TextureFilterType::SizeZoneVariance : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SizeZoneVariance Means",featureMeans->ElementAt(i))); break; case TextureFilterType::ZoneEntropy : featureList.push_back(std::make_pair("SizeZone ("+ strRange+") ZoneEntropy Means",featureMeans->ElementAt(i))); break; default: break; } } } mitk::GIFGrayLevelSizeZone::GIFGrayLevelSizeZone(): m_Range(1.0), m_UseCtRange(false), m_Direction(0) { + SetShortName("glsz"); + SetLongName("greylevelsizezone"); } mitk::GIFGrayLevelSizeZone::FeatureListType mitk::GIFGrayLevelSizeZone::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; ParameterStruct params; params.m_UseCtRange=m_UseCtRange; params.m_Range = m_Range; params.m_Direction = m_Direction; AccessByItk_3(image, CalculateGrayLevelSizeZoneFeatures, mask, featureList,params); return featureList; } mitk::GIFGrayLevelSizeZone::FeatureNameListType mitk::GIFGrayLevelSizeZone::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("SizeZone SmallZoneEmphasis Means"); featureList.push_back("SizeZone SmallZoneEmphasis Std."); featureList.push_back("SizeZone LargeZoneEmphasis Means"); featureList.push_back("SizeZone LargeZoneEmphasis Std."); featureList.push_back("SizeZone GreyLevelNonuniformity Means"); featureList.push_back("SizeZone GreyLevelNonuniformity Std."); featureList.push_back("SizeZone SizeZoneNonuniformity Means"); featureList.push_back("SizeZone SizeZoneNonuniformity Std."); featureList.push_back("SizeZone LowGreyLevelZoneEmphasis Means"); featureList.push_back("SizeZone LowGreyLevelZoneEmphasis Std."); featureList.push_back("SizeZone HighGreyLevelZoneEmphasis Means"); featureList.push_back("SizeZone HighGreyLevelZoneEmphasis Std."); featureList.push_back("SizeZone SmallZoneLowGreyLevelEmphasis Means"); featureList.push_back("SizeZone SmallZoneLowGreyLevelEmphasis Std."); featureList.push_back("SizeZone SmallZoneHighGreyLevelEmphasis Means"); featureList.push_back("SizeZone SmallZoneHighGreyLevelEmphasis Std."); featureList.push_back("SizeZone LargeZoneHighGreyLevelEmphasis Means"); featureList.push_back("SizeZone LargeZoneHighGreyLevelEmphasis Std."); return featureList; } + + + +void mitk::GIFGrayLevelSizeZone::AddArguments(mitkCommandLineParser &parser) +{ + std::string name = GetOptionPrefix(); + + parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features (new implementation)", us::Any()); + parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::String, "Cooc 2 Range", "Define the range that is used (Semicolon-separated)", us::Any()); + parser.addArgument(name + "::direction", name + "::dir", mitkCommandLineParser::String, "Int", "Allows to specify the direction for Cooc and RL. 0: All directions, 1: Only single direction (Test purpose), 2,3,4... without dimension 0,1,2... ", us::Any()); + parser.addArgument(name + "::ctrange", name + "::ctrange", mitkCommandLineParser::String, "Int", "Turn on if CT images are used", us::Any()); +} + +void +mitk::GIFGrayLevelSizeZone::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())) + { + int direction = 0; + if (parsedArgs.count(name + "::direction")) + { + direction = SplitDouble(parsedArgs[name + "::direction"].ToString(), ';')[0]; + } + std::vector ranges; + if (parsedArgs.count(name + "::range")) + { + ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); + } + else + { + ranges.push_back(1); + } + + if (parsedArgs.count(name + "::ctrange")) + { + SetUseCtRange(true); + } + + + for (std::size_t i = 0; i < ranges.size(); ++i) + { + MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; + this->SetRange(ranges[i]); + this->SetDirection(direction); + auto localResults = this->CalculateFeatures(feature, maskNoNAN); + featureList.insert(featureList.end(), localResults.begin(), localResults.end()); + MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; + } + } + +} + + diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp index 2f7511c0af..c3d887bfc9 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp @@ -1,174 +1,224 @@ /*=================================================================== 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 // ITK #include #include // STL #include template void CalculateGrayLevelNeighbourhoodGreyLevelDifferenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFNeighbourhoodGreyLevelDifference::FeatureListType & featureList, mitk::GIFNeighbourhoodGreyLevelDifference::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::Statistics::EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter FilterType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef typename FilterType::NeighbourhoodGreyLevelDifferenceFeaturesFilterType TextureFilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename FilterType::Pointer filter = FilterType::New(); typename FilterType::OffsetVector::Pointer newOffset = FilterType::OffsetVector::New(); auto oldOffsets = filter->GetOffsets(); auto oldOffsetsIterator = oldOffsets->Begin(); while (oldOffsetsIterator != oldOffsets->End()) { bool continueOuterLoop = false; typename FilterType::OffsetType offset = oldOffsetsIterator->Value(); for (unsigned int i = 0; i < VImageDimension; ++i) { if (params.m_Direction == i + 2 && offset[i] != 0) { continueOuterLoop = true; } } if (params.m_Direction == 1) { offset[0] = 0; offset[1] = 0; offset[2] = 1; newOffset->push_back(offset); break; } oldOffsetsIterator++; if (continueOuterLoop) continue; newOffset->push_back(offset); } filter->SetOffsets(newOffset); // All features are required typename FilterType::FeatureNameVectorPointer requestedFeatures = FilterType::FeatureNameVector::New(); requestedFeatures->push_back(TextureFilterType::Coarseness); requestedFeatures->push_back(TextureFilterType::Contrast); requestedFeatures->push_back(TextureFilterType::Busyness); requestedFeatures->push_back(TextureFilterType::Complexity); requestedFeatures->push_back(TextureFilterType::Strength); typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); filter->SetInput(itkImage); filter->SetMaskImage(maskImage); filter->SetRequestedFeatures(requestedFeatures); int rangeOfPixels = params.m_Range; if (rangeOfPixels < 2) rangeOfPixels = 256; if (params.m_UseCtRange) { filter->SetPixelValueMinMax((TPixel)(-1024.5),(TPixel)(3096.5)); filter->SetNumberOfBinsPerAxis(3096.5+1024.5); } else { filter->SetPixelValueMinMax(minMaxComputer->GetMinimum(),minMaxComputer->GetMaximum()); filter->SetNumberOfBinsPerAxis(rangeOfPixels); } filter->SetDistanceValueMinMax(0,rangeOfPixels); filter->Update(); auto featureMeans = filter->GetFeatureMeans (); auto featureStd = filter->GetFeatureStandardDeviations(); std::ostringstream ss; ss << rangeOfPixels; std::string strRange = ss.str(); for (std::size_t i = 0; i < featureMeans->size(); ++i) { switch (i) { case TextureFilterType::Coarseness : featureList.push_back(std::make_pair("NeighbourhoodGreyLevelDifference ("+ strRange+") Coarseness Means",featureMeans->ElementAt(i))); break; case TextureFilterType::Contrast : featureList.push_back(std::make_pair("NeighbourhoodGreyLevelDifference ("+ strRange+") Contrast Means",featureMeans->ElementAt(i))); break; case TextureFilterType::Busyness : featureList.push_back(std::make_pair("NeighbourhoodGreyLevelDifference ("+ strRange+") Busyness Means",featureMeans->ElementAt(i))); break; case TextureFilterType::Complexity : featureList.push_back(std::make_pair("NeighbourhoodGreyLevelDifference ("+ strRange+") Complexity Means",featureMeans->ElementAt(i))); break; case TextureFilterType::Strength : featureList.push_back(std::make_pair("NeighbourhoodGreyLevelDifference ("+ strRange+") Strength Means",featureMeans->ElementAt(i))); break; default: break; } } } mitk::GIFNeighbourhoodGreyLevelDifference::GIFNeighbourhoodGreyLevelDifference(): m_Range(1.0), m_UseCtRange(false), m_Direction(0) { + SetShortName("ngld"); + SetLongName("NeighbourhoodGreyLevelDifference"); } mitk::GIFNeighbourhoodGreyLevelDifference::FeatureListType mitk::GIFNeighbourhoodGreyLevelDifference::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; ParameterStruct params; params.m_UseCtRange=m_UseCtRange; params.m_Range = m_Range; params.m_Direction = m_Direction; AccessByItk_3(image, CalculateGrayLevelNeighbourhoodGreyLevelDifferenceFeatures, mask, featureList,params); return featureList; } mitk::GIFNeighbourhoodGreyLevelDifference::FeatureNameListType mitk::GIFNeighbourhoodGreyLevelDifference::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("NeighbourhoodGreyLevelDifference. Coarseness Means"); featureList.push_back("NeighbourhoodGreyLevelDifference. Coarseness Std."); featureList.push_back("NeighbourhoodGreyLevelDifference. Contrast Means"); featureList.push_back("NeighbourhoodGreyLevelDifference. Contrast Std."); featureList.push_back("NeighbourhoodGreyLevelDifference. Busyness Means"); featureList.push_back("NeighbourhoodGreyLevelDifference. Busyness Std."); featureList.push_back("NeighbourhoodGreyLevelDifference. Complexity Means"); featureList.push_back("NeighbourhoodGreyLevelDifference. Complexity Std."); featureList.push_back("NeighbourhoodGreyLevelDifference. Strength Means"); featureList.push_back("NeighbourhoodGreyLevelDifference. Strength Std."); return featureList; } + + + +void mitk::GIFNeighbourhoodGreyLevelDifference::AddArguments(mitkCommandLineParser &parser) +{ + std::string name = GetOptionPrefix(); + + parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features (new implementation)", us::Any()); + parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::String, "Cooc 2 Range", "Define the range that is used (Semicolon-separated)", us::Any()); + parser.addArgument(name + "::direction", name + "::dir", mitkCommandLineParser::String, "Int", "Allows to specify the direction for Cooc and RL. 0: All directions, 1: Only single direction (Test purpose), 2,3,4... without dimension 0,1,2... ", us::Any()); +} + +void +mitk::GIFNeighbourhoodGreyLevelDifference::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())) + { + int direction = 0; + if (parsedArgs.count(name + "::direction")) + { + direction = SplitDouble(parsedArgs[name + "::direction"].ToString(), ';')[0]; + } + std::vector ranges; + if (parsedArgs.count(name + "::range")) + { + ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); + } + else + { + ranges.push_back(1); + } + + for (std::size_t i = 0; i < ranges.size(); ++i) + { + MITK_INFO << "Start calculating Neighbourhood Grey Level Difference with range " << ranges[i] << "...."; + this->SetRange(ranges[i]); + this->SetDirection(direction); + auto localResults = this->CalculateFeatures(feature, maskNoNAN); + featureList.insert(featureList.end(), localResults.begin(), localResults.end()); + MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; + } + } + +} + diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp new file mode 100644 index 0000000000..6cd3977cef --- /dev/null +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp @@ -0,0 +1,394 @@ +#include + +// MITK +#include +#include +#include + +// ITK +#include +#include +#include +#include + +// STL +#include + +static +void MatrixFeaturesTo(mitk::NGLDMMatrixFeatures features, + std::string prefix, + mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType &featureList); +static +void CalculateMeanAndStdDevFeatures(std::vector featureList, + mitk::NGLDMMatrixFeatures &meanFeature, + mitk::NGLDMMatrixFeatures &stdFeature); +static +void NormalizeMatrixFeature(mitk::NGLDMMatrixFeatures &features, + std::size_t number); + + + + +mitk::NGLDMMatrixHolder::NGLDMMatrixHolder(double min, double max, int number, int depenence) : + m_MinimumRange(min), + m_MaximumRange(max), + m_NumberOfBins(number), + m_NumberOfDependences(depenence) +{ + m_Matrix.resize(number, depenence); + m_Matrix.fill(0); + m_Stepsize = (max - min) / (number); +} + +int mitk::NGLDMMatrixHolder::IntensityToIndex(double intensity) +{ + return std::floor((intensity - m_MinimumRange) / m_Stepsize); +} + +double mitk::NGLDMMatrixHolder::IndexToMinIntensity(int index) +{ + return m_MinimumRange + index * m_Stepsize; +} +double mitk::NGLDMMatrixHolder::IndexToMeanIntensity(int index) +{ + return m_MinimumRange + (index+0.5) * m_Stepsize; +} +double mitk::NGLDMMatrixHolder::IndexToMaxIntensity(int index) +{ + return m_MinimumRange + (index + 1) * m_Stepsize; +} + +template +void +CalculateNGLDMMatrix(itk::Image* itkImage, + itk::Image* mask, + int alpha, + int range, + mitk::NGLDMMatrixHolder &holder) +{ + typedef itk::Image ImageType; + typedef itk::Image MaskImageType; + typedef itk::NeighborhoodIterator ShapeIterType; + typedef itk::NeighborhoodIterator ShapeMaskIterType; + + holder.m_NumberOfCompleteNeighbourhoods = 0; + holder.m_NumberOfNeighbourhoods = 0; + holder.m_NumberOfNeighbourVoxels = 0; + holder.m_NumberOfDependenceNeighbourVoxels = 0; + + itk::Size radius; + radius.Fill(range); + ShapeIterType imageIter(radius, itkImage, itkImage->GetLargestPossibleRegion()); + ShapeMaskIterType maskIter(radius, mask, mask->GetLargestPossibleRegion()); + + auto region = mask->GetLargestPossibleRegion(); + + auto center = imageIter.Size() / 2; + auto iterSize = imageIter.Size(); + holder.m_NeighbourhoodSize = iterSize-1; + while (!maskIter.IsAtEnd()) + { + int sameValues = 0; + bool completeNeighbourhood = true; + + int i = holder.IntensityToIndex(imageIter.GetCenterPixel()); + + if (imageIter.GetCenterPixel() != imageIter.GetCenterPixel()) + { + ++imageIter; + ++maskIter; + continue; + } + + for (int position = 0; position < iterSize; ++position) + { + if ((position == center) || + ( ! region.IsInside(maskIter.GetIndex(position)))) + { + completeNeighbourhood = false; + continue; + } + bool isInBounds; + auto jIntensity = imageIter.GetPixel(position, isInBounds); + auto jMask = maskIter.GetPixel(position, isInBounds); + if (jMask < 0 || (jIntensity != jIntensity) || ( ! isInBounds)) + { + completeNeighbourhood = false; + continue; + } + + int j = holder.IntensityToIndex(jIntensity); + holder.m_NumberOfNeighbourVoxels += 1; + if (std::abs(i - j) <= alpha) + { + holder.m_NumberOfDependenceNeighbourVoxels += 1; + ++sameValues; + } + } + holder.m_Matrix(i, sameValues) += 1; + holder.m_NumberOfNeighbourhoods += 1; + if (completeNeighbourhood) + { + holder.m_NumberOfCompleteNeighbourhoods += 1; + } + + ++imageIter; + ++maskIter; + } + +} + +void CalculateFeatures( + mitk::NGLDMMatrixHolder &holder, + mitk::NGLDMMatrixFeatures & results + ) +{ + auto sijMatrix = holder.m_Matrix; + auto piMatrix = holder.m_Matrix; + auto pjMatrix = holder.m_Matrix; + double Ng = holder.m_NumberOfBins; + int NgSize = holder.m_NumberOfBins; + double Ns = sijMatrix.sum(); + piMatrix.rowwise().normalize(); + pjMatrix.colwise().normalize(); + + for (int i = 0; i < holder.m_NumberOfBins; ++i) + { + double sj = 0; + for (int j = 0; j < holder.m_NumberOfDependences; ++j) + { + double iInt = holder.IndexToMeanIntensity(i); + double sij = sijMatrix(i, j); + double k = j+1; + double pij = sij / Ns; + + results.LowDependenceEmphasis += sij / k / k; + results.HighDependenceEmphasis += sij * k*k; + if (iInt != 0) + { + results.LowGreyLevelCountEmphasis += sij / iInt / iInt; + } + results.HighGreyLevelCountEmphasis += sij * iInt*iInt; + if (iInt != 0) + { + results.LowDependenceLowGreyLevelEmphasis += sij / k / k / iInt / iInt; + } + results.LowDependenceHighGreyLevelEmphasis += sij * iInt*iInt / k / k; + if (iInt != 0) + { + results.HighDependenceLowGreyLevelEmphasis += sij *k * k / iInt / iInt; + } + results.HighDependenceHighGreyLevelEmphasis += sij * k*k*iInt*iInt; + + + results.MeanGreyLevelCount += iInt * pij; + results.MeanDependenceCount += k * pij; + if (pij > 0) + { + results.DependenceCountEntropy -= pij * std::log(pij) / std::log(2); + } + sj += sij; + } + results.GreyLevelNonUniformity += sj*sj; + results.GreyLevelNonUniformityNormalised += sj*sj; + } + + for (int j = 0; j < holder.m_NumberOfDependences; ++j) + { + double si = 0; + for (int i = 0; i < holder.m_NumberOfBins; ++i) + { + double sij = sijMatrix(i, j); + si += sij; + } + results.DependenceCountNonUniformity += si*si; + results.DependenceCountNonUniformityNormalised += si*si; + } + for (int i = 0; i < holder.m_NumberOfBins; ++i) + { + for (int j = 0; j < holder.m_NumberOfDependences; ++j) + { + double iInt = holder.IndexToMeanIntensity(i); + double sij = sijMatrix(i, j); + double k = j + 1; + double pij = sij / Ns; + + results.GreyLevelVariance += (iInt - results.MeanGreyLevelCount)* (iInt - results.MeanGreyLevelCount) * pij; + results.DependenceCountVariance += (k - results.MeanDependenceCount)* (k - results.MeanDependenceCount) * pij; + } + } + results.LowDependenceEmphasis /= Ns; + results.HighDependenceEmphasis /= Ns; + results.LowGreyLevelCountEmphasis /= Ns; + results.HighGreyLevelCountEmphasis /= Ns; + results.LowDependenceLowGreyLevelEmphasis /= Ns; + results.LowDependenceHighGreyLevelEmphasis /= Ns; + results.HighDependenceLowGreyLevelEmphasis /= Ns; + results.HighDependenceHighGreyLevelEmphasis /= Ns; + + results.GreyLevelNonUniformity /= Ns; + results.GreyLevelNonUniformityNormalised /= (Ns*Ns); + results.DependenceCountNonUniformity /= Ns; + results.DependenceCountNonUniformityNormalised /= (Ns*Ns); + + results.DependenceCountPercentage = 1; + + results.ExpectedNeighbourhoodSize = holder.m_NeighbourhoodSize; + results.AverageNeighbourhoodSize = holder.m_NumberOfNeighbourVoxels / (1.0 * holder.m_NumberOfNeighbourhoods); + results.AverageIncompleteNeighbourhoodSize = (holder.m_NumberOfNeighbourVoxels - holder.m_NumberOfCompleteNeighbourhoods* holder.m_NeighbourhoodSize) / (1.0 * (holder.m_NumberOfNeighbourhoods - holder.m_NumberOfCompleteNeighbourhoods)); + results.PercentageOfCompleteNeighbourhoods = holder.m_NumberOfCompleteNeighbourhoods / (1.0 * holder.m_NumberOfNeighbourhoods); + results.PercentageOfDependenceNeighbours = holder.m_NumberOfDependenceNeighbourVoxels / (1.0 * holder.m_NumberOfNeighbourVoxels); +} + +template +void +CalculateCoocurenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType & featureList, mitk::GIFNeighbouringGreyLevelDependenceFeature::GIFNeighbouringGreyLevelDependenceFeatureConfiguration config) +{ + typedef itk::Image ImageType; + typedef itk::Image MaskType; + typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; + typedef itk::Neighborhood NeighborhoodType; + typedef itk::Offset OffsetType; + + + /////////////////////////////////////////////////////////////////////////////////////////////// + + + typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); + minMaxComputer->SetImage(itkImage); + minMaxComputer->Compute(); + + double rangeMin = -0.5 + minMaxComputer->GetMinimum(); + double rangeMax = 0.5 + minMaxComputer->GetMaximum(); + + // Define Range + int numberOfBins = 6; + + typename MaskType::Pointer maskImage = MaskType::New(); + mitk::CastToItkImage(mask, maskImage); + + std::vector resultVector; + mitk::NGLDMMatrixHolder holderOverall(rangeMin, rangeMax, numberOfBins,37); + mitk::NGLDMMatrixFeatures overallFeature; + CalculateNGLDMMatrix(itkImage, maskImage, config.alpha, config.range, holderOverall); + CalculateFeatures(holderOverall, overallFeature); + //NormalizeMatrixFeature(overallFeature, offsetVector.size()); + + std::ostringstream ss; + ss << config.range; + std::string strRange = ss.str(); + + MatrixFeaturesTo(overallFeature, "NGLDM (" + strRange + ") overall", featureList); +} + + +static +void MatrixFeaturesTo(mitk::NGLDMMatrixFeatures features, + std::string prefix, + mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType &featureList) +{ + featureList.push_back(std::make_pair(prefix + " Low Dependence Emphasis", features.LowDependenceEmphasis)); + featureList.push_back(std::make_pair(prefix + " High Dependence Emphasis", features.HighDependenceEmphasis)); + featureList.push_back(std::make_pair(prefix + " Low Grey Level Count Emphasis", features.LowGreyLevelCountEmphasis)); + featureList.push_back(std::make_pair(prefix + " High Grey Level Count Emphasis", features.HighGreyLevelCountEmphasis)); + featureList.push_back(std::make_pair(prefix + " Low Dependence Low Grey Level Emphasis", features.LowDependenceLowGreyLevelEmphasis)); + featureList.push_back(std::make_pair(prefix + " Low Dependence High Grey Level Emphasis", features.LowDependenceHighGreyLevelEmphasis)); + featureList.push_back(std::make_pair(prefix + " High Dependence Low Grey Level Emphasis", features.HighDependenceLowGreyLevelEmphasis)); + featureList.push_back(std::make_pair(prefix + " High Dependence Low Grey Level Emphasis", features.HighDependenceHighGreyLevelEmphasis)); + + featureList.push_back(std::make_pair(prefix + " Grey Level Non-Uniformity", features.GreyLevelNonUniformity)); + featureList.push_back(std::make_pair(prefix + " Grey Level Non-Uniformity Normalised", features.GreyLevelNonUniformityNormalised)); + featureList.push_back(std::make_pair(prefix + " Dependence Count Non-Uniformity", features.DependenceCountNonUniformity)); + featureList.push_back(std::make_pair(prefix + " Dependence Count Non-Uniformtiy Normalised", features.DependenceCountNonUniformityNormalised)); + + featureList.push_back(std::make_pair(prefix + " Dependence Count Percentage", features.DependenceCountPercentage)); + featureList.push_back(std::make_pair(prefix + " Grey Level Mean", features.MeanGreyLevelCount)); + featureList.push_back(std::make_pair(prefix + " Grey Level Variance", features.GreyLevelVariance)); + featureList.push_back(std::make_pair(prefix + " Dependence Count Mean", features.MeanDependenceCount)); + featureList.push_back(std::make_pair(prefix + " Dependence Count Variance", features.DependenceCountVariance)); + featureList.push_back(std::make_pair(prefix + " Dependence Count Entropy", features.DependenceCountEntropy)); + + featureList.push_back(std::make_pair(prefix + " Expected Neighbourhood Size", features.ExpectedNeighbourhoodSize)); + featureList.push_back(std::make_pair(prefix + " Average Neighbourhood Size", features.AverageNeighbourhoodSize)); + featureList.push_back(std::make_pair(prefix + " Average Incomplete Neighbourhood Size", features.AverageIncompleteNeighbourhoodSize)); + featureList.push_back(std::make_pair(prefix + " Percentage of complete Neighbourhoods", features.PercentageOfCompleteNeighbourhoods)); + featureList.push_back(std::make_pair(prefix + " Percentage of Dependence Neighbour Voxels", features.PercentageOfDependenceNeighbours)); + +} + +mitk::GIFNeighbouringGreyLevelDependenceFeature::GIFNeighbouringGreyLevelDependenceFeature(): +m_Range(1.0), m_Direction(0) +{ + SetShortName("ngldm"); + SetLongName("ngldm"); +} + +mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType mitk::GIFNeighbouringGreyLevelDependenceFeature::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) +{ + FeatureListType featureList; + + GIFNeighbouringGreyLevelDependenceFeatureConfiguration config; + config.direction = m_Direction; + config.range = m_Range; + config.alpha = 0; + + AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config); + + return featureList; +} + +mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureNameListType mitk::GIFNeighbouringGreyLevelDependenceFeature::GetFeatureNames() +{ + FeatureNameListType featureList; + featureList.push_back("co-occ. Energy Means"); + featureList.push_back("co-occ. Energy Std."); + + return featureList; +} + + + +void mitk::GIFNeighbouringGreyLevelDependenceFeature::AddArguments(mitkCommandLineParser &parser) +{ + std::string name = GetOptionPrefix(); + + parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Calculate Neighbouring grey level dependence based features", "Calculate Neighbouring grey level dependence based features", us::Any()); + parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::String, "Cooc 2 Range", "Define the range that is used (Semicolon-separated)", us::Any()); + parser.addArgument(name + "::direction", name + "::dir", mitkCommandLineParser::String, "Int", "Allows to specify the direction for Cooc and RL. 0: All directions, 1: Only single direction (Test purpose), 2,3,4... without dimension 0,1,2... ", us::Any()); +} + +void +mitk::GIFNeighbouringGreyLevelDependenceFeature::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())) + { + int direction = 0; + if (parsedArgs.count(name + "::direction")) + { + direction = SplitDouble(parsedArgs[name + "::direction"].ToString(), ';')[0]; + } + std::vector ranges; + if (parsedArgs.count(name + "::range")) + { + ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); + } + else + { + ranges.push_back(1); + } + + for (std::size_t i = 0; i < ranges.size(); ++i) + { + MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; + this->SetRange(ranges[i]); + this->SetDirection(direction); + auto localResults = this->CalculateFeatures(feature, maskNoNAN); + featureList.insert(featureList.end(), localResults.begin(), localResults.end()); + MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; + } + } +} + diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp index efa7a5bebb..0d10c146bc 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp @@ -1,402 +1,426 @@ /*=================================================================== 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 // VTK #include #include #include // STL #include #include // OpenCV #include template void CalculateVolumeStatistic(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFVolumetricStatistics::FeatureListType & featureList) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::LabelStatisticsImageFilter FilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename FilterType::Pointer labelStatisticsImageFilter = FilterType::New(); labelStatisticsImageFilter->SetInput( itkImage ); labelStatisticsImageFilter->SetLabelInput(maskImage); labelStatisticsImageFilter->Update(); double volume = labelStatisticsImageFilter->GetCount(1); double voxelVolume = 1; for (int i = 0; i < (int)(VImageDimension); ++i) { volume *= itkImage->GetSpacing()[i]; voxelVolume *= itkImage->GetSpacing()[i]; } featureList.push_back(std::make_pair("Volumetric Features Volume (pixel based)", volume)); featureList.push_back(std::make_pair("Volumetric Features Voxel Volume", voxelVolume)); } template void CalculateLargestDiameter(itk::Image* mask, mitk::Image::Pointer valueImage, mitk::GIFVolumetricStatistics::FeatureListType & featureList) { typedef itk::Image ValueImageType; typename ValueImageType::Pointer itkValueImage = ValueImageType::New(); mitk::CastToItkImage(valueImage, itkValueImage); typedef itk::Image ImageType; typedef typename ImageType::PointType PointType; typename ImageType::SizeType radius; for (int i=0; i < (int)VImageDimension; ++i) radius[i] = 1; itk::NeighborhoodIterator iterator(radius, mask, mask->GetRequestedRegion()); itk::NeighborhoodIterator valueIter(radius, itkValueImage, itkValueImage->GetRequestedRegion()); std::vector borderPoints; // // Calculate surface in different directions // double surface = 0; std::vector directionSurface; for (int i = 0; i < (int)(iterator.Size()); ++i) { auto offset = iterator.GetOffset(i); double deltaS = 1; int nonZeros = 0; for (int j = 0; j < VImageDimension; ++j) { if (offset[j] != 0 && nonZeros == 0) { for (int k = 0; k < VImageDimension; ++k) { if (k != j) deltaS *= mask->GetSpacing()[k]; } nonZeros++; } else if (offset[j] != 0) { deltaS = 0; } } if (nonZeros < 1) deltaS = 0; directionSurface.push_back(deltaS); } // // Prepare calulation of Centre of mass shift // PointType normalCenter(0); PointType normalCenterUncorrected(0); PointType weightedCenter(0); PointType currentPoint; int numberOfPoints = 0; int numberOfPointsUncorrected = 0; double sumOfPoints = 0; while(!iterator.IsAtEnd()) { if (iterator.GetCenterPixel() == 0) { ++iterator; ++valueIter; continue; } mask->TransformIndexToPhysicalPoint(iterator.GetIndex(), currentPoint); normalCenterUncorrected += currentPoint.GetVectorFromOrigin(); ++numberOfPointsUncorrected; double intensityValue = valueIter.GetCenterPixel(); if (intensityValue == intensityValue) { normalCenter += currentPoint.GetVectorFromOrigin(); weightedCenter += currentPoint.GetVectorFromOrigin() * intensityValue; sumOfPoints += intensityValue; ++numberOfPoints; } bool border = false; for (int i = 0; i < (int)(iterator.Size()); ++i) { if (iterator.GetPixel(i) == 0 || ( ! iterator.IndexInBounds(i))) { border = true; surface += directionSurface[i]; //break; } } if (border) { auto centerIndex = iterator.GetIndex(); PointType centerPoint; mask->TransformIndexToPhysicalPoint(centerIndex, centerPoint ); borderPoints.push_back(centerPoint); } ++iterator; ++valueIter; } auto normalCenterVector = normalCenter.GetVectorFromOrigin() / numberOfPoints; auto normalCenterVectorUncorrected = normalCenter.GetVectorFromOrigin() / numberOfPointsUncorrected; auto weightedCenterVector = weightedCenter.GetVectorFromOrigin() / sumOfPoints; auto differenceOfCentersUncorrected = (normalCenterVectorUncorrected - weightedCenterVector).GetNorm(); auto differenceOfCenters = (normalCenterVector - weightedCenterVector).GetNorm(); double longestDiameter = 0; unsigned long numberOfBorderPoints = borderPoints.size(); for (int i = 0; i < (int)numberOfBorderPoints; ++i) { auto point = borderPoints[i]; for (int j = i; j < (int)numberOfBorderPoints; ++j) { double newDiameter=point.EuclideanDistanceTo(borderPoints[j]); if (newDiameter > longestDiameter) longestDiameter = newDiameter; } } featureList.push_back(std::make_pair("Volumetric Features Maximum 3D diameter", longestDiameter)); featureList.push_back(std::make_pair("Volumetric Features Surface (Voxel based)", surface)); featureList.push_back(std::make_pair("Volumetric Features Centre of mass shift", differenceOfCenters)); featureList.push_back(std::make_pair("Volumetric Features Centre of mass shift (Uncorrected)", differenceOfCentersUncorrected)); } mitk::GIFVolumetricStatistics::GIFVolumetricStatistics() { + SetLongName("volume"); + SetShortName("vol"); } mitk::GIFVolumetricStatistics::FeatureListType mitk::GIFVolumetricStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; AccessByItk_2(image, CalculateVolumeStatistic, mask, featureList); AccessByItk_2(mask, CalculateLargestDiameter, image, featureList); vtkSmartPointer mesher = vtkSmartPointer::New(); vtkSmartPointer stats = vtkSmartPointer::New(); mesher->SetInputData(mask->GetVtkImageData()); stats->SetInputConnection(mesher->GetOutputPort()); stats->Update(); double pi = vnl_math::pi; double meshVolume = stats->GetVolume(); double meshSurf = stats->GetSurfaceArea(); double pixelVolume = featureList[0].second; double pixelSurface = featureList[3].second; MITK_INFO << "Surf: " << pixelSurface << " Vol " << pixelVolume; double compactness1 = pixelVolume / (std::sqrt(pi) * std::pow(meshSurf, 2.0 / 3.0)); double compactness1Pixel = pixelVolume / (std::sqrt(pi) * std::pow(pixelSurface, 2.0 / 3.0)); //This is the definition used by Aertz. However, due to 2/3 this feature is not demensionless. Use compactness3 instead. double compactness2 = 36 * pi*pixelVolume*pixelVolume / meshSurf / meshSurf / meshSurf; double compactness2Pixel = 36 * pi*pixelVolume*pixelVolume / pixelSurface / pixelSurface / pixelSurface; double compactness3 = pixelVolume / (std::sqrt(pi) * std::pow(meshSurf, 3.0 / 2.0)); double compactness3Pixel = pixelVolume / (std::sqrt(pi) * std::pow(pixelSurface, 3.0 / 2.0)); double sphericity = std::pow(pi, 1 / 3.0) *std::pow(6 * pixelVolume, 2.0 / 3.0) / meshSurf; double sphericityPixel = std::pow(pi, 1 / 3.0) *std::pow(6 * pixelVolume, 2.0 / 3.0) / pixelSurface; double surfaceToVolume = meshSurf / pixelVolume; double surfaceToVolumePixel = pixelSurface / pixelVolume; double sphericalDisproportion = meshSurf / 4 / pi / std::pow(3.0 / 4.0 / pi * pixelVolume, 2.0 / 3.0); double sphericalDisproportionPixel = pixelSurface / 4 / pi / std::pow(3.0 / 4.0 / pi * pixelVolume, 2.0 / 3.0); double asphericity = std::pow(1.0/compactness2, (1.0 / 3.0)) - 1; double asphericityPixel = std::pow(1.0/compactness2Pixel, (1.0 / 3.0)) - 1; //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]; std::vector pointsForPCA; std::vector pointsForPCAUncorrected; 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; mitkPixelTypeMultiplex5( mitk::FastSinglePixelAccess, image->GetChannelDescriptor().GetPixelType(), image, image->GetVolumeData(), index, pxImage, 0); mitkPixelTypeMultiplex5( mitk::FastSinglePixelAccess, mask->GetChannelDescriptor().GetPixelType(), mask, mask->GetVolumeData(), index, pxMask, 0); //Check if voxel is contained in segmentation if (pxMask > 0) { cv::Point3d tmp; tmp.x = x * xd; tmp.y = y * yd; tmp.z = z * zd; pointsForPCAUncorrected.push_back(tmp); if (pxImage == pxImage) { pointsForPCA.push_back(tmp); } } } } } //Calculate PCA Features int sz = pointsForPCA.size(); cv::Mat data_pts = cv::Mat(sz, 3, CV_64FC1); for (int i = 0; i < data_pts.rows; ++i) { data_pts.at(i, 0) = pointsForPCA[i].x; data_pts.at(i, 1) = pointsForPCA[i].y; data_pts.at(i, 2) = pointsForPCA[i].z; } //Calculate PCA Features int szUC = pointsForPCAUncorrected.size(); cv::Mat data_ptsUC = cv::Mat(szUC, 3, CV_64FC1); for (int i = 0; i < data_ptsUC.rows; ++i) { data_ptsUC.at(i, 0) = pointsForPCAUncorrected[i].x; data_ptsUC.at(i, 1) = pointsForPCAUncorrected[i].y; data_ptsUC.at(i, 2) = pointsForPCAUncorrected[i].z; } //Perform PCA analysis cv::PCA pca_analysis(data_pts, cv::Mat(), CV_PCA_DATA_AS_ROW); cv::PCA pca_analysisUC(data_ptsUC, cv::Mat(), CV_PCA_DATA_AS_ROW); //Store the eigenvalues std::vector eigen_val(3); std::vector eigen_valUC(3); for (int i = 0; i < 3; ++i) { eigen_val[i] = pca_analysis.eigenvalues.at(0, i); eigen_valUC[i] = pca_analysisUC.eigenvalues.at(0, i); } std::sort(eigen_val.begin(), eigen_val.end()); std::sort(eigen_valUC.begin(), eigen_valUC.end()); double major = eigen_val[2]; double minor = eigen_val[1]; double least = eigen_val[0]; double elongation = major == 0 ? 0 : minor/major; double flatness = major == 0 ? 0 : least / major; double majorUC = eigen_valUC[2]; double minorUC = eigen_valUC[1]; double leastUC = eigen_valUC[0]; double elongationUC = majorUC == 0 ? 0 : minorUC / majorUC; double flatnessUC = majorUC == 0 ? 0 : leastUC / majorUC; featureList.push_back(std::make_pair("Volumetric Features Volume (mesh based)",meshVolume)); featureList.push_back(std::make_pair("Volumetric Features Surface area",meshSurf)); featureList.push_back(std::make_pair("Volumetric Features Surface to volume ratio",surfaceToVolume)); featureList.push_back(std::make_pair("Volumetric Features Sphericity",sphericity)); featureList.push_back(std::make_pair("Volumetric Features Asphericity",asphericity)); featureList.push_back(std::make_pair("Volumetric Features Compactness 1",compactness1)); featureList.push_back(std::make_pair("Volumetric Features Compactness 2",compactness2)); featureList.push_back(std::make_pair("Volumetric Features Compactness 3",compactness3)); featureList.push_back(std::make_pair("Volumetric Features Spherical disproportion", sphericalDisproportion)); featureList.push_back(std::make_pair("Volumetric Features Surface to volume ratio (Voxel based)", surfaceToVolumePixel)); featureList.push_back(std::make_pair("Volumetric Features Sphericity (Voxel based)", sphericityPixel)); featureList.push_back(std::make_pair("Volumetric Features Asphericity (Voxel based)", asphericityPixel)); featureList.push_back(std::make_pair("Volumetric Features Compactness 1 (Voxel based)", compactness1Pixel)); featureList.push_back(std::make_pair("Volumetric Features Compactness 2 (Voxel based)", compactness2Pixel)); featureList.push_back(std::make_pair("Volumetric Features Compactness 3 (Voxel based)", compactness3Pixel)); featureList.push_back(std::make_pair("Volumetric Features Spherical disproportion (Voxel based)", sphericalDisproportionPixel)); featureList.push_back(std::make_pair("Volumetric Features PCA Major Axis",major)); featureList.push_back(std::make_pair("Volumetric Features PCA Minor Axis",minor)); featureList.push_back(std::make_pair("Volumetric Features PCA Least Axis",least)); featureList.push_back(std::make_pair("Volumetric Features PCA Elongation",elongation)); featureList.push_back(std::make_pair("Volumetric Features PCA Flatness",flatness)); featureList.push_back(std::make_pair("Volumetric Features PCA Major Axis (Uncorrected)", majorUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Minor Axis (Uncorrected)", minorUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Least Axis (Uncorrected)", leastUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Elongation (Uncorrected)", elongationUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Flatness (Uncorrected)", flatnessUC)); return featureList; } mitk::GIFVolumetricStatistics::FeatureNameListType mitk::GIFVolumetricStatistics::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("Volumetric Features Compactness 1"); featureList.push_back("Volumetric Features Compactness 2"); featureList.push_back("Volumetric Features Compactness 3"); featureList.push_back("Volumetric Features Sphericity"); featureList.push_back("Volumetric Features Asphericity"); featureList.push_back("Volumetric Features Surface to volume ratio"); featureList.push_back("Volumetric Features Surface area"); featureList.push_back("Volumetric Features Volume (mesh based)"); featureList.push_back("Volumetric Features Volume (pixel based)"); featureList.push_back("Volumetric Features Spherical disproportion"); featureList.push_back("Volumetric Features Maximum 3D diameter"); return featureList; } + + +void mitk::GIFVolumetricStatistics::AddArguments(mitkCommandLineParser &parser) +{ + std::string name = GetOptionPrefix(); + + parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Volume-Statistic", "calculates volume based features", us::Any()); +} + +void +mitk::GIFVolumetricStatistics::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList) +{ + auto parsedArgs = GetParameter(); + if (parsedArgs.count(GetLongName())) + { + MITK_INFO << "Start calculating volumetric features ...."; + auto localResults = this->CalculateFeatures(feature, mask); + featureList.insert(featureList.end(), localResults.begin(), localResults.end()); + MITK_INFO << "Finished calculating volumetric features...."; + } +} +