diff --git a/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp b/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp index c48da79eb9..9c3ff8700f 100644 --- a/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp +++ b/Modules/Classification/CLMiniApps/CLGlobalImageFeatures.cpp @@ -1,353 +1,353 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkCLPolyToNrrd_cpp #define mitkCLPolyToNrrd_cpp #include "time.h" #include #include #include #include "mitkCommandLineParser.h" -#include +#include #include #include #include #include #include #include #include #include #include "itkNearestNeighborInterpolateImageFunction.h" #include "itkResampleImageFilter.h" typedef itk::Image< double, 3 > FloatImageType; typedef itk::Image< unsigned char, 3 > MaskImageType; static std::vector splitDouble(std::string str, char delimiter) { std::vector internal; std::stringstream ss(str); // Turn the string into a stream. std::string tok; double val; while (std::getline(ss, tok, delimiter)) { std::stringstream s2(tok); s2 >> val; internal.push_back(val); } return internal; } template void ResampleImage(itk::Image* itkImage, float resolution, mitk::Image::Pointer& newImage) { typedef itk::Image ImageType; typedef itk::ResampleImageFilter ResampleFilterType; typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); auto spacing = itkImage->GetSpacing(); auto size = itkImage->GetLargestPossibleRegion().GetSize(); for (int i = 0; i < VImageDimension; ++i) { size[i] = size[i] / (1.0*resolution)*(1.0*spacing[i])+1.0; } spacing.Fill(resolution); resampler->SetInput(itkImage); resampler->SetSize(size); resampler->SetOutputSpacing(spacing); resampler->SetOutputOrigin(itkImage->GetOrigin()); resampler->SetOutputDirection(itkImage->GetDirection()); resampler->Update(); newImage->InitializeByItk(resampler->GetOutput()); mitk::GrabItkImageMemory(resampler->GetOutput(), newImage); } static void CreateNoNaNMask(mitk::Image::Pointer mask, mitk::Image::Pointer ref, mitk::Image::Pointer& newMask) { MaskImageType::Pointer itkMask = MaskImageType::New(); FloatImageType::Pointer itkValue = FloatImageType::New(); mitk::CastToItkImage(mask, itkMask); mitk::CastToItkImage(ref, itkValue); typedef itk::ImageDuplicator< MaskImageType > DuplicatorType; DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage(itkMask); duplicator->Update(); auto tmpMask = duplicator->GetOutput(); itk::ImageRegionIterator mask1Iter(itkMask, itkMask->GetLargestPossibleRegion()); itk::ImageRegionIterator mask2Iter(tmpMask, tmpMask->GetLargestPossibleRegion()); itk::ImageRegionIterator imageIter(itkValue, itkValue->GetLargestPossibleRegion()); while (!mask1Iter.IsAtEnd()) { mask2Iter.Set(0); if (mask1Iter.Value() > 0) { // Is not NaN if (imageIter.Value() == imageIter.Value()) { mask2Iter.Set(1); } } ++mask1Iter; ++mask2Iter; ++imageIter; } newMask->InitializeByItk(tmpMask); mitk::GrabItkImageMemory(tmpMask, newMask); } static void ResampleMask(mitk::Image::Pointer mask, mitk::Image::Pointer ref, mitk::Image::Pointer& newMask) { typedef itk::NearestNeighborInterpolateImageFunction< MaskImageType> NearestNeighborInterpolateImageFunctionType; typedef itk::ResampleImageFilter ResampleFilterType; NearestNeighborInterpolateImageFunctionType::Pointer nn_interpolator = NearestNeighborInterpolateImageFunctionType::New(); MaskImageType::Pointer itkMoving = MaskImageType::New(); MaskImageType::Pointer itkRef = MaskImageType::New(); mitk::CastToItkImage(mask, itkMoving); mitk::CastToItkImage(ref, itkRef); ResampleFilterType::Pointer resampler = ResampleFilterType::New(); resampler->SetInput(itkMoving); resampler->SetReferenceImage(itkRef); resampler->UseReferenceImageOn(); resampler->SetInterpolator(nn_interpolator); resampler->Update(); newMask->InitializeByItk(resampler->GetOutput()); mitk::GrabItkImageMemory(resampler->GetOutput(), newMask); } int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); // required params parser.addArgument("image", "i", mitkCommandLineParser::InputImage, "Input Image", "Path to the input VTK polydata", us::Any(), false); parser.addArgument("mask", "m", mitkCommandLineParser::InputImage, "Input Mask", "Mask Image that specifies the area over for the statistic, (Values = 1)", us::Any(), false); parser.addArgument("output", "o", mitkCommandLineParser::OutputFile, "Output text file", "Target file. The output statistic is appended to this file.", us::Any(), false); parser.addArgument("cooccurence","cooc",mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features",us::Any()); parser.addArgument("volume","vol",mitkCommandLineParser::String, "Use Volume-Statistic", "calculates volume based features",us::Any()); parser.addArgument("run-length","rl",mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features",us::Any()); parser.addArgument("first-order","fo",mitkCommandLineParser::String, "Use First Order Features", "calculates First order based features",us::Any()); parser.addArgument("header","head",mitkCommandLineParser::String,"Add Header (Labels) to output","",us::Any()); parser.addArgument("description","d",mitkCommandLineParser::String,"Text","Description that is added to the output",us::Any()); parser.addArgument("same-space", "sp", mitkCommandLineParser::String, "Bool", "Set the spacing of all images to equal. Otherwise an error will be thrown. ", us::Any()); parser.addArgument("resample-mask", "rm", mitkCommandLineParser::Bool, "Bool", "Resamples the mask to the resolution of the input image ", us::Any()); parser.addArgument("save-resample-mask", "srm", mitkCommandLineParser::String, "String", "If specified the resampled mask is saved to this path (if -rm is 1)", us::Any()); parser.addArgument("fixed-isotropic", "fi", mitkCommandLineParser::Float, "Float", "Input image resampled to fixed isotropic resolution given in mm. Should be used with resample-mask ", us::Any()); parser.addArgument("direction", "dir", mitkCommandLineParser::String, "Int", "Allows to specify the direction for Cooc and RL. 0: All directions, 1: Only single direction (Test purpose), 2,3,4... Without dimension 0,1,2... ", us::Any()); // Miniapp Infos parser.setCategory("Classification Tools"); parser.setTitle("Global Image Feature calculator"); parser.setDescription("Calculates different global statistics for a given segmentation / image combination"); parser.setContributor("MBI"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) { return EXIT_FAILURE; } if ( parsedArgs.count("help") || parsedArgs.count("h")) { return EXIT_SUCCESS; } MITK_INFO << "Version: "<< 1.7; //bool useCooc = parsedArgs.count("cooccurence"); bool resampleMask = false; if (parsedArgs.count("resample-mask")) { resampleMask = us::any_cast(parsedArgs["resample-mask"]); } mitk::Image::Pointer image = mitk::IOUtil::LoadImage(parsedArgs["image"].ToString()); mitk::Image::Pointer mask = mitk::IOUtil::LoadImage(parsedArgs["mask"].ToString()); if (parsedArgs.count("fixed-isotropic")) { mitk::Image::Pointer newImage = mitk::Image::New(); float resolution = us::any_cast(parsedArgs["fixed-isotropic"]); AccessByItk_2(image, ResampleImage, resolution, newImage); image = newImage; } if (resampleMask) { mitk::Image::Pointer newMaskImage = mitk::Image::New(); ResampleMask(mask, image, newMaskImage); mask = newMaskImage; if (parsedArgs.count("save-resample-mask")) { mitk::IOUtil::SaveImage(mask, parsedArgs["save-resample-mask"].ToString()); } } bool fixDifferentSpaces = parsedArgs.count("same-space"); if ( ! mitk::Equal(mask->GetGeometry(0)->GetOrigin(), image->GetGeometry(0)->GetOrigin())) { MITK_INFO << "Not equal Origins"; if (fixDifferentSpaces) { image->GetGeometry(0)->SetOrigin(mask->GetGeometry(0)->GetOrigin()); } else { return -1; } } if ( ! mitk::Equal(mask->GetGeometry(0)->GetSpacing(), image->GetGeometry(0)->GetSpacing())) { MITK_INFO << "Not equal Sapcings"; if (fixDifferentSpaces) { image->GetGeometry(0)->SetSpacing(mask->GetGeometry(0)->GetSpacing()); } else { return -1; } } int direction = 0; if (parsedArgs.count("direction")) { direction = splitDouble(parsedArgs["direction"].ToString(), ';')[0]; } mitk::Image::Pointer maskNoNaN = mitk::Image::New(); CreateNoNaNMask(mask, image, maskNoNaN); mitk::AbstractGlobalImageFeature::FeatureListType stats; //////////////////////////////////////////////////////////////// // Calculate First Order Features //////////////////////////////////////////////////////////////// if (parsedArgs.count("first-order")) { MITK_INFO << "Start calculating first order statistics...."; mitk::GIFFirstOrderStatistics::Pointer firstOrderCalculator = mitk::GIFFirstOrderStatistics::New(); auto localResults = firstOrderCalculator->CalculateFeatures(image, maskNoNaN); stats.insert(stats.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating first order statistics...."; } //////////////////////////////////////////////////////////////// // Calculate Volume based Features //////////////////////////////////////////////////////////////// if (parsedArgs.count("volume")) { MITK_INFO << "Start calculating volumetric ...."; mitk::GIFVolumetricStatistics::Pointer volCalculator = mitk::GIFVolumetricStatistics::New(); auto localResults = volCalculator->CalculateFeatures(image, mask); stats.insert(stats.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating volumetric...."; } //////////////////////////////////////////////////////////////// // Calculate Co-occurence Features //////////////////////////////////////////////////////////////// if (parsedArgs.count("cooccurence")) { auto ranges = splitDouble(parsedArgs["cooccurence"].ToString(),';'); for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; - mitk::GIFCooccurenceMatrix::Pointer coocCalculator = mitk::GIFCooccurenceMatrix::New(); + mitk::GIFCooccurenceMatrix2::Pointer coocCalculator = mitk::GIFCooccurenceMatrix2::New(); coocCalculator->SetRange(ranges[i]); coocCalculator->SetDirection(direction); - auto localResults = coocCalculator->CalculateFeatures(image, mask); + auto localResults = coocCalculator->CalculateFeatures(image, maskNoNaN); stats.insert(stats.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; } } //////////////////////////////////////////////////////////////// // Calculate Run-Length Features //////////////////////////////////////////////////////////////// if (parsedArgs.count("run-length")) { auto ranges = splitDouble(parsedArgs["run-length"].ToString(),';'); for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating run-length with number of bins " << ranges[i] << "...."; mitk::GIFGrayLevelRunLength::Pointer calculator = mitk::GIFGrayLevelRunLength::New(); calculator->SetRange(ranges[i]); auto localResults = calculator->CalculateFeatures(image, mask); stats.insert(stats.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating run-length with number of bins " << ranges[i] << "...."; } } for (std::size_t i = 0; i < stats.size(); ++i) { std::cout << stats[i].first << " - " << stats[i].second < \ const \ typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * \ EnhancedHistogramToTextureFeaturesFilter< THistogram > \ ::Get##name##Output() const \ { \ return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(id) ); \ } \ \ template< typename THistogram > \ typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType \ EnhancedHistogramToTextureFeaturesFilter< THistogram > \ ::Get##name() const \ { \ return this->Get##name##Output()->Get(); \ } namespace itk { namespace Statistics { //constructor template< typename THistogram > EnhancedHistogramToTextureFeaturesFilter< THistogram >::EnhancedHistogramToTextureFeaturesFilter(void) { this->ProcessObject::SetNumberOfRequiredInputs(1); // allocate the data objects for the outputs which are // just decorators real types for ( int i = 0; i < 28; ++i ) { this->ProcessObject::SetNthOutput( i, this->MakeOutput(i) ); } } template< typename THistogram > void EnhancedHistogramToTextureFeaturesFilter< THistogram > ::SetInput(const HistogramType *histogram) { this->ProcessObject::SetNthInput( 0, const_cast< HistogramType * >( histogram ) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::HistogramType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetInput() const { return itkDynamicCastInDebugMode< const HistogramType * >( this->GetPrimaryInput() ); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::DataObjectPointer EnhancedHistogramToTextureFeaturesFilter< THistogram > ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed(idx) ) { return MeasurementObjectType::New().GetPointer(); } template< typename THistogram > void EnhancedHistogramToTextureFeaturesFilter< THistogram >::GenerateData(void) { typedef typename HistogramType::ConstIterator HistogramIterator; const HistogramType *inputHistogram = this->GetInput(); //Normalize the absolute frequencies and populate the relative frequency //container TotalRelativeFrequencyType totalFrequency = static_cast< TotalRelativeFrequencyType >( inputHistogram->GetTotalFrequency() ); m_RelativeFrequencyContainer.clear(); typename HistogramType::SizeValueType binsPerAxis = inputHistogram->GetSize(0); std::vector sumP; std::vector diffP; sumP.resize(2*binsPerAxis,0.0); diffP.resize(binsPerAxis,0.0); double numberOfPixels = 0; for ( HistogramIterator hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { AbsoluteFrequencyType frequency = hit.GetFrequency(); RelativeFrequencyType relativeFrequency = (totalFrequency > 0) ? frequency / totalFrequency : 0; m_RelativeFrequencyContainer.push_back(relativeFrequency); IndexType index = inputHistogram->GetIndex( hit.GetInstanceIdentifier() ); sumP[index[0] + index[1]] += relativeFrequency; diffP[std::abs(index[0] - index[1])] += relativeFrequency; //if (index[1] == 0) numberOfPixels += frequency; } // Now get the various means and variances. This is takes two passes // through the histogram. double pixelMean; double marginalMean; double marginalDevSquared; double pixelVariance; this->ComputeMeansAndVariances(pixelMean, marginalMean, marginalDevSquared, pixelVariance); // Finally compute the texture features. Another one pass. MeasurementType energy = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType entropy = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType correlation = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inverseDifferenceMoment = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inertia = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType clusterShade = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType clusterProminence = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType haralickCorrelation = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType autocorrelation = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType contrast = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType dissimilarity = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType maximumProbability = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inverseVariance = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType homogeneity1 = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType clusterTendency = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType variance = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType sumAverage = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType sumEntropy = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType sumVariance = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType diffAverage = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType diffEntropy = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType diffVariance = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inverseDifferenceMomentNormalized = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inverseDifferenceNormalized = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inverseDifference = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType averageProbability = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType firstMeasureOfInformationCorrelation = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType secondMeasureOfInformationCorrelation = NumericTraits< MeasurementType >::ZeroValue(); double pixelVarianceSquared = pixelVariance * pixelVariance; // Variance is only used in correlation. If variance is 0, then // (index[0] - pixelMean) * (index[1] - pixelMean) // should be zero as well. In this case, set the variance to 1. in // order to avoid NaN correlation. if( Math::FloatAlmostEqual( pixelVarianceSquared, 0.0, 4, 2*NumericTraits::epsilon() ) ) { pixelVarianceSquared = 1.; } const double log2 = std::log(2.0); typename RelativeFrequencyContainerType::const_iterator rFreqIterator = m_RelativeFrequencyContainer.begin(); IndexType globalIndex(2); for ( HistogramIterator hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { RelativeFrequencyType frequency = *rFreqIterator; ++rFreqIterator; if ( frequency == 0 ) { continue; // no use doing these calculations if we're just multiplying by // zero. } IndexType index = inputHistogram->GetIndex( hit.GetInstanceIdentifier() ); globalIndex = index; energy += frequency * frequency; entropy -= ( frequency > 0.0001 ) ? frequency *std::log(frequency) / log2:0; correlation += ( ( index[0] - pixelMean ) * ( index[1] - pixelMean ) * frequency ) / pixelVarianceSquared; inverseDifferenceMoment += frequency / ( 1.0 + ( index[0] - index[1] ) * ( index[0] - index[1] ) ); inertia += ( index[0] - index[1] ) * ( index[0] - index[1] ) * frequency; clusterShade += std::pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 3 ) * frequency; clusterProminence += std::pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 4 ) * frequency; haralickCorrelation += index[0] * index[1] * frequency; // New Features added for Aerts compatibility autocorrelation +=index[0] * index[1] * frequency; contrast += (index[0] - index[1]) * (index[0] - index[1]) * frequency; dissimilarity += std::abs(index[0] - index[1]) * frequency; - maximumProbability = std::max(maximumProbability, frequency); + maximumProbability +=std::max(maximumProbability, frequency); averageProbability += frequency * index[0]; if (index[0] != index[1]) inverseVariance += frequency / ((index[0] - index[1])*(index[0] - index[1])); homogeneity1 +=frequency / ( 1.0 + std::abs( index[0] - index[1] )); clusterTendency += std::pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 2 ) * frequency; variance += std::pow( ( index[0] - pixelMean ), 2) * frequency; if (numberOfPixels > 0) { inverseDifferenceMomentNormalized += frequency / ( 1.0 + ( index[0] - index[1] ) * ( index[0] - index[1] ) / numberOfPixels / numberOfPixels); inverseDifferenceNormalized += frequency / ( 1.0 + std::abs( index[0] - index[1] ) / numberOfPixels ); } else { inverseDifferenceMomentNormalized = 0; inverseDifferenceNormalized = 0; } inverseDifference += frequency / ( 1.0 + std::abs( index[0] - index[1] ) ); } for (int i = 0; i < (int)sumP.size();++i) { double frequency = sumP[i]; sumAverage += i * frequency; sumEntropy -= ( frequency > 0.0001 ) ? frequency *std::log(frequency) / log2:0; } for (int i = 0; i < (int)sumP.size();++i) { double frequency = sumP[i]; sumVariance += (i-sumAverage)*(i-sumAverage) * frequency; } for (int i = 0; i < (int)diffP.size();++i) { double frequency = diffP[i]; diffAverage += i * frequency; diffEntropy -= ( frequency > 0.0001 ) ? frequency *std::log(frequency) / log2:0; } for (int i = 0; i < (int)diffP.size();++i) { double frequency = diffP[i]; sumVariance += (i-diffAverage)*(i-diffAverage) * frequency; } if (marginalDevSquared > 0) { haralickCorrelation = ( haralickCorrelation - marginalMean * marginalMean ) / marginalDevSquared; } else { haralickCorrelation =0; } //Calculate the margin probs std::vector pi_margins; std::vector pj_margins; //pi. for ( int i = 1; i <= inputHistogram->GetSize(0); i++ ) { double pi_tmp= 0.0; for( int j = 1; j <= inputHistogram->GetSize(1); j++ ) { globalIndex[0] = i; globalIndex[1] = j; pi_tmp += inputHistogram->GetFrequency(globalIndex); } pi_margins.push_back(pi_tmp); } //pj. for ( int j = 1; j <= inputHistogram->GetSize(1); j++ ) { double pj_tmp= 0.0; for( int i = 1; i <= inputHistogram->GetSize(0); i++ ) { globalIndex[0] = i; globalIndex[1] = j; pj_tmp += inputHistogram->GetFrequency(globalIndex); } pj_margins.push_back(pj_tmp); } //Calculate HX double hx = 0.0; for ( int i = 0; i < inputHistogram->GetSize(0); i++ ) { double pi_margin = pi_margins[i]; hx -= ( pi_margin > 0.0001 ) ? pi_margin *std::log(pi_margin) / log2:0; } //Calculate HXY1 double hxy1 = 0.0; for ( int i = 0; i < inputHistogram->GetSize(0); i++ ) { for ( int j = 0; j < inputHistogram->GetSize(1); j++ ) { globalIndex[0] = i; globalIndex[1] = j; double pi_margin = pi_margins[i]; double pj_margin = pj_margins[j]; double p_ij = inputHistogram->GetFrequency(globalIndex); hxy1 -= ( (pi_margin * pj_margin) > 0.0001 ) ? p_ij *std::log(pi_margin * pj_margin) / log2:0; } } //Calculate HXY2 double hxy2 = 0.0; for ( int i = 0; i < inputHistogram->GetSize(0); i++ ) { for ( int j = 0; j < inputHistogram->GetSize(1); j++ ) { globalIndex[0] = i; globalIndex[1] = j; double pi_margin = pi_margins[i]; double pj_margin = pj_margins[j]; hxy1 -= ( (pi_margin * pj_margin) > 0.0001 ) ? (pi_margin * pj_margin) *std::log(pi_margin * pj_margin) / log2:0; } } firstMeasureOfInformationCorrelation = (entropy - hxy1) / hx; secondMeasureOfInformationCorrelation = (entropy > hxy2) ? 0 : std::sqrt(1 - std::exp(-2*(hxy2 - entropy))); MeasurementObjectType *energyOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(0) ); energyOutputObject->Set(energy); MeasurementObjectType *entropyOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(1) ); entropyOutputObject->Set(entropy); MeasurementObjectType *correlationOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(2) ); correlationOutputObject->Set(correlation); MeasurementObjectType *inverseDifferenceMomentOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(3) ); inverseDifferenceMomentOutputObject->Set(inverseDifferenceMoment); MeasurementObjectType *inertiaOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(4) ); inertiaOutputObject->Set(inertia); MeasurementObjectType *clusterShadeOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(5) ); clusterShadeOutputObject->Set(clusterShade); MeasurementObjectType *clusterProminenceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(6) ); clusterProminenceOutputObject->Set(clusterProminence); MeasurementObjectType *haralickCorrelationOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(7) ); haralickCorrelationOutputObject->Set(haralickCorrelation); MeasurementObjectType *autocorrelationOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(8) ); autocorrelationOutputObject->Set(autocorrelation); MeasurementObjectType *contrastOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(9) ); contrastOutputObject->Set(contrast); MeasurementObjectType *dissimilarityOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(10) ); dissimilarityOutputObject->Set(dissimilarity); MeasurementObjectType *maximumProbabilityOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(11) ); maximumProbabilityOutputObject->Set(maximumProbability); MeasurementObjectType *inverseVarianceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(12) ); inverseVarianceOutputObject->Set(inverseVariance); MeasurementObjectType *homogeneity1OutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(13) ); homogeneity1OutputObject->Set(homogeneity1); MeasurementObjectType *clusterTendencyOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(14) ); clusterTendencyOutputObject->Set(clusterTendency); MeasurementObjectType *varianceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(15) ); varianceOutputObject->Set(variance); MeasurementObjectType *sumAverageOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(16) ); sumAverageOutputObject->Set(sumAverage); MeasurementObjectType *sumEntropyOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(17) ); sumEntropyOutputObject->Set(sumEntropy); MeasurementObjectType *sumVarianceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(18) ); sumVarianceOutputObject->Set(sumVariance); MeasurementObjectType *diffAverageOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(19) ); diffAverageOutputObject->Set(diffAverage); MeasurementObjectType *diffEntropyOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(20) ); diffEntropyOutputObject->Set(diffEntropy); MeasurementObjectType *diffVarianceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(21) ); diffVarianceOutputObject->Set(diffVariance); MeasurementObjectType *inverseDifferenceMomentNormalizedOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(22) ); inverseDifferenceMomentNormalizedOutputObject->Set(inverseDifferenceMomentNormalized); MeasurementObjectType *inverseDifferenceNormalizedOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(23) ); inverseDifferenceNormalizedOutputObject->Set(inverseDifferenceNormalized); MeasurementObjectType *inverseDifferenceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(24) ); inverseDifferenceOutputObject->Set(inverseDifference); MeasurementObjectType *jointAverageOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(25) ); jointAverageOutputObject->Set(averageProbability); MeasurementObjectType *firstMeasureOfInformationCorrelationOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(26) ); firstMeasureOfInformationCorrelationOutputObject->Set(firstMeasureOfInformationCorrelation); MeasurementObjectType *secondMeasureOfInformationCorrelationOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(27) ); secondMeasureOfInformationCorrelationOutputObject->Set(secondMeasureOfInformationCorrelation); } template< typename THistogram > void EnhancedHistogramToTextureFeaturesFilter< THistogram >::ComputeMeansAndVariances(double & pixelMean, double & marginalMean, double & marginalDevSquared, double & pixelVariance) { // This function takes two passes through the histogram and two passes through // an array of the same length as a histogram axis. This could probably be // cleverly compressed to one pass, but it's not clear that that's necessary. typedef typename HistogramType::ConstIterator HistogramIterator; const HistogramType *inputHistogram = this->GetInput(); // Initialize everything typename HistogramType::SizeValueType binsPerAxis = inputHistogram->GetSize(0); double *marginalSums = new double[binsPerAxis]; for ( double *ms_It = marginalSums; ms_It < marginalSums + binsPerAxis; ms_It++ ) { *ms_It = 0; } pixelMean = 0; typename RelativeFrequencyContainerType::const_iterator rFreqIterator = m_RelativeFrequencyContainer.begin(); // Ok, now do the first pass through the histogram to get the marginal sums // and compute the pixel mean HistogramIterator hit = inputHistogram->Begin(); while ( hit != inputHistogram->End() ) { RelativeFrequencyType frequency = *rFreqIterator; IndexType index = inputHistogram->GetIndex( hit.GetInstanceIdentifier() ); pixelMean += index[0] * frequency; marginalSums[index[0]] += frequency; ++hit; ++rFreqIterator; } /* Now get the mean and deviaton of the marginal sums. Compute incremental mean and SD, a la Knuth, "The Art of Computer Programming, Volume 2: Seminumerical Algorithms", section 4.2.2. Compute mean and standard deviation using the recurrence relation: M(1) = x(1), M(k) = M(k-1) + (x(k) - M(k-1) ) / k S(1) = 0, S(k) = S(k-1) + (x(k) - M(k-1)) * (x(k) - M(k)) for 2 <= k <= n, then sigma = std::sqrt(S(n) / n) (or divide by n-1 for sample SD instead of population SD). */ marginalMean = marginalSums[0]; marginalDevSquared = 0; for ( unsigned int arrayIndex = 1; arrayIndex < binsPerAxis; arrayIndex++ ) { int k = arrayIndex + 1; double M_k_minus_1 = marginalMean; double S_k_minus_1 = marginalDevSquared; double x_k = marginalSums[arrayIndex]; double M_k = M_k_minus_1 + ( x_k - M_k_minus_1 ) / k; double S_k = S_k_minus_1 + ( x_k - M_k_minus_1 ) * ( x_k - M_k ); marginalMean = M_k; marginalDevSquared = S_k; } marginalDevSquared = marginalDevSquared / binsPerAxis; rFreqIterator = m_RelativeFrequencyContainer.begin(); // OK, now compute the pixel variances. pixelVariance = 0; for ( hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { RelativeFrequencyType frequency = *rFreqIterator; IndexType index = inputHistogram->GetIndex( hit.GetInstanceIdentifier() ); pixelVariance += ( index[0] - pixelMean ) * ( index[0] - pixelMean ) * frequency; ++rFreqIterator; } delete[] marginalSums; } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetEnergyOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(0) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetEntropyOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(1) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetCorrelationOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(2) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetInverseDifferenceMomentOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(3) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetInertiaOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(4) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetClusterShadeOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(5) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetClusterProminenceOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(6) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetHaralickCorrelationOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(7) ); } itkMacroGLCMFeature(Autocorrelation,8) itkMacroGLCMFeature(Contrast,9) itkMacroGLCMFeature(Dissimilarity,10) itkMacroGLCMFeature(MaximumProbability,11) itkMacroGLCMFeature(InverseVariance,12) itkMacroGLCMFeature(Homogeneity1,13) itkMacroGLCMFeature(ClusterTendency,14) itkMacroGLCMFeature(Variance,15) itkMacroGLCMFeature(SumAverage,16) itkMacroGLCMFeature(SumEntropy,17) itkMacroGLCMFeature(SumVariance,18) itkMacroGLCMFeature(DifferenceAverage,19) itkMacroGLCMFeature(DifferenceEntropy,20) itkMacroGLCMFeature(DifferenceVariance,21) itkMacroGLCMFeature(InverseDifferenceMomentNormalized,22) itkMacroGLCMFeature(InverseDifferenceNormalized,23) itkMacroGLCMFeature(InverseDifference,24) itkMacroGLCMFeature(JointAverage,25) itkMacroGLCMFeature(FirstMeasureOfInformationCorrelation,26) itkMacroGLCMFeature(SecondMeasureOfInformationCorrelation,27) template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetEnergy() const { return this->GetEnergyOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetEntropy() const { return this->GetEntropyOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetCorrelation() const { return this->GetCorrelationOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetInverseDifferenceMoment() const { return this->GetInverseDifferenceMomentOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetInertia() const { return this->GetInertiaOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetClusterShade() const { return this->GetClusterShadeOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetClusterProminence() const { return this->GetClusterProminenceOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetHaralickCorrelation() const { return this->GetHaralickCorrelationOutput()->Get(); } #define itkMacroGLCMFeatureSwitch(name) \ case name : \ return this->Get##name() template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetFeature(TextureFeatureName feature) { switch ( feature ) { itkMacroGLCMFeatureSwitch(Energy); itkMacroGLCMFeatureSwitch(Entropy); itkMacroGLCMFeatureSwitch(Correlation); itkMacroGLCMFeatureSwitch(InverseDifferenceMoment); itkMacroGLCMFeatureSwitch(Inertia); itkMacroGLCMFeatureSwitch(ClusterShade); itkMacroGLCMFeatureSwitch(ClusterProminence); itkMacroGLCMFeatureSwitch(HaralickCorrelation); itkMacroGLCMFeatureSwitch(Autocorrelation); itkMacroGLCMFeatureSwitch(Contrast); itkMacroGLCMFeatureSwitch(Dissimilarity); itkMacroGLCMFeatureSwitch(MaximumProbability); itkMacroGLCMFeatureSwitch(InverseVariance); itkMacroGLCMFeatureSwitch(Homogeneity1); itkMacroGLCMFeatureSwitch(ClusterTendency); itkMacroGLCMFeatureSwitch(Variance); itkMacroGLCMFeatureSwitch(SumAverage); itkMacroGLCMFeatureSwitch(SumEntropy); itkMacroGLCMFeatureSwitch(SumVariance); itkMacroGLCMFeatureSwitch(DifferenceAverage); itkMacroGLCMFeatureSwitch(DifferenceEntropy); itkMacroGLCMFeatureSwitch(DifferenceVariance); itkMacroGLCMFeatureSwitch(InverseDifferenceMomentNormalized); itkMacroGLCMFeatureSwitch(InverseDifferenceNormalized); itkMacroGLCMFeatureSwitch(InverseDifference); itkMacroGLCMFeatureSwitch(JointAverage); itkMacroGLCMFeatureSwitch(FirstMeasureOfInformationCorrelation); itkMacroGLCMFeatureSwitch(SecondMeasureOfInformationCorrelation); default: return 0; } } #undef itkMacroGLCMFeatureSwitch template< typename THistogram > void EnhancedHistogramToTextureFeaturesFilter< THistogram > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); } } // end of namespace Statistics } // end of namespace itk #endif diff --git a/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h b/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h new file mode 100644 index 0000000000..fc3dad39f1 --- /dev/null +++ b/Modules/Classification/CLUtilities/include/mitkGIFCooccurenceMatrix2.h @@ -0,0 +1,74 @@ +#ifndef mitkGIFCooccurenceMatrix2_h +#define mitkGIFCooccurenceMatrix2_h + +#include +#include +#include + +#include + +namespace mitk +{ + struct CoocurenceMatrixHolder + { + public: + CoocurenceMatrixHolder(double min, double max, int number); + + int IntensityToIndex(double intensity); + double IndexToMinIntensity(int index); + double IndexToMeanIntensity(int index); + double IndexToMaxIntensity(int index); + + double m_MinimumRange; + double m_MaximumRange; + double m_Stepsize; + int m_NumberOfBins; + Eigen::MatrixXd m_Matrix; + + }; + + struct CoocurenceMatrixFeatures + { + public: + double JointMaximum; + double JointAverage; + }; + + + class MITKCLUTILITIES_EXPORT GIFCooccurenceMatrix2 : public AbstractGlobalImageFeature + { + public: + mitkClassMacro(GIFCooccurenceMatrix2,AbstractGlobalImageFeature) + itkFactorylessNewMacro(Self) + itkCloneMacro(Self) + + GIFCooccurenceMatrix2(); + + /** + * \brief Calculates the Cooccurence-Matrix based features for this class. + */ + virtual FeatureListType CalculateFeatures(const Image::Pointer & image, const Image::Pointer &feature) override; + + /** + * \brief Returns a list of the names of all features that are calculated from this class + */ + virtual FeatureNameListType GetFeatureNames() override; + + itkGetConstMacro(Range,double); + itkSetMacro(Range, double); + itkGetConstMacro(Direction, unsigned int); + itkSetMacro(Direction, unsigned int); + + struct GIFCooccurenceMatrix2Configuration + { + double range; + unsigned int direction; + }; + + private: + double m_Range; + unsigned int m_Direction; + }; + +} +#endif //mitkGIFCooccurenceMatrix2_h diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp new file mode 100644 index 0000000000..7a1d6f4b36 --- /dev/null +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp @@ -0,0 +1,452 @@ +#include + +// MITK +#include +#include +#include + +// ITK +#include +#include +#include +#include + +// STL +#include + +mitk::CoocurenceMatrixHolder::CoocurenceMatrixHolder(double min, double max, int number) : +m_MinimumRange(min), +m_MaximumRange(max), +m_NumberOfBins(number) +{ + m_Matrix.resize(number, number); + m_Matrix.fill(0); + m_Stepsize = (max - min) / (number); +} + +int mitk::CoocurenceMatrixHolder::IntensityToIndex(double intensity) +{ + return std::floor((intensity - m_MinimumRange) / m_Stepsize); +} + +double mitk::CoocurenceMatrixHolder::IndexToMinIntensity(int index) +{ + return m_MinimumRange + index * m_Stepsize; +} +double mitk::CoocurenceMatrixHolder::IndexToMeanIntensity(int index) +{ + return m_MinimumRange + (index+0.5) * m_Stepsize; +} +double mitk::CoocurenceMatrixHolder::IndexToMaxIntensity(int index) +{ + return m_MinimumRange + (index + 1) * m_Stepsize; +} + +template +void +CalculateCoOcMatrix(itk::Image* itkImage, + itk::Image* mask, + itk::Offset offset, + int range, + mitk::CoocurenceMatrixHolder &holder) +{ + typedef itk::Image ImageType; + typedef itk::Image MaskImageType; + typedef itk::ShapedNeighborhoodIterator ShapeIterType; + typedef itk::ShapedNeighborhoodIterator ShapeMaskIterType; + typedef itk::ImageRegionConstIterator ConstIterType; + typedef itk::ImageRegionConstIterator ConstMaskIterType; + + + itk::Size radius; + radius.Fill(range+1); + ShapeIterType imageOffsetIter(radius, itkImage, itkImage->GetLargestPossibleRegion()); + ShapeMaskIterType maskOffsetIter(radius, mask, mask->GetLargestPossibleRegion()); + imageOffsetIter.ActivateOffset(offset); + maskOffsetIter.ActivateOffset(offset); + ConstIterType imageIter(itkImage, itkImage->GetLargestPossibleRegion()); + ConstMaskIterType maskIter(mask, mask->GetLargestPossibleRegion()); + // iterator.GetIndex() + ci.GetNeighborhoodOffset() + auto region = mask->GetLargestPossibleRegion(); + + + while (!maskIter.IsAtEnd()) + { + auto ciMask = maskOffsetIter.Begin(); + auto ciValue = imageOffsetIter.Begin(); + if (maskIter.Value() > 0 && + ciMask.Get() > 0 && + imageIter.Get() == imageIter.Get() && + ciValue.Get() == ciValue.Get() && + region.IsInside(maskOffsetIter.GetIndex() + ciMask.GetNeighborhoodOffset())) + { + int i = holder.IntensityToIndex(imageIter.Get()); + int j = holder.IntensityToIndex(ciValue.Get()); + holder.m_Matrix(i, j) += 1; + holder.m_Matrix(j, i) += 1; + } + ++imageOffsetIter; + ++maskOffsetIter; + ++imageIter; + ++maskIter; + } +} + +void CalculateFeatures( + mitk::CoocurenceMatrixHolder &holder, + mitk::CoocurenceMatrixFeatures & results + ) +{ + auto pijMatrix = holder.m_Matrix; + auto piMatrix = holder.m_Matrix; + auto pjMatrix = holder.m_Matrix; + 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; + } + + results.JointMaximum += pijMatrix.maxCoeff(); + + + //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); + //} + + +} + +template +void +CalculateCoocurenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFCooccurenceMatrix2::FeatureListType & featureList, mitk::GIFCooccurenceMatrix2::GIFCooccurenceMatrix2Configuration config) +{ + typedef itk::Image ImageType; + typedef itk::Image MaskType; + typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; + typedef itk::Neighborhood NeighborhoodType; + typedef itk::Offset OffsetType; + + typedef itk::Statistics::EnhancedScalarImageToTextureFeaturesFilter FilterType; + + typedef typename FilterType::TextureFeaturesFilterType TextureFilterType; + + + + /////////////////////////////////////////////////////////////////////////////////////////////// + + typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); + minMaxComputer->SetImage(itkImage); + minMaxComputer->Compute(); + + double rangeMin = -0.5 + minMaxComputer->GetMinimum(); + double rangeMax = 0.5 + minMaxComputer->GetMaximum(); + + // Define Range + int numberOfBins = 6; + + typename MaskType::Pointer maskImage = MaskType::New(); + mitk::CastToItkImage(mask, maskImage); + + + //Find possible directions + std::vector < itk::Offset > offsetVector; + NeighborhoodType hood; + hood.SetRadius(1); + unsigned int centerIndex = hood.GetCenterNeighborhoodIndex(); + OffsetType offset; + for (unsigned int d = 0; d < centerIndex; d++) + { + offset = hood.GetOffset(d); + bool useOffset = true; + for (int i = 0; i < VImageDimension; ++i) + { + offset[i] *= config.range; + if (config.direction == i + 2 && offset[i] != 0) + { + useOffset = false; + } + } + if (useOffset) + { + offsetVector.push_back(offset); + } + } + + mitk::CoocurenceMatrixFeatures coocResults; + coocResults.JointAverage = 0; + coocResults.JointMaximum = 0; + + mitk::CoocurenceMatrixHolder holder(rangeMin, rangeMax, numberOfBins); + for (int i = 0; i < offsetVector.size(); ++i) + { + offset = offsetVector[i]; + MITK_INFO << offset; + CalculateCoOcMatrix(itkImage, maskImage, offset, config.range, holder); + } + CalculateFeatures(holder, coocResults); + //coocResults.JointMaximum /= offsetVector.size(); + MITK_INFO << coocResults.JointMaximum; + + /////////////////////////////////////////////////////////////////////////////////////////////// + /* + + typename FilterType::Pointer filter = FilterType::New(); + + typename FilterType::OffsetVector::Pointer newOffset = FilterType::OffsetVector::New(); + auto oldOffsets = filter->GetOffsets(); + auto oldOffsetsIterator = oldOffsets->Begin(); + while(oldOffsetsIterator != oldOffsets->End()) + { + bool continueOuterLoop = false; + typename FilterType::OffsetType offset = oldOffsetsIterator->Value(); + for (unsigned int i = 0; i < VImageDimension; ++i) + { + offset[i] *= config.range; + if (config.direction == i + 2 && offset[i] != 0) + { + continueOuterLoop = true; + } + } + if (config.direction == 1) + { + offset[0] = 0; + offset[1] = 0; + offset[2] = 1; + newOffset->push_back(offset); + break; + } + + + oldOffsetsIterator++; + if (continueOuterLoop) + continue; + newOffset->push_back(offset); + + } + filter->SetOffsets(newOffset); + + + filter->Update(); + + auto featureMeans = filter->GetFeatureMeans (); + auto featureStd = filter->GetFeatureStandardDeviations(); + + std::ostringstream ss; + ss << config.range; + std::string strRange = ss.str(); + for (std::size_t i = 0; i < featureMeans->size(); ++i) + { + switch (i) + { + case TextureFilterType::Energy : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Energy Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Energy Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::Entropy : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Entropy Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Entropy Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::Correlation : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Correlation Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Correlation Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::InverseDifferenceMoment : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMoment Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMoment Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::Inertia : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Inertia Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Inertia Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::ClusterShade : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterShade Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterShade Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::ClusterProminence : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterProminence Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterProminence Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::HaralickCorrelation : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") HaralickCorrelation Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") HaralickCorrelation Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::Autocorrelation : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Autocorrelation Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Autocorrelation Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::Contrast : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Contrast Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Contrast Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::Dissimilarity : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Dissimilarity Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Dissimilarity Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::MaximumProbability : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") MaximumProbability Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") MaximumProbability Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::InverseVariance : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseVariance Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseVariance Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::Homogeneity1 : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Homogeneity1 Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Homogeneity1 Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::ClusterTendency : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterTendency Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterTendency Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::Variance : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Variance Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Variance Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::SumAverage : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumAverage Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumAverage Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::SumEntropy : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumEntropy Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumEntropy Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::SumVariance : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumVariance Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumVariance Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::DifferenceAverage : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceAverage Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceAverage Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::DifferenceEntropy : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceEntropy Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceEntropy Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::DifferenceVariance : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceVariance Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceVariance Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::InverseDifferenceMomentNormalized : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMomentNormalized Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMomentNormalized Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::InverseDifferenceNormalized : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceNormalized Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceNormalized Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::InverseDifference : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifference Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifference Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::JointAverage : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") JointAverage Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") JointAverage Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::FirstMeasureOfInformationCorrelation : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") FirstMeasureOfInformationCorrelation Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") FirstMeasureOfInformationCorrelation Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::SecondMeasureOfInformationCorrelation : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SecondMeasureOfInformationCorrelation Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SecondMeasureOfInformationCorrelation Std.",featureStd->ElementAt(i))); + break; + default: + break; + } + }*/ +} + +mitk::GIFCooccurenceMatrix2::GIFCooccurenceMatrix2(): +m_Range(1.0), m_Direction(0) +{ +} + +mitk::GIFCooccurenceMatrix2::FeatureListType mitk::GIFCooccurenceMatrix2::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) +{ + FeatureListType featureList; + + GIFCooccurenceMatrix2Configuration config; + config.direction = m_Direction; + config.range = m_Range; + + AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config); + + return featureList; +} + +mitk::GIFCooccurenceMatrix2::FeatureNameListType mitk::GIFCooccurenceMatrix2::GetFeatureNames() +{ + FeatureNameListType featureList; + featureList.push_back("co-occ. Energy Means"); + featureList.push_back("co-occ. Energy Std."); + featureList.push_back("co-occ. Entropy Means"); + featureList.push_back("co-occ. Entropy Std."); + featureList.push_back("co-occ. Correlation Means"); + featureList.push_back("co-occ. Correlation Std."); + featureList.push_back("co-occ. InverseDifferenceMoment Means"); + featureList.push_back("co-occ. InverseDifferenceMoment Std."); + featureList.push_back("co-occ. Inertia Means"); + featureList.push_back("co-occ. Inertia Std."); + featureList.push_back("co-occ. ClusterShade Means"); + featureList.push_back("co-occ. ClusterShade Std."); + featureList.push_back("co-occ. ClusterProminence Means"); + featureList.push_back("co-occ. ClusterProminence Std."); + featureList.push_back("co-occ. HaralickCorrelation Means"); + featureList.push_back("co-occ. HaralickCorrelation Std."); + featureList.push_back("co-occ. Autocorrelation Means"); + featureList.push_back("co-occ. Autocorrelation Std."); + featureList.push_back("co-occ. Contrast Means"); + featureList.push_back("co-occ. Contrast Std."); + featureList.push_back("co-occ. Dissimilarity Means"); + featureList.push_back("co-occ. Dissimilarity Std."); + featureList.push_back("co-occ. MaximumProbability Means"); + featureList.push_back("co-occ. MaximumProbability Std."); + featureList.push_back("co-occ. InverseVariance Means"); + featureList.push_back("co-occ. InverseVariance Std."); + featureList.push_back("co-occ. Homogeneity1 Means"); + featureList.push_back("co-occ. Homogeneity1 Std."); + featureList.push_back("co-occ. ClusterTendency Means"); + featureList.push_back("co-occ. ClusterTendency Std."); + featureList.push_back("co-occ. Variance Means"); + featureList.push_back("co-occ. Variance Std."); + featureList.push_back("co-occ. SumAverage Means"); + featureList.push_back("co-occ. SumAverage Std."); + featureList.push_back("co-occ. SumEntropy Means"); + featureList.push_back("co-occ. SumEntropy Std."); + featureList.push_back("co-occ. SumVariance Means"); + featureList.push_back("co-occ. SumVariance Std."); + featureList.push_back("co-occ. DifferenceAverage Means"); + featureList.push_back("co-occ. DifferenceAverage Std."); + featureList.push_back("co-occ. DifferenceEntropy Means"); + featureList.push_back("co-occ. DifferenceEntropy Std."); + featureList.push_back("co-occ. DifferenceVariance Means"); + featureList.push_back("co-occ. DifferenceVariance Std."); + featureList.push_back("co-occ. InverseDifferenceMomentNormalized Means"); + featureList.push_back("co-occ. InverseDifferenceMomentNormalized Std."); + featureList.push_back("co-occ. InverseDifferenceNormalized Means"); + featureList.push_back("co-occ. InverseDifferenceNormalized Std."); + featureList.push_back("co-occ. InverseDifference Means"); + featureList.push_back("co-occ. InverseDifference Std."); + featureList.push_back("co-occ. JointAverage Means"); + featureList.push_back("co-occ. JointAverage Std."); + featureList.push_back("co-occ. FirstMeasurementOfInformationCorrelation Means"); + featureList.push_back("co-occ. FirstMeasurementOfInformationCorrelation Std."); + featureList.push_back("co-occ. SecondMeasurementOfInformationCorrelation Means"); + featureList.push_back("co-occ. SecondMeasurementOfInformationCorrelation Std."); + return featureList; +}