diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx index 9d5a963f4c..ff154fcf33 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx @@ -1,760 +1,760 @@ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ #ifndef __itkEnhancedHistogramToTextureFeaturesFilter_hxx #define __itkEnhancedHistogramToTextureFeaturesFilter_hxx #include "itkEnhancedHistogramToTextureFeaturesFilter.h" #include "itkNumericTraits.h" #include "vnl/vnl_math.h" #include "itkMath.h" #define itkMacroGLCMFeature(name, id) \ template< typename THistogram > \ 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); 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++ ) + for ( std::size_t i = 1; i <= inputHistogram->GetSize(0); i++ ) { double pi_tmp= 0.0; - for( int j = 1; j <= inputHistogram->GetSize(1); j++ ) + for( std::size_t 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++ ) + for ( std::size_t j = 1; j <= inputHistogram->GetSize(1); j++ ) { double pj_tmp= 0.0; - for( int i = 1; i <= inputHistogram->GetSize(0); i++ ) + for( std::size_t 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++ ) + for ( std::size_t 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 ( std::size_t i = 0; i < inputHistogram->GetSize(0); i++ ) { - for ( int j = 0; j < inputHistogram->GetSize(1); j++ ) + for ( std::size_t 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 ( std::size_t i = 0; i < inputHistogram->GetSize(0); i++ ) { - for ( int j = 0; j < inputHistogram->GetSize(1); j++ ) + for ( std::size_t 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/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx index 2f1343d659..8b3a423693 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx @@ -1,400 +1,400 @@ /*=================================================================== 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. ===================================================================*/ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ #ifndef __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter_hxx #define __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter_hxx #include "itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.h" #include "itkConstNeighborhoodIterator.h" #include "itkNeighborhood.h" #include "vnl/vnl_math.h" #include "itkMacro.h" #include "itkRescaleIntensityImageFilter.h" #include "itkMaskImageFilter.h" #include "itkLabelStatisticsImageFilter.h" #include "itkScalarConnectedComponentImageFilter.h" #include "itkRelabelComponentImageFilter.h" #include "itkCastImageFilter.h" #include namespace itk { namespace Statistics { template EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter() : m_NumberOfBinsPerAxis( itkGetStaticConstMacro( DefaultBinsPerAxis ) ), m_Min( NumericTraits::NonpositiveMin() ), m_Max( NumericTraits::max() ), m_MinDistance( NumericTraits::ZeroValue() ), m_MaxDistance( NumericTraits::max() ), m_InsidePixelValue( NumericTraits::OneValue() ) { this->SetNumberOfRequiredInputs( 1 ); this->SetNumberOfRequiredOutputs( 1 ); const unsigned int measurementVectorSize = 1; this->ProcessObject::SetNthOutput( 0, this->MakeOutput( 0 ) ); this->ProcessObject::SetNthOutput( 1, this->MakeOutput( 1 ) ); HistogramType *output = const_cast( this->GetOutput() ); output->SetMeasurementVectorSize( measurementVectorSize ); m_siMatrix = new double[m_NumberOfBinsPerAxis]; - for(int i = 0; i < m_NumberOfBinsPerAxis; i++) + for(unsigned int i = 0; i < m_NumberOfBinsPerAxis; i++) { m_siMatrix[i] = 0; } this->m_LowerBound.SetSize( measurementVectorSize ); this->m_UpperBound.SetSize( measurementVectorSize ); this->m_LowerBound[0] = this->m_Min; this->m_LowerBound[1] = this->m_MinDistance; this->m_UpperBound[0] = this->m_Max; this->m_UpperBound[1] = this->m_MaxDistance; } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::SetOffset( const OffsetType offset ) { OffsetVectorPointer offsetVector = OffsetVector::New(); offsetVector->push_back( offset ); this->SetOffsets( offsetVector ); MITK_WARN << "We now have " << this->GetOffsets()->size() << " offsets in matrixgenerator"; } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::AddOffsets( const std::vector _offsets ) { OffsetVectorPointer offsetVector = OffsetVector::New(); typename OffsetVector::ConstIterator offsets; //MITK_WARN << "We have " << this->GetOffsets()->size() << " offsets!"; - for( int i = 0; i < _offsets.size(); i++) + for( std::size_t i = 0; i < _offsets.size(); i++) { offsetVector->push_back(_offsets[i]); auto k = _offsets[i]; this->NormalizeOffsetDirection(k); offsetVector->push_back(k); } this->SetOffsets( offsetVector ); MITK_WARN << "We now have " << this->GetOffsets()->size() << " offsets in matrixgenerator"; } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::SetInput( const ImageType *image ) { // Process object is not const-correct so the const_cast is required here this->ProcessObject::SetNthInput( 0, const_cast( image ) ); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::SetMaskImage( const ImageType *image ) { // Process object is not const-correct so the const_cast is required here this->ProcessObject::SetNthInput( 1, const_cast( image ) ); } template const TImageType * EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetInput() const { if( this->GetNumberOfInputs() < 1 ) { return ITK_NULLPTR; } return static_cast( this->ProcessObject::GetInput( 0 ) ); } template const TImageType * EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetMaskImage() const { if( this->GetNumberOfInputs() < 2 ) { return ITK_NULLPTR; } return static_cast( this->ProcessObject::GetInput( 1 ) ); } template const typename EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter::HistogramType * EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetOutput() const { const HistogramType *output = static_cast( this->ProcessObject::GetOutput( 0 ) ); return output; } template double* EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetSiMatrix() const { return m_siMatrix; } template typename EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter::DataObjectPointer EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed( idx ) ) { return HistogramType::New().GetPointer(); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GenerateData() { HistogramType *output = static_cast( this->ProcessObject::GetOutput( 0 ) ); const ImageType * inputImage = this->GetInput(); const ImageType * maskImage = this->GetMaskImage(); // First, create an appropriate histogram with the right number of bins // and mins and maxes correct for the image type. typename HistogramType::SizeType size( output->GetMeasurementVectorSize() ); size.Fill( this->m_NumberOfBinsPerAxis ); this->m_LowerBound[0] = this->m_Min; this->m_LowerBound[1] = this->m_MinDistance; this->m_UpperBound[0] = this->m_Max; this->m_UpperBound[1] = this->m_MaxDistance; output->Initialize( size, this->m_LowerBound, this->m_UpperBound ); MeasurementVectorType run( output->GetMeasurementVectorSize() ); typename HistogramType::IndexType hIndex; //Cast the image to a float image - with no respect to the incoming image //to prevent some non-templated itk issues typedef itk::Image FloatImageType; typedef itk::CastImageFilter CastFilterType; typename CastFilterType::Pointer caster = CastFilterType::New(); caster->SetInput(inputImage); caster->Update(); typename FloatImageType::Pointer floatImage = caster->GetOutput(); //Cast the mask to an unsigned short image - with no respect to the incomimg maskimage //to prevent some non-templated itk issues typedef unsigned short LabelPixelType; typedef itk::Image LabelImageType; typedef itk::CastImageFilter MaskCastFilterType; typename MaskCastFilterType::Pointer maskCaster = MaskCastFilterType::New(); maskCaster->SetInput(maskImage); maskCaster->Update(); //Set all values out of the mask to nans. typedef itk::MaskImageFilter< FloatImageType, LabelImageType, FloatImageType > MaskFilterType; typename MaskFilterType::Pointer maskFilter = MaskFilterType::New(); maskFilter->SetInput(floatImage); maskFilter->SetMaskImage(maskCaster->GetOutput()); maskFilter->SetOutsideValue( std::numeric_limits::quiet_NaN()); maskFilter->Update(); FloatImageType::Pointer floatImageMasked = maskFilter->GetOutput(); typedef ConstNeighborhoodIterator NeighborhoodIteratorType; typename NeighborhoodIteratorType::RadiusType radius; radius.Fill( 1 ); NeighborhoodIteratorType neighborIt( radius, inputImage, inputImage->GetRequestedRegion() ); for( neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt ) { const PixelType centerPixelIntensity = neighborIt.GetCenterPixel(); IndexType centerIndex = neighborIt.GetIndex(); FloatImageType::IndexType cIndex; cIndex[0] = centerIndex[0]; cIndex[1] = centerIndex[1]; cIndex[2] = centerIndex[2]; float centerFloatPixel = floatImageMasked->GetPixel(cIndex); int px = 0; PixelType sum = 0.0; bool canCalculate = true; typename OffsetVector::ConstIterator offsets; for( offsets = this->GetOffsets()->Begin(); offsets != this->GetOffsets()->End(); offsets++ ) { OffsetType offset = offsets.Value(); IndexType index; index = centerIndex + offset; if(!inputImage->GetRequestedRegion().IsInside(index)) { canCalculate = false; break; } PixelType offsetPixel = inputImage->GetPixel(index); FloatImageType::IndexType fIndex; fIndex[0] = index[0]; fIndex[1] = index[1]; fIndex[2] = index[2]; float floatPixel = floatImageMasked->GetPixel(fIndex); //We have a nan here if(floatPixel != floatPixel || centerFloatPixel!= centerFloatPixel) { canCalculate = false; break; } sum += offsetPixel; px++; } //If we have a nan in the neighbourhood, continue if(!canCalculate) continue; PixelType mean = sum / px; double si = std::abs(mean-centerPixelIntensity); run[0] = centerPixelIntensity; //Check for NaN and inf if(run[0] == run[0] && !std::isinf(std::abs(run[0]))) { output->GetIndex( run, hIndex ); output->IncreaseFrequencyOfIndex( hIndex, 1 ); m_siMatrix[hIndex[0]] += si; } //MITK_WARN << " -> In this round we added: " << centerIndex << " with value " << centerPixelIntensity << " and si = " << si; //MITK_WARN << " -> Values are now siMatrix["< Values are now niMatrix: " << output->GetFrequency(hIndex) << "/" << run; } } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::SetPixelValueMinMax( PixelType min, PixelType max ) { if( this->m_Min != min || this->m_Max != max ) { itkDebugMacro( "setting Min to " << min << "and Max to " << max ); this->m_Min = min; this->m_Max = max; this->Modified(); } } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::SetDistanceValueMinMax( RealType min, RealType max ) { if( this->m_MinDistance != min || this->m_MaxDistance != max ) { itkDebugMacro( "setting MinDistance to " << min << "and MaxDistance to " << max ); this->m_MinDistance = min; this->m_MaxDistance = max; this->Modified(); } } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::PrintSelf( std::ostream& os, Indent indent ) const { Superclass::PrintSelf( os,indent ); os << indent << "Offsets: " << this->GetOffsets() << std::endl; os << indent << "Min: " << this->m_Min << std::endl; os << indent << "Max: " << this->m_Max << std::endl; os << indent << "Min distance: " << this->m_MinDistance << std::endl; os << indent << "Max distance: " << this->m_MaxDistance << std::endl; os << indent << "NumberOfBinsPerAxis: " << this->m_NumberOfBinsPerAxis << std::endl; os << indent << "InsidePixelValue: " << this->m_InsidePixelValue << std::endl; } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::NormalizeOffsetDirection(OffsetType &offset) { //MITK_WARN <<" -> NGTDM old offset = " << offset; itkDebugMacro("old offset = " << offset << std::endl); int sign = 1; bool metLastNonZero = false; for (int i = offset.GetOffsetDimension()-1; i>=0; i--) { if (metLastNonZero) { offset[i] *= sign; } else if (offset[i] != 0) { sign = (offset[i] > 0 ) ? 1 : -1; metLastNonZero = true; offset[i] *= sign; } } //MITK_WARN << " ->NGTDM new offset = " << offset; itkDebugMacro("new offset = " << offset << std::endl); } } // end of namespace Statistics } // end of namespace itk #endif diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp index 7752f3e433..fa8778bd9a 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix2.cpp @@ -1,604 +1,604 @@ #include // MITK #include #include #include // ITK #include #include #include #include // STL #include static void MatrixFeaturesTo(mitk::CoocurenceMatrixFeatures features, std::string prefix, mitk::GIFCooccurenceMatrix2::FeatureListType &featureList); static void CalculateMeanAndStdDevFeatures(std::vector featureList, mitk::CoocurenceMatrixFeatures &meanFeature, mitk::CoocurenceMatrixFeatures &stdFeature); static void NormalizeMatrixFeature(mitk::CoocurenceMatrixFeatures &features, std::size_t number); mitk::CoocurenceMatrixHolder::CoocurenceMatrixHolder(double min, double max, int number) : m_MinimumRange(min), m_MaximumRange(max), m_NumberOfBins(number) { m_Matrix.resize(number, number); m_Matrix.fill(0); m_Stepsize = (max - min) / (number); } int mitk::CoocurenceMatrixHolder::IntensityToIndex(double intensity) { return std::floor((intensity - m_MinimumRange) / m_Stepsize); } double mitk::CoocurenceMatrixHolder::IndexToMinIntensity(int index) { return m_MinimumRange + index * m_Stepsize; } double mitk::CoocurenceMatrixHolder::IndexToMeanIntensity(int index) { return m_MinimumRange + (index+0.5) * m_Stepsize; } double mitk::CoocurenceMatrixHolder::IndexToMaxIntensity(int index) { return m_MinimumRange + (index + 1) * m_Stepsize; } template void CalculateCoOcMatrix(itk::Image* itkImage, itk::Image* mask, itk::Offset offset, int range, mitk::CoocurenceMatrixHolder &holder) { typedef itk::Image ImageType; typedef itk::Image MaskImageType; typedef itk::ShapedNeighborhoodIterator ShapeIterType; typedef itk::ShapedNeighborhoodIterator ShapeMaskIterType; typedef itk::ImageRegionConstIterator ConstIterType; typedef itk::ImageRegionConstIterator ConstMaskIterType; itk::Size radius; radius.Fill(range+1); ShapeIterType imageOffsetIter(radius, itkImage, itkImage->GetLargestPossibleRegion()); ShapeMaskIterType maskOffsetIter(radius, mask, mask->GetLargestPossibleRegion()); imageOffsetIter.ActivateOffset(offset); maskOffsetIter.ActivateOffset(offset); ConstIterType imageIter(itkImage, itkImage->GetLargestPossibleRegion()); ConstMaskIterType maskIter(mask, mask->GetLargestPossibleRegion()); // iterator.GetIndex() + ci.GetNeighborhoodOffset() auto region = mask->GetLargestPossibleRegion(); while (!maskIter.IsAtEnd()) { auto ciMask = maskOffsetIter.Begin(); auto ciValue = imageOffsetIter.Begin(); if (maskIter.Value() > 0 && ciMask.Get() > 0 && imageIter.Get() == imageIter.Get() && ciValue.Get() == ciValue.Get() && region.IsInside(maskOffsetIter.GetIndex() + ciMask.GetNeighborhoodOffset())) { int i = holder.IntensityToIndex(imageIter.Get()); int j = holder.IntensityToIndex(ciValue.Get()); holder.m_Matrix(i, j) += 1; holder.m_Matrix(j, i) += 1; } ++imageOffsetIter; ++maskOffsetIter; ++imageIter; ++maskIter; } } void CalculateFeatures( mitk::CoocurenceMatrixHolder &holder, mitk::CoocurenceMatrixFeatures & results ) { auto pijMatrix = holder.m_Matrix; auto piMatrix = holder.m_Matrix; auto pjMatrix = holder.m_Matrix; double Ng = holder.m_NumberOfBins; int NgSize = holder.m_NumberOfBins; pijMatrix /= pijMatrix.sum(); piMatrix.rowwise().normalize(); pjMatrix.colwise().normalize(); for (int i = 0; i < holder.m_NumberOfBins; ++i) for (int j = 0; j < holder.m_NumberOfBins; ++j) { if (pijMatrix(i, j) != pijMatrix(i, j)) pijMatrix(i, j) = 0; if (piMatrix(i, j) != piMatrix(i, j)) piMatrix(i, j) = 0; if (pjMatrix(i, j) != pjMatrix(i, j)) pjMatrix(i, j) = 0; } Eigen::VectorXd piVector = pijMatrix.colwise().sum(); Eigen::VectorXd pjVector = pijMatrix.rowwise().sum(); double sigmai = 0;; for (int i = 0; i < holder.m_NumberOfBins; ++i) { double iInt = holder.IndexToMeanIntensity(i); results.RowAverage += iInt * piVector(i); if (piVector(i) > 0) { results.RowEntropy -= piVector(i) * std::log(piVector(i)) / std::log(2); } } for (int i = 0; i < holder.m_NumberOfBins; ++i) { double iInt = holder.IndexToMeanIntensity(i); results.RowVariance += (iInt - results.RowAverage)*(iInt - results.RowAverage) * piVector(i); } results.RowMaximum = piVector.maxCoeff(); sigmai = std::sqrt(results.RowVariance); Eigen::VectorXd pimj(NgSize); pimj.fill(0); Eigen::VectorXd pipj(2*NgSize); pipj.fill(0); results.JointMaximum += pijMatrix.maxCoeff(); for (int i = 0; i < holder.m_NumberOfBins; ++i) { for (int j = 0; j < holder.m_NumberOfBins; ++j) { double iInt = holder.IndexToMeanIntensity(i); double jInt = holder.IndexToMeanIntensity(j); double pij = pijMatrix(i, j); int deltaK = (i - j)>0?(i-j) : (j-i); pimj(deltaK) += pij; pipj(i + j) += pij; results.JointAverage += iInt * pij; if (pij > 0) { results.JointEntropy -= pij * std::log(pij) / std::log(2); results.FirstRowColumnEntropy -= pij * std::log(piVector(i)*pjVector(j)) / std::log(2); } if (piVector(i) > 0 && pjVector(j) > 0 ) { results.SecondRowColumnEntropy -= piVector(i)*pjVector(j) * std::log(piVector(i)*pjVector(j)) / std::log(2); } results.AngularSecondMoment += pij*pij; results.Contrast += (iInt - jInt)* (iInt - jInt) * pij; results.Dissimilarity += std::abs(iInt - jInt) * pij; results.InverseDifference += pij / (1 + (std::abs(iInt - jInt))); results.InverseDifferenceNormalised += pij / (1 + (std::abs(iInt - jInt) / Ng)); results.InverseDifferenceMoment += pij / (1 + (iInt - jInt)*(iInt - jInt)); results.InverseDifferenceMomentNormalised += pij / (1 + (iInt - jInt)*(iInt - jInt)/Ng/Ng); results.Autocorrelation += iInt*jInt * pij; double cluster = (iInt + jInt - 2 * results.RowAverage); results.ClusterTendency += cluster*cluster * pij; results.ClusterShade += cluster*cluster*cluster * pij; results.ClusterProminence += cluster*cluster*cluster*cluster * pij; if (iInt != jInt) { results.InverseVariance += pij / (iInt - jInt) / (iInt - jInt); } } } results.Correlation = 1 / sigmai / sigmai * (-results.RowAverage*results.RowAverage+ results.Autocorrelation); results.FirstMeasureOfInformationCorrelation = (results.JointEntropy - results.FirstRowColumnEntropy) / results.RowEntropy; if (results.JointEntropy < results.SecondRowColumnEntropy) { results.SecondMeasureOfInformationCorrelation = sqrt(1 - exp(-2 * (results.SecondRowColumnEntropy - results.JointEntropy))); } else { results.SecondMeasureOfInformationCorrelation = 0; } for (int i = 0; i < holder.m_NumberOfBins; ++i) { for (int j = 0; j < holder.m_NumberOfBins; ++j) { double iInt = holder.IndexToMeanIntensity(i); //double jInt = holder.IndexToMeanIntensity(j); double pij = pijMatrix(i, j); results.JointVariance += (iInt - results.JointAverage)* (iInt - results.JointAverage)*pij; } } for (int k = 0; k < NgSize; ++k) { results.DifferenceAverage += k* pimj(k); if (pimj(k) > 0) { results.DifferenceEntropy -= pimj(k) * log(pimj(k)) / std::log(2); } } for (int k = 0; k < NgSize; ++k) { results.DifferenceVariance += (results.DifferenceAverage-k)* (results.DifferenceAverage-k)*pimj(k); } for (int k = 0; k <2* NgSize ; ++k) { results.SumAverage += (2+k)* pipj(k); if (pipj(k) > 0) { results.SumEntropy -= pipj(k) * log(pipj(k)) / std::log(2); } } for (int k = 0; k < 2*NgSize; ++k) { results.SumVariance += (2+k - results.SumAverage)* (2+k - results.SumAverage)*pipj(k); } //MITK_INFO << std::endl << holder.m_Matrix; //MITK_INFO << std::endl << pijMatrix; //MITK_INFO << std::endl << piMatrix; //MITK_INFO << std::endl << pjMatrix; //for (int i = 0; i < holder.m_NumberOfBins; ++i) //{ // MITK_INFO << "Bin " << i << " Min: " << holder.IndexToMinIntensity(i) << " Max: " << holder.IndexToMaxIntensity(i); //} //MITK_INFO << pimj; //MITK_INFO << pipj; } template void CalculateCoocurenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFCooccurenceMatrix2::FeatureListType & featureList, mitk::GIFCooccurenceMatrix2::GIFCooccurenceMatrix2Configuration config) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef itk::Neighborhood NeighborhoodType; typedef itk::Offset OffsetType; /////////////////////////////////////////////////////////////////////////////////////////////// typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); double rangeMin = -0.5 + minMaxComputer->GetMinimum(); double rangeMax = 0.5 + minMaxComputer->GetMaximum(); if (config.UseMinimumIntensity) rangeMin = config.MinimumIntensity; if (config.UseMaximumIntensity) rangeMax = config.MaximumIntensity; // Define Range int numberOfBins = config.Bins; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); //Find possible directions std::vector < itk::Offset > offsetVector; NeighborhoodType hood; hood.SetRadius(1); unsigned int centerIndex = hood.GetCenterNeighborhoodIndex(); OffsetType offset; for (unsigned int d = 0; d < centerIndex; d++) { offset = hood.GetOffset(d); bool useOffset = true; - for (int i = 0; i < VImageDimension; ++i) + 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 (int i = 0; i < offsetVector.size(); ++i) + 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-occ. 2 (" + strRange + ") overall", featureList); MatrixFeaturesTo(featureMean, "co-occ. 2 (" + strRange + ") mean", featureList); MatrixFeaturesTo(featureStd, "co-occ. 2 (" + strRange + ") std.dev.", featureList); } static void MatrixFeaturesTo(mitk::CoocurenceMatrixFeatures features, std::string prefix, mitk::GIFCooccurenceMatrix2::FeatureListType &featureList) { featureList.push_back(std::make_pair(prefix + " Joint Maximum", features.JointMaximum)); featureList.push_back(std::make_pair(prefix + " Joint Average", features.JointAverage)); featureList.push_back(std::make_pair(prefix + " Joint Variance", features.JointVariance)); featureList.push_back(std::make_pair(prefix + " Joint Entropy", features.JointEntropy)); featureList.push_back(std::make_pair(prefix + " Row Maximum", features.RowMaximum)); featureList.push_back(std::make_pair(prefix + " Row Average", features.RowAverage)); featureList.push_back(std::make_pair(prefix + " Row Variance", features.RowVariance)); featureList.push_back(std::make_pair(prefix + " Row Entropy", features.RowEntropy)); featureList.push_back(std::make_pair(prefix + " First Row-Column Entropy", features.FirstRowColumnEntropy)); featureList.push_back(std::make_pair(prefix + " Second Row-Column Entropy", features.SecondRowColumnEntropy)); featureList.push_back(std::make_pair(prefix + " Difference Average", features.DifferenceAverage)); featureList.push_back(std::make_pair(prefix + " Difference Variance", features.DifferenceVariance)); featureList.push_back(std::make_pair(prefix + " Difference Entropy", features.DifferenceEntropy)); featureList.push_back(std::make_pair(prefix + " Sum Average", features.SumAverage)); featureList.push_back(std::make_pair(prefix + " Sum Variance", features.SumVariance)); featureList.push_back(std::make_pair(prefix + " Sum Entropy", features.SumEntropy)); featureList.push_back(std::make_pair(prefix + " Angular Second Moment", features.AngularSecondMoment)); featureList.push_back(std::make_pair(prefix + " Contrast", features.Contrast)); featureList.push_back(std::make_pair(prefix + " Dissimilarity", features.Dissimilarity)); featureList.push_back(std::make_pair(prefix + " Inverse Difference", features.InverseDifference)); featureList.push_back(std::make_pair(prefix + " Inverse Difference Normalized", features.InverseDifferenceNormalised)); featureList.push_back(std::make_pair(prefix + " Inverse Difference Moment", features.InverseDifferenceMoment)); featureList.push_back(std::make_pair(prefix + " Inverse Difference Moment Normalized", features.InverseDifferenceMomentNormalised)); featureList.push_back(std::make_pair(prefix + " Inverse Variance", features.InverseVariance)); featureList.push_back(std::make_pair(prefix + " Correlation", features.Correlation)); featureList.push_back(std::make_pair(prefix + " Autocorrleation", features.Autocorrelation)); featureList.push_back(std::make_pair(prefix + " Cluster Tendency", features.ClusterTendency)); featureList.push_back(std::make_pair(prefix + " Cluster Shade", features.ClusterShade)); featureList.push_back(std::make_pair(prefix + " Cluster Prominence", features.ClusterProminence)); featureList.push_back(std::make_pair(prefix + " First Measure of Information Correlation", features.FirstMeasureOfInformationCorrelation)); featureList.push_back(std::make_pair(prefix + " Second Measure of Information Correlation", features.SecondMeasureOfInformationCorrelation)); } static void CalculateMeanAndStdDevFeatures(std::vector featureList, mitk::CoocurenceMatrixFeatures &meanFeature, mitk::CoocurenceMatrixFeatures &stdFeature) { #define ADDFEATURE(a) meanFeature.a += featureList[i].a;stdFeature.a += featureList[i].a*featureList[i].a #define CALCVARIANCE(a) stdFeature.a =sqrt(stdFeature.a - meanFeature.a*meanFeature.a) for (std::size_t i = 0; i < featureList.size(); ++i) { ADDFEATURE(JointMaximum); ADDFEATURE(JointAverage); ADDFEATURE(JointVariance); ADDFEATURE(JointEntropy); ADDFEATURE(RowMaximum); ADDFEATURE(RowAverage); ADDFEATURE(RowVariance); ADDFEATURE(RowEntropy); ADDFEATURE(FirstRowColumnEntropy); ADDFEATURE(SecondRowColumnEntropy); ADDFEATURE(DifferenceAverage); ADDFEATURE(DifferenceVariance); ADDFEATURE(DifferenceEntropy); ADDFEATURE(SumAverage); ADDFEATURE(SumVariance); ADDFEATURE(SumEntropy); ADDFEATURE(AngularSecondMoment); ADDFEATURE(Contrast); ADDFEATURE(Dissimilarity); ADDFEATURE(InverseDifference); ADDFEATURE(InverseDifferenceNormalised); ADDFEATURE(InverseDifferenceMoment); ADDFEATURE(InverseDifferenceMomentNormalised); ADDFEATURE(InverseVariance); ADDFEATURE(Correlation); ADDFEATURE(Autocorrelation); ADDFEATURE(ClusterShade); ADDFEATURE(ClusterTendency); ADDFEATURE(ClusterProminence); ADDFEATURE(FirstMeasureOfInformationCorrelation); ADDFEATURE(SecondMeasureOfInformationCorrelation); } NormalizeMatrixFeature(meanFeature, featureList.size()); NormalizeMatrixFeature(stdFeature, featureList.size()); CALCVARIANCE(JointMaximum); CALCVARIANCE(JointAverage); CALCVARIANCE(JointVariance); CALCVARIANCE(JointEntropy); CALCVARIANCE(RowMaximum); CALCVARIANCE(RowAverage); CALCVARIANCE(RowVariance); CALCVARIANCE(RowEntropy); CALCVARIANCE(FirstRowColumnEntropy); CALCVARIANCE(SecondRowColumnEntropy); CALCVARIANCE(DifferenceAverage); CALCVARIANCE(DifferenceVariance); CALCVARIANCE(DifferenceEntropy); CALCVARIANCE(SumAverage); CALCVARIANCE(SumVariance); CALCVARIANCE(SumEntropy); CALCVARIANCE(AngularSecondMoment); CALCVARIANCE(Contrast); CALCVARIANCE(Dissimilarity); CALCVARIANCE(InverseDifference); CALCVARIANCE(InverseDifferenceNormalised); CALCVARIANCE(InverseDifferenceMoment); CALCVARIANCE(InverseDifferenceMomentNormalised); CALCVARIANCE(InverseVariance); CALCVARIANCE(Correlation); CALCVARIANCE(Autocorrelation); CALCVARIANCE(ClusterShade); CALCVARIANCE(ClusterTendency); CALCVARIANCE(ClusterProminence); CALCVARIANCE(FirstMeasureOfInformationCorrelation); CALCVARIANCE(SecondMeasureOfInformationCorrelation); #undef ADDFEATURE #undef CALCVARIANCE } static void NormalizeMatrixFeature(mitk::CoocurenceMatrixFeatures &features, std::size_t number) { features.JointMaximum = features.JointMaximum / number; features.JointAverage = features.JointAverage / number; features.JointVariance = features.JointVariance / number; features.JointEntropy = features.JointEntropy / number; features.RowMaximum = features.RowMaximum / number; features.RowAverage = features.RowAverage / number; features.RowVariance = features.RowVariance / number; features.RowEntropy = features.RowEntropy / number; features.FirstRowColumnEntropy = features.FirstRowColumnEntropy / number; features.SecondRowColumnEntropy = features.SecondRowColumnEntropy / number; features.DifferenceAverage = features.DifferenceAverage / number; features.DifferenceVariance = features.DifferenceVariance / number; features.DifferenceEntropy = features.DifferenceEntropy / number; features.SumAverage = features.SumAverage / number; features.SumVariance = features.SumVariance / number; features.SumEntropy = features.SumEntropy / number; features.AngularSecondMoment = features.AngularSecondMoment / number; features.Contrast = features.Contrast / number; features.Dissimilarity = features.Dissimilarity / number; features.InverseDifference = features.InverseDifference / number; features.InverseDifferenceNormalised = features.InverseDifferenceNormalised / number; features.InverseDifferenceMoment = features.InverseDifferenceMoment / number; features.InverseDifferenceMomentNormalised = features.InverseDifferenceMomentNormalised / number; features.InverseVariance = features.InverseVariance / number; features.Correlation = features.Correlation / number; features.Autocorrelation = features.Autocorrelation / number; features.ClusterShade = features.ClusterShade / number; features.ClusterTendency = features.ClusterTendency / number; features.ClusterProminence = features.ClusterProminence / number; features.FirstMeasureOfInformationCorrelation = features.FirstMeasureOfInformationCorrelation / number; features.SecondMeasureOfInformationCorrelation = features.SecondMeasureOfInformationCorrelation / number; } mitk::GIFCooccurenceMatrix2::GIFCooccurenceMatrix2(): m_Range(1.0), m_Bins(128) { SetShortName("cooc2"); SetLongName("cooccurence2"); } mitk::GIFCooccurenceMatrix2::FeatureListType mitk::GIFCooccurenceMatrix2::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; GIFCooccurenceMatrix2Configuration config; config.direction = GetDirection(); config.range = m_Range; config.MinimumIntensity = GetMinimumIntensity(); config.MaximumIntensity = GetMaximumIntensity(); config.UseMinimumIntensity = GetUseMinimumIntensity(); config.UseMaximumIntensity = GetUseMaximumIntensity(); config.Bins = GetBins(); AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config); return featureList; } mitk::GIFCooccurenceMatrix2::FeatureNameListType mitk::GIFCooccurenceMatrix2::GetFeatureNames() { FeatureNameListType featureList; return featureList; } void mitk::GIFCooccurenceMatrix2::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Co-occurence matrix", "calculates Co-occurence based features (new implementation)", us::Any()); parser.addArgument(name+"::range", name+"::range", mitkCommandLineParser::String, "Cooc 2 Range", "Define the range that is used (Semicolon-separated)", us::Any()); parser.addArgument(name + "::bins", name + "::bins", mitkCommandLineParser::String, "Cooc 2 Number of Bins", "Define the number of bins that is used ", us::Any()); } void mitk::GIFCooccurenceMatrix2::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); std::string name = GetOptionPrefix(); if (parsedArgs.count(GetLongName())) { std::vector ranges; if (parsedArgs.count(name + "::range")) { ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); } else { ranges.push_back(1); } if (parsedArgs.count(name + "::bins")) { auto bins = SplitDouble(parsedArgs[name + "::bins"].ToString(), ';')[0]; this->SetBins(bins); } for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; this->SetRange(ranges[i]); auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; } } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp index 464734c7d3..9c39868704 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp @@ -1,404 +1,404 @@ #include // MITK #include #include #include // ITK #include #include #include #include // STL #include static void MatrixFeaturesTo(mitk::NGLDMMatrixFeatures features, std::string prefix, mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType &featureList); mitk::NGLDMMatrixHolder::NGLDMMatrixHolder(double min, double max, int number, int depenence) : m_MinimumRange(min), m_MaximumRange(max), m_Stepsize(0), m_NumberOfDependences(depenence), m_NumberOfBins(number), m_NeighbourhoodSize(1), m_NumberOfNeighbourVoxels(0), m_NumberOfDependenceNeighbourVoxels(0), m_NumberOfNeighbourhoods(0), m_NumberOfCompleteNeighbourhoods(0) { m_Matrix.resize(number, depenence); m_Matrix.fill(0); m_Stepsize = (max - min) / (number); } int mitk::NGLDMMatrixHolder::IntensityToIndex(double intensity) { return std::floor((intensity - m_MinimumRange) / m_Stepsize); } double mitk::NGLDMMatrixHolder::IndexToMinIntensity(int index) { return m_MinimumRange + index * m_Stepsize; } double mitk::NGLDMMatrixHolder::IndexToMeanIntensity(int index) { return m_MinimumRange + (index+0.5) * m_Stepsize; } double mitk::NGLDMMatrixHolder::IndexToMaxIntensity(int index) { return m_MinimumRange + (index + 1) * m_Stepsize; } template void CalculateNGLDMMatrix(itk::Image* itkImage, itk::Image* mask, int alpha, int range, - int direction, + unsigned int direction, mitk::NGLDMMatrixHolder &holder) { typedef itk::Image ImageType; typedef itk::Image MaskImageType; typedef itk::NeighborhoodIterator ShapeIterType; typedef itk::NeighborhoodIterator ShapeMaskIterType; holder.m_NumberOfCompleteNeighbourhoods = 0; holder.m_NumberOfNeighbourhoods = 0; holder.m_NumberOfNeighbourVoxels = 0; holder.m_NumberOfDependenceNeighbourVoxels = 0; itk::Size radius; radius.Fill(range); if ((direction > 1) && (direction +2 GetLargestPossibleRegion()); ShapeMaskIterType maskIter(radius, mask, mask->GetLargestPossibleRegion()); auto region = mask->GetLargestPossibleRegion(); auto center = imageIter.Size() / 2; auto iterSize = imageIter.Size(); holder.m_NeighbourhoodSize = iterSize-1; while (!maskIter.IsAtEnd()) { int sameValues = 0; bool completeNeighbourhood = true; int i = holder.IntensityToIndex(imageIter.GetCenterPixel()); if (imageIter.GetCenterPixel() != imageIter.GetCenterPixel()) { ++imageIter; ++maskIter; continue; } - for (int position = 0; position < iterSize; ++position) + for (unsigned int position = 0; position < iterSize; ++position) { if ((position == center) || ( ! region.IsInside(maskIter.GetIndex(position)))) { completeNeighbourhood = false; continue; } bool isInBounds; auto jIntensity = imageIter.GetPixel(position, isInBounds); auto jMask = maskIter.GetPixel(position, isInBounds); if (jMask < 0 || (jIntensity != jIntensity) || ( ! isInBounds)) { completeNeighbourhood = false; continue; } int j = holder.IntensityToIndex(jIntensity); holder.m_NumberOfNeighbourVoxels += 1; if (std::abs(i - j) <= alpha) { holder.m_NumberOfDependenceNeighbourVoxels += 1; ++sameValues; } } holder.m_Matrix(i, sameValues) += 1; holder.m_NumberOfNeighbourhoods += 1; if (completeNeighbourhood) { holder.m_NumberOfCompleteNeighbourhoods += 1; } ++imageIter; ++maskIter; } } void LocalCalculateFeatures( mitk::NGLDMMatrixHolder &holder, mitk::NGLDMMatrixFeatures & results ) { auto sijMatrix = holder.m_Matrix; auto piMatrix = holder.m_Matrix; auto pjMatrix = holder.m_Matrix; // double Ng = holder.m_NumberOfBins; // int NgSize = holder.m_NumberOfBins; double Ns = sijMatrix.sum(); piMatrix.rowwise().normalize(); pjMatrix.colwise().normalize(); for (int i = 0; i < holder.m_NumberOfBins; ++i) { double sj = 0; for (int j = 0; j < holder.m_NumberOfDependences; ++j) { double iInt = holder.IndexToMeanIntensity(i); double sij = sijMatrix(i, j); double k = j+1; double pij = sij / Ns; results.LowDependenceEmphasis += sij / k / k; results.HighDependenceEmphasis += sij * k*k; if (iInt != 0) { results.LowGreyLevelCountEmphasis += sij / iInt / iInt; } results.HighGreyLevelCountEmphasis += sij * iInt*iInt; if (iInt != 0) { results.LowDependenceLowGreyLevelEmphasis += sij / k / k / iInt / iInt; } results.LowDependenceHighGreyLevelEmphasis += sij * iInt*iInt / k / k; if (iInt != 0) { results.HighDependenceLowGreyLevelEmphasis += sij *k * k / iInt / iInt; } results.HighDependenceHighGreyLevelEmphasis += sij * k*k*iInt*iInt; results.MeanGreyLevelCount += iInt * pij; results.MeanDependenceCount += k * pij; if (pij > 0) { results.DependenceCountEntropy -= pij * std::log(pij) / std::log(2); } sj += sij; } results.GreyLevelNonUniformity += sj*sj; results.GreyLevelNonUniformityNormalised += sj*sj; } for (int j = 0; j < holder.m_NumberOfDependences; ++j) { double si = 0; for (int i = 0; i < holder.m_NumberOfBins; ++i) { double sij = sijMatrix(i, j); si += sij; } results.DependenceCountNonUniformity += si*si; results.DependenceCountNonUniformityNormalised += si*si; } for (int i = 0; i < holder.m_NumberOfBins; ++i) { for (int j = 0; j < holder.m_NumberOfDependences; ++j) { double iInt = holder.IndexToMeanIntensity(i); double sij = sijMatrix(i, j); double k = j + 1; double pij = sij / Ns; results.GreyLevelVariance += (iInt - results.MeanGreyLevelCount)* (iInt - results.MeanGreyLevelCount) * pij; results.DependenceCountVariance += (k - results.MeanDependenceCount)* (k - results.MeanDependenceCount) * pij; } } results.LowDependenceEmphasis /= Ns; results.HighDependenceEmphasis /= Ns; results.LowGreyLevelCountEmphasis /= Ns; results.HighGreyLevelCountEmphasis /= Ns; results.LowDependenceLowGreyLevelEmphasis /= Ns; results.LowDependenceHighGreyLevelEmphasis /= Ns; results.HighDependenceLowGreyLevelEmphasis /= Ns; results.HighDependenceHighGreyLevelEmphasis /= Ns; results.GreyLevelNonUniformity /= Ns; results.GreyLevelNonUniformityNormalised /= (Ns*Ns); results.DependenceCountNonUniformity /= Ns; results.DependenceCountNonUniformityNormalised /= (Ns*Ns); results.DependenceCountPercentage = 1; results.ExpectedNeighbourhoodSize = holder.m_NeighbourhoodSize; results.AverageNeighbourhoodSize = holder.m_NumberOfNeighbourVoxels / (1.0 * holder.m_NumberOfNeighbourhoods); results.AverageIncompleteNeighbourhoodSize = (holder.m_NumberOfNeighbourVoxels - holder.m_NumberOfCompleteNeighbourhoods* holder.m_NeighbourhoodSize) / (1.0 * (holder.m_NumberOfNeighbourhoods - holder.m_NumberOfCompleteNeighbourhoods)); results.PercentageOfCompleteNeighbourhoods = holder.m_NumberOfCompleteNeighbourhoods / (1.0 * holder.m_NumberOfNeighbourhoods); results.PercentageOfDependenceNeighbours = holder.m_NumberOfDependenceNeighbourVoxels / (1.0 * holder.m_NumberOfNeighbourVoxels); } template void CalculateCoocurenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType & featureList, mitk::GIFNeighbouringGreyLevelDependenceFeature::GIFNeighbouringGreyLevelDependenceFeatureConfiguration config) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; /////////////////////////////////////////////////////////////////////////////////////////////// typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); double rangeMin = -0.5 + minMaxComputer->GetMinimum(); double rangeMax = 0.5 + minMaxComputer->GetMaximum(); if (config.UseMinimumIntensity) rangeMin = config.MinimumIntensity; if (config.UseMaximumIntensity) rangeMax = config.MaximumIntensity; // Define Range int numberOfBins = config.Bins; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); std::vector resultVector; mitk::NGLDMMatrixHolder holderOverall(rangeMin, rangeMax, numberOfBins,37); mitk::NGLDMMatrixFeatures overallFeature; CalculateNGLDMMatrix(itkImage, maskImage, config.alpha, config.range, config.direction, holderOverall); LocalCalculateFeatures(holderOverall, overallFeature); std::ostringstream ss; ss << config.range; std::string strRange = ss.str(); MatrixFeaturesTo(overallFeature, "NGLDM (" + strRange + ") overall", featureList); } static void MatrixFeaturesTo(mitk::NGLDMMatrixFeatures features, std::string prefix, mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType &featureList) { featureList.push_back(std::make_pair(prefix + " Low Dependence Emphasis", features.LowDependenceEmphasis)); featureList.push_back(std::make_pair(prefix + " High Dependence Emphasis", features.HighDependenceEmphasis)); featureList.push_back(std::make_pair(prefix + " Low Grey Level Count Emphasis", features.LowGreyLevelCountEmphasis)); featureList.push_back(std::make_pair(prefix + " High Grey Level Count Emphasis", features.HighGreyLevelCountEmphasis)); featureList.push_back(std::make_pair(prefix + " Low Dependence Low Grey Level Emphasis", features.LowDependenceLowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + " Low Dependence High Grey Level Emphasis", features.LowDependenceHighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + " High Dependence Low Grey Level Emphasis", features.HighDependenceLowGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + " High Dependence High Grey Level Emphasis", features.HighDependenceHighGreyLevelEmphasis)); featureList.push_back(std::make_pair(prefix + " Grey Level Non-Uniformity", features.GreyLevelNonUniformity)); featureList.push_back(std::make_pair(prefix + " Grey Level Non-Uniformity Normalised", features.GreyLevelNonUniformityNormalised)); featureList.push_back(std::make_pair(prefix + " Dependence Count Non-Uniformity", features.DependenceCountNonUniformity)); featureList.push_back(std::make_pair(prefix + " Dependence Count Non-Uniformtiy Normalised", features.DependenceCountNonUniformityNormalised)); featureList.push_back(std::make_pair(prefix + " Dependence Count Percentage", features.DependenceCountPercentage)); featureList.push_back(std::make_pair(prefix + " Grey Level Mean", features.MeanGreyLevelCount)); featureList.push_back(std::make_pair(prefix + " Grey Level Variance", features.GreyLevelVariance)); featureList.push_back(std::make_pair(prefix + " Dependence Count Mean", features.MeanDependenceCount)); featureList.push_back(std::make_pair(prefix + " Dependence Count Variance", features.DependenceCountVariance)); featureList.push_back(std::make_pair(prefix + " Dependence Count Entropy", features.DependenceCountEntropy)); featureList.push_back(std::make_pair(prefix + " Expected Neighbourhood Size", features.ExpectedNeighbourhoodSize)); featureList.push_back(std::make_pair(prefix + " Average Neighbourhood Size", features.AverageNeighbourhoodSize)); featureList.push_back(std::make_pair(prefix + " Average Incomplete Neighbourhood Size", features.AverageIncompleteNeighbourhoodSize)); featureList.push_back(std::make_pair(prefix + " Percentage of complete Neighbourhoods", features.PercentageOfCompleteNeighbourhoods)); featureList.push_back(std::make_pair(prefix + " Percentage of Dependence Neighbour Voxels", features.PercentageOfDependenceNeighbours)); } mitk::GIFNeighbouringGreyLevelDependenceFeature::GIFNeighbouringGreyLevelDependenceFeature() : m_Range(1.0), m_Bins(6) { SetShortName("ngldm"); SetLongName("ngldm"); } mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType mitk::GIFNeighbouringGreyLevelDependenceFeature::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; GIFNeighbouringGreyLevelDependenceFeatureConfiguration config; config.direction = GetDirection(); config.range = m_Range; config.alpha = 0; config.MinimumIntensity = GetMinimumIntensity(); config.MaximumIntensity = GetMaximumIntensity(); config.UseMinimumIntensity = GetUseMinimumIntensity(); config.UseMaximumIntensity = GetUseMaximumIntensity(); config.Bins = GetBins(); AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config); return featureList; } mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureNameListType mitk::GIFNeighbouringGreyLevelDependenceFeature::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("co-occ. Energy Means"); featureList.push_back("co-occ. Energy Std."); return featureList; } void mitk::GIFNeighbouringGreyLevelDependenceFeature::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Calculate Neighbouring grey level dependence based features", "Calculate Neighbouring grey level dependence based features", us::Any()); parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::String, "NGLD Range", "Define the range that is used (Semicolon-separated)", us::Any()); parser.addArgument(name + "::bins", name + "::bins", mitkCommandLineParser::String, "NGLD Number of Bins", "Define the number of bins that is used ", us::Any()); } void mitk::GIFNeighbouringGreyLevelDependenceFeature::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &, const Image::Pointer &maskNoNAN, FeatureListType &featureList) { auto parsedArgs = GetParameter(); std::string name = GetOptionPrefix(); if (parsedArgs.count(GetLongName())) { std::vector ranges; if (parsedArgs.count(name + "::range")) { ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';'); } else { ranges.push_back(1); } if (parsedArgs.count(name + "::bins")) { auto bins = SplitDouble(parsedArgs[name + "::bins"].ToString(), ';')[0]; this->SetBins(bins); } for (std::size_t i = 0; i < ranges.size(); ++i) { MITK_INFO << "Start calculating coocurence with range " << ranges[i] << "...."; this->SetRange(ranges[i]); auto localResults = this->CalculateFeatures(feature, maskNoNAN); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating coocurence with range " << ranges[i] << "...."; } } } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp index 027758b5db..0734a02c5e 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp @@ -1,431 +1,431 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include #include #include // ITK #include #include // VTK #include #include #include // STL #include #include // OpenCV #include template void CalculateVolumeStatistic(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFVolumetricStatistics::FeatureListType & featureList) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::LabelStatisticsImageFilter FilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename FilterType::Pointer labelStatisticsImageFilter = FilterType::New(); labelStatisticsImageFilter->SetInput( itkImage ); labelStatisticsImageFilter->SetLabelInput(maskImage); labelStatisticsImageFilter->Update(); double volume = labelStatisticsImageFilter->GetCount(1); double voxelVolume = 1; for (int i = 0; i < (int)(VImageDimension); ++i) { volume *= itkImage->GetSpacing()[i]; voxelVolume *= itkImage->GetSpacing()[i]; } featureList.push_back(std::make_pair("Volumetric Features Volume (pixel based)", volume)); featureList.push_back(std::make_pair("Volumetric Features Voxel Volume", voxelVolume)); } template void CalculateLargestDiameter(itk::Image* mask, mitk::Image::Pointer valueImage, mitk::GIFVolumetricStatistics::FeatureListType & featureList) { typedef itk::Image ValueImageType; typename ValueImageType::Pointer itkValueImage = ValueImageType::New(); mitk::CastToItkImage(valueImage, itkValueImage); typedef itk::Image ImageType; typedef typename ImageType::PointType PointType; typename ImageType::SizeType radius; for (int i=0; i < (int)VImageDimension; ++i) radius[i] = 1; itk::NeighborhoodIterator iterator(radius, mask, mask->GetRequestedRegion()); itk::NeighborhoodIterator valueIter(radius, itkValueImage, itkValueImage->GetRequestedRegion()); std::vector borderPoints; // // Calculate surface in different directions // double surface = 0; std::vector directionSurface; for (int i = 0; i < (int)(iterator.Size()); ++i) { auto offset = iterator.GetOffset(i); double deltaS = 1; int nonZeros = 0; - for (int j = 0; j < VImageDimension; ++j) + for (unsigned int j = 0; j < VImageDimension; ++j) { if (offset[j] != 0 && nonZeros == 0) { - for (int k = 0; k < VImageDimension; ++k) + for (unsigned int k = 0; k < VImageDimension; ++k) { if (k != j) deltaS *= mask->GetSpacing()[k]; } nonZeros++; } else if (offset[j] != 0) { deltaS = 0; } } if (nonZeros < 1) deltaS = 0; directionSurface.push_back(deltaS); } // // Prepare calulation of Centre of mass shift // PointType normalCenter(0); PointType normalCenterUncorrected(0); PointType weightedCenter(0); PointType currentPoint; int numberOfPoints = 0; int numberOfPointsUncorrected = 0; double sumOfPoints = 0; while(!iterator.IsAtEnd()) { if (iterator.GetCenterPixel() == 0) { ++iterator; ++valueIter; continue; } mask->TransformIndexToPhysicalPoint(iterator.GetIndex(), currentPoint); normalCenterUncorrected += currentPoint.GetVectorFromOrigin(); ++numberOfPointsUncorrected; double intensityValue = valueIter.GetCenterPixel(); if (intensityValue == intensityValue) { normalCenter += currentPoint.GetVectorFromOrigin(); weightedCenter += currentPoint.GetVectorFromOrigin() * intensityValue; sumOfPoints += intensityValue; ++numberOfPoints; } bool border = false; for (int i = 0; i < (int)(iterator.Size()); ++i) { if (iterator.GetPixel(i) == 0 || ( ! iterator.IndexInBounds(i))) { border = true; surface += directionSurface[i]; //break; } } if (border) { auto centerIndex = iterator.GetIndex(); PointType centerPoint; mask->TransformIndexToPhysicalPoint(centerIndex, centerPoint ); borderPoints.push_back(centerPoint); } ++iterator; ++valueIter; } auto normalCenterVector = normalCenter.GetVectorFromOrigin() / numberOfPoints; auto normalCenterVectorUncorrected = normalCenter.GetVectorFromOrigin() / numberOfPointsUncorrected; auto weightedCenterVector = weightedCenter.GetVectorFromOrigin() / sumOfPoints; auto differenceOfCentersUncorrected = (normalCenterVectorUncorrected - weightedCenterVector).GetNorm(); auto differenceOfCenters = (normalCenterVector - weightedCenterVector).GetNorm(); double longestDiameter = 0; unsigned long numberOfBorderPoints = borderPoints.size(); for (int i = 0; i < (int)numberOfBorderPoints; ++i) { auto point = borderPoints[i]; for (int j = i; j < (int)numberOfBorderPoints; ++j) { double newDiameter=point.EuclideanDistanceTo(borderPoints[j]); if (newDiameter > longestDiameter) longestDiameter = newDiameter; } } featureList.push_back(std::make_pair("Volumetric Features Maximum 3D diameter", longestDiameter)); featureList.push_back(std::make_pair("Volumetric Features Surface (Voxel based)", surface)); featureList.push_back(std::make_pair("Volumetric Features Centre of mass shift", differenceOfCenters)); featureList.push_back(std::make_pair("Volumetric Features Centre of mass shift (Uncorrected)", differenceOfCentersUncorrected)); } mitk::GIFVolumetricStatistics::GIFVolumetricStatistics() { SetLongName("volume"); SetShortName("vol"); } mitk::GIFVolumetricStatistics::FeatureListType mitk::GIFVolumetricStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; if (image->GetDimension() < 3) { return featureList; } AccessByItk_2(image, CalculateVolumeStatistic, mask, featureList); AccessByItk_2(mask, CalculateLargestDiameter, image, featureList); vtkSmartPointer mesher = vtkSmartPointer::New(); vtkSmartPointer stats = vtkSmartPointer::New(); mesher->SetInputData(mask->GetVtkImageData()); stats->SetInputConnection(mesher->GetOutputPort()); stats->Update(); double pi = vnl_math::pi; double meshVolume = stats->GetVolume(); double meshSurf = stats->GetSurfaceArea(); double pixelVolume = featureList[0].second; double pixelSurface = featureList[3].second; MITK_INFO << "Surf: " << pixelSurface << " Vol " << pixelVolume; double compactness1 = pixelVolume / (std::sqrt(pi) * std::pow(meshSurf, 2.0 / 3.0)); double compactness1Pixel = pixelVolume / (std::sqrt(pi) * std::pow(pixelSurface, 2.0 / 3.0)); //This is the definition used by Aertz. However, due to 2/3 this feature is not demensionless. Use compactness3 instead. double compactness2 = 36 * pi*pixelVolume*pixelVolume / meshSurf / meshSurf / meshSurf; double compactness2Pixel = 36 * pi*pixelVolume*pixelVolume / pixelSurface / pixelSurface / pixelSurface; double compactness3 = pixelVolume / (std::sqrt(pi) * std::pow(meshSurf, 3.0 / 2.0)); double compactness3Pixel = pixelVolume / (std::sqrt(pi) * std::pow(pixelSurface, 3.0 / 2.0)); double sphericity = std::pow(pi, 1 / 3.0) *std::pow(6 * pixelVolume, 2.0 / 3.0) / meshSurf; double sphericityPixel = std::pow(pi, 1 / 3.0) *std::pow(6 * pixelVolume, 2.0 / 3.0) / pixelSurface; double surfaceToVolume = meshSurf / pixelVolume; double surfaceToVolumePixel = pixelSurface / pixelVolume; double sphericalDisproportion = meshSurf / 4 / pi / std::pow(3.0 / 4.0 / pi * pixelVolume, 2.0 / 3.0); double sphericalDisproportionPixel = pixelSurface / 4 / pi / std::pow(3.0 / 4.0 / pi * pixelVolume, 2.0 / 3.0); double asphericity = std::pow(1.0/compactness2, (1.0 / 3.0)) - 1; double asphericityPixel = std::pow(1.0/compactness2Pixel, (1.0 / 3.0)) - 1; //Calculate center of mass shift int xx = mask->GetDimensions()[0]; int yy = mask->GetDimensions()[1]; int zz = mask->GetDimensions()[2]; double xd = mask->GetGeometry()->GetSpacing()[0]; double yd = mask->GetGeometry()->GetSpacing()[1]; double zd = mask->GetGeometry()->GetSpacing()[2]; std::vector pointsForPCA; std::vector pointsForPCAUncorrected; for (int x = 0; x < xx; x++) { for (int y = 0; y < yy; y++) { for (int z = 0; z < zz; z++) { itk::Image::IndexType index; index[0] = x; index[1] = y; index[2] = z; mitk::ScalarType pxImage; mitk::ScalarType pxMask; mitkPixelTypeMultiplex5( mitk::FastSinglePixelAccess, image->GetChannelDescriptor().GetPixelType(), image, image->GetVolumeData(), index, pxImage, 0); mitkPixelTypeMultiplex5( mitk::FastSinglePixelAccess, mask->GetChannelDescriptor().GetPixelType(), mask, mask->GetVolumeData(), index, pxMask, 0); //Check if voxel is contained in segmentation if (pxMask > 0) { cv::Point3d tmp; tmp.x = x * xd; tmp.y = y * yd; tmp.z = z * zd; pointsForPCAUncorrected.push_back(tmp); if (pxImage == pxImage) { pointsForPCA.push_back(tmp); } } } } } //Calculate PCA Features int sz = pointsForPCA.size(); cv::Mat data_pts = cv::Mat(sz, 3, CV_64FC1); for (int i = 0; i < data_pts.rows; ++i) { data_pts.at(i, 0) = pointsForPCA[i].x; data_pts.at(i, 1) = pointsForPCA[i].y; data_pts.at(i, 2) = pointsForPCA[i].z; } //Calculate PCA Features int szUC = pointsForPCAUncorrected.size(); cv::Mat data_ptsUC = cv::Mat(szUC, 3, CV_64FC1); for (int i = 0; i < data_ptsUC.rows; ++i) { data_ptsUC.at(i, 0) = pointsForPCAUncorrected[i].x; data_ptsUC.at(i, 1) = pointsForPCAUncorrected[i].y; data_ptsUC.at(i, 2) = pointsForPCAUncorrected[i].z; } //Perform PCA analysis cv::PCA pca_analysis(data_pts, cv::Mat(), CV_PCA_DATA_AS_ROW); cv::PCA pca_analysisUC(data_ptsUC, cv::Mat(), CV_PCA_DATA_AS_ROW); //Store the eigenvalues std::vector eigen_val(3); std::vector eigen_valUC(3); for (int i = 0; i < 3; ++i) { eigen_val[i] = pca_analysis.eigenvalues.at(0, i); eigen_valUC[i] = pca_analysisUC.eigenvalues.at(0, i); } std::sort(eigen_val.begin(), eigen_val.end()); std::sort(eigen_valUC.begin(), eigen_valUC.end()); double major = eigen_val[2]; double minor = eigen_val[1]; double least = eigen_val[0]; double elongation = major == 0 ? 0 : minor/major; double flatness = major == 0 ? 0 : least / major; double majorUC = eigen_valUC[2]; double minorUC = eigen_valUC[1]; double leastUC = eigen_valUC[0]; double elongationUC = majorUC == 0 ? 0 : minorUC / majorUC; double flatnessUC = majorUC == 0 ? 0 : leastUC / majorUC; featureList.push_back(std::make_pair("Volumetric Features Volume (mesh based)",meshVolume)); featureList.push_back(std::make_pair("Volumetric Features Surface area",meshSurf)); featureList.push_back(std::make_pair("Volumetric Features Surface to volume ratio",surfaceToVolume)); featureList.push_back(std::make_pair("Volumetric Features Sphericity",sphericity)); featureList.push_back(std::make_pair("Volumetric Features Asphericity",asphericity)); featureList.push_back(std::make_pair("Volumetric Features Compactness 1",compactness1)); featureList.push_back(std::make_pair("Volumetric Features Compactness 2",compactness2)); featureList.push_back(std::make_pair("Volumetric Features Compactness 3",compactness3)); featureList.push_back(std::make_pair("Volumetric Features Spherical disproportion", sphericalDisproportion)); featureList.push_back(std::make_pair("Volumetric Features Surface to volume ratio (Voxel based)", surfaceToVolumePixel)); featureList.push_back(std::make_pair("Volumetric Features Sphericity (Voxel based)", sphericityPixel)); featureList.push_back(std::make_pair("Volumetric Features Asphericity (Voxel based)", asphericityPixel)); featureList.push_back(std::make_pair("Volumetric Features Compactness 1 (Voxel based)", compactness1Pixel)); featureList.push_back(std::make_pair("Volumetric Features Compactness 2 (Voxel based)", compactness2Pixel)); featureList.push_back(std::make_pair("Volumetric Features Compactness 3 (Voxel based)", compactness3Pixel)); featureList.push_back(std::make_pair("Volumetric Features Spherical disproportion (Voxel based)", sphericalDisproportionPixel)); featureList.push_back(std::make_pair("Volumetric Features PCA Major Axis",major)); featureList.push_back(std::make_pair("Volumetric Features PCA Minor Axis",minor)); featureList.push_back(std::make_pair("Volumetric Features PCA Least Axis",least)); featureList.push_back(std::make_pair("Volumetric Features PCA Elongation",elongation)); featureList.push_back(std::make_pair("Volumetric Features PCA Flatness",flatness)); featureList.push_back(std::make_pair("Volumetric Features PCA Major Axis (Uncorrected)", majorUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Minor Axis (Uncorrected)", minorUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Least Axis (Uncorrected)", leastUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Elongation (Uncorrected)", elongationUC)); featureList.push_back(std::make_pair("Volumetric Features PCA Flatness (Uncorrected)", flatnessUC)); return featureList; } mitk::GIFVolumetricStatistics::FeatureNameListType mitk::GIFVolumetricStatistics::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("Volumetric Features Compactness 1"); featureList.push_back("Volumetric Features Compactness 2"); featureList.push_back("Volumetric Features Compactness 3"); featureList.push_back("Volumetric Features Sphericity"); featureList.push_back("Volumetric Features Asphericity"); featureList.push_back("Volumetric Features Surface to volume ratio"); featureList.push_back("Volumetric Features Surface area"); featureList.push_back("Volumetric Features Volume (mesh based)"); featureList.push_back("Volumetric Features Volume (pixel based)"); featureList.push_back("Volumetric Features Spherical disproportion"); featureList.push_back("Volumetric Features Maximum 3D diameter"); return featureList; } void mitk::GIFVolumetricStatistics::AddArguments(mitkCommandLineParser &parser) { std::string name = GetOptionPrefix(); parser.addArgument(GetLongName(), name, mitkCommandLineParser::String, "Use Volume-Statistic", "calculates volume based features", us::Any()); } void mitk::GIFVolumetricStatistics::CalculateFeaturesUsingParameters(const Image::Pointer & feature, const Image::Pointer &mask, const Image::Pointer &, FeatureListType &featureList) { auto parsedArgs = GetParameter(); if (parsedArgs.count(GetLongName())) { MITK_INFO << "Start calculating volumetric features ...."; auto localResults = this->CalculateFeatures(feature, mask); featureList.insert(featureList.end(), localResults.begin(), localResults.end()); MITK_INFO << "Finished calculating volumetric features...."; } }