diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.h b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.h index 0d68bb9fbd..56586f46ae 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.h +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.h @@ -1,274 +1,278 @@ /*========================================================================= * * 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_h #define __itkEnhancedHistogramToTextureFeaturesFilter_h #include "itkHistogram.h" #include "itkMacro.h" #include "itkProcessObject.h" #include "itkSimpleDataObjectDecorator.h" /** Get built-in type. Creates member Get"name"() (e.g., GetVisibility()); */ #define itkMacroGLCMFeatureGetter(name) \ const MeasurementObjectType * Get##name##Output() const; \ \ MeasurementType Get##name() const; namespace itk { namespace Statistics { /** \class EnhancedHistogramToTextureFeaturesFilter * \brief This class computes texture feature coefficients from a grey level * co-occurrence matrix. * * This class computes features that summarize image texture, given a grey level * co-occurrence matrix (generated by a ScalarImageToCooccurrenceMatrixFilter * or related class). * * The features calculated are as follows (where \f$ g(i, j) \f$ is the element in * cell i, j of a a normalized GLCM): * * "Energy" \f$ = f_1 = \sum_{i,j}g(i, j)^2 \f$ * * "Entropy" \f$ = f_2 = -\sum_{i,j}g(i, j) \log_2 g(i, j)\f$, or 0 if \f$g(i, j) = 0\f$ * * "Correlation" \f$ = f_3 = \sum_{i,j}\frac{(i - \mu)(j - \mu)g(i, j)}{\sigma^2} \f$ * * "Difference Moment" \f$= f_4 = \sum_{i,j}\frac{1}{1 + (i - j)^2}g(i, j) \f$ * * "Inertia" \f$ = f_5 = \sum_{i,j}(i - j)^2g(i, j) \f$ (sometimes called "contrast.") * * "Cluster Shade" \f$ = f_6 = \sum_{i,j}((i - \mu) + (j - \mu))^3 g(i, j) \f$ * * "Cluster Prominence" \f$ = f_7 = \sum_{i,j}((i - \mu) + (j - \mu))^4 g(i, j) \f$ * * "Haralick's Correlation" \f$ = f_8 = \frac{\sum_{i,j}(i, j) g(i, j) -\mu_t^2}{\sigma_t^2} \f$ * where \f$\mu_t\f$ and \f$\sigma_t\f$ are the mean and standard deviation of the row * (or column, due to symmetry) sums. * * Above, \f$ \mu = \f$ (weighted pixel average) \f$ = \sum_{i,j}i \cdot g(i, j) = * \sum_{i,j}j \cdot g(i, j) \f$ (due to matrix summetry), and * * \f$ \sigma = \f$ (weighted pixel variance) \f$ = \sum_{i,j}(i - \mu)^2 \cdot g(i, j) = * \sum_{i,j}(j - \mu)^2 \cdot g(i, j) \f$ (due to matrix summetry) * * A good texture feature set to use is the Conners, Trivedi and Harlow set: * features 1, 2, 4, 5, 6, and 7. There is some correlation between the various * features, so using all of them at the same time is not necessarialy a good idea. * * NOTA BENE: The input histogram will be forcably normalized! * This algorithm takes three passes through the input * histogram if the histogram was already normalized, and four if not. * * Web references: * * http://www.cssip.uq.edu.au/meastex/www/algs/algs/algs.html * http://www.ucalgary.ca/~mhallbey/texture/texture_tutorial.html * * Print references: * * Haralick, R.M., K. Shanmugam and I. Dinstein. 1973. Textural Features for * Image Classification. IEEE Transactions on Systems, Man and Cybernetics. * SMC-3(6):610-620. * * Haralick, R.M. 1979. Statistical and Structural Approaches to Texture. * Proceedings of the IEEE, 67:786-804. * * R.W. Conners and C.A. Harlow. A Theoretical Comaprison of Texture Algorithms. * IEEE Transactions on Pattern Analysis and Machine Intelligence, 2:204-222, 1980. * * R.W. Conners, M.M. Trivedi, and C.A. Harlow. Segmentation of a High-Resolution * Urban Scene using Texture Operators. Computer Vision, Graphics and Image * Processing, 25:273-310, 1984. * * \sa ScalarImageToCooccurrenceMatrixFilter * \sa ScalarImageToTextureFeaturesFilter * * Author: Zachary Pincus * \ingroup ITKStatistics */ template< typename THistogram > class EnhancedHistogramToTextureFeaturesFilter:public ProcessObject { public: /** Standard typedefs */ typedef EnhancedHistogramToTextureFeaturesFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; /** Run-time type information (and related methods). */ itkTypeMacro(EnhancedHistogramToTextureFeaturesFilter, ProcessObject); /** standard New() method support */ itkNewMacro(Self); typedef THistogram HistogramType; typedef typename HistogramType::Pointer HistogramPointer; typedef typename HistogramType::ConstPointer HistogramConstPointer; typedef typename HistogramType::MeasurementType MeasurementType; typedef typename HistogramType::MeasurementVectorType MeasurementVectorType; typedef typename HistogramType::IndexType IndexType; typedef typename HistogramType::AbsoluteFrequencyType AbsoluteFrequencyType; typedef typename HistogramType::RelativeFrequencyType RelativeFrequencyType; typedef typename HistogramType::TotalAbsoluteFrequencyType TotalAbsoluteFrequencyType; typedef typename HistogramType::TotalRelativeFrequencyType TotalRelativeFrequencyType; /** Container to hold relative frequencies of the histogram */ typedef std::vector< RelativeFrequencyType > RelativeFrequencyContainerType; /** Method to Set/Get the input Histogram */ using Superclass::SetInput; void SetInput(const HistogramType *histogram); const HistogramType * GetInput() const; /** Smart Pointer type to a DataObject. */ typedef DataObject::Pointer DataObjectPointer; /** Type of DataObjects used for scalar outputs */ typedef SimpleDataObjectDecorator< MeasurementType > MeasurementObjectType; /** Return energy texture value. */ MeasurementType GetEnergy() const; const MeasurementObjectType * GetEnergyOutput() const; /** Return entropy texture value. */ MeasurementType GetEntropy() const; const MeasurementObjectType * GetEntropyOutput() const; /** return correlation texture value. */ MeasurementType GetCorrelation() const; const MeasurementObjectType * GetCorrelationOutput() const; /** Return inverse difference moment texture value. */ MeasurementType GetInverseDifferenceMoment() const; const MeasurementObjectType * GetInverseDifferenceMomentOutput() const; /** Return inertia texture value. */ MeasurementType GetInertia() const; const MeasurementObjectType * GetInertiaOutput() const; /** Return cluster shade texture value. */ MeasurementType GetClusterShade() const; const MeasurementObjectType * GetClusterShadeOutput() const; /** Return cluster prominence texture value. */ MeasurementType GetClusterProminence() const; const MeasurementObjectType * GetClusterProminenceOutput() const; /** Return Haralick correlation texture value. */ MeasurementType GetHaralickCorrelation() const; const MeasurementObjectType * GetHaralickCorrelationOutput() const; itkMacroGLCMFeatureGetter(Autocorrelation); itkMacroGLCMFeatureGetter(Contrast); itkMacroGLCMFeatureGetter(Dissimilarity); itkMacroGLCMFeatureGetter(MaximumProbability); itkMacroGLCMFeatureGetter(InverseVariance); itkMacroGLCMFeatureGetter(Homogeneity1); itkMacroGLCMFeatureGetter(ClusterTendency); itkMacroGLCMFeatureGetter(Variance); itkMacroGLCMFeatureGetter(SumAverage); itkMacroGLCMFeatureGetter(SumEntropy); itkMacroGLCMFeatureGetter(SumVariance); itkMacroGLCMFeatureGetter(DifferenceAverage); itkMacroGLCMFeatureGetter(DifferenceEntropy); itkMacroGLCMFeatureGetter(DifferenceVariance); itkMacroGLCMFeatureGetter(InverseDifferenceMomentNormalized); itkMacroGLCMFeatureGetter(InverseDifferenceNormalized); itkMacroGLCMFeatureGetter(InverseDifference); itkMacroGLCMFeatureGetter(JointAverage); + itkMacroGLCMFeatureGetter(FirstMeasureOfInformationCorrelation); + itkMacroGLCMFeatureGetter(SecondMeasureOfInformationCorrelation); /** Texture feature types */ typedef enum { Energy, Entropy, Correlation, InverseDifferenceMoment, Inertia, ClusterShade, ClusterProminence, HaralickCorrelation, Autocorrelation, Contrast, Dissimilarity, MaximumProbability, InverseVariance, Homogeneity1, ClusterTendency, Variance, SumAverage, SumEntropy, SumVariance, DifferenceAverage, DifferenceEntropy, DifferenceVariance, InverseDifferenceMomentNormalized, InverseDifferenceNormalized, InverseDifference, JointAverage, + FirstMeasureOfInformationCorrelation, + SecondMeasureOfInformationCorrelation, InvalidFeatureName } TextureFeatureName; /** convenience method to access the texture values */ MeasurementType GetFeature(TextureFeatureName name); protected: EnhancedHistogramToTextureFeaturesFilter(); ~EnhancedHistogramToTextureFeaturesFilter() {} virtual void PrintSelf(std::ostream & os, Indent indent) const ITK_OVERRIDE; /** Make a DataObject to be used for output output. */ typedef ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType; using Superclass::MakeOutput; virtual DataObjectPointer MakeOutput(DataObjectPointerArraySizeType) ITK_OVERRIDE; virtual void GenerateData() ITK_OVERRIDE; private: EnhancedHistogramToTextureFeaturesFilter(const Self &); //purposely not implemented void operator=(const Self &); //purposely not implemented void ComputeMeansAndVariances(double & pixelMean, double & marginalMean, double & marginalDevSquared, double & pixelVariance); RelativeFrequencyContainerType m_RelativeFrequencyContainer; }; } // end of namespace Statistics } // end of namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkEnhancedHistogramToTextureFeaturesFilter.hxx" #endif #endif diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx index de729c8b03..fa37529483 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx @@ -1,665 +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 +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 ) { - //constructor - template< typename THistogram > - EnhancedHistogramToTextureFeaturesFilter< THistogram >::EnhancedHistogramToTextureFeaturesFilter(void) + 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; + + for ( HistogramIterator hit = inputHistogram->Begin(); + hit != inputHistogram->End(); ++hit ) + { + RelativeFrequencyType frequency = *rFreqIterator; + ++rFreqIterator; + if ( frequency == 0 ) { - this->ProcessObject::SetNumberOfRequiredInputs(1); - - // allocate the data objects for the outputs which are - // just decorators real types - for ( int i = 0; i < 26; ++i ) - { - this->ProcessObject::SetNthOutput( i, this->MakeOutput(i) ); - } + continue; // no use doing these calculations if we're just multiplying by + // zero. } - template< typename THistogram > - void - EnhancedHistogramToTextureFeaturesFilter< THistogram > - ::SetInput(const HistogramType *histogram) + 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) { - this->ProcessObject::SetNthInput( 0, const_cast< HistogramType * >( histogram ) ); + 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 ); } - - template< typename THistogram > - const typename - EnhancedHistogramToTextureFeaturesFilter< THistogram >::HistogramType * - EnhancedHistogramToTextureFeaturesFilter< THistogram > - ::GetInput() const + else { - return itkDynamicCastInDebugMode< const HistogramType * >( this->GetPrimaryInput() ); + inverseDifferenceMomentNormalized = 0; + inverseDifferenceNormalized = 0; } + inverseDifference += frequency / ( 1.0 + std::abs( index[0] - index[1] ) ); + } - template< typename THistogram > - typename - EnhancedHistogramToTextureFeaturesFilter< THistogram >::DataObjectPointer - EnhancedHistogramToTextureFeaturesFilter< THistogram > - ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed(idx) ) + 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++ ) { - return MeasurementObjectType::New().GetPointer(); + globalIndex[0] = i; + globalIndex[1] = j; + pi_tmp += inputHistogram->GetFrequency(globalIndex); } + pi_margins.push_back(pi_tmp); + } - template< typename THistogram > - void - EnhancedHistogramToTextureFeaturesFilter< THistogram >::GenerateData(void) + //pj. + for ( int j = 1; j <= inputHistogram->GetSize(1); j++ ) + { + double pj_tmp= 0.0; + for( int i = 1; i <= inputHistogram->GetSize(0); i++ ) { - typedef typename HistogramType::ConstIterator HistogramIterator; + globalIndex[0] = i; + globalIndex[1] = j; + pj_tmp += inputHistogram->GetFrequency(globalIndex); + } + pj_margins.push_back(pj_tmp); + } - const HistogramType *inputHistogram = this->GetInput(); + //Calculate HX + double hx = 0.0; - //Normalize the absolute frequencies and populate the relative frequency - //container - TotalRelativeFrequencyType totalFrequency = - static_cast< TotalRelativeFrequencyType >( inputHistogram->GetTotalFrequency() ); + 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; + } - m_RelativeFrequencyContainer.clear(); + //Calculate HXY1 + double hxy1 = 0.0; - typename HistogramType::SizeValueType binsPerAxis = inputHistogram->GetSize(0); - std::vector sumP; - std::vector diffP; + for ( int i = 0; i < inputHistogram->GetSize(0); i++ ) + { + for ( int j = 0; j < inputHistogram->GetSize(1); j++ ) + { + globalIndex[0] = i; + globalIndex[1] = j; - sumP.resize(2*binsPerAxis,0.0); - diffP.resize(binsPerAxis,0.0); + double pi_margin = pi_margins[i]; + double pj_margin = pj_margins[j]; - double numberOfPixels = 0; + double p_ij = inputHistogram->GetFrequency(globalIndex); - 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(); - - - 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(); - - 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() ); - 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; - } - - 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); + hxy1 -= ( (pi_margin * pj_margin) > 0.0001 ) ? p_ij *std::log(pi_margin * pj_margin) / log2:0; } + } + + //Calculate HXY2 + double hxy2 = 0.0; - template< typename THistogram > - void - EnhancedHistogramToTextureFeaturesFilter< THistogram >::ComputeMeansAndVariances(double & pixelMean, - double & marginalMean, - double & marginalDevSquared, - double & pixelVariance) + for ( int i = 0; i < inputHistogram->GetSize(0); i++ ) + { + for ( int j = 0; j < inputHistogram->GetSize(1); j++ ) { - // 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. + 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; + typedef typename HistogramType::ConstIterator HistogramIterator; - const HistogramType *inputHistogram = this->GetInput(); + const HistogramType *inputHistogram = this->GetInput(); - // Initialize everything - typename HistogramType::SizeValueType binsPerAxis = inputHistogram->GetSize(0); - double *marginalSums = new double[binsPerAxis]; - for ( double *ms_It = marginalSums; + // 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. + { + *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; - } + 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; + } - template< typename THistogram > - const - typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * - EnhancedHistogramToTextureFeaturesFilter< THistogram > - ::GetEnergyOutput() const - { - return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(0) ); - } + delete[] marginalSums; +} - 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 > +::GetEnergyOutput() const +{ + return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(0) ); +} - 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 > +::GetEntropyOutput() const +{ + return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(1) ); +} - 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 > +::GetCorrelationOutput() const +{ + return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(2) ); +} - 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 > +::GetInverseDifferenceMomentOutput() const +{ + return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(3) ); +} - 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 > +::GetInertiaOutput() const +{ + return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(4) ); +} - 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 > +::GetClusterShadeOutput() const +{ + return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(5) ); +} - template< typename THistogram > - const - typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * - EnhancedHistogramToTextureFeaturesFilter< THistogram > - ::GetHaralickCorrelationOutput() const - { - return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(7) ); - } +template< typename THistogram > +const +typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * +EnhancedHistogramToTextureFeaturesFilter< THistogram > +::GetClusterProminenceOutput() const +{ + return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(6) ); +} - 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) - - template< typename THistogram > - typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType - EnhancedHistogramToTextureFeaturesFilter< THistogram > - ::GetEnergy() const - { - return this->GetEnergyOutput()->Get(); - } +template< typename THistogram > +const +typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * +EnhancedHistogramToTextureFeaturesFilter< THistogram > +::GetHaralickCorrelationOutput() const +{ + return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(7) ); +} - template< typename THistogram > - typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType - EnhancedHistogramToTextureFeaturesFilter< THistogram > - ::GetEntropy() const - { - return this->GetEntropyOutput()->Get(); - } +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 > - ::GetCorrelation() const - { - return this->GetCorrelationOutput()->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 > - ::GetInverseDifferenceMoment() const - { - return this->GetInverseDifferenceMomentOutput()->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 > - ::GetInertia() const - { - return this->GetInertiaOutput()->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 > - ::GetClusterShade() const - { - return this->GetClusterShadeOutput()->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 > - ::GetClusterProminence() const - { - return this->GetClusterProminenceOutput()->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 > - ::GetHaralickCorrelation() const - { - return this->GetHaralickCorrelationOutput()->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); - default: - return 0; - } - } + 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 +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/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp index e2839b1336..f086dbc940 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp @@ -1,296 +1,312 @@ #include // MITK #include #include #include // ITK #include #include // STL #include template void CalculateCoocurenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFCooccurenceMatrix::FeatureListType & featureList, mitk::GIFCooccurenceMatrix::GIFCooccurenceMatrixConfiguration config) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::Statistics::EnhancedScalarImageToTextureFeaturesFilter FilterType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef typename FilterType::TextureFeaturesFilterType TextureFilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename FilterType::Pointer filter = FilterType::New(); typename FilterType::OffsetVector::Pointer newOffset = FilterType::OffsetVector::New(); auto oldOffsets = filter->GetOffsets(); auto oldOffsetsIterator = oldOffsets->Begin(); while(oldOffsetsIterator != oldOffsets->End()) { bool continueOuterLoop = false; typename FilterType::OffsetType offset = oldOffsetsIterator->Value(); for (unsigned int i = 0; i < VImageDimension; ++i) { offset[i] *= config.range; if (config.direction == i + 2 && offset[i] != 0) { continueOuterLoop = true; } } if (config.direction == 1) { offset[0] = 0; offset[1] = 0; offset[2] = 1; newOffset->push_back(offset); break; } oldOffsetsIterator++; if (continueOuterLoop) continue; newOffset->push_back(offset); } filter->SetOffsets(newOffset); // All features are required typename FilterType::FeatureNameVectorPointer requestedFeatures = FilterType::FeatureNameVector::New(); requestedFeatures->push_back(TextureFilterType::Energy); requestedFeatures->push_back(TextureFilterType::Entropy); requestedFeatures->push_back(TextureFilterType::Correlation); requestedFeatures->push_back(TextureFilterType::InverseDifferenceMoment); requestedFeatures->push_back(TextureFilterType::Inertia); requestedFeatures->push_back(TextureFilterType::ClusterShade); requestedFeatures->push_back(TextureFilterType::ClusterProminence); requestedFeatures->push_back(TextureFilterType::HaralickCorrelation); requestedFeatures->push_back(TextureFilterType::Autocorrelation); requestedFeatures->push_back(TextureFilterType::Contrast); requestedFeatures->push_back(TextureFilterType::Dissimilarity); requestedFeatures->push_back(TextureFilterType::MaximumProbability); requestedFeatures->push_back(TextureFilterType::InverseVariance); requestedFeatures->push_back(TextureFilterType::Homogeneity1); requestedFeatures->push_back(TextureFilterType::ClusterTendency); requestedFeatures->push_back(TextureFilterType::Variance); requestedFeatures->push_back(TextureFilterType::SumAverage); requestedFeatures->push_back(TextureFilterType::SumEntropy); requestedFeatures->push_back(TextureFilterType::SumVariance); requestedFeatures->push_back(TextureFilterType::DifferenceAverage); requestedFeatures->push_back(TextureFilterType::DifferenceEntropy); requestedFeatures->push_back(TextureFilterType::DifferenceVariance); requestedFeatures->push_back(TextureFilterType::InverseDifferenceMomentNormalized); requestedFeatures->push_back(TextureFilterType::InverseDifferenceNormalized); requestedFeatures->push_back(TextureFilterType::InverseDifference); requestedFeatures->push_back(TextureFilterType::JointAverage); + requestedFeatures->push_back(TextureFilterType::FirstMeasureOfInformationCorrelation); + requestedFeatures->push_back(TextureFilterType::SecondMeasureOfInformationCorrelation); typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); filter->SetInput(itkImage); filter->SetMaskImage(maskImage); filter->SetRequestedFeatures(requestedFeatures); filter->SetPixelValueMinMax(minMaxComputer->GetMinimum()-0.5,minMaxComputer->GetMaximum()+0.5); //filter->SetPixelValueMinMax(-1024,3096); //filter->SetNumberOfBinsPerAxis(5); filter->Update(); auto featureMeans = filter->GetFeatureMeans (); auto featureStd = filter->GetFeatureStandardDeviations(); std::ostringstream ss; ss << config.range; std::string strRange = ss.str(); for (std::size_t i = 0; i < featureMeans->size(); ++i) { switch (i) { case TextureFilterType::Energy : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Energy Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Energy Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Entropy : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Entropy Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Entropy Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Correlation : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Correlation Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Correlation Std.",featureStd->ElementAt(i))); break; case TextureFilterType::InverseDifferenceMoment : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMoment Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMoment Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Inertia : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Inertia Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Inertia Std.",featureStd->ElementAt(i))); break; case TextureFilterType::ClusterShade : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterShade Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterShade Std.",featureStd->ElementAt(i))); break; case TextureFilterType::ClusterProminence : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterProminence Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterProminence Std.",featureStd->ElementAt(i))); break; case TextureFilterType::HaralickCorrelation : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") HaralickCorrelation Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") HaralickCorrelation Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Autocorrelation : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Autocorrelation Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Autocorrelation Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Contrast : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Contrast Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Contrast Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Dissimilarity : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Dissimilarity Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Dissimilarity Std.",featureStd->ElementAt(i))); break; case TextureFilterType::MaximumProbability : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") MaximumProbability Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") MaximumProbability Std.",featureStd->ElementAt(i))); break; case TextureFilterType::InverseVariance : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseVariance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseVariance Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Homogeneity1 : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Homogeneity1 Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Homogeneity1 Std.",featureStd->ElementAt(i))); break; case TextureFilterType::ClusterTendency : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterTendency Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") ClusterTendency Std.",featureStd->ElementAt(i))); break; case TextureFilterType::Variance : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Variance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") Variance Std.",featureStd->ElementAt(i))); break; case TextureFilterType::SumAverage : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumAverage Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumAverage Std.",featureStd->ElementAt(i))); break; case TextureFilterType::SumEntropy : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumEntropy Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumEntropy Std.",featureStd->ElementAt(i))); break; case TextureFilterType::SumVariance : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumVariance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SumVariance Std.",featureStd->ElementAt(i))); break; case TextureFilterType::DifferenceAverage : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceAverage Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceAverage Std.",featureStd->ElementAt(i))); break; case TextureFilterType::DifferenceEntropy : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceEntropy Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceEntropy Std.",featureStd->ElementAt(i))); break; case TextureFilterType::DifferenceVariance : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceVariance Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") DifferenceVariance Std.",featureStd->ElementAt(i))); break; case TextureFilterType::InverseDifferenceMomentNormalized : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMomentNormalized Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceMomentNormalized Std.",featureStd->ElementAt(i))); break; case TextureFilterType::InverseDifferenceNormalized : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceNormalized Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifferenceNormalized Std.",featureStd->ElementAt(i))); break; case TextureFilterType::InverseDifference : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifference Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") InverseDifference Std.",featureStd->ElementAt(i))); break; case TextureFilterType::JointAverage : featureList.push_back(std::make_pair("co-occ. ("+ strRange+") JointAverage Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("co-occ. ("+ strRange+") JointAverage Std.",featureStd->ElementAt(i))); break; + case TextureFilterType::FirstMeasureOfInformationCorrelation : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") FirstMeasureOfInformationCorrelation Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") FirstMeasureOfInformationCorrelation Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::SecondMeasureOfInformationCorrelation : + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SecondMeasureOfInformationCorrelation Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("co-occ. ("+ strRange+") SecondMeasureOfInformationCorrelation Std.",featureStd->ElementAt(i))); + break; default: break; } } } mitk::GIFCooccurenceMatrix::GIFCooccurenceMatrix(): m_Range(1.0), m_Direction(0) { } mitk::GIFCooccurenceMatrix::FeatureListType mitk::GIFCooccurenceMatrix::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; GIFCooccurenceMatrixConfiguration config; config.direction = m_Direction; config.range = m_Range; AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config); return featureList; } mitk::GIFCooccurenceMatrix::FeatureNameListType mitk::GIFCooccurenceMatrix::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("co-occ. Energy Means"); featureList.push_back("co-occ. Energy Std."); featureList.push_back("co-occ. Entropy Means"); featureList.push_back("co-occ. Entropy Std."); featureList.push_back("co-occ. Correlation Means"); featureList.push_back("co-occ. Correlation Std."); featureList.push_back("co-occ. InverseDifferenceMoment Means"); featureList.push_back("co-occ. InverseDifferenceMoment Std."); featureList.push_back("co-occ. Inertia Means"); featureList.push_back("co-occ. Inertia Std."); featureList.push_back("co-occ. ClusterShade Means"); featureList.push_back("co-occ. ClusterShade Std."); featureList.push_back("co-occ. ClusterProminence Means"); featureList.push_back("co-occ. ClusterProminence Std."); featureList.push_back("co-occ. HaralickCorrelation Means"); featureList.push_back("co-occ. HaralickCorrelation Std."); featureList.push_back("co-occ. Autocorrelation Means"); featureList.push_back("co-occ. Autocorrelation Std."); featureList.push_back("co-occ. Contrast Means"); featureList.push_back("co-occ. Contrast Std."); featureList.push_back("co-occ. Dissimilarity Means"); featureList.push_back("co-occ. Dissimilarity Std."); featureList.push_back("co-occ. MaximumProbability Means"); featureList.push_back("co-occ. MaximumProbability Std."); featureList.push_back("co-occ. InverseVariance Means"); featureList.push_back("co-occ. InverseVariance Std."); featureList.push_back("co-occ. Homogeneity1 Means"); featureList.push_back("co-occ. Homogeneity1 Std."); featureList.push_back("co-occ. ClusterTendency Means"); featureList.push_back("co-occ. ClusterTendency Std."); featureList.push_back("co-occ. Variance Means"); featureList.push_back("co-occ. Variance Std."); featureList.push_back("co-occ. SumAverage Means"); featureList.push_back("co-occ. SumAverage Std."); featureList.push_back("co-occ. SumEntropy Means"); featureList.push_back("co-occ. SumEntropy Std."); featureList.push_back("co-occ. SumVariance Means"); featureList.push_back("co-occ. SumVariance Std."); featureList.push_back("co-occ. DifferenceAverage Means"); featureList.push_back("co-occ. DifferenceAverage Std."); featureList.push_back("co-occ. DifferenceEntropy Means"); featureList.push_back("co-occ. DifferenceEntropy Std."); featureList.push_back("co-occ. DifferenceVariance Means"); featureList.push_back("co-occ. DifferenceVariance Std."); featureList.push_back("co-occ. InverseDifferenceMomentNormalized Means"); featureList.push_back("co-occ. InverseDifferenceMomentNormalized Std."); featureList.push_back("co-occ. InverseDifferenceNormalized Means"); featureList.push_back("co-occ. InverseDifferenceNormalized Std."); featureList.push_back("co-occ. InverseDifference Means"); featureList.push_back("co-occ. InverseDifference Std."); + featureList.push_back("co-occ. JointAverage Means"); + featureList.push_back("co-occ. JointAverage Std."); + featureList.push_back("co-occ. FirstMeasurementOfInformationCorrelation Means"); + featureList.push_back("co-occ. FirstMeasurementOfInformationCorrelation Std."); + featureList.push_back("co-occ. SecondMeasurementOfInformationCorrelation Means"); + featureList.push_back("co-occ. SecondMeasurementOfInformationCorrelation Std."); return featureList; }