diff --git a/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp b/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp index 20c3037159..b0f40a3232 100644 --- a/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp +++ b/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp @@ -1,593 +1,593 @@ /*=================================================================== 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 #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); } template static void CreateNoNaNMask(itk::Image* itkValue, mitk::Image::Pointer mask, mitk::Image::Pointer& newMask) { typedef itk::Image< TPixel, VImageDimension> LFloatImageType; typedef itk::Image< unsigned char, VImageDimension> LMaskImageType; - LMaskImageType::Pointer itkMask = LMaskImageType::New(); + typename LMaskImageType::Pointer itkMask = LMaskImageType::New(); mitk::CastToItkImage(mask, itkMask); typedef itk::ImageDuplicator< LMaskImageType > DuplicatorType; - DuplicatorType::Pointer duplicator = DuplicatorType::New(); + typename 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); } template static void ResampleMask(itk::Image* itkMoving, mitk::Image::Pointer ref, mitk::Image::Pointer& newMask) { typedef itk::Image< TPixel, VImageDimension> LMaskImageType; typedef itk::NearestNeighborInterpolateImageFunction< LMaskImageType> NearestNeighborInterpolateImageFunctionType; typedef itk::ResampleImageFilter ResampleFilterType; - NearestNeighborInterpolateImageFunctionType::Pointer nn_interpolator = NearestNeighborInterpolateImageFunctionType::New(); - LMaskImageType::Pointer itkRef = LMaskImageType::New(); + typename NearestNeighborInterpolateImageFunctionType::Pointer nn_interpolator = NearestNeighborInterpolateImageFunctionType::New(); + typename LMaskImageType::Pointer itkRef = LMaskImageType::New(); mitk::CastToItkImage(ref, itkRef); - ResampleFilterType::Pointer resampler = ResampleFilterType::New(); + typename 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); } static void ExtractSlicesFromImages(mitk::Image::Pointer image, mitk::Image::Pointer mask, mitk::Image::Pointer maskNoNaN, int direction, std::vector &imageVector, std::vector &maskVector, std::vector &maskNoNaNVector) { typedef itk::Image< double, 2 > FloatImage2DType; typedef itk::Image< unsigned char, 2 > MaskImage2DType; FloatImageType::Pointer itkFloat = FloatImageType::New(); MaskImageType::Pointer itkMask = MaskImageType::New(); MaskImageType::Pointer itkMaskNoNaN = MaskImageType::New(); mitk::CastToItkImage(mask, itkMask); mitk::CastToItkImage(maskNoNaN, itkMaskNoNaN); mitk::CastToItkImage(image, itkFloat); int idxA, idxB, idxC; switch (direction) { case 0: idxA = 1; idxB = 2; idxC = 0; break; case 1: idxA = 0; idxB = 2; idxC = 1; break; case 2: idxA = 0; idxB = 1; idxC = 2; break; default: idxA = 1; idxB = 2; idxC = 0; break; } auto imageSize = image->GetLargestPossibleRegion().GetSize(); FloatImageType::IndexType index3D; FloatImage2DType::IndexType index2D; FloatImage2DType::SpacingType spacing2D; spacing2D[0] = itkFloat->GetSpacing()[idxA]; spacing2D[1] = itkFloat->GetSpacing()[idxB]; for (int i = 0; i < imageSize[idxC]; ++i) { FloatImage2DType::RegionType region; FloatImage2DType::IndexType start; FloatImage2DType::SizeType size; start[0] = 0; start[1] = 0; size[0] = imageSize[idxA]; size[1] = imageSize[idxB]; region.SetIndex(start); region.SetSize(size); FloatImage2DType::Pointer image2D = FloatImage2DType::New(); image2D->SetRegions(region); image2D->Allocate(); MaskImage2DType::Pointer mask2D = MaskImage2DType::New(); mask2D->SetRegions(region); mask2D->Allocate(); MaskImage2DType::Pointer masnNoNaN2D = MaskImage2DType::New(); masnNoNaN2D->SetRegions(region); masnNoNaN2D->Allocate(); for (int a = 0; a < imageSize[idxA]; ++a) { for (int b = 0; b < imageSize[idxB]; ++b) { index3D[idxA] = a; index3D[idxB] = b; index3D[idxC] = i; index2D[0] = a; index2D[1] = b; image2D->SetPixel(index2D, itkFloat->GetPixel(index3D)); mask2D->SetPixel(index2D, itkMask->GetPixel(index3D)); masnNoNaN2D->SetPixel(index2D, itkMaskNoNaN->GetPixel(index3D)); } } image2D->SetSpacing(spacing2D); mask2D->SetSpacing(spacing2D); masnNoNaN2D->SetSpacing(spacing2D); mitk::Image::Pointer tmpFloatImage = mitk::Image::New(); tmpFloatImage->InitializeByItk(image2D.GetPointer()); mitk::GrabItkImageMemory(image2D, tmpFloatImage); mitk::Image::Pointer tmpMaskImage = mitk::Image::New(); tmpMaskImage->InitializeByItk(mask2D.GetPointer()); mitk::GrabItkImageMemory(mask2D, tmpMaskImage); mitk::Image::Pointer tmpMaskNoNaNImage = mitk::Image::New(); tmpMaskNoNaNImage->InitializeByItk(masnNoNaN2D.GetPointer()); mitk::GrabItkImageMemory(masnNoNaN2D, tmpMaskNoNaNImage); imageVector.push_back(tmpFloatImage); maskVector.push_back(tmpMaskImage); maskNoNaNVector.push_back(tmpMaskNoNaNImage); } } int main(int argc, char* argv[]) { mitk::GIFFirstOrderStatistics::Pointer firstOrderCalculator = mitk::GIFFirstOrderStatistics::New(); mitk::GIFVolumetricStatistics::Pointer volCalculator = mitk::GIFVolumetricStatistics::New(); mitk::GIFCooccurenceMatrix::Pointer coocCalculator = mitk::GIFCooccurenceMatrix::New(); mitk::GIFCooccurenceMatrix2::Pointer cooc2Calculator = mitk::GIFCooccurenceMatrix2::New(); mitk::GIFNeighbouringGreyLevelDependenceFeature::Pointer ngldCalculator = mitk::GIFNeighbouringGreyLevelDependenceFeature::New(); mitk::GIFGrayLevelRunLength::Pointer rlCalculator = mitk::GIFGrayLevelRunLength::New(); mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); // required params parser.addArgument("image", "i", mitkCommandLineParser::InputImage, "Input Image", "Path to the input image file", us::Any(), false); parser.addArgument("mask", "m", mitkCommandLineParser::InputImage, "Input Mask", "Path to the mask Image that specifies the area over for the statistic (Values = 1)", us::Any(), false); parser.addArgument("output", "o", mitkCommandLineParser::OutputFile, "Output text file", "Path to output file. The output statistic is appended to this file.", us::Any(), false); parser.addArgument("--","-", mitkCommandLineParser::String, "---", "---", us::Any(),true); firstOrderCalculator->AddArguments(parser); // Does not support single direction volCalculator->AddArguments(parser); // Does not support single direction coocCalculator->AddArguments(parser); // cooc2Calculator->AddArguments(parser); // Needs parameter fixing ngldCalculator->AddArguments(parser); // Needs parameter fixing rlCalculator->AddArguments(parser); // Does not support single direction parser.addArgument("--", "-", mitkCommandLineParser::String, "---", "---", us::Any(), true); //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()); parser.addArgument("slice-wise", "slice", mitkCommandLineParser::String, "Int", "Allows to specify if the image is processed slice-wise (number giving direction) ", us::Any()); parser.addArgument("minimum-intensity", "minimum", mitkCommandLineParser::String, "Float", "Minimum intensity. If set, it is overwritten by more specific intensity minima", us::Any()); parser.addArgument("maximum-intensity", "maximum", mitkCommandLineParser::String, "Float", "Maximum intensity. If set, it is overwritten by more specific intensity maxima", us::Any()); parser.addArgument("bins", "bins", mitkCommandLineParser::String, "Int", "Number of bins if bins are used. If set, it is overwritten by more specific bin count", 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.9; //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::Image::Pointer mask; mitk::Image::Pointer tmpImage = mitk::IOUtil::LoadImage(parsedArgs["image"].ToString()); mitk::Image::Pointer tmpMask = mitk::IOUtil::LoadImage(parsedArgs["mask"].ToString()); image = tmpImage; mask = tmpMask; if ((image->GetDimension() != mask->GetDimension())) { MITK_INFO << "Dimension of image does not match. "; MITK_INFO << "Correct one image, may affect the result"; if (image->GetDimension() == 2) { mitk::Convert2Dto3DImageFilter::Pointer multiFilter2 = mitk::Convert2Dto3DImageFilter::New(); multiFilter2->SetInput(tmpImage); multiFilter2->Update(); image = multiFilter2->GetOutput(); } if (mask->GetDimension() == 2) { mitk::Convert2Dto3DImageFilter::Pointer multiFilter3 = mitk::Convert2Dto3DImageFilter::New(); multiFilter3->SetInput(tmpMask); multiFilter3->Update(); mask = multiFilter3->GetOutput(); } } 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(); AccessByItk_2(mask, ResampleMask, 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) { MITK_INFO << "Warning!"; MITK_INFO << "The origin of the input image and the mask do not match. They are"; MITK_INFO << "now corrected. Please check to make sure that the images still match"; 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) { MITK_INFO << "Warning!"; MITK_INFO << "The spacing of the mask was set to match the spacing of the input image."; MITK_INFO << "This might cause unintended spacing of the mask image"; image->GetGeometry(0)->SetSpacing(mask->GetGeometry(0)->GetSpacing()); } else { MITK_INFO << "The spacing of the mask and the input images is not equal."; MITK_INFO << "Terminating the programm. You may use the '-fi' option"; return -1; } } int direction = 0; if (parsedArgs.count("direction")) { direction = splitDouble(parsedArgs["direction"].ToString(), ';')[0]; } MITK_INFO << "Start creating Mask without NaN"; mitk::Image::Pointer maskNoNaN = mitk::Image::New(); AccessByItk_2(image, CreateNoNaNMask, mask, maskNoNaN); //CreateNoNaNMask(mask, image, maskNoNaN); bool sliceWise = false; int sliceDirection = 0; int currentSlice = 0; bool imageToProcess = true; std::vector floatVector; std::vector maskVector; std::vector maskNoNaNVector; if (parsedArgs.count("slice-wise")) { MITK_INFO << "Enabled slice-wise"; sliceWise = true; sliceDirection = splitDouble(parsedArgs["slice-wise"].ToString(), ';')[0]; MITK_INFO << sliceDirection; ExtractSlicesFromImages(image, mask, maskNoNaN, sliceDirection, floatVector, maskVector, maskNoNaNVector); MITK_INFO << "Slice"; } if (parsedArgs.count("minimum-intensity")) { float minimum = splitDouble(parsedArgs["minimum-intensity"].ToString(), ';')[0]; firstOrderCalculator->SetMinimumIntensity(minimum); firstOrderCalculator->SetUseMinimumIntensity(true); volCalculator->SetMinimumIntensity(minimum); volCalculator->SetUseMinimumIntensity(true); coocCalculator->SetMinimumIntensity(minimum); coocCalculator->SetUseMinimumIntensity(true); cooc2Calculator->SetMinimumIntensity(minimum); cooc2Calculator->SetUseMinimumIntensity(true); ngldCalculator->SetMinimumIntensity(minimum); ngldCalculator->SetUseMinimumIntensity(true); rlCalculator->SetMinimumIntensity(minimum); rlCalculator->SetUseMinimumIntensity(true); } if (parsedArgs.count("maximum-intensity")) { float minimum = splitDouble(parsedArgs["maximum-intensity"].ToString(), ';')[0]; firstOrderCalculator->SetMaximumIntensity(minimum); firstOrderCalculator->SetUseMaximumIntensity(true); volCalculator->SetMaximumIntensity(minimum); volCalculator->SetUseMaximumIntensity(true); coocCalculator->SetMaximumIntensity(minimum); coocCalculator->SetUseMaximumIntensity(true); cooc2Calculator->SetMaximumIntensity(minimum); cooc2Calculator->SetUseMaximumIntensity(true); ngldCalculator->SetMaximumIntensity(minimum); ngldCalculator->SetUseMaximumIntensity(true); rlCalculator->SetMaximumIntensity(minimum); rlCalculator->SetUseMaximumIntensity(true); } if (parsedArgs.count("bins")) { int minimum = splitDouble(parsedArgs["bins"].ToString(), ';')[0]; firstOrderCalculator->SetBins(minimum); volCalculator->SetBins(minimum); coocCalculator->SetBins(minimum); cooc2Calculator->SetBins(minimum); ngldCalculator->SetBins(minimum); rlCalculator->SetBins(minimum); } firstOrderCalculator->SetParameter(parsedArgs); firstOrderCalculator->SetDirection(direction); volCalculator->SetParameter(parsedArgs); volCalculator->SetDirection(direction); coocCalculator->SetParameter(parsedArgs); coocCalculator->SetDirection(direction); cooc2Calculator->SetParameter(parsedArgs); cooc2Calculator->SetDirection(direction); ngldCalculator->SetParameter(parsedArgs); ngldCalculator->SetDirection(direction); rlCalculator->SetParameter(parsedArgs); rlCalculator->SetDirection(direction); std::ofstream output(parsedArgs["output"].ToString(), std::ios::app); mitk::Image::Pointer cImage = image; mitk::Image::Pointer cMask = mask; mitk::Image::Pointer cMaskNoNaN = maskNoNaN; while (imageToProcess) { if (sliceWise) { cImage = floatVector[currentSlice]; cMask = maskVector[currentSlice]; cMaskNoNaN = maskNoNaNVector[currentSlice]; imageToProcess = (floatVector.size()-1 > (currentSlice)) ? true : false ; } else { imageToProcess = false; } mitk::AbstractGlobalImageFeature::FeatureListType stats; firstOrderCalculator->CalculateFeaturesUsingParameters(cImage, cMask, cMaskNoNaN, stats); volCalculator->CalculateFeaturesUsingParameters(cImage, cMask, cMaskNoNaN, stats); coocCalculator->CalculateFeaturesUsingParameters(cImage, cMask, cMaskNoNaN, stats); cooc2Calculator->CalculateFeaturesUsingParameters(cImage, cMask, cMaskNoNaN, stats); ngldCalculator->CalculateFeaturesUsingParameters(cImage, cMask, cMaskNoNaN, stats); rlCalculator->CalculateFeaturesUsingParameters(cImage, cMask, cMaskNoNaN, stats); for (std::size_t i = 0; i < stats.size(); ++i) { std::cout << stats[i].first << " - " << stats[i].second << std::endl; } if (parsedArgs.count("header") && (currentSlice == 0)) { if (parsedArgs.count("description")) { output << "Description" << ";"; } if (sliceWise) { output << "SliceNumber" << ";"; } for (std::size_t i = 0; i < stats.size(); ++i) { output << stats[i].first << ";"; } output << "EndOfMeasurement;"; output << std::endl; } if (parsedArgs.count("description")) { output << parsedArgs["description"].ToString() << ";"; } if (sliceWise) { output << currentSlice << ";"; } for (std::size_t i = 0; i < stats.size(); ++i) { output << stats[i].second << ";"; } output << "EndOfMeasurement;"; output << std::endl; ++currentSlice; } if (sliceWise) { output << "EndOfFile" << std::endl; } output.close(); return 0; } #endif diff --git a/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h b/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h index cb16ec35f7..2665167adc 100644 --- a/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h +++ b/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h @@ -1,149 +1,147 @@ #ifndef mitkGIFCooccurenceMatrix2_h #define mitkGIFCooccurenceMatrix2_h #include #include #include -//#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); 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(Bins, int); itkSetMacro(Bins, int); struct GIFCooccurenceMatrix2Configuration { double range; unsigned int direction; double MinimumIntensity; bool UseMinimumIntensity; double MaximumIntensity; bool UseMaximumIntensity; int Bins; }; private: double m_Range; int m_Bins; }; } #endif //mitkGIFCooccurenceMatrix2_h diff --git a/Modules/Classification/CLUtilities/include/mitkGIFNeighbouringGreyLevelDependenceFeatures.h b/Modules/Classification/CLUtilities/include/mitkGIFNeighbouringGreyLevelDependenceFeatures.h index 47834e8f17..cdc38d25c8 100644 --- a/Modules/Classification/CLUtilities/include/mitkGIFNeighbouringGreyLevelDependenceFeatures.h +++ b/Modules/Classification/CLUtilities/include/mitkGIFNeighbouringGreyLevelDependenceFeatures.h @@ -1,146 +1,146 @@ #ifndef mitkGIFNeighbouringGreyLevelDependenceFeatures_h #define mitkGIFNeighbouringGreyLevelDependenceFeatures_h #include #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(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; double MinimumIntensity; bool UseMinimumIntensity; double MaximumIntensity; bool UseMaximumIntensity; int Bins; }; private: double m_Range; int m_Alpha; int m_Bins; }; } #endif //mitkGIFNeighbouringGreyLevelDependenceFeature_h diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp index 6e9100079b..7752f3e433 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp @@ -1,610 +1,604 @@ #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 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(); if (config.UseMinimumIntensity) rangeMin = config.MinimumIntensity; if (config.UseMaximumIntensity) rangeMax = config.MaximumIntensity; // Define Range int numberOfBins = config.Bins; 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); } } if (config.direction == 1) { offsetVector.clear(); offset[0] = 0; offset[1] = 0; offset[2] = 1; } std::vector resultVector; mitk::CoocurenceMatrixHolder holderOverall(rangeMin, rangeMax, numberOfBins); mitk::CoocurenceMatrixFeatures overallFeature; for (int i = 0; i < offsetVector.size(); ++i) { if (config.direction > 1) { if (offsetVector[i][config.direction - 2] != 0) { continue; } } offset = offsetVector[i]; 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) + for (std::size_t 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_Bins(128) { SetShortName("cooc2"); SetLongName("cooccurence2"); } mitk::GIFCooccurenceMatrix2::FeatureListType mitk::GIFCooccurenceMatrix2::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; GIFCooccurenceMatrix2Configuration config; config.direction = GetDirection(); config.range = m_Range; config.MinimumIntensity = GetMinimumIntensity(); config.MaximumIntensity = GetMaximumIntensity(); config.UseMinimumIntensity = GetUseMinimumIntensity(); config.UseMaximumIntensity = GetUseMaximumIntensity(); config.Bins = GetBins(); AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config); return featureList; } mitk::GIFCooccurenceMatrix2::FeatureNameListType mitk::GIFCooccurenceMatrix2::GetFeatureNames() { FeatureNameListType featureList; 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 + "::bins", name + "::bins", mitkCommandLineParser::String, "Cooc 2 Number of Bins", "Define the number of bins that is used ", us::Any()); } void mitk::GIFCooccurenceMatrix2::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); std::string name = GetOptionPrefix(); if (parsedArgs.count(GetLongName())) { std::vector ranges; if (parsedArgs.count(name + "::range")) { ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); } else { ranges.push_back(1); } if (parsedArgs.count(name + "::bins")) { auto bins = SplitDouble(parsedArgs[name + "::bins"].ToString(), ';')[0]; this->SetBins(bins); } for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; this->SetRange(ranges[i]); 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/mitkGIFGrayLevelSizeZone.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp index ce16f9d139..e0e80f4be9 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp @@ -1,284 +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) { 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 = GetDirection(); 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) +mitk::GIFGrayLevelSizeZone::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, 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 b845130954..47fe2d6e2a 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp @@ -1,224 +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) { 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 = GetDirection(); 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) +mitk::GIFNeighbourhoodGreyLevelDifference::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, 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 index c9a537408e..464734c7d3 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp @@ -1,410 +1,404 @@ #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_Stepsize(0), + m_NumberOfDependences(depenence), m_NumberOfBins(number), - m_NumberOfDependences(depenence) + m_NeighbourhoodSize(1), + m_NumberOfNeighbourVoxels(0), + m_NumberOfDependenceNeighbourVoxels(0), + m_NumberOfNeighbourhoods(0), + m_NumberOfCompleteNeighbourhoods(0) { 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, int direction, 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); if ((direction > 1) && (direction +2 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 LocalCalculateFeatures( 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 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(); if (config.UseMinimumIntensity) rangeMin = config.MinimumIntensity; if (config.UseMaximumIntensity) rangeMax = config.MaximumIntensity; // Define Range int numberOfBins = config.Bins; 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, config.direction, holderOverall); LocalCalculateFeatures(holderOverall, overallFeature); 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 High 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_Bins(6) { SetShortName("ngldm"); SetLongName("ngldm"); } mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType mitk::GIFNeighbouringGreyLevelDependenceFeature::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; GIFNeighbouringGreyLevelDependenceFeatureConfiguration config; config.direction = GetDirection(); config.range = m_Range; config.alpha = 0; config.MinimumIntensity = GetMinimumIntensity(); config.MaximumIntensity = GetMaximumIntensity(); config.UseMinimumIntensity = GetUseMinimumIntensity(); config.UseMaximumIntensity = GetUseMaximumIntensity(); config.Bins = GetBins(); 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, "NGLD Range", "Define the range that is used (Semicolon-separated)", us::Any()); parser.addArgument(name + "::bins", name + "::bins", mitkCommandLineParser::String, "NGLD Number of Bins", "Define the number of bins that is used ", us::Any()); } void mitk::GIFNeighbouringGreyLevelDependenceFeature::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); std::string name = GetOptionPrefix(); if (parsedArgs.count(GetLongName())) { std::vector ranges; if (parsedArgs.count(name + "::range")) { ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); } else { ranges.push_back(1); } if (parsedArgs.count(name + "::bins")) { auto bins = SplitDouble(parsedArgs[name + "::bins"].ToString(), ';')[0]; this->SetBins(bins); } for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; this->SetRange(ranges[i]); auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; } } }