diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp index 0652ded536..ee96579555 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp @@ -1,591 +1,590 @@ #include // MITK #include #include #include // ITK #include #include #include // STL #include #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) { int index = std::floor((intensity - m_MinimumRange) / m_Stepsize); return std::max(0, std::min(index, m_NumberOfBins - 1)); } 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 = i + 1;// 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 = i + 1; // 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 iInt = i + 1;// holder.IndexToMeanIntensity(i); double jInt = j + 1;// 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 iInt = i + 1; 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::Neighborhood NeighborhoodType; typedef itk::Offset OffsetType; /////////////////////////////////////////////////////////////////////////////////////////////// double rangeMin = config.MinimumIntensity; double rangeMax = config.MaximumIntensity; 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 (unsigned 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 (std::size_t 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-occurenced Based Features (" + strRange + ")::Overall", featureList); MatrixFeaturesTo(featureMean, "Co-occurenced Based Features (" + strRange + ")::Mean", featureList); MatrixFeaturesTo(featureStd, "Co-occurenced Based Features (" + 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 (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) { SetShortName("cooc2"); SetLongName("cooccurence2"); } mitk::GIFCooccurenceMatrix2::FeatureListType mitk::GIFCooccurenceMatrix2::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { InitializeQuantifier(image, mask); FeatureListType featureList; GIFCooccurenceMatrix2Configuration config; config.direction = GetDirection(); config.range = m_Range; config.MinimumIntensity = GetQuantifier()->GetMinimum(); config.MaximumIntensity = GetQuantifier()->GetMaximum(); config.Bins = GetQuantifier()->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::Bool, "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()); AddQuantifierArguments(parser); } 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())) { InitializeQuantifierFromParameters(feature, maskNoNAN); 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]); 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/mitkGIFFirstOrderHistogramStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderHistogramStatistics.cpp index 526772545c..574e5c3470 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderHistogramStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderHistogramStatistics.cpp @@ -1,332 +1,331 @@ /*=================================================================== 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 #include // STL #include #include #include #define GET_VARIABLE_INDEX(value) \ mv[0] = value; \ histogram->GetIndex(mv, resultingIndex); template void CalculateFirstOrderHistogramStatistics(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFFirstOrderHistogramStatistics::FeatureListType & featureList, mitk::GIFFirstOrderHistogramStatistics::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 FilterType::Pointer labelStatisticsImageFilter = FilterType::New(); labelStatisticsImageFilter->SetInput(itkImage); labelStatisticsImageFilter->SetLabelInput(maskImage); labelStatisticsImageFilter->SetUseHistograms(true); labelStatisticsImageFilter->SetHistogramParameters(params.Bins, params.MinimumIntensity, params.MaximumIntensity); labelStatisticsImageFilter->Update(); - HistogramType::MeasurementVectorType mv(1); + typename HistogramType::MeasurementVectorType mv(1); mv[0] = 4.1; - HistogramType::IndexType resultingIndex; + typename HistogramType::IndexType resultingIndex; auto histogram = labelStatisticsImageFilter->GetHistogram(1); double meanValue = 0; // labelStatisticsImageFilter->GetMean(1); GET_VARIABLE_INDEX(meanValue); double meanIndex = 0; // resultingIndex[0]; double medianValue = labelStatisticsImageFilter->GetMedian(1); GET_VARIABLE_INDEX(medianValue); double medianIndex = resultingIndex[0]; double minimumValue = labelStatisticsImageFilter->GetMinimum(1); GET_VARIABLE_INDEX(minimumValue); double minimumIndex = resultingIndex[0]; double p10Value = histogram->Quantile(0, 0.10); GET_VARIABLE_INDEX(p10Value); double p10Index = resultingIndex[0]; double p25Value = histogram->Quantile(0, 0.25); GET_VARIABLE_INDEX(p25Value); double p25Index = resultingIndex[0]; double p75Value = histogram->Quantile(0, 0.75); GET_VARIABLE_INDEX(p75Value); double p75Index = resultingIndex[0]; double p90Value = histogram->Quantile(0, 0.90); GET_VARIABLE_INDEX(p90Value); double p90Index = resultingIndex[0]; double maximumValue = labelStatisticsImageFilter->GetMaximum(1); GET_VARIABLE_INDEX(maximumValue); double maximumIndex = resultingIndex[0]; double Log2 = log(2); HIndexType index; HIndexType index2; index.SetSize(1); index2.SetSize(1); double binWidth = histogram->GetBinMax(0, 0) - histogram->GetBinMin(0, 0); double count = labelStatisticsImageFilter->GetCount(1); double robustMeanValue = 0; double robustMeanIndex = 0; double robustCount = 0; for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; double frequence = histogram->GetFrequency(index); double probability = frequence / count; double voxelValue = histogram->GetBinMin(0, i) + binWidth * 0.5; meanValue += probability * voxelValue; meanIndex += probability * (i); if ((i >= p10Index) && (i <= p90Index)) { robustMeanValue += frequence * voxelValue; robustMeanIndex += frequence * i; robustCount += frequence; } } robustMeanValue /= robustCount; robustMeanIndex /= robustCount; double varianceValue = 0; double varianceIndex = 0; double skewnessValue = 0; double skewnessIndex = 0; double kurtosisValue = 0; double kurtosisIndex = 0; double modeValue = 0; double modeIndex = 0; double modeFrequence = 0; double meanAbsoluteDeviationValue = 0; double meanAbsoluteDeviationIndex = 0; double robustMeanAbsoluteDeviationValue = 0; double robustMeanAbsoluteDeivationIndex = 0; double medianAbsoluteDeviationValue = 0; double medianAbsoluteDeviationIndex = 0; double coefficientOfVariationValue = 0; double coefficientOfVariationIndex = 0; double quantileCoefficientOfDispersionValue = 0; double quantileCoefficientOfDispersionIndex = 0; double entropyValue = 0; double entropyIndex = 0; double uniformityValue = 0; double uniformityIndex = 0; double maximumGradientValue = std::numeric_limits::min(); double maximumGradientIndex = 0; double minimumGradientValue = std::numeric_limits::max(); double minimumGradientIndex = 0; double gradient = 0; for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; double frequence = histogram->GetFrequency(index); double probability = frequence / count; double voxelValue = histogram->GetBinMin(0, i) + binWidth * 0.5; double deltaValue = (voxelValue - meanValue); double deltaIndex = (i - meanIndex); varianceValue += probability * deltaValue * deltaValue; varianceIndex += probability * deltaIndex * deltaIndex; skewnessValue += probability * deltaValue * deltaValue * deltaValue; skewnessIndex += probability * deltaIndex * deltaIndex * deltaIndex; kurtosisValue += probability * deltaValue * deltaValue * deltaValue * deltaValue; kurtosisIndex += probability * deltaIndex * deltaIndex * deltaIndex * deltaIndex; if (modeFrequence < frequence) { modeFrequence = frequence; modeValue = voxelValue; modeIndex = i; } meanAbsoluteDeviationValue += probability * std::abs(deltaValue); meanAbsoluteDeviationIndex += probability * std::abs(deltaIndex); if ((i >= p10Index) && (i <= p90Index)) { robustMeanAbsoluteDeviationValue += frequence * std::abs(voxelValue - robustMeanValue); robustMeanAbsoluteDeivationIndex += frequence * std::abs(i*1.0 - robustMeanIndex*1.0); } medianAbsoluteDeviationValue += probability * std::abs(voxelValue - medianValue); medianAbsoluteDeviationIndex += probability * std::abs(i*1.0 - medianIndex); if (probability > 0.0000001) { entropyValue -= probability * std::log(probability) / Log2; entropyIndex = entropyValue; } uniformityValue += probability*probability; uniformityIndex = uniformityValue; if (i == 0) { index[0] = 1; index2[0] = 0; gradient = histogram->GetFrequency(index)*1.0 - histogram->GetFrequency(index2)*1.0; } else if (i == (int)(histogram->GetSize(0)) - 1) { index[0] = i; index2[0] = i - 1; gradient = histogram->GetFrequency(index)*1.0 - histogram->GetFrequency(index2)*1.0; } else { index[0] = i+1; index2[0] = i - 1; gradient = (histogram->GetFrequency(index)*1.0 - histogram->GetFrequency(index2)*1.0) / 2.0; } if (gradient > maximumGradientValue) { maximumGradientValue = gradient; maximumGradientIndex = i + 1; } if (gradient < minimumGradientValue) { minimumGradientValue = gradient; minimumGradientIndex = i + 1; } } skewnessValue = skewnessValue / (varianceValue * std::sqrt(varianceValue)); skewnessIndex = skewnessIndex / (varianceIndex * std::sqrt(varianceIndex)); kurtosisValue = kurtosisValue / (varianceValue * varianceValue) - 3; // Excess Kurtosis kurtosisIndex = kurtosisIndex / (varianceIndex * varianceIndex) - 3; // Excess Kurtosis coefficientOfVariationValue = std::sqrt(varianceValue) / meanValue; coefficientOfVariationIndex = std::sqrt(varianceIndex) / (meanIndex+1); quantileCoefficientOfDispersionValue = (p75Value - p25Value) / (p75Value + p25Value); quantileCoefficientOfDispersionIndex = (p75Index - p25Index) / (p75Index + p25Index); robustMeanAbsoluteDeviationValue /= robustCount; robustMeanAbsoluteDeivationIndex /= robustCount; featureList.push_back(std::make_pair("FirstOrderHistogram::Mean Value", meanValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Variance Value", varianceValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Skewness Value", skewnessValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Excess Kurtosis Value", kurtosisValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Median Value", medianValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Minimum Value", minimumValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Percentile 10 Value", p10Value)); featureList.push_back(std::make_pair("FirstOrderHistogram::Percentile 90 Value", p90Value)); featureList.push_back(std::make_pair("FirstOrderHistogram::Maximum Value", maximumValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Mode Value", modeValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Interquantile Range Value", p75Value - p25Value)); featureList.push_back(std::make_pair("FirstOrderHistogram::Range Value", maximumValue - minimumValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Mean Absolute Deviation Value", meanAbsoluteDeviationValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Robust Mean Absolute Deviation Value", robustMeanAbsoluteDeviationValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Median Absolute Deviation Value", medianAbsoluteDeviationValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Coefficient of Variation Value", coefficientOfVariationValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Quantile coefficient of Dispersion Value", quantileCoefficientOfDispersionValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Entropy Value", entropyValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Uniformity Value", uniformityValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Robust Mean Value", robustMeanValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Mean Index", meanIndex + 1 )); featureList.push_back(std::make_pair("FirstOrderHistogram::Variance Index", varianceIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Skewness Index", skewnessIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Excess Kurtosis Index", kurtosisIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Median Index", medianIndex + 1)); featureList.push_back(std::make_pair("FirstOrderHistogram::Minimum Index", minimumIndex + 1)); featureList.push_back(std::make_pair("FirstOrderHistogram::Percentile 10 Index", p10Index + 1)); featureList.push_back(std::make_pair("FirstOrderHistogram::Percentile 90 Index", p90Index + 1)); featureList.push_back(std::make_pair("FirstOrderHistogram::Maximum Index", maximumIndex + 1)); featureList.push_back(std::make_pair("FirstOrderHistogram::Mode Index", modeIndex + 1)); featureList.push_back(std::make_pair("FirstOrderHistogram::Interquantile Range Index", p75Index - p25Index)); featureList.push_back(std::make_pair("FirstOrderHistogram::Range Index", maximumIndex - minimumIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Mean Absolute Deviation Index", meanAbsoluteDeviationIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Robust Mean Absolute Deviation Index", robustMeanAbsoluteDeivationIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Median Absolute Deviation Index", medianAbsoluteDeviationIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Coefficient of Variation Index", coefficientOfVariationIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Quantile coefficient of Dispersion Index", quantileCoefficientOfDispersionIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Entropy Index", entropyIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Uniformity Index", uniformityIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Maximum Gradient", maximumGradientValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Maximum Gradient Index", maximumGradientIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Minimum Gradient", minimumGradientValue)); featureList.push_back(std::make_pair("FirstOrderHistogram::Minimum Gradient Index", minimumGradientIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Robust Mean Index", robustMeanIndex)); featureList.push_back(std::make_pair("FirstOrderHistogram::Number of Bins", histogram->GetSize(0))); featureList.push_back(std::make_pair("FirstOrderHistogram::Bin Size", binWidth)); } mitk::GIFFirstOrderHistogramStatistics::GIFFirstOrderHistogramStatistics() : m_HistogramSize(256), m_UseCtRange(false), m_BinSize(-1) { SetShortName("foh"); SetLongName("first-order-histogram"); } mitk::GIFFirstOrderHistogramStatistics::FeatureListType mitk::GIFFirstOrderHistogramStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { InitializeQuantifier(image, mask); FeatureListType featureList; ParameterStruct params; params.MinimumIntensity = GetQuantifier()->GetMinimum(); params.MaximumIntensity = GetQuantifier()->GetMaximum(); params.Bins = GetQuantifier()->GetBins(); AccessByItk_3(image, CalculateFirstOrderHistogramStatistics, mask, featureList, params); return featureList; } mitk::GIFFirstOrderHistogramStatistics::FeatureNameListType mitk::GIFFirstOrderHistogramStatistics::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFFirstOrderHistogramStatistics::AddArguments(mitkCommandLineParser &parser) { AddQuantifierArguments(parser); std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Histogram based First order features", "calculates first order features based on a histogram", us::Any()); } void mitk::GIFFirstOrderHistogramStatistics::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { std::string name = GetOptionPrefix(); auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { InitializeQuantifierFromParameters(feature, mask); MITK_INFO << "Start calculating first order histogram features ...."; auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating first order histogram features...."; } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp index c171967bc9..6f41e925c0 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp @@ -1,322 +1,322 @@ /*=================================================================== 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); double voxelVolume = 1; for (unsigned int i = 0; i < std::min(3, VImageDimension); ++i) voxelVolume *= itkImage->GetSpacing()[i]; double voxelSpace = 1; for (unsigned int i = 0; i < VImageDimension; ++i) voxelSpace *= itkImage->GetSpacing()[i]; 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); double min = params.MinimumIntensity; double max = params.MaximumIntensity; labelStatisticsImageFilter->SetHistogramParameters(params.Bins, min,max); 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 mean = labelStatisticsImageFilter->GetMean(1); double median = labelStatisticsImageFilter->GetMedian(1); auto histogram = labelStatisticsImageFilter->GetHistogram(1); bool histogramIsCalculated = histogram; HIndexType index; index.SetSize(1); double uniformity = 0; double entropy = 0; double squared_sum = 0; double kurtosis = 0; double mean_absolut_deviation = 0; double median_absolut_deviation = 0; double skewness = 0; double sum_prob = 0; double binWidth = 0; double p05th = 0, p10th = 0, p15th = 0, p20th = 0, p25th = 0, p30th = 0, p35th = 0, p40th = 0, p45th = 0, p50th = 0; double p55th = 0, p60th = 0, p65th = 0, p70th = 0, p75th = 0, p80th = 0, p85th = 0, p90th = 0, p95th = 0; double voxelValue = 0; if (histogramIsCalculated) { binWidth = histogram->GetBinMax(0, 0) - histogram->GetBinMin(0, 0); p05th = histogram->Quantile(0, 0.05); p10th = histogram->Quantile(0, 0.10); p15th = histogram->Quantile(0, 0.15); p20th = histogram->Quantile(0, 0.20); p25th = histogram->Quantile(0, 0.25); p30th = histogram->Quantile(0, 0.30); p35th = histogram->Quantile(0, 0.35); p40th = histogram->Quantile(0, 0.40); p45th = histogram->Quantile(0, 0.45); p50th = histogram->Quantile(0, 0.50); p55th = histogram->Quantile(0, 0.55); p60th = histogram->Quantile(0, 0.60); p65th = histogram->Quantile(0, 0.65); p70th = histogram->Quantile(0, 0.70); p75th = histogram->Quantile(0, 0.75); p80th = histogram->Quantile(0, 0.80); p85th = histogram->Quantile(0, 0.85); p90th = histogram->Quantile(0, 0.90); p95th = histogram->Quantile(0, 0.95); } double Log2=log(2); double mode_bin; double mode_value = 0; double variance = 0; if (histogramIsCalculated) { for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; double prob = histogram->GetFrequency(index); if (prob < 0.00000001) continue; voxelValue = histogram->GetBinMin(0, i) + binWidth * 0.5; if (prob > mode_value) { mode_value = prob; mode_bin = voxelValue; } sum_prob += prob; squared_sum += prob * voxelValue*voxelValue; prob /= count; mean_absolut_deviation += prob* std::abs(voxelValue - mean); median_absolut_deviation += prob* std::abs(voxelValue - median); variance += prob * (voxelValue - mean) * (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 uncorrected_std_dev = std::sqrt(variance); double rms = std::sqrt(squared_sum / count); kurtosis = kurtosis / (variance * variance); skewness = skewness / (variance * uncorrected_std_dev); double coveredGrayValueRange = range / imageRange; double coefficient_of_variation = (mean == 0) ? 0 : std::sqrt(variance) / mean; double quantile_coefficient_of_dispersion = (p75th - p25th) / (p75th + p25th); //Calculate the robust mean absolute deviation //First, set all frequencies to 0 that are <10th or >90th percentile double meanRobust = 0.0; double robustMeanAbsoluteDeviation = 0.0; if (histogramIsCalculated) { 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 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(); 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("First Order::Mean", labelStatisticsImageFilter->GetMean(1))); featureList.push_back(std::make_pair("First Order::Unbiased Variance", labelStatisticsImageFilter->GetVariance(1))); //Siehe Definition von Unbiased Variance estimation. (Wird nicht durch n sondern durch n-1 normalisiert) featureList.push_back(std::make_pair("First Order::Biased Variance", variance)); featureList.push_back(std::make_pair("First Order::Skewness", skewness)); featureList.push_back(std::make_pair("First Order::Kurtosis", kurtosis)); featureList.push_back(std::make_pair("First Order::Median", labelStatisticsImageFilter->GetMedian(1))); featureList.push_back(std::make_pair("First Order::Minimum", labelStatisticsImageFilter->GetMinimum(1))); featureList.push_back(std::make_pair("First Order::Maximum", labelStatisticsImageFilter->GetMaximum(1))); featureList.push_back(std::make_pair("First Order::Range", range)); featureList.push_back(std::make_pair("First Order::Mean Absolute Deviation", mean_absolut_deviation)); featureList.push_back(std::make_pair("First Order::Robust Mean Absolute Deviation", robustMeanAbsoluteDeviation)); featureList.push_back(std::make_pair("First Order::Median Absolute Deviation", median_absolut_deviation)); featureList.push_back(std::make_pair("First Order::Coefficient Of Variation", coefficient_of_variation)); featureList.push_back(std::make_pair("First Order::Quantile Coefficient Of Dispersion", quantile_coefficient_of_dispersion)); featureList.push_back(std::make_pair("First Order::Energy", squared_sum)); featureList.push_back(std::make_pair("First Order::Root Mean Square", rms)); - HistogramType::MeasurementVectorType mv(1); + typename HistogramType::MeasurementVectorType mv(1); mv[0] = 0; - HistogramType::IndexType resultingIndex; + typename HistogramType::IndexType resultingIndex; histogram->GetIndex(mv, resultingIndex); featureList.push_back(std::make_pair("First Order::Robust Mean", meanRobust)); featureList.push_back(std::make_pair("First Order::Uniformity", uniformity)); featureList.push_back(std::make_pair("First Order::Entropy", entropy)); featureList.push_back(std::make_pair("First Order::Excess Kurtosis", kurtosis - 3)); featureList.push_back(std::make_pair("First Order::Covered Image Intensity Range", coveredGrayValueRange)); featureList.push_back(std::make_pair("First Order::Sum", labelStatisticsImageFilter->GetSum(1))); featureList.push_back(std::make_pair("First Order::Mode", mode_bin)); featureList.push_back(std::make_pair("First Order::Mode Probability", mode_value)); featureList.push_back(std::make_pair("First Order::Unbiased Standard deviation", labelStatisticsImageFilter->GetSigma(1))); featureList.push_back(std::make_pair("First Order::Biased Standard deviation", sqrt(variance))); featureList.push_back(std::make_pair("First Order::Number Of Voxels", labelStatisticsImageFilter->GetCount(1))); featureList.push_back(std::make_pair("First Order::05th Percentile", p05th)); featureList.push_back(std::make_pair("First Order::10th Percentile", p10th)); featureList.push_back(std::make_pair("First Order::15th Percentile", p15th)); featureList.push_back(std::make_pair("First Order::20th Percentile", p20th)); featureList.push_back(std::make_pair("First Order::25th Percentile", p25th)); featureList.push_back(std::make_pair("First Order::30th Percentile", p30th)); featureList.push_back(std::make_pair("First Order::35th Percentile", p35th)); featureList.push_back(std::make_pair("First Order::40th Percentile", p40th)); featureList.push_back(std::make_pair("First Order::45th Percentile", p45th)); featureList.push_back(std::make_pair("First Order::50th Percentile", p50th)); featureList.push_back(std::make_pair("First Order::55th Percentile", p55th)); featureList.push_back(std::make_pair("First Order::60th Percentile", p60th)); featureList.push_back(std::make_pair("First Order::65th Percentile", p65th)); featureList.push_back(std::make_pair("First Order::70th Percentile", p70th)); featureList.push_back(std::make_pair("First Order::75th Percentile", p75th)); featureList.push_back(std::make_pair("First Order::80th Percentile", p80th)); featureList.push_back(std::make_pair("First Order::85th Percentile", p85th)); featureList.push_back(std::make_pair("First Order::90th Percentile", p90th)); featureList.push_back(std::make_pair("First Order::95th Percentile", p95th)); featureList.push_back(std::make_pair("First Order::Interquartile Range", (p75th - p25th))); featureList.push_back(std::make_pair("First Order::Image Dimension", VImageDimension)); featureList.push_back(std::make_pair("First Order::Voxel Space", voxelSpace)); featureList.push_back(std::make_pair("First Order::Voxel Volume", voxelVolume)); } mitk::GIFFirstOrderStatistics::GIFFirstOrderStatistics() { SetShortName("fo"); SetLongName("first-order"); } mitk::GIFFirstOrderStatistics::FeatureListType mitk::GIFFirstOrderStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { InitializeQuantifier(image, mask); FeatureListType featureList; ParameterStruct params; params.MinimumIntensity = GetQuantifier()->GetMinimum(); params.MaximumIntensity = GetQuantifier()->GetMaximum(); params.Bins = GetQuantifier()->GetBins(); AccessByItk_3(image, CalculateFirstOrderStatistics, mask, featureList, params); return featureList; } mitk::GIFFirstOrderStatistics::FeatureNameListType mitk::GIFFirstOrderStatistics::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("First Order::Minimum"); featureList.push_back("First Order::Maximum"); featureList.push_back("First Order::Mean"); featureList.push_back("First Order::Variance"); featureList.push_back("First Order::Sum"); featureList.push_back("First Order::Median"); featureList.push_back("First Order::Standard deviation"); featureList.push_back("First Order::No. of Voxel"); return featureList; } void mitk::GIFFirstOrderStatistics::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Volume-Statistic", "calculates volume based features", us::Any()); AddQuantifierArguments(parser); } void mitk::GIFFirstOrderStatistics::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { InitializeQuantifierFromParameters(feature, maskNoNAN); MITK_INFO << "Start calculating first order features ...."; auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating first order features...."; } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp index 3ef0c99140..27f28103be 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelDistanceZone.cpp @@ -1,471 +1,469 @@ #include // MITK #include #include #include #include #include // ITK #include #include #include #include #include static void MatrixFeaturesTo(mitk::GreyLevelDistanceZoneFeatures features, std::string prefix, mitk::GIFGreyLevelDistanceZone::FeatureListType &featureList); mitk::GreyLevelDistanceZoneMatrixHolder::GreyLevelDistanceZoneMatrixHolder(double min, double max, int number, int maxSize) : m_MinimumRange(min), m_MaximumRange(max), m_NumberOfBins(number), m_MaximumSize(maxSize), m_NumerOfVoxels(0) { m_Matrix.resize(number, maxSize); m_Matrix.fill(0); m_Stepsize = (max - min) / (number); } int mitk::GreyLevelDistanceZoneMatrixHolder::IntensityToIndex(double intensity) { return std::floor((intensity - m_MinimumRange) / m_Stepsize); } double mitk::GreyLevelDistanceZoneMatrixHolder::IndexToMinIntensity(int index) { return m_MinimumRange + index * m_Stepsize; } double mitk::GreyLevelDistanceZoneMatrixHolder::IndexToMeanIntensity(int index) { return m_MinimumRange + (index+0.5) * m_Stepsize; } double mitk::GreyLevelDistanceZoneMatrixHolder::IndexToMaxIntensity(int index) { return m_MinimumRange + (index + 1) * m_Stepsize; } template int CalculateGlSZMatrix(itk::Image* itkImage, itk::Image* mask, itk::Image* distanceImage, std::vector > offsets, bool estimateLargestRegion, mitk::GreyLevelDistanceZoneMatrixHolder &holder) { typedef itk::Image ImageType; typedef itk::Image MaskImageType; typedef typename ImageType::IndexType IndexType; typedef itk::ImageRegionIteratorWithIndex ConstIterType; typedef itk::ImageRegionIteratorWithIndex ConstMaskIterType; auto region = mask->GetLargestPossibleRegion(); typename MaskImageType::RegionType newRegion; newRegion.SetSize(region.GetSize()); newRegion.SetIndex(region.GetIndex()); ConstIterType imageIter(itkImage, itkImage->GetLargestPossibleRegion()); ConstMaskIterType maskIter(mask, mask->GetLargestPossibleRegion()); typename MaskImageType::Pointer visitedImage = MaskImageType::New(); visitedImage->SetRegions(newRegion); visitedImage->Allocate(); visitedImage->FillBuffer(0); int largestRegion = 0; holder.m_NumberOfBins = 0; while (!maskIter.IsAtEnd()) { if (maskIter.Value() > 0 ) { auto startIntensityIndex = holder.IntensityToIndex(imageIter.Value()); std::vector indices; indices.push_back(maskIter.GetIndex()); unsigned int steps = 0; int smallestDistance = 500; while (indices.size() > 0) { auto currentIndex = indices.back(); indices.pop_back(); if (!region.IsInside(currentIndex)) { continue; } auto wasVisited = visitedImage->GetPixel(currentIndex); auto newIntensityIndex = holder.IntensityToIndex(itkImage->GetPixel(currentIndex)); auto isInMask = mask->GetPixel(currentIndex); if ((isInMask > 0) && (newIntensityIndex == startIntensityIndex) && (wasVisited < 1)) { ++(holder.m_NumerOfVoxels); smallestDistance = (smallestDistance > distanceImage->GetPixel(currentIndex)) ? distanceImage->GetPixel(currentIndex) : smallestDistance; ++steps; visitedImage->SetPixel(currentIndex, 1); for (auto offset : offsets) { auto newIndex = currentIndex + offset; indices.push_back(newIndex); newIndex = currentIndex - offset; indices.push_back(newIndex); } } } if (steps > 0) { largestRegion = std::max(steps, largestRegion); steps = std::min(steps, holder.m_MaximumSize); if (!estimateLargestRegion) { holder.m_Matrix(startIntensityIndex, smallestDistance-1) += 1; } } } ++imageIter; ++maskIter; } return largestRegion; } template void itkErode( itk::Image *sourceImage, mitk::Image::Pointer &resultImage, bool &segmentationRemaining) { typedef itk::Image ImageType; typedef itk::BinaryCrossStructuringElement CrossType; typedef typename itk::BinaryErodeImageFilter CrossErodeFilterType; CrossType cross; typename CrossType::SizeType size; size.Fill(1); cross.SetRadius(size); cross.CreateStructuringElement(); typename CrossErodeFilterType::Pointer erodeFilter = CrossErodeFilterType::New(); erodeFilter->SetKernel(cross); erodeFilter->SetInput(sourceImage); erodeFilter->SetErodeValue(1); erodeFilter->UpdateLargestPossibleRegion(); segmentationRemaining = false; itk::ImageRegionConstIterator iter(erodeFilter->GetOutput(), erodeFilter->GetOutput()->GetLargestPossibleRegion()); while ((!iter.IsAtEnd()) && (!segmentationRemaining)) { if (iter.Get() > 0) segmentationRemaining = true; ++iter; } mitk::CastToMitkImage(erodeFilter->GetOutput(), resultImage); } template void itkAdd( itk::Image *imageA, mitk::Image::Pointer mitkImageB, mitk::Image::Pointer &resultImage) { typedef itk::Image ImageType; typedef itk::AddImageFilter AddFilterType; typename ImageType::Pointer imageB; mitk::CastToItkImage(mitkImageB, imageB); typename AddFilterType::Pointer addFilter = AddFilterType::New(); addFilter->SetInput1(imageA); addFilter->SetInput2(imageB); addFilter->UpdateLargestPossibleRegion(); mitk::CastToMitkImage(addFilter->GetOutput(), resultImage); } void erode(mitk::Image::Pointer input, mitk::Image::Pointer &output, bool &segmentationRemaining) { AccessByItk_2(input, itkErode, output, segmentationRemaining); } void add(mitk::Image::Pointer inputA, mitk::Image::Pointer inputB, mitk::Image::Pointer &output) { AccessByItk_2(inputA, itkAdd, inputB, output); } void erodeAndAdd(mitk::Image::Pointer input, mitk::Image::Pointer& finalOutput, int &maxDistance) { mitk::Image::Pointer workingImage = input; mitk::Image::Pointer output = input; maxDistance = 0; bool segmentationRemaining = true; while (segmentationRemaining) { mitk::Image::Pointer erodedImage = mitk::Image::New(); mitk::Image::Pointer result = mitk::Image::New(); erode(workingImage, erodedImage, segmentationRemaining); workingImage = erodedImage; add(output, erodedImage, result); output = result; ++maxDistance; } finalOutput = output; } void static CalculateFeatures( mitk::GreyLevelDistanceZoneMatrixHolder &holder, mitk::GreyLevelDistanceZoneFeatures & results ) { auto SgzMatrix = holder.m_Matrix; auto pgzMatrix = holder.m_Matrix; auto pgMatrix = holder.m_Matrix; auto pzMatrix = holder.m_Matrix; double Ns = pgzMatrix.sum(); pgzMatrix /= Ns; pgMatrix.rowwise().normalize(); pzMatrix.colwise().normalize(); for (int i = 0; i < holder.m_NumberOfBins; ++i) for (int j = 0; j < holder.m_NumberOfBins; ++j) { if (pgzMatrix(i, j) != pgzMatrix(i, j)) pgzMatrix(i, j) = 0; if (pgMatrix(i, j) != pgMatrix(i, j)) pgMatrix(i, j) = 0; if (pzMatrix(i, j) != pzMatrix(i, j)) pzMatrix(i, j) = 0; } Eigen::VectorXd SgVector = SgzMatrix.rowwise().sum(); Eigen::VectorXd SzVector = SgzMatrix.colwise().sum(); for (int j = 0; j < SzVector.size(); ++j) { results.SmallDistanceEmphasis += SzVector(j) / (j+1) / (j+1); results.LargeDistanceEmphasis += SzVector(j) * (j + 1.0) * (j + 1.0); results.ZoneDistanceNonUniformity += SzVector(j) * SzVector(j); results.ZoneDistanceNoneUniformityNormalized += SzVector(j) * SzVector(j); } for (int i = 0; i < SgVector.size(); ++i) { results.LowGreyLevelEmphasis += SgVector(i) / (i + 1) / (i + 1); results.HighGreyLevelEmphasis += SgVector(i) * (i + 1) * (i + 1); results.GreyLevelNonUniformity += SgVector(i)*SgVector(i); results.GreyLevelNonUniformityNormalized += SgVector(i)*SgVector(i); } for (int i = 0; i < SgzMatrix.rows(); ++i) { for (int j = 0; j < SgzMatrix.cols(); ++j) { results.SmallDistanceLowGreyLevelEmphasis += SgzMatrix(i, j) / (i + 1) / (i + 1) / (j + 1) / (j + 1); results.SmallDistanceHighGreyLevelEmphasis += SgzMatrix(i, j) * (i + 1) * (i + 1) / (j + 1) / (j + 1); results.LargeDistanceLowGreyLevelEmphasis += SgzMatrix(i, j) / (i + 1) / (i + 1) * (j + 1.0) * (j + 1.0); results.LargeDistanceHighGreyLevelEmphasis += SgzMatrix(i, j) * (i + 1) * (i + 1) * (j + 1.0) * (j + 1.0); results.ZonePercentage += SgzMatrix(i, j); results.GreyLevelMean += (i + 1)*pgzMatrix(i, j); results.ZoneDistanceMean += (j + 1)*pgzMatrix(i, j); if (pgzMatrix(i, j) > 0) results.ZoneDistanceEntropy -= pgzMatrix(i, j) * std::log(pgzMatrix(i, j)) / std::log(2); } } for (int i = 0; i < SgzMatrix.rows(); ++i) { for (int j = 0; j < SgzMatrix.cols(); ++j) { results.GreyLevelVariance += (i + 1 - results.GreyLevelMean)*(i + 1 - results.GreyLevelMean)*pgzMatrix(i, j); results.ZoneDistanceVariance += (j + 1 - results.ZoneDistanceMean)*(j + 1 - results.ZoneDistanceMean)*pgzMatrix(i, j); } } results.SmallDistanceEmphasis /= Ns; results.LargeDistanceEmphasis /= Ns; results.LowGreyLevelEmphasis /= Ns; results.HighGreyLevelEmphasis /= Ns; results.SmallDistanceLowGreyLevelEmphasis /= Ns; results.SmallDistanceHighGreyLevelEmphasis /= Ns; results.LargeDistanceLowGreyLevelEmphasis /= Ns; results.LargeDistanceHighGreyLevelEmphasis /= Ns; results.GreyLevelNonUniformity /= Ns; results.GreyLevelNonUniformityNormalized /= Ns*Ns; results.ZoneDistanceNonUniformity /= Ns; results.ZoneDistanceNoneUniformityNormalized /= Ns*Ns; results.ZonePercentage = Ns / holder.m_NumerOfVoxels;// results.ZonePercentage; } template static void CalculateGreyLevelDistanceZoneFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFGreyLevelDistanceZone::FeatureListType & featureList, mitk::GIFGreyLevelDistanceZone::GIFGreyLevelDistanceZoneConfiguration config) { - typedef itk::Image ImageType; typedef itk::Image MaskType; - typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef itk::Neighborhood NeighborhoodType; typedef itk::Offset OffsetType; /////////////////////////////////////////////////////////////////////////////////////////////// double rangeMin = config.MinimumIntensity; double rangeMax = config.MaximumIntensity; int numberOfBins = config.Bins; int maximumDistance = 0; mitk::Image::Pointer mitkDistanceImage = mitk::Image::New(); erodeAndAdd(config.distanceMask, mitkDistanceImage, maximumDistance); typename MaskType::Pointer distanceImage = MaskType::New(); mitk::CastToItkImage(mitkDistanceImage, distanceImage); typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); //Find possible directions std::vector < itk::Offset > offsetVector; NeighborhoodType hood; hood.SetRadius(1); unsigned int centerIndex = hood.GetCenterNeighborhoodIndex(); OffsetType offset; for (unsigned int d = 0; d < centerIndex; d++) { offset = hood.GetOffset(d); bool useOffset = true; for (unsigned int i = 0; i < VImageDimension; ++i) { if ((config.direction == i + 2) && offset[i] != 0) { useOffset = false; } } if (useOffset) { offsetVector.push_back(offset); } } if (config.direction == 1) { offsetVector.clear(); offset[0] = 0; offset[1] = 0; offset[2] = 1; offsetVector.push_back(offset); } MITK_INFO << "Maximum Distance: " << maximumDistance; std::vector resultVector; mitk::GreyLevelDistanceZoneMatrixHolder holderOverall(rangeMin, rangeMax, numberOfBins, maximumDistance+1); mitk::GreyLevelDistanceZoneFeatures overallFeature; CalculateGlSZMatrix(itkImage, maskImage, distanceImage, offsetVector, false, holderOverall); CalculateFeatures(holderOverall, overallFeature); MatrixFeaturesTo(overallFeature, "Grey Level Distance Zone::", featureList); } static void MatrixFeaturesTo(mitk::GreyLevelDistanceZoneFeatures features, std::string prefix, mitk::GIFGreyLevelDistanceZone::FeatureListType &featureList) { featureList.push_back(std::make_pair(prefix + "Small Distance Emphasis", features.SmallDistanceEmphasis)); featureList.push_back(std::make_pair(prefix + "Large Distance Emphasis", features.LargeDistanceEmphasis)); featureList.push_back(std::make_pair(prefix + "Low Grey Level Emphasis", features.LowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "High Grey Level Emphasis", features.HighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Small Distance Low Grey Level Emphasis", features.SmallDistanceLowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Small Distance High Grey Level Emphasis", features.SmallDistanceHighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Large Distance Low Grey Level Emphasis", features.LargeDistanceLowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Large Distance High Grey Level Emphasis", features.LargeDistanceHighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity", features.GreyLevelNonUniformity)); featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity Normalized", features.GreyLevelNonUniformityNormalized)); featureList.push_back(std::make_pair(prefix + "Distance Size Non-Uniformity", features.ZoneDistanceNonUniformity)); featureList.push_back(std::make_pair(prefix + "Distance Size Non-Uniformity Normalized", features.ZoneDistanceNoneUniformityNormalized)); featureList.push_back(std::make_pair(prefix + "Zone Percentage", features.ZonePercentage)); featureList.push_back(std::make_pair(prefix + "Grey Level Mean", features.GreyLevelMean)); featureList.push_back(std::make_pair(prefix + "Grey Level Variance", features.GreyLevelVariance)); featureList.push_back(std::make_pair(prefix + "Zone Distance Mean", features.ZoneDistanceMean)); featureList.push_back(std::make_pair(prefix + "Zone Distance Variance", features.ZoneDistanceVariance)); featureList.push_back(std::make_pair(prefix + "Zone Distance Entropy", features.ZoneDistanceEntropy)); featureList.push_back(std::make_pair(prefix + "Grey Level Entropy", features.ZoneDistanceEntropy)); } mitk::GIFGreyLevelDistanceZone::GIFGreyLevelDistanceZone() { SetShortName("gldz"); SetLongName("distance-zone"); } mitk::GIFGreyLevelDistanceZone::FeatureListType mitk::GIFGreyLevelDistanceZone::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { InitializeQuantifier(image, mask); FeatureListType featureList; GIFGreyLevelDistanceZoneConfiguration config; config.direction = GetDirection(); if (GetMorphMask().IsNull()) { config.distanceMask = mask->Clone(); } else { config.distanceMask = GetMorphMask(); } config.MinimumIntensity = GetQuantifier()->GetMinimum(); config.MaximumIntensity = GetQuantifier()->GetMaximum(); config.Bins = GetQuantifier()->GetBins(); AccessByItk_3(image, CalculateGreyLevelDistanceZoneFeatures, mask, featureList, config); return featureList; } mitk::GIFGreyLevelDistanceZone::FeatureNameListType mitk::GIFGreyLevelDistanceZone::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFGreyLevelDistanceZone::AddArguments(mitkCommandLineParser &parser) { AddQuantifierArguments(parser); std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Grey Level Distance Zone", "Calculates the size zone based features.", us::Any()); } void mitk::GIFGreyLevelDistanceZone::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); std::string name = GetOptionPrefix(); if (parsedArgs.count(GetLongName())) { InitializeQuantifierFromParameters(feature, mask); MITK_INFO << "Start calculating Grey Level Distance Zone ...."; auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating Grey Level Distance Zone."; } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelSizeZone.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelSizeZone.cpp index 5a6e7ae016..4187b8ffbf 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelSizeZone.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGreyLevelSizeZone.cpp @@ -1,365 +1,364 @@ #include // MITK #include #include #include // ITK #include // STL static void MatrixFeaturesTo(mitk::GreyLevelSizeZoneFeatures features, std::string prefix, mitk::GIFGreyLevelSizeZone::FeatureListType &featureList); mitk::GreyLevelSizeZoneMatrixHolder::GreyLevelSizeZoneMatrixHolder(double min, double max, int number, int maxSize) : m_MinimumRange(min), m_MaximumRange(max), m_NumberOfBins(number), m_MaximumSize(maxSize) { m_Matrix.resize(number, maxSize); m_Matrix.fill(0); m_Stepsize = (max - min) / (number); } int mitk::GreyLevelSizeZoneMatrixHolder::IntensityToIndex(double intensity) { return std::floor((intensity - m_MinimumRange) / m_Stepsize); } double mitk::GreyLevelSizeZoneMatrixHolder::IndexToMinIntensity(int index) { return m_MinimumRange + index * m_Stepsize; } double mitk::GreyLevelSizeZoneMatrixHolder::IndexToMeanIntensity(int index) { return m_MinimumRange + (index+0.5) * m_Stepsize; } double mitk::GreyLevelSizeZoneMatrixHolder::IndexToMaxIntensity(int index) { return m_MinimumRange + (index + 1) * m_Stepsize; } template static int CalculateGlSZMatrix(itk::Image* itkImage, itk::Image* mask, std::vector > offsets, bool estimateLargestRegion, mitk::GreyLevelSizeZoneMatrixHolder &holder) { typedef itk::Image ImageType; typedef itk::Image MaskImageType; typedef typename ImageType::IndexType IndexType; typedef itk::ImageRegionIteratorWithIndex ConstIterType; typedef itk::ImageRegionIteratorWithIndex ConstMaskIterType; auto region = mask->GetLargestPossibleRegion(); typename MaskImageType::RegionType newRegion; newRegion.SetSize(region.GetSize()); newRegion.SetIndex(region.GetIndex()); ConstIterType imageIter(itkImage, itkImage->GetLargestPossibleRegion()); ConstMaskIterType maskIter(mask, mask->GetLargestPossibleRegion()); typename MaskImageType::Pointer visitedImage = MaskImageType::New(); visitedImage->SetRegions(newRegion); visitedImage->Allocate(); visitedImage->FillBuffer(0); int largestRegion = 0; while (!maskIter.IsAtEnd()) { if (maskIter.Value() > 0 ) { auto startIntensityIndex = holder.IntensityToIndex(imageIter.Value()); std::vector indices; indices.push_back(maskIter.GetIndex()); unsigned int steps = 0; while (indices.size() > 0) { auto currentIndex = indices.back(); indices.pop_back(); if (!region.IsInside(currentIndex)) { continue; } auto wasVisited = visitedImage->GetPixel(currentIndex); auto newIntensityIndex = holder.IntensityToIndex(itkImage->GetPixel(currentIndex)); auto isInMask = mask->GetPixel(currentIndex); if ((isInMask > 0) && (newIntensityIndex == startIntensityIndex) && (wasVisited < 1)) { ++steps; visitedImage->SetPixel(currentIndex, 1); for (auto offset : offsets) { auto newIndex = currentIndex + offset; indices.push_back(newIndex); newIndex = currentIndex - offset; indices.push_back(newIndex); } } } if (steps > 0) { largestRegion = std::max(steps, largestRegion); steps = std::min(steps, holder.m_MaximumSize); if (!estimateLargestRegion) { holder.m_Matrix(startIntensityIndex, steps - 1) += 1; } } } ++imageIter; ++maskIter; } return largestRegion; } static void CalculateFeatures( mitk::GreyLevelSizeZoneMatrixHolder &holder, mitk::GreyLevelSizeZoneFeatures & results ) { auto SgzMatrix = holder.m_Matrix; auto pgzMatrix = holder.m_Matrix; auto pgMatrix = holder.m_Matrix; auto pzMatrix = holder.m_Matrix; double Ns = pgzMatrix.sum(); pgzMatrix /= Ns; pgMatrix.rowwise().normalize(); pzMatrix.colwise().normalize(); for (int i = 0; i < holder.m_NumberOfBins; ++i) for (int j = 0; j < holder.m_NumberOfBins; ++j) { if (pgzMatrix(i, j) != pgzMatrix(i, j)) pgzMatrix(i, j) = 0; if (pgMatrix(i, j) != pgMatrix(i, j)) pgMatrix(i, j) = 0; if (pzMatrix(i, j) != pzMatrix(i, j)) pzMatrix(i, j) = 0; } Eigen::VectorXd SgVector = SgzMatrix.rowwise().sum(); Eigen::VectorXd SzVector = SgzMatrix.colwise().sum(); for (int j = 0; j < SzVector.size(); ++j) { results.SmallZoneEmphasis += SzVector(j) / (j + 1) / (j + 1); results.LargeZoneEmphasis += SzVector(j) * (j + 1.0) * (j + 1.0); results.ZoneSizeNonUniformity += SzVector(j) * SzVector(j); results.ZoneSizeNoneUniformityNormalized += SzVector(j) * SzVector(j); } for (int i = 0; i < SgVector.size(); ++i) { results.LowGreyLevelEmphasis += SgVector(i) / (i + 1) / (i + 1); results.HighGreyLevelEmphasis += SgVector(i) * (i + 1) * (i + 1); results.GreyLevelNonUniformity += SgVector(i)*SgVector(i); results.GreyLevelNonUniformityNormalized += SgVector(i)*SgVector(i); } for (int i = 0; i < SgzMatrix.rows(); ++i) { for (int j = 0; j < SgzMatrix.cols(); ++j) { results.SmallZoneLowGreyLevelEmphasis += SgzMatrix(i, j) / (i + 1) / (i + 1) / (j + 1) / (j + 1); results.SmallZoneHighGreyLevelEmphasis += SgzMatrix(i, j) * (i + 1) * (i + 1) / (j + 1) / (j + 1); results.LargeZoneLowGreyLevelEmphasis += SgzMatrix(i, j) / (i + 1) / (i + 1) * (j + 1.0) * (j + 1.0); results.LargeZoneHighGreyLevelEmphasis += SgzMatrix(i, j) * (i + 1) * (i + 1) * (j + 1.0) * (j + 1.0); results.ZonePercentage += SgzMatrix(i, j)*(j + 1); results.GreyLevelMean += (i + 1)*pgzMatrix(i, j); results.ZoneSizeMean += (j + 1)*pgzMatrix(i, j); if (pgzMatrix(i, j) > 0) results.ZoneSizeEntropy -= pgzMatrix(i, j) * std::log(pgzMatrix(i, j)) / std::log(2); } } for (int i = 0; i < SgzMatrix.rows(); ++i) { for (int j = 0; j < SgzMatrix.cols(); ++j) { results.GreyLevelVariance += (i + 1 - results.GreyLevelMean)*(i + 1 - results.GreyLevelMean)*pgzMatrix(i, j); results.ZoneSizeVariance += (j + 1 - results.ZoneSizeMean)*(j + 1 - results.ZoneSizeMean)*pgzMatrix(i, j); } } results.SmallZoneEmphasis /= Ns; results.LargeZoneEmphasis /= Ns; results.LowGreyLevelEmphasis /= Ns; results.HighGreyLevelEmphasis /= Ns; results.SmallZoneLowGreyLevelEmphasis /= Ns; results.SmallZoneHighGreyLevelEmphasis /= Ns; results.LargeZoneLowGreyLevelEmphasis /= Ns; results.LargeZoneHighGreyLevelEmphasis /= Ns; results.GreyLevelNonUniformity /= Ns; results.GreyLevelNonUniformityNormalized /= Ns*Ns; results.ZoneSizeNonUniformity /= Ns; results.ZoneSizeNoneUniformityNormalized /= Ns*Ns; results.ZonePercentage = Ns / results.ZonePercentage; } template static void CalculateGreyLevelSizeZoneFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFGreyLevelSizeZone::FeatureListType & featureList, mitk::GIFGreyLevelSizeZone::GIFGreyLevelSizeZoneConfiguration config) { - typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::Neighborhood NeighborhoodType; typedef itk::Offset OffsetType; /////////////////////////////////////////////////////////////////////////////////////////////// double rangeMin = config.MinimumIntensity; double rangeMax = config.MaximumIntensity; 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 (unsigned int i = 0; i < VImageDimension; ++i) { if ((config.direction == i + 2) && offset[i] != 0) { useOffset = false; } } if (useOffset) { offsetVector.push_back(offset); MITK_INFO << offset; } } if (config.direction == 1) { offsetVector.clear(); offset[0] = 0; offset[1] = 0; offset[2] = 1; offsetVector.push_back(offset); } std::vector resultVector; mitk::GreyLevelSizeZoneMatrixHolder tmpHolder(rangeMin, rangeMax, numberOfBins, 3); int largestRegion = CalculateGlSZMatrix(itkImage, maskImage, offsetVector, true, tmpHolder); mitk::GreyLevelSizeZoneMatrixHolder holderOverall(rangeMin, rangeMax, numberOfBins,largestRegion); mitk::GreyLevelSizeZoneFeatures overallFeature; CalculateGlSZMatrix(itkImage, maskImage, offsetVector, false, holderOverall); CalculateFeatures(holderOverall, overallFeature); MatrixFeaturesTo(overallFeature, "Grey Level Size Zone::", featureList); } static void MatrixFeaturesTo(mitk::GreyLevelSizeZoneFeatures features, std::string prefix, mitk::GIFGreyLevelSizeZone::FeatureListType &featureList) { featureList.push_back(std::make_pair(prefix + "Small Zone Emphasis", features.SmallZoneEmphasis)); featureList.push_back(std::make_pair(prefix + "Large Zone Emphasis", features.LargeZoneEmphasis)); featureList.push_back(std::make_pair(prefix + "Low Grey Level Emphasis", features.LowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "High Grey Level Emphasis", features.HighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Small Zone Low Grey Level Emphasis", features.SmallZoneLowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Small Zone High Grey Level Emphasis", features.SmallZoneHighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Large Zone Low Grey Level Emphasis", features.LargeZoneLowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Large Zone High Grey Level Emphasis", features.LargeZoneHighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity", features.GreyLevelNonUniformity)); featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity Normalized", features.GreyLevelNonUniformityNormalized)); featureList.push_back(std::make_pair(prefix + "Zone Size Non-Uniformity", features.ZoneSizeNonUniformity)); featureList.push_back(std::make_pair(prefix + "Zone Size Non-Uniformity Normalized", features.ZoneSizeNoneUniformityNormalized)); featureList.push_back(std::make_pair(prefix + "Zone Percentage", features.ZonePercentage)); featureList.push_back(std::make_pair(prefix + "Grey Level Mean", features.GreyLevelMean)); featureList.push_back(std::make_pair(prefix + "Grey Level Variance", features.GreyLevelVariance)); featureList.push_back(std::make_pair(prefix + "Zone Size Mean", features.ZoneSizeMean)); featureList.push_back(std::make_pair(prefix + "Zone Size Variance", features.ZoneSizeVariance)); featureList.push_back(std::make_pair(prefix + "Zone Size Entropy", features.ZoneSizeEntropy)); } mitk::GIFGreyLevelSizeZone::GIFGreyLevelSizeZone() { SetShortName("glsz"); SetLongName("grey-level-sizezone"); } mitk::GIFGreyLevelSizeZone::FeatureListType mitk::GIFGreyLevelSizeZone::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { InitializeQuantifier(image, mask); FeatureListType featureList; GIFGreyLevelSizeZoneConfiguration config; config.direction = GetDirection(); config.MinimumIntensity = GetQuantifier()->GetMinimum(); config.MaximumIntensity = GetQuantifier()->GetMaximum(); config.Bins = GetQuantifier()->GetBins(); AccessByItk_3(image, CalculateGreyLevelSizeZoneFeatures, mask, featureList, config); return featureList; } mitk::GIFGreyLevelSizeZone::FeatureNameListType mitk::GIFGreyLevelSizeZone::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFGreyLevelSizeZone::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Grey Level Size Zone", "Calculates the size zone based features.", us::Any()); AddQuantifierArguments(parser); } void mitk::GIFGreyLevelSizeZone::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())) { InitializeQuantifierFromParameters(feature, maskNoNAN); MITK_INFO << "Start calculating Grey leve size zone ..."; auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating Grey level size zone ..."; } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFImageDescriptionFeatures.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFImageDescriptionFeatures.cpp index 88a07ee7c0..f1ac63c92f 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFImageDescriptionFeatures.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFImageDescriptionFeatures.cpp @@ -1,178 +1,175 @@ /*=================================================================== 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 #include // STL #include template static void CalculateFirstOrderStatistics(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFImageDescriptionFeatures::FeatureListType & featureList) { 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); unsigned int imageDimensionX = itkImage->GetLargestPossibleRegion().GetSize()[0]; unsigned int imageDimensionY = itkImage->GetLargestPossibleRegion().GetSize()[1]; unsigned int imageDimensionZ = itkImage->GetLargestPossibleRegion().GetSize()[2]; double imageVoxelSpacingX = itkImage->GetSpacing()[0]; double imageVoxelSpacingY = itkImage->GetSpacing()[1]; double imageVoxelSpacingZ = itkImage->GetSpacing()[2]; typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); double imageMinimum = minMaxComputer->GetMinimum(); double imageMaximum = minMaxComputer->GetMaximum(); unsigned int maskDimensionX = maskImage->GetLargestPossibleRegion().GetSize()[0]; unsigned int maskDimensionY = maskImage->GetLargestPossibleRegion().GetSize()[1]; unsigned int maskDimensionZ = maskImage->GetLargestPossibleRegion().GetSize()[2]; double maskVoxelSpacingX = maskImage->GetSpacing()[0]; double maskVoxelSpacingY = maskImage->GetSpacing()[1]; double maskVoxelSpacingZ = maskImage->GetSpacing()[2]; unsigned int voxelCount = 0; unsigned int maskVoxelCount = 0; double maskMinimum = imageMaximum; double maskMaximum = imageMinimum; double imageMean = 0; double maskMean = 0; unsigned int maskMinimumX = maskDimensionX; unsigned int maskMaximumX = 0; unsigned int maskMinimumY = maskDimensionY; unsigned int maskMaximumY = 0; unsigned int maskMinimumZ = maskDimensionZ; unsigned int maskMaximumZ = 0; itk::ImageRegionConstIteratorWithIndex imIter(itkImage, itkImage->GetLargestPossibleRegion()); itk::ImageRegionConstIteratorWithIndex maIter(maskImage, maskImage->GetLargestPossibleRegion()); while (!imIter.IsAtEnd()) { auto pixelValue = imIter.Get(); if (maIter.Get() > 0) { ++maskVoxelCount; maskMean += pixelValue; maskMinimum = (maskMinimum > pixelValue) ? pixelValue : maskMinimum; maskMaximum = (maskMaximum < pixelValue) ? pixelValue : maskMaximum; maskMinimumX = (maskMinimumX > imIter.GetIndex()[0]) ? imIter.GetIndex()[0] : maskMinimumX; maskMaximumX = (maskMaximumX < imIter.GetIndex()[0]) ? imIter.GetIndex()[0] : maskMaximumX; maskMinimumY = (maskMinimumY > imIter.GetIndex()[1]) ? imIter.GetIndex()[1] : maskMinimumY; maskMaximumY = (maskMaximumY < imIter.GetIndex()[1]) ? imIter.GetIndex()[1] : maskMaximumY; maskMinimumZ = (maskMinimumZ > imIter.GetIndex()[2]) ? imIter.GetIndex()[2] : maskMinimumZ; maskMaximumZ = (maskMaximumZ < imIter.GetIndex()[2]) ? imIter.GetIndex()[2] : maskMaximumZ; } ++voxelCount; imageMean += pixelValue; ++imIter; ++maIter; } imageMean /= voxelCount; maskMean /= maskVoxelCount; featureList.push_back(std::make_pair("Diagnostic (Image)::Dimension X", imageDimensionX)); featureList.push_back(std::make_pair("Diagnostic (Image)::Dimension Y", imageDimensionY)); featureList.push_back(std::make_pair("Diagnostic (Image)::Dimension Z", imageDimensionZ)); featureList.push_back(std::make_pair("Diagnostic (Image)::Spacing X", imageVoxelSpacingX)); featureList.push_back(std::make_pair("Diagnostic (Image)::Spacing Y", imageVoxelSpacingY)); featureList.push_back(std::make_pair("Diagnostic (Image)::Spacing Z", imageVoxelSpacingZ)); featureList.push_back(std::make_pair("Diagnostic (Image)::Mean intensity", imageMean)); featureList.push_back(std::make_pair("Diagnostic (Image)::Minimum intensity", imageMinimum)); featureList.push_back(std::make_pair("Diagnostic (Image)::Maximum intensity", imageMaximum)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Dimension X", maskDimensionX)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Dimension Y", maskDimensionY)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Dimension Z", maskDimensionZ)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Mask bounding box X", maskMaximumX - maskMinimumX + 1)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Mask bounding box Y", maskMaximumY - maskMinimumY + 1)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Mask bounding box Z", maskMaximumZ - maskMinimumZ + 1)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Spacing X", maskVoxelSpacingX)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Spacing Y", maskVoxelSpacingY)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Spacing Z", maskVoxelSpacingZ)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Voxel Count ", maskVoxelCount)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Mean intensity", maskMean)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Minimum intensity", maskMinimum)); featureList.push_back(std::make_pair("Diagnostic (Mask)::Maximum intensity", maskMaximum)); } mitk::GIFImageDescriptionFeatures::GIFImageDescriptionFeatures() { SetShortName("id"); SetLongName("image-diagnostic"); } mitk::GIFImageDescriptionFeatures::FeatureListType mitk::GIFImageDescriptionFeatures::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; AccessByItk_2(image, CalculateFirstOrderStatistics, mask, featureList); return featureList; } mitk::GIFImageDescriptionFeatures::FeatureNameListType mitk::GIFImageDescriptionFeatures::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFImageDescriptionFeatures::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Image Description", "calculates image description features", us::Any()); } void mitk::GIFImageDescriptionFeatures::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { MITK_INFO << "Start calculating image description features...."; auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating image description features...."; } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFIntensityVolumeHistogramFeatures.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFIntensityVolumeHistogramFeatures.cpp index c94b39dfc2..0badfacbbb 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFIntensityVolumeHistogramFeatures.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFIntensityVolumeHistogramFeatures.cpp @@ -1,158 +1,157 @@ /*=================================================================== 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 // ITK #include #include // STL #include template static void CalculateIntensityPeak(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::IntensityQuantifier::Pointer quantifier, mitk::GIFIntensityVolumeHistogramFeatures::FeatureListType & featureList) { - int bins = 1000; typedef itk::Image ImageType; typedef itk::Image MaskType; typename MaskType::Pointer itkMask = MaskType::New(); mitk::CastToItkImage(mask, itkMask); itk::ImageRegionConstIterator iter(itkImage, itkImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator iterMask(itkMask, itkMask->GetLargestPossibleRegion()); MITK_INFO << "Quantification: " << quantifier->GetMinimum() << " to " << quantifier->GetMaximum() << " with " << quantifier->GetBins()<< " bins"; iter.GoToBegin(); iterMask.GoToBegin(); std::vector hist; hist.resize(quantifier->GetBins() , 0); int count = 0; while (!iter.IsAtEnd()) { if (iterMask.Get() > 0) { double value = iter.Get(); //std::size_t index = std::floor((value - minimum) / (maximum - minimum) * (bins-1)); std::size_t index = quantifier->IntensityToIndex(value); ++count; hist[index] += 1.0;// / count; } ++iterMask; ++iter; } bool notFoundIntenstiy010 = true; bool notFoundIntenstiy090 = true; double intensity010 = -1; double intensity090 = -1; double fraction = 0; double auc = 0; bool firstRound = true; for (int i = quantifier->GetBins()-1; i >= 0; --i) { hist[i] /= count; hist[i] += fraction; fraction = hist[i]; if (!firstRound) { auc += 0.5 * (hist[i] + hist[i+1]) / (quantifier->GetBins()-1); } firstRound = false; if (notFoundIntenstiy010 && fraction > 0.1) { intensity010 = quantifier->IndexToMeanIntensity(i + 1); notFoundIntenstiy010 = false; } if (notFoundIntenstiy090 && fraction > 0.9) { intensity090 = quantifier->IndexToMeanIntensity(i + 1); notFoundIntenstiy090 = false; } } unsigned int index010 = std::ceil(quantifier->GetBins() * 0.1); unsigned int index090 = std::floor(quantifier->GetBins() * 0.9); featureList.push_back(std::make_pair("Intensity Volume Histogram::Volume fration at 0.10 intensity", hist[index010])); featureList.push_back(std::make_pair("Intensity Volume Histogram::Volume fration at 0.90 intensity", hist[index090])); featureList.push_back(std::make_pair("Intensity Volume Histogram::Intensity at 0.10 volume", intensity010)); featureList.push_back(std::make_pair("Intensity Volume Histogram::Intensity at 0.90 volume", intensity090)); featureList.push_back(std::make_pair("Intensity Volume Histogram::Difference volume fration at 0.10 and 0.90 intensity", std::abs(hist[index010] - hist[index090]))); featureList.push_back(std::make_pair("Intensity Volume Histogram::Difference intensity at 0.10 and 0.90 volume", std::abs(intensity090 - intensity010))); featureList.push_back(std::make_pair("Intensity Volume Histogram::Area under IVH curve", auc)); //featureList.push_back(std::make_pair("Local Intensity Global Intensity Peak", globalPeakValue)); } mitk::GIFIntensityVolumeHistogramFeatures::GIFIntensityVolumeHistogramFeatures() { SetLongName("intensity-volume-histogram"); SetShortName("ivoh"); } mitk::GIFIntensityVolumeHistogramFeatures::FeatureListType mitk::GIFIntensityVolumeHistogramFeatures::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { InitializeQuantifier(image, mask, 1000); FeatureListType featureList; AccessByItk_3(image, CalculateIntensityPeak, mask, GetQuantifier(), featureList); return featureList; } mitk::GIFIntensityVolumeHistogramFeatures::FeatureNameListType mitk::GIFIntensityVolumeHistogramFeatures::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFIntensityVolumeHistogramFeatures::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Local Intensity", "calculates local intensity based features", us::Any()); AddQuantifierArguments(parser); } void -mitk::GIFIntensityVolumeHistogramFeatures::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNaN, FeatureListType &featureList) +mitk::GIFIntensityVolumeHistogramFeatures::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &, FeatureListType &featureList) { InitializeQuantifierFromParameters(feature, mask, 1000); auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { MITK_INFO << "Start calculating local intensity features ...."; auto localResults = this->CalculateFeatures(feature, mask); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating local intensity features...."; } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFLocalIntensity.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFLocalIntensity.cpp index ac5e25cb30..c943a3e828 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFLocalIntensity.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFLocalIntensity.cpp @@ -1,159 +1,159 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include #include #include // ITK #include #include // STL #include template static void CalculateIntensityPeak(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFLocalIntensity::FeatureListType & featureList, double range) { typedef itk::Image ImageType; typedef itk::Image MaskType; typename MaskType::Pointer itkMask = MaskType::New(); mitk::CastToItkImage(mask, itkMask); double minimumSpacing = std::numeric_limits::max(); itkImage->GetSpacing(); - for (int i = 0; i < VImageDimension; ++i) + for (unsigned int i = 0; i < VImageDimension; ++i) { minimumSpacing = (minimumSpacing < itkImage->GetSpacing()[i]) ? minimumSpacing : itkImage->GetSpacing()[i]; } - ImageType::SizeType regionSize; + typename ImageType::SizeType regionSize; int offset = std::ceil(range / minimumSpacing); regionSize.Fill(offset); itk::NeighborhoodIterator iter(regionSize, itkImage, itkImage->GetLargestPossibleRegion()); itk::NeighborhoodIterator iterMask(regionSize, itkMask, itkMask->GetLargestPossibleRegion()); - ImageType::PointType origin; - ImageType::PointType localPoint; + typename ImageType::PointType origin; + typename ImageType::PointType localPoint; itk::Index index; double tmpPeakValue; double globalPeakValue = 0; double localPeakValue = 0; TPixel localMaximum = 0; int count = 0; while (!iter.IsAtEnd()) { if (iterMask.GetCenterPixel() > 0) { tmpPeakValue = 0; count = 0; index = iter.GetIndex(); itkImage->TransformIndexToPhysicalPoint(index, origin); - for (int i = 0; i < iter.Size(); ++i) + for (itk::SizeValueType i = 0; i < iter.Size(); ++i) { itkImage->TransformIndexToPhysicalPoint(iter.GetIndex(i), localPoint); double dist = origin.EuclideanDistanceTo(localPoint); if (dist < 6.2) { if (iter.IndexInBounds(i)) { tmpPeakValue += iter.GetPixel(i); ++count; } } } tmpPeakValue /= count; globalPeakValue = std::max(tmpPeakValue, globalPeakValue); if (localMaximum == iter.GetCenterPixel()) { localPeakValue = std::max(tmpPeakValue,localPeakValue); } else if (localMaximum < iter.GetCenterPixel()) { localMaximum = iter.GetCenterPixel(); localPeakValue = tmpPeakValue; } } ++iterMask; ++iter; } featureList.push_back(std::make_pair("Local Intensity::Local Intensity Peak", localPeakValue)); featureList.push_back(std::make_pair("Local Intensity::Global Intensity Peak", globalPeakValue)); } mitk::GIFLocalIntensity::GIFLocalIntensity() : m_Range(6.2) { SetLongName("local-intensity"); SetShortName("loci"); } mitk::GIFLocalIntensity::FeatureListType mitk::GIFLocalIntensity::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; if (image->GetDimension() < 3) { return featureList; } double range = GetRange(); AccessByItk_3(image, CalculateIntensityPeak, mask, featureList, range); return featureList; } mitk::GIFLocalIntensity::FeatureNameListType mitk::GIFLocalIntensity::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFLocalIntensity::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Local Intensity", "calculates local intensity based features", us::Any()); parser.addArgument(name + "::range", name+"::range", mitkCommandLineParser::Float, "Range for the local intensity", "Give the range that should be used for the local intensity in mm", us::Any()); } void mitk::GIFLocalIntensity::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &, FeatureListType &featureList) { std::string name = GetOptionPrefix(); auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { if (parsedArgs.count(name + "::range")) { double range = us::any_cast(parsedArgs[name + "::range"]); this->SetRange(range); } MITK_INFO << "Start calculating local intensity features ...."; auto localResults = this->CalculateFeatures(feature, mask); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating local intensity features...."; } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyToneDifferenceFeatures.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyToneDifferenceFeatures.cpp index 4a19d344b5..ca63758f61 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyToneDifferenceFeatures.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyToneDifferenceFeatures.cpp @@ -1,202 +1,202 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include #include #include // ITK #include #include // STL #include template static void CalculateIntensityPeak(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::NGTDFeatures params, mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::FeatureListType & featureList) { typedef itk::Image ImageType; typedef itk::Image MaskType; typename MaskType::Pointer itkMask = MaskType::New(); mitk::CastToItkImage(mask, itkMask); - ImageType::SizeType regionSize; + typename ImageType::SizeType regionSize; regionSize.Fill(params.Range); itk::NeighborhoodIterator iter(regionSize, itkImage, itkImage->GetLargestPossibleRegion()); itk::NeighborhoodIterator iterMask(regionSize, itkMask, itkMask->GetLargestPossibleRegion()); std::vector pVector; std::vector sVector; pVector.resize(params.quantifier->GetBins(), 0); sVector.resize(params.quantifier->GetBins(), 0); int count = 0; while (!iter.IsAtEnd()) { if (iterMask.GetCenterPixel() > 0) { int localCount = 0; double localMean = 0; unsigned int localIndex = params.quantifier->IntensityToIndex(iter.GetCenterPixel()); - for (int i = 0; i < iter.Size(); ++i) + for (itk::SizeValueType i = 0; i < iter.Size(); ++i) { if (i == (iter.Size() / 2)) continue; if (iterMask.GetPixel(i) > 0) { ++localCount; localMean += params.quantifier->IntensityToIndex(iter.GetPixel(i)) + 1; } } if (localCount > 0) { localMean /= localCount; } localMean = std::abs(localIndex + 1 - localMean); pVector[localIndex] += 1; sVector[localIndex] += localMean; ++count; } ++iterMask; ++iter; } unsigned int Ngp = 0; for (unsigned int i = 0; i < params.quantifier->GetBins(); ++i) { if (pVector[i] > 0.1) { ++Ngp; } pVector[i] /= count; } double sumS = 0; double sumStimesP = 0; double contrastA = 0; double busynessA = 0; double complexity = 0; double strengthA = 0; for (unsigned int i = 0; i < params.quantifier->GetBins(); ++i) { sumS += sVector[i]; sumStimesP += pVector[i] * sVector[i]; for (unsigned int j = 0; j < params.quantifier->GetBins(); ++j) { double iMinusj = 1.0*i - 1.0*j; contrastA += pVector[i] * pVector[j] * iMinusj*iMinusj; if ((pVector[i] > 0) && (pVector[j] > 0)) { busynessA += std::abs((i + 1.0)*pVector[i] - (j + 1.0)*pVector[j]); complexity += std::abs(iMinusj)*(pVector[i] * sVector[i] + pVector[j] * sVector[j]) / (pVector[i] + pVector[j]); strengthA += (pVector[i] + pVector[j])*iMinusj*iMinusj; } } } double coarsness = 1.0 / std::min(sumStimesP, 1000000); double contrast = 0; double busyness = 0; if (Ngp > 1) { contrast = contrastA / Ngp / (Ngp - 1) / count * sumS; busyness = sumStimesP / busynessA; } complexity /= count; double strength = 0; if (sumS > 0) { strength = strengthA / sumS; } featureList.push_back(std::make_pair("Neighbourhood Grey Tone Difference::Coarsness", coarsness)); featureList.push_back(std::make_pair("Neighbourhood Grey Tone Difference::Contrast", contrast)); featureList.push_back(std::make_pair("Neighbourhood Grey Tone Difference::Busyness", busyness)); featureList.push_back(std::make_pair("Neighbourhood Grey Tone Difference::Complexity", complexity)); featureList.push_back(std::make_pair("Neighbourhood Grey Tone Difference::Strength", strength)); } mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::GIFNeighbourhoodGreyToneDifferenceFeatures() : m_Range(1) { SetLongName("neighbourhood-grey-tone-difference"); SetShortName("ngtd"); } mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::FeatureListType mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; InitializeQuantifierFromParameters(image, mask); NGTDFeatures params; params.Range = GetRange(); params.quantifier = GetQuantifier(); AccessByItk_3(image, CalculateIntensityPeak, mask, params, featureList); return featureList; } mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::FeatureNameListType mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::AddArguments(mitkCommandLineParser &parser) { AddQuantifierArguments(parser); std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Neighbourhood Grey Tone Difference", "calculates Neighborhood Grey Tone based features", us::Any()); parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::Int, "Range for the local intensity", "Give the range that should be used for the local intensity in mm", us::Any()); } void -mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &maskNoNaN, FeatureListType &featureList) +mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &, FeatureListType &featureList) { InitializeQuantifierFromParameters(feature, mask); std::string name = GetOptionPrefix(); auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { MITK_INFO << "Start calculating Neighbourhood Grey Tone Difference features ...."; if (parsedArgs.count(name + "::range")) { int range = us::any_cast(parsedArgs[name + "::range"]); this->SetRange(range); } auto localResults = this->CalculateFeatures(feature, mask); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating Neighbourhood Grey Tone Difference features...."; } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricDensityStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricDensityStatistics.cpp index a0151b1f42..42b314f95e 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricDensityStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricDensityStatistics.cpp @@ -1,515 +1,511 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include #include #include // ITK #include #include #include #include // VTK #include #include #include #include #include #include #include #include // STL #include #include // Eigen #include template void CalculateVolumeDensityStatistic(itk::Image* itkImage, mitk::Image::Pointer mask, double volume, mitk::GIFVolumetricDensityStatistics::FeatureListType & featureList) { typedef itk::Image ImageType; typedef itk::Image MaskType; - typedef itk::LabelGeometryImageFilter FilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); itk::ImageRegionConstIteratorWithIndex imgA(itkImage, itkImage->GetLargestPossibleRegion()); itk::ImageRegionConstIteratorWithIndex imgB(itkImage, itkImage->GetLargestPossibleRegion()); itk::ImageRegionConstIteratorWithIndex maskA(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIteratorWithIndex maskB(maskImage, maskImage->GetLargestPossibleRegion()); double moranA = 0; double moranB = 0; double geary = 0; double Nv = 0; double w_ij = 0; double mean = 0; - ImageType::PointType pointA; - ImageType::PointType pointB; + typename ImageType::PointType pointA; + typename ImageType::PointType pointB; while (!imgA.IsAtEnd()) { if (maskA.Get() > 0) { Nv += 1; mean += imgA.Get(); } ++imgA; ++maskA; } mean /= Nv; imgA.GoToBegin(); maskA.GoToBegin(); while (!imgA.IsAtEnd()) { if (maskA.Get() > 0) { imgB.GoToBegin(); maskB.GoToBegin(); while (!imgB.IsAtEnd()) { if ((imgA.GetIndex() == imgB.GetIndex()) || (maskB.Get() < 1)) { ++imgB; ++maskB; continue; } itkImage->TransformIndexToPhysicalPoint(maskA.GetIndex(), pointA); itkImage->TransformIndexToPhysicalPoint(maskB.GetIndex(), pointB); double w = 1 / pointA.EuclideanDistanceTo(pointB); moranA += w*(imgA.Get() - mean)* (imgB.Get() - mean); geary += w * (imgA.Get() - imgB.Get()) * (imgA.Get() - imgB.Get()); w_ij += w; ++imgB; ++maskB; } moranB += (imgA.Get() - mean)* (imgA.Get() - mean); } ++imgA; ++maskA; } featureList.push_back(std::make_pair("Morphological Density::Volume integrated intensity", volume* mean)); featureList.push_back(std::make_pair("Morphological Density::Volume Moran's I index", Nv / w_ij * moranA / moranB)); featureList.push_back(std::make_pair("Morphological Density::Volume Geary's C measure", ( Nv -1 ) / 2 / w_ij * geary/ moranB)); } void calculateMOBB(vtkPointSet *pointset, double &volume, double &surface) { - auto cell = pointset->GetCell(0); volume = std::numeric_limits::max(); for (int cellID = 0; cellID < pointset->GetNumberOfCells(); ++cellID) { auto cell = pointset->GetCell(cellID); for (int edgeID = 0; edgeID < 3; ++edgeID) { auto edge = cell->GetEdge(edgeID); double pA[3], pB[3]; double pAA[3], pBB[3]; vtkSmartPointer transform = vtkSmartPointer::New(); transform->PostMultiply(); pointset->GetPoint(edge->GetPointId(0), pA); pointset->GetPoint(edge->GetPointId(1), pB); double angleZ = std::atan2((- pA[2] + pB[2]) ,(pA[1] - pB[1])); angleZ *= 180 / vnl_math::pi; if (pA[2] == pB[2]) angleZ = 0; transform->RotateX(angleZ); transform->TransformPoint(pA, pAA); transform->TransformPoint(pB, pBB); double angleY = std::atan2((pAA[1] -pBB[1]) ,-(pAA[0] - pBB[0])); angleY *= 180 / vnl_math::pi; if (pAA[1] == pBB[1]) angleY = 0; transform->RotateZ(angleY); double p0[3]; pointset->GetPoint(edge->GetPointId(0), p0); double curMinX = std::numeric_limits::max(); double curMaxX = std::numeric_limits::lowest(); double curMinY = std::numeric_limits::max(); double curMaxY = std::numeric_limits::lowest(); double curMinZ = std::numeric_limits::max(); double curMaxZ = std::numeric_limits::lowest(); for (int pointID = 0; pointID < pointset->GetNumberOfPoints(); ++pointID) { double p[3]; double p2[3]; pointset->GetPoint(pointID, p); p[0] -= p0[0]; p[1] -= p0[1]; p[2] -= p0[2]; transform->TransformPoint(p, p2); curMinX = std::min(p2[0], curMinX); curMaxX = std::max(p2[0], curMaxX); curMinY = std::min(p2[1], curMinY); curMaxY = std::max(p2[1], curMaxY); curMinZ = std::min(p2[2], curMinZ); curMaxZ = std::max(p2[2], curMaxZ); //std::cout << pointID << " (" << p[0] << "|" << p[1] << "|" << p[2] << ") (" << p2[0] << "|" << p2[1] << "|" << p2[2] << ")" << std::endl; } if ((curMaxX - curMinX)*(curMaxY - curMinY)*(curMaxZ - curMinZ) < volume) { volume = (curMaxX - curMinX)*(curMaxY - curMinY)*(curMaxZ - curMinZ); surface = (curMaxX - curMinX)*(curMaxX - curMinX) + (curMaxY - curMinY)*(curMaxY - curMinY) + (curMaxZ - curMinZ)*(curMaxZ - curMinZ); surface *= 2; } } } } void calculateMEE(vtkPointSet *pointset, double &vol, double &surf, double tolerance=0.0000001) { // Inspired by https://github.com/smdabdoub/ProkaryMetrics/blob/master/calc/fitting.py int numberOfPoints = pointset->GetNumberOfPoints(); int dimension = 3; Eigen::MatrixXd points(3, numberOfPoints); Eigen::MatrixXd Q(3+1, numberOfPoints); double p[3]; for (int i = 0; i < numberOfPoints; ++i) { pointset->GetPoint(i, p); points(0, i) = p[0]; points(1, i) = p[1]; points(2, i) = p[2]; Q(0, i) = p[0]; Q(1, i) = p[1]; Q(2, i) = p[2]; Q(3, i) = p[3]; } int count = 1; double error = 1; Eigen::VectorXd u_vector(numberOfPoints); u_vector.fill(1.0 / numberOfPoints); Eigen::DiagonalMatrix u = u_vector.asDiagonal(); Eigen::VectorXd ones(dimension + 1); ones.fill(1); Eigen::MatrixXd Ones = ones.asDiagonal(); // Khachiyan Algorithm while (error > tolerance) { auto Qt = Q.transpose(); Eigen::MatrixXd X = Q*u*Qt; Eigen::FullPivHouseholderQR qr(X); Eigen::MatrixXd Xi = qr.solve(Ones); Eigen::MatrixXd M = Qt * Xi * Q; double maximumValue = M(0, 0); int maximumPosition = 0; for (int i = 0; i < numberOfPoints; ++i) { if (maximumValue < M(i, i)) { maximumValue = M(i, i); maximumPosition = i; } } double stepsize = (maximumValue - dimension - 1) / ((dimension + 1) * (maximumValue - 1)); Eigen::DiagonalMatrix new_u = (1.0 - stepsize) * u; new_u.diagonal()[maximumPosition] = (new_u.diagonal())(maximumPosition) + stepsize; ++count; error = (new_u.diagonal() - u.diagonal()).norm(); u.diagonal() = new_u.diagonal(); } // U = u Eigen::MatrixXd Ai = points * u * points.transpose() - points * u *(points * u).transpose(); Eigen::FullPivHouseholderQR qr(Ai); Eigen::VectorXd ones2(dimension); ones2.fill(1); Eigen::MatrixXd Ones2 = ones2.asDiagonal(); Eigen::MatrixXd A = qr.solve(Ones2)*1.0/dimension; Eigen::JacobiSVD svd(A); double c = 1 / sqrt(svd.singularValues()[0]); double b = 1 / sqrt(svd.singularValues()[1]); double a = 1 / sqrt(svd.singularValues()[2]); double V = 4 * vnl_math::pi*a*b*c / 3; double ad_mvee= 0; double alpha = std::sqrt(1 - b*b / a / a); double beta = std::sqrt(1 - c*c / a / a); for (int i = 0; i < 20; ++i) { ad_mvee += 4 * vnl_math::pi*a*b*(alpha*alpha + beta*beta) / (2 * alpha*beta) * (std::pow(alpha*beta, i)) / (1 - 4 * i*i); } vol = V; surf = ad_mvee; } mitk::GIFVolumetricDensityStatistics::FeatureListType mitk::GIFVolumetricDensityStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; if (image->GetDimension() < 3) { return featureList; } vtkSmartPointer mesher = vtkSmartPointer::New(); vtkSmartPointer stats = vtkSmartPointer::New(); vtkSmartPointer stats2 = vtkSmartPointer::New(); mesher->SetInputData(mask->GetVtkImageData()); mesher->SetValue(0, 0.5); stats->SetInputConnection(mesher->GetOutputPort()); stats->Update(); vtkSmartPointer delaunay = vtkSmartPointer< vtkDelaunay3D >::New(); delaunay->SetInputConnection(mesher->GetOutputPort()); delaunay->SetAlpha(0); delaunay->Update(); vtkSmartPointer geometryFilter = vtkSmartPointer::New(); geometryFilter->SetInputConnection(delaunay->GetOutputPort()); geometryFilter->Update(); stats2->SetInputConnection(geometryFilter->GetOutputPort()); stats2->Update(); - auto data = mesher->GetOutput()->GetPointData(); - double vol_mvee; double surf_mvee; calculateMEE(mesher->GetOutput(), vol_mvee, surf_mvee); double vol_mobb; double surf_mobb; calculateMOBB(geometryFilter->GetOutput(), vol_mobb, surf_mobb); double pi = vnl_math::pi; double meshVolume = stats->GetVolume(); double meshSurf = stats->GetSurfaceArea(); AccessByItk_3(image, CalculateVolumeDensityStatistic, mask, meshVolume, featureList); //Calculate center of mass shift int xx = mask->GetDimensions()[0]; int yy = mask->GetDimensions()[1]; int zz = mask->GetDimensions()[2]; double xd = mask->GetGeometry()->GetSpacing()[0]; double yd = mask->GetGeometry()->GetSpacing()[1]; double zd = mask->GetGeometry()->GetSpacing()[2]; int minimumX=xx; int maximumX=0; int minimumY=yy; int maximumY=0; int minimumZ=zz; int maximumZ=0; vtkSmartPointer dataset1Arr = vtkSmartPointer::New(); vtkSmartPointer dataset2Arr = vtkSmartPointer::New(); vtkSmartPointer dataset3Arr = vtkSmartPointer::New(); dataset1Arr->SetNumberOfComponents(1); dataset2Arr->SetNumberOfComponents(1); dataset3Arr->SetNumberOfComponents(1); dataset1Arr->SetName("M1"); dataset2Arr->SetName("M2"); dataset3Arr->SetName("M3"); vtkSmartPointer dataset1ArrU = vtkSmartPointer::New(); vtkSmartPointer dataset2ArrU = vtkSmartPointer::New(); vtkSmartPointer dataset3ArrU = vtkSmartPointer::New(); dataset1ArrU->SetNumberOfComponents(1); dataset2ArrU->SetNumberOfComponents(1); dataset3ArrU->SetNumberOfComponents(1); dataset1ArrU->SetName("M1"); dataset2ArrU->SetName("M2"); dataset3ArrU->SetName("M3"); vtkSmartPointer points = vtkSmartPointer< vtkPoints >::New(); for (int x = 0; x < xx; x++) { for (int y = 0; y < yy; y++) { for (int z = 0; z < zz; z++) { itk::Image::IndexType index; index[0] = x; index[1] = y; index[2] = z; mitk::ScalarType pxImage; mitk::ScalarType pxMask; 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) { minimumX = std::min(x, minimumX); minimumY = std::min(y, minimumY); minimumZ = std::min(z, minimumZ); maximumX = std::max(x, maximumX); maximumY = std::max(y, maximumY); maximumZ = std::max(z, maximumZ); points->InsertNextPoint(x*xd, y*yd, z*zd); if (pxImage == pxImage) { dataset1Arr->InsertNextValue(x*xd); dataset2Arr->InsertNextValue(y*yd); dataset3Arr->InsertNextValue(z*zd); } } } } } vtkSmartPointer datasetTable = vtkSmartPointer::New(); datasetTable->AddColumn(dataset1Arr); datasetTable->AddColumn(dataset2Arr); datasetTable->AddColumn(dataset3Arr); vtkSmartPointer pcaStatistics = vtkSmartPointer::New(); pcaStatistics->SetInputData(vtkStatisticsAlgorithm::INPUT_DATA, datasetTable); pcaStatistics->SetColumnStatus("M1", 1); pcaStatistics->SetColumnStatus("M2", 1); pcaStatistics->SetColumnStatus("M3", 1); pcaStatistics->RequestSelectedColumns(); pcaStatistics->SetDeriveOption(true); pcaStatistics->Update(); vtkSmartPointer eigenvalues = vtkSmartPointer::New(); pcaStatistics->GetEigenvalues(eigenvalues); std::vector eigen_val(3); eigen_val[2] = eigenvalues->GetValue(0); eigen_val[1] = eigenvalues->GetValue(1); eigen_val[0] = eigenvalues->GetValue(2); double major = 2*sqrt(eigen_val[2]); double minor = 2*sqrt(eigen_val[1]); double least = 2*sqrt(eigen_val[0]); double alpha = std::sqrt(1 - minor*minor / major / major); double beta = std::sqrt(1 - least*least / major / major); double a = (maximumX - minimumX+1) * xd; double b = (maximumY - minimumY+1) * yd; double c = (maximumZ - minimumZ+1) * zd; double vd_aabb = meshVolume / (a*b*c); double ad_aabb = meshSurf / (2 * a*b + 2 * a*c + 2 * b*c); double vd_aee = 3 * meshVolume / (4.0*pi*major*minor*least); double ad_aee = 0; for (int i = 0; i < 20; ++i) { ad_aee += 4 * pi*major*minor*(alpha*alpha + beta*beta) / (2 * alpha*beta) * (std::pow(alpha*beta, i)) / (1 - 4 * i*i); } ad_aee = meshSurf / ad_aee; double vd_ch = meshVolume / stats2->GetVolume(); double ad_ch = meshSurf / stats2->GetSurfaceArea(); featureList.push_back(std::make_pair("Morphological Density::Volume density axis-aligned bounding box", vd_aabb)); featureList.push_back(std::make_pair("Morphological Density::Surface density axis-aligned bounding box", ad_aabb)); featureList.push_back(std::make_pair("Morphological Density::Volume density oriented minimum bounding box", meshVolume / vol_mobb)); featureList.push_back(std::make_pair("Morphological Density::Surface density oriented minimum bounding box", meshSurf / surf_mobb)); featureList.push_back(std::make_pair("Morphological Density::Volume density approx. enclosing ellipsoid", vd_aee)); featureList.push_back(std::make_pair("Morphological Density::Surface density approx. enclosing ellipsoid", ad_aee)); featureList.push_back(std::make_pair("Morphological Density::Volume density approx. minimum volume enclosing ellipsoid", meshVolume / vol_mvee)); featureList.push_back(std::make_pair("Morphological Density::Surface density approx. minimum volume enclosing ellipsoid", meshSurf / surf_mvee)); featureList.push_back(std::make_pair("Morphological Density::Volume density convex hull", vd_ch)); featureList.push_back(std::make_pair("Morphological Density::Surface density convex hull", ad_ch)); return featureList; } mitk::GIFVolumetricDensityStatistics::GIFVolumetricDensityStatistics() { SetLongName("volume-density"); SetShortName("volden"); } mitk::GIFVolumetricDensityStatistics::FeatureNameListType mitk::GIFVolumetricDensityStatistics::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFVolumetricDensityStatistics::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Volume-Density Statistic", "calculates volume density based features", us::Any()); } void mitk::GIFVolumetricDensityStatistics::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &, FeatureListType &featureList) { auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { MITK_INFO << "Start calculating volumetric density features ...."; auto localResults = this->CalculateFeatures(feature, mask); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating volumetric density features...."; } }