diff --git a/Modules/Classification/CLUtilities/CMakeLists.txt b/Modules/Classification/CLUtilities/CMakeLists.txt index 7e4e3a5f17..c3e7ffa662 100644 --- a/Modules/Classification/CLUtilities/CMakeLists.txt +++ b/Modules/Classification/CLUtilities/CMakeLists.txt @@ -1,7 +1,7 @@ MITK_CREATE_MODULE( DEPENDS MitkCore MitkCLCore - PACKAGE_DEPENDS PUBLIC Eigen + PACKAGE_DEPENDS PUBLIC Eigen OpenCV WARNINGS_AS_ERRORS ) add_subdirectory(test) diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.h b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.h index 692c47f1a8..44d67da15c 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.h +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.h @@ -1,216 +1,231 @@ /*=================================================================== 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 __itkEnhancedHistogramToRunLengthFeaturesFilter_h #define __itkEnhancedHistogramToRunLengthFeaturesFilter_h #include "itkHistogram.h" #include "itkMacro.h" #include "itkProcessObject.h" #include "itkSimpleDataObjectDecorator.h" namespace itk { namespace Statistics { /** \class EnhancedHistogramToRunLengthFeaturesFilter * \brief This class computes texture feature coefficients from a grey level * run-length matrix. * * By default, run length features are computed for each spatial * direction and then averaged afterward, so it is possible to access the * standard deviations of the texture features. These values give a clue as * to texture anisotropy. However, doing this is much more work, because it * involved computing one for each offset given. To compute a single matrix * using the first offset, call FastCalculationsOn(). If this is called, * then the texture standard deviations will not be computed (and will be set * to zero), but texture computation will be much faster. * * This class is templated over the input histogram type. * * Print references: * M. M. Galloway. Texture analysis using gray level run lengths. Computer * Graphics and Image Processing, 4:172-179, 1975. * * A. Chu, C. M. Sehgal, and J. F. Greenleaf. Use of gray value distribution of * run lengths for texture analysis. Pattern Recognition Letters, 11:415-420, * 1990. * * B. R. Dasarathy and E. B. Holder. Image characterizations based on joint * gray-level run-length distributions. Pattern Recognition Letters, 12:490-502, * 1991. * * IJ article: http://hdl.handle.net/1926/1374 * * \sa ScalarImageToRunLengthFeaturesFilter * \sa ScalarImageToRunLengthMatrixFilter * \sa EnhancedHistogramToRunLengthFeaturesFilter * * \author: Nick Tustison * \ingroup ITKStatistics */ template< typename THistogram > class EnhancedHistogramToRunLengthFeaturesFilter : public ProcessObject { public: /** Standard typedefs */ typedef EnhancedHistogramToRunLengthFeaturesFilter Self; typedef ProcessObject Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; /** Run-time type information (and related methods). */ itkTypeMacro( EnhancedHistogramToRunLengthFeaturesFilter, 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:: TotalAbsoluteFrequencyType FrequencyType; /** 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 MeasurementObjectType; /** Methods to return the short run emphasis. */ MeasurementType GetShortRunEmphasis() const; const MeasurementObjectType* GetShortRunEmphasisOutput() const; /** Methods to return the long run emphasis. */ MeasurementType GetLongRunEmphasis() const; const MeasurementObjectType* GetLongRunEmphasisOutput() const; /** Methods to return the grey level nonuniformity. */ MeasurementType GetGreyLevelNonuniformity() const; const MeasurementObjectType* GetGreyLevelNonuniformityOutput() const; /** Methods to return the run length nonuniformity. */ MeasurementType GetRunLengthNonuniformity() const; const MeasurementObjectType* GetRunLengthNonuniformityOutput() const; /** Methods to return the low grey level run emphasis. */ MeasurementType GetLowGreyLevelRunEmphasis() const; const MeasurementObjectType* GetLowGreyLevelRunEmphasisOutput() const; /** Methods to return the high grey level run emphasis. */ MeasurementType GetHighGreyLevelRunEmphasis() const; const MeasurementObjectType* GetHighGreyLevelRunEmphasisOutput() const; /** Methods to return the short run low grey level run emphasis. */ MeasurementType GetShortRunLowGreyLevelEmphasis() const; const MeasurementObjectType* GetShortRunLowGreyLevelEmphasisOutput() const; /** Methods to return the short run high grey level run emphasis. */ MeasurementType GetShortRunHighGreyLevelEmphasis() const; const MeasurementObjectType* GetShortRunHighGreyLevelEmphasisOutput() const; /** Methods to return the long run low grey level run emphasis. */ MeasurementType GetLongRunLowGreyLevelEmphasis() const; const MeasurementObjectType* GetLongRunLowGreyLevelEmphasisOutput() const; /** Methods to return the long run high grey level run emphasis. */ MeasurementType GetLongRunHighGreyLevelEmphasis() const; const MeasurementObjectType* GetLongRunHighGreyLevelEmphasisOutput() const; /** Methods to return the long run high grey level run emphasis. */ MeasurementType GetRunPercentage() const; const MeasurementObjectType* GetRunPercentageOutput() const; /** Methods to return the long run high grey level run emphasis. */ MeasurementType GetNumberOfRuns() const; const MeasurementObjectType* GetNumberOfRunsOutput() const; + /** Methods to return the grey level variance. */ + MeasurementType GetGreyLevelVariance() const; + const MeasurementObjectType* GetGreyLevelVarianceOutput() const; + + /** Methods to return the run length variance. */ + MeasurementType GetRunLengthVariance() const; + const MeasurementObjectType* GetRunLengthVarianceOutput() const; + + /** Methods to return the run entropy. */ + MeasurementType GetRunEntropy() const; + const MeasurementObjectType* GetRunEntropyOutput() const; + itkGetMacro( TotalNumberOfRuns, unsigned long ); itkGetConstMacro(NumberOfVoxels, unsigned long); itkSetMacro(NumberOfVoxels, unsigned long); /** Run-length feature types */ typedef enum { ShortRunEmphasis, LongRunEmphasis, GreyLevelNonuniformity, RunLengthNonuniformity, LowGreyLevelRunEmphasis, HighGreyLevelRunEmphasis, ShortRunLowGreyLevelEmphasis, ShortRunHighGreyLevelEmphasis, LongRunLowGreyLevelEmphasis, LongRunHighGreyLevelEmphasis, RunPercentage, - NumberOfRuns + NumberOfRuns, + GreyLevelVariance, + RunLengthVariance, + RunEntropy } RunLengthFeatureName; /** convenience method to access the run length values */ MeasurementType GetFeature( RunLengthFeatureName name ); protected: EnhancedHistogramToRunLengthFeaturesFilter(); ~EnhancedHistogramToRunLengthFeaturesFilter() {}; 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: EnhancedHistogramToRunLengthFeaturesFilter(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented unsigned long m_TotalNumberOfRuns; unsigned long m_NumberOfVoxels; }; } // end of namespace Statistics } // end of namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkEnhancedHistogramToRunLengthFeaturesFilter.hxx" #endif #endif diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.hxx index 39c7acf636..b405d40f8d 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToRunLengthFeaturesFilter.hxx @@ -1,494 +1,602 @@ /*=================================================================== 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 __itkEnhancedHistogramToRunLengthFeaturesFilter_hxx #define __itkEnhancedHistogramToRunLengthFeaturesFilter_hxx #include "itkEnhancedHistogramToRunLengthFeaturesFilter.h" #include "itkNumericTraits.h" #include "vnl/vnl_math.h" namespace itk { namespace Statistics { //constructor template EnhancedHistogramToRunLengthFeaturesFilter ::EnhancedHistogramToRunLengthFeaturesFilter() : m_NumberOfVoxels(1) { this->ProcessObject::SetNumberOfRequiredInputs( 1 ); // allocate the data objects for the outputs which are // just decorators real types - for( unsigned int i = 0; i < 12; i++ ) + for( unsigned int i = 0; i < 15; i++ ) { this->ProcessObject::SetNthOutput( i, this->MakeOutput( i ) ); } } template void EnhancedHistogramToRunLengthFeaturesFilter< THistogram> ::SetInput( const HistogramType *histogram ) { this->ProcessObject::SetNthInput( 0, const_cast( histogram ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::HistogramType * EnhancedHistogramToRunLengthFeaturesFilter< THistogram> ::GetInput() const { if ( this->GetNumberOfInputs() < 1 ) { return ITK_NULLPTR; } return itkDynamicCastInDebugMode(this->ProcessObject::GetInput( 0 ) ); } template typename EnhancedHistogramToRunLengthFeaturesFilter::DataObjectPointer EnhancedHistogramToRunLengthFeaturesFilter ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed( idx ) ) { return MeasurementObjectType::New().GetPointer(); } template void EnhancedHistogramToRunLengthFeaturesFilter< THistogram>:: GenerateData( void ) { const HistogramType * inputHistogram = this->GetInput(); this->m_TotalNumberOfRuns = static_cast ( inputHistogram->GetTotalFrequency() ); MeasurementType shortRunEmphasis = NumericTraits::ZeroValue(); MeasurementType longRunEmphasis = NumericTraits::ZeroValue(); MeasurementType greyLevelNonuniformity = NumericTraits::ZeroValue(); MeasurementType runLengthNonuniformity = NumericTraits::ZeroValue(); MeasurementType lowGreyLevelRunEmphasis = NumericTraits::ZeroValue(); MeasurementType highGreyLevelRunEmphasis = NumericTraits::ZeroValue(); MeasurementType shortRunLowGreyLevelEmphasis = NumericTraits::ZeroValue(); MeasurementType shortRunHighGreyLevelEmphasis = NumericTraits::ZeroValue(); MeasurementType longRunLowGreyLevelEmphasis = NumericTraits::ZeroValue(); MeasurementType longRunHighGreyLevelEmphasis = NumericTraits::ZeroValue(); MeasurementType runPercentage = NumericTraits::ZeroValue(); MeasurementType numberOfRuns = NumericTraits::ZeroValue(); + //Added 15.07.2016 + MeasurementType greyLevelVariance = NumericTraits::ZeroValue(); + MeasurementType runLengthVariance = NumericTraits::ZeroValue(); + MeasurementType runEntropy = NumericTraits::ZeroValue(); vnl_vector greyLevelNonuniformityVector( inputHistogram->GetSize()[0], 0.0 ); vnl_vector runLengthNonuniformityVector( inputHistogram->GetSize()[1], 0.0 ); typedef typename HistogramType::ConstIterator HistogramIterator; + + double mu_i = 0.0; + double mu_j = 0.0; + + //Calculate the means. + for ( HistogramIterator hit = inputHistogram->Begin(); + hit != inputHistogram->End(); ++hit ) + { + MeasurementType frequency = hit.GetFrequency(); + if ( frequency == 0 ) + { + continue; + } + MeasurementVectorType measurement = hit.GetMeasurementVector(); + IndexType index = hit.GetIndex(); + + double i = index[0] + 1; + double j = index[1] + 1; + + double p_ij = frequency / m_TotalNumberOfRuns; + + mu_i += i * p_ij; + mu_j += j * p_ij; + } + + //Calculate the other features. + const double log2 = std::log(2.0); + for ( HistogramIterator hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { MeasurementType frequency = hit.GetFrequency(); if ( frequency == 0 ) { continue; } MeasurementVectorType measurement = hit.GetMeasurementVector(); IndexType index = hit.GetIndex(); // inputHistogram->GetIndex( hit.GetInstanceIdentifier() ); double i2 = static_cast( ( index[0] + 1 ) * ( index[0] + 1 ) ); double j2 = static_cast( ( index[1] + 1 ) * ( index[1] + 1 ) ); + double i = index[0] + 1; + double j = index[1] + 1; + + double p_ij = frequency / m_TotalNumberOfRuns; + + greyLevelVariance += ((i - mu_i) * (i - mu_i) * p_ij); + runLengthVariance += ((j - mu_j) * (j - mu_j) * p_ij); + runEntropy -= ( p_ij > 0.0001 ) ? p_ij *std::log(p_ij) / log2 : 0; + // Traditional measures shortRunEmphasis += ( frequency / j2 ); longRunEmphasis += ( frequency * j2 ); greyLevelNonuniformityVector[index[0]] += frequency; runLengthNonuniformityVector[index[1]] += frequency; // measures from Chu et al. lowGreyLevelRunEmphasis += ( frequency / i2 ); highGreyLevelRunEmphasis += ( frequency * i2 ); // measures from Dasarathy and Holder shortRunLowGreyLevelEmphasis += ( frequency / ( i2 * j2 ) ); shortRunHighGreyLevelEmphasis += ( frequency * i2 / j2 ); longRunLowGreyLevelEmphasis += ( frequency * j2 / i2 ); longRunHighGreyLevelEmphasis += ( frequency * i2 * j2 ); } greyLevelNonuniformity = greyLevelNonuniformityVector.squared_magnitude(); runLengthNonuniformity = runLengthNonuniformityVector.squared_magnitude(); // Normalize all measures by the total number of runs if (this->m_TotalNumberOfRuns > 0) { shortRunEmphasis /= static_cast( this->m_TotalNumberOfRuns ); longRunEmphasis /= static_cast( this->m_TotalNumberOfRuns ); greyLevelNonuniformity /= static_cast( this->m_TotalNumberOfRuns ); runLengthNonuniformity /= static_cast( this->m_TotalNumberOfRuns ); lowGreyLevelRunEmphasis /= static_cast( this->m_TotalNumberOfRuns ); highGreyLevelRunEmphasis /= static_cast( this->m_TotalNumberOfRuns ); shortRunLowGreyLevelEmphasis /= static_cast( this->m_TotalNumberOfRuns ); shortRunHighGreyLevelEmphasis /= static_cast( this->m_TotalNumberOfRuns ); longRunLowGreyLevelEmphasis /= static_cast( this->m_TotalNumberOfRuns ); longRunHighGreyLevelEmphasis /= static_cast( this->m_TotalNumberOfRuns ); runPercentage = static_cast( this->m_TotalNumberOfRuns ) / static_cast( this->m_NumberOfVoxels ); numberOfRuns = static_cast( this->m_TotalNumberOfRuns ) ; } else { shortRunEmphasis = 0; longRunEmphasis = 0; greyLevelNonuniformity = 0; runLengthNonuniformity= 0; lowGreyLevelRunEmphasis = 0; highGreyLevelRunEmphasis = 0; shortRunLowGreyLevelEmphasis = 0; shortRunHighGreyLevelEmphasis= 0; longRunLowGreyLevelEmphasis = 0; longRunHighGreyLevelEmphasis = 0; runPercentage = 0; numberOfRuns = static_cast( this->m_TotalNumberOfRuns ) ; } MeasurementObjectType* shortRunEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 0 ) ); shortRunEmphasisOutputObject->Set( shortRunEmphasis ); MeasurementObjectType* longRunEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 1 ) ); longRunEmphasisOutputObject->Set( longRunEmphasis ); MeasurementObjectType* greyLevelNonuniformityOutputObject = static_cast( this->ProcessObject::GetOutput( 2 ) ); greyLevelNonuniformityOutputObject->Set( greyLevelNonuniformity ); MeasurementObjectType* runLengthNonuniformityOutputObject = static_cast( this->ProcessObject::GetOutput( 3 ) ); runLengthNonuniformityOutputObject->Set( runLengthNonuniformity ); MeasurementObjectType* lowGreyLevelRunEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 4 ) ); lowGreyLevelRunEmphasisOutputObject->Set( lowGreyLevelRunEmphasis ); MeasurementObjectType* highGreyLevelRunEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 5 ) ); highGreyLevelRunEmphasisOutputObject->Set( highGreyLevelRunEmphasis ); MeasurementObjectType* shortRunLowGreyLevelEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 6 ) ); shortRunLowGreyLevelEmphasisOutputObject->Set( shortRunLowGreyLevelEmphasis ); MeasurementObjectType* shortRunHighGreyLevelEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 7 ) ); shortRunHighGreyLevelEmphasisOutputObject->Set( shortRunHighGreyLevelEmphasis ); MeasurementObjectType* longRunLowGreyLevelEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 8 ) ); longRunLowGreyLevelEmphasisOutputObject->Set( longRunLowGreyLevelEmphasis ); MeasurementObjectType* longRunHighGreyLevelEmphasisOutputObject = static_cast( this->ProcessObject::GetOutput( 9 ) ); longRunHighGreyLevelEmphasisOutputObject->Set( longRunHighGreyLevelEmphasis ); MeasurementObjectType* runPercentagesOutputObject = static_cast( this->ProcessObject::GetOutput( 10 ) ); runPercentagesOutputObject->Set( runPercentage ); MeasurementObjectType* numberOfRunsOutputObject = static_cast( this->ProcessObject::GetOutput( 11 ) ); numberOfRunsOutputObject->Set( numberOfRuns ); + + MeasurementObjectType* greyLevelVarianceOutputObject = + static_cast( this->ProcessObject::GetOutput( 12 ) ); + greyLevelVarianceOutputObject->Set( greyLevelVariance ); + + MeasurementObjectType* runLengthVarianceOutputObject = + static_cast( this->ProcessObject::GetOutput( 13 ) ); + runLengthVarianceOutputObject->Set( runLengthVariance ); + + MeasurementObjectType* runEntropyOutputObject = + static_cast( this->ProcessObject::GetOutput( 14 ) ); + runEntropyOutputObject->Set( runEntropy ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunEmphasisOutput() const { return itkDynamicCastInDebugMode(this->ProcessObject::GetOutput( 0 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 1 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetGreyLevelNonuniformityOutput() const { return itkDynamicCastInDebugMode(this->ProcessObject::GetOutput( 2 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetRunLengthNonuniformityOutput() const { return itkDynamicCastInDebugMode(this->ProcessObject::GetOutput( 3 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetLowGreyLevelRunEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 4 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetHighGreyLevelRunEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 5 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunLowGreyLevelEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 6 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunHighGreyLevelEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 7 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunLowGreyLevelEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 8 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunHighGreyLevelEmphasisOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 9 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetRunPercentageOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 10 ) ); } template const typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* EnhancedHistogramToRunLengthFeaturesFilter ::GetNumberOfRunsOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 11 ) ); } + template + const + typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* + EnhancedHistogramToRunLengthFeaturesFilter + ::GetGreyLevelVarianceOutput() const + { + return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 12 ) ); + } + + template + const + typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* + EnhancedHistogramToRunLengthFeaturesFilter + ::GetRunLengthVarianceOutput() const + { + return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 13 ) ); + } + + template + const + typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementObjectType* + EnhancedHistogramToRunLengthFeaturesFilter + ::GetRunEntropyOutput() const + { + return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 14 ) ); + } + template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunEmphasis() const { return this->GetShortRunEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunEmphasis() const { return this->GetLongRunEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetGreyLevelNonuniformity() const { return this->GetGreyLevelNonuniformityOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetRunLengthNonuniformity() const { return this->GetRunLengthNonuniformityOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetLowGreyLevelRunEmphasis() const { return this->GetLowGreyLevelRunEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetHighGreyLevelRunEmphasis() const { return this->GetHighGreyLevelRunEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunLowGreyLevelEmphasis() const { return this->GetShortRunLowGreyLevelEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetShortRunHighGreyLevelEmphasis() const { return this->GetShortRunHighGreyLevelEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunLowGreyLevelEmphasis() const { return this->GetLongRunLowGreyLevelEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetLongRunHighGreyLevelEmphasis() const { return this->GetLongRunHighGreyLevelEmphasisOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetRunPercentage() const { return this->GetRunPercentageOutput()->Get(); } template typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetNumberOfRuns() const { return this->GetNumberOfRunsOutput()->Get(); } + template + typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType + EnhancedHistogramToRunLengthFeaturesFilter + ::GetGreyLevelVariance() const + { + return this->GetGreyLevelVarianceOutput()->Get(); + } + template + typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType + EnhancedHistogramToRunLengthFeaturesFilter + ::GetRunLengthVariance() const + { + return this->GetRunLengthVarianceOutput()->Get(); + } + template + typename EnhancedHistogramToRunLengthFeaturesFilter::MeasurementType + EnhancedHistogramToRunLengthFeaturesFilter + ::GetRunEntropy() const + { + return this->GetRunEntropyOutput()->Get(); + } + template typename EnhancedHistogramToRunLengthFeaturesFilter< THistogram>::MeasurementType EnhancedHistogramToRunLengthFeaturesFilter ::GetFeature( RunLengthFeatureName feature ) { switch( feature ) { case ShortRunEmphasis: return this->GetShortRunEmphasis(); case LongRunEmphasis: return this->GetLongRunEmphasis(); case GreyLevelNonuniformity: return this->GetGreyLevelNonuniformity(); case RunLengthNonuniformity: return this->GetRunLengthNonuniformity(); case LowGreyLevelRunEmphasis: return this->GetLowGreyLevelRunEmphasis(); case HighGreyLevelRunEmphasis: return this->GetHighGreyLevelRunEmphasis(); case ShortRunLowGreyLevelEmphasis: return this->GetShortRunLowGreyLevelEmphasis(); case ShortRunHighGreyLevelEmphasis: return this->GetShortRunHighGreyLevelEmphasis(); case LongRunLowGreyLevelEmphasis: return this->GetLongRunLowGreyLevelEmphasis(); case LongRunHighGreyLevelEmphasis: return this->GetLongRunHighGreyLevelEmphasis(); case RunPercentage: return this->GetRunPercentage(); case NumberOfRuns: return this->GetNumberOfRuns(); + case GreyLevelVariance: + return this->GetGreyLevelVariance(); + case RunLengthVariance: + return this->GetRunLengthVariance(); + case RunEntropy: + return this->GetRunEntropy(); default: return 0; } } template< typename THistogram> void EnhancedHistogramToRunLengthFeaturesFilter< THistogram>:: PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf( os,indent ); } } // end of namespace Statistics } // end of namespace itk -#endif \ No newline at end of file +#endif diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.h b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.h index ed12263e6e..0d68bb9fbd 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.h +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.h @@ -1,272 +1,274 @@ /*========================================================================= * * 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); /** 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, 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 a257e6735d..de729c8b03 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx @@ -1,656 +1,665 @@ /*========================================================================= * * 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 < 25; ++i ) + for ( int i = 0; i < 26; ++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(); + 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); } 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) 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); 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 \ No newline at end of file +#endif diff --git a/Modules/Classification/CLUtilities/include/mitkRandomImageSampler.h b/Modules/Classification/CLUtilities/include/mitkRandomImageSampler.h index b737d78330..f2f11a6bbd 100644 --- a/Modules/Classification/CLUtilities/include/mitkRandomImageSampler.h +++ b/Modules/Classification/CLUtilities/include/mitkRandomImageSampler.h @@ -1,111 +1,121 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkRandomImageSampler_h #define __mitkRandomImageSampler_h #include "MitkCLUtilitiesExports.h" //MITK #include #include "mitkImageToImageFilter.h" #include namespace mitk { enum RandomImageSamplerMode { SINGLE_ACCEPTANCE_RATE, CLASS_DEPENDEND_ACCEPTANCE_RATE, SINGLE_NUMBER_OF_ACCEPTANCE, CLASS_DEPENDEND_NUMBER_OF_ACCEPTANCE }; class MITKCLUTILITIES_EXPORT RandomImageSampler : public ImageToImageFilter { public: mitkClassMacro( RandomImageSampler , ImageToImageFilter ); itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkSetMacro(SamplingMode, RandomImageSamplerMode); itkGetConstMacro(SamplingMode, RandomImageSamplerMode); itkSetMacro(AcceptanceRate, double); itkGetConstMacro(AcceptanceRate, double); - itkSetMacro(AcceptanceRateVector, std::vector); + //itkSetMacro(AcceptanceRateVector, std::vector); + void SetAcceptanceRateVector(std::vector arg) + { + m_AcceptanceRateVector = arg; + } + itkGetConstMacro(AcceptanceRateVector, std::vector); itkSetMacro(NumberOfSamples, unsigned int); itkGetConstMacro(NumberOfSamples, unsigned int); - itkSetMacro(NumberOfSamplesVector, std::vector); + //itkSetMacro(NumberOfSamplesVector, std::vector); + void SetNumberOfSamplesVector(std::vector arg) + { + m_NumberOfSamplesVector = arg; + } + itkGetConstMacro(NumberOfSamplesVector, std::vector); private: /*! \brief standard constructor */ RandomImageSampler(); /*! \brief standard destructor */ ~RandomImageSampler(); /*! \brief Method generating the output information of this filter (e.g. image dimension, image type, etc.). The interface ImageToImageFilter requires this implementation. Everything is taken from the input image. */ virtual void GenerateOutputInformation() override; /*! \brief Method generating the output of this filter. Called in the updated process of the pipeline. This method generates the smoothed output image. */ virtual void GenerateData() override; /*! \brief Internal templated method calling the ITK bilteral filter. Here the actual filtering is performed. */ template void ItkImageProcessing(const itk::Image* itkImage); /*! \brief Internal templated method calling the ITK bilteral filter. Here the actual filtering is performed. */ template void ItkImageProcessingClassDependendSampling(const itk::Image* itkImage); /*! \brief Internal templated method calling the ITK bilteral filter. Here the actual filtering is performed. */ template void ItkImageProcessingFixedNumberSampling(const itk::Image* itkImage); /*! \brief Internal templated method calling the ITK bilteral filter. Here the actual filtering is performed. */ template void ItkImageProcessingClassDependendNumberSampling(const itk::Image* itkImage); double m_AcceptanceRate; std::vector m_AcceptanceRateVector; unsigned int m_NumberOfSamples; std::vector m_NumberOfSamplesVector; RandomImageSamplerMode m_SamplingMode; }; } //END mitk namespace #endif diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp index 70896f696c..e2839b1336 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp @@ -1,291 +1,296 @@ #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); 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; 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."); return featureList; -} \ No newline at end of file +} diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp index b35e6afa7f..3f1c7d6236 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp @@ -1,165 +1,211 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include // ITK #include #include // STL #include template void - CalculateFirstOrderStatistics(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFFirstOrderStatistics::FeatureListType & featureList, mitk::GIFFirstOrderStatistics::ParameterStruct params) +CalculateFirstOrderStatistics(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFFirstOrderStatistics::FeatureListType & featureList, mitk::GIFFirstOrderStatistics::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::LabelStatisticsImageFilter FilterType; typedef typename FilterType::HistogramType HistogramType; typedef typename HistogramType::IndexType HIndexType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); double imageRange = minMaxComputer->GetMaximum() - minMaxComputer->GetMinimum(); typename FilterType::Pointer labelStatisticsImageFilter = FilterType::New(); labelStatisticsImageFilter->SetInput( itkImage ); labelStatisticsImageFilter->SetLabelInput(maskImage); labelStatisticsImageFilter->SetUseHistograms(true); if (params.m_UseCtRange) { labelStatisticsImageFilter->SetHistogramParameters(1024.5+3096.5, -1024.5,3096.5); } else { labelStatisticsImageFilter->SetHistogramParameters(params.m_HistogramSize, minMaxComputer->GetMinimum(),minMaxComputer->GetMaximum()); } labelStatisticsImageFilter->Update(); // --------------- Range -------------------- double range = labelStatisticsImageFilter->GetMaximum(1) - labelStatisticsImageFilter->GetMinimum(1); // --------------- Uniformity, Entropy -------------------- double count = labelStatisticsImageFilter->GetCount(1); //double std_dev = labelStatisticsImageFilter->GetSigma(1); double uncorrected_std_dev = std::sqrt((count - 1) / count * labelStatisticsImageFilter->GetVariance(1)); double mean = labelStatisticsImageFilter->GetMean(1); auto histogram = labelStatisticsImageFilter->GetHistogram(1); HIndexType index; index.SetSize(1); double binWidth = histogram->GetBinMax(0, 0) - histogram->GetBinMin(0, 0); double uniformity = 0; double entropy = 0; double squared_sum = 0; double kurtosis = 0; double mean_absolut_deviation = 0; double skewness = 0; double sum_prob = 0; + double p10th = histogram->Quantile(0,0.10); + double p25th = histogram->Quantile(0,0.25); + double p75th = histogram->Quantile(0,0.75); + double p90th = histogram->Quantile(0,0.90); + double Log2=log(2); for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) { index[0] = i; double prob = histogram->GetFrequency(index); if (prob < 0.1) continue; double voxelValue = histogram->GetBinMin(0, i) +binWidth * 0.5; sum_prob += prob; squared_sum += prob * voxelValue*voxelValue; prob /= count; mean_absolut_deviation += prob* std::abs(voxelValue - mean); kurtosis +=prob* (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean); skewness += prob* (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean); uniformity += prob*prob; if (prob > 0) { entropy += prob * std::log(prob) / Log2; } } double rms = std::sqrt(squared_sum / count); kurtosis = kurtosis / (uncorrected_std_dev*uncorrected_std_dev * uncorrected_std_dev*uncorrected_std_dev); skewness = skewness / (uncorrected_std_dev*uncorrected_std_dev * uncorrected_std_dev); //mean_absolut_deviation = mean_absolut_deviation; double coveredGrayValueRange = range / imageRange; featureList.push_back(std::make_pair("FirstOrder Range",range)); featureList.push_back(std::make_pair("FirstOrder Uniformity",uniformity)); featureList.push_back(std::make_pair("FirstOrder Entropy",entropy)); featureList.push_back(std::make_pair("FirstOrder Energy",squared_sum)); featureList.push_back(std::make_pair("FirstOrder RMS",rms)); featureList.push_back(std::make_pair("FirstOrder Kurtosis",kurtosis)); featureList.push_back(std::make_pair("FirstOrder Skewness",skewness)); featureList.push_back(std::make_pair("FirstOrder Mean absolute deviation",mean_absolut_deviation)); featureList.push_back(std::make_pair("FirstOrder Covered Image Intensity Range",coveredGrayValueRange)); featureList.push_back(std::make_pair("FirstOrder Minimum",labelStatisticsImageFilter->GetMinimum(1))); featureList.push_back(std::make_pair("FirstOrder Maximum",labelStatisticsImageFilter->GetMaximum(1))); featureList.push_back(std::make_pair("FirstOrder Mean",labelStatisticsImageFilter->GetMean(1))); featureList.push_back(std::make_pair("FirstOrder Variance",labelStatisticsImageFilter->GetVariance(1))); featureList.push_back(std::make_pair("FirstOrder Sum",labelStatisticsImageFilter->GetSum(1))); featureList.push_back(std::make_pair("FirstOrder Median",labelStatisticsImageFilter->GetMedian(1))); featureList.push_back(std::make_pair("FirstOrder Standard deviation",labelStatisticsImageFilter->GetSigma(1))); featureList.push_back(std::make_pair("FirstOrder No. of Voxel",labelStatisticsImageFilter->GetCount(1))); + + featureList.push_back(std::make_pair("FirstOrder 10th Percentile",p10th)); + featureList.push_back(std::make_pair("FirstOrder 90th Percentile",p90th)); + featureList.push_back(std::make_pair("FirstOrder Interquartile Range",(p75th - p25th))); + + //Calculate the robus mean absolute deviation + //First, set all frequencies to 0 that are <10th or >90th percentile + for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) + { + index[0] = i; + MITK_WARN << "Old frequency of " << histogram->GetBinMin(0,i) << ": " << histogram->GetFrequency(index); + if(histogram->GetBinMin(0,i) < p10th) + histogram->SetFrequencyOfIndex(index, 0); + else if(histogram->GetBinMin(0,i) > p90th) + histogram->SetFrequencyOfIndex(index, 0); + MITK_WARN << "New frequency of " << histogram->GetBinMin(0,i) << ": " << histogram->GetFrequency(index); + + } + + //Calculate the mean + double meanRobust = 0.0; + for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) + { + index[0] = i; + meanRobust += histogram->GetFrequency(index) * 0.5 * (histogram->GetBinMin(0,i) + histogram->GetBinMax(0,i)); + } + meanRobust = meanRobust / histogram->GetTotalFrequency(); + + double robustMeanAbsoluteDeviation = 0.0; + for (int i = 0; i < (int)(histogram->GetSize(0)); ++i) + { + index[0] = i; + robustMeanAbsoluteDeviation += std::abs(histogram->GetFrequency(index) * + ( (0.5 * (histogram->GetBinMin(0,i) + histogram->GetBinMax(0,i))) + - meanRobust + )); + } + robustMeanAbsoluteDeviation = robustMeanAbsoluteDeviation / histogram->GetTotalFrequency(); + + featureList.push_back(std::make_pair("FirstOrder Robust Mean Absolute Deviation",robustMeanAbsoluteDeviation)); + } mitk::GIFFirstOrderStatistics::GIFFirstOrderStatistics() : m_HistogramSize(256), m_UseCtRange(false) { } mitk::GIFFirstOrderStatistics::FeatureListType mitk::GIFFirstOrderStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; ParameterStruct params; params.m_HistogramSize = this->m_HistogramSize; params.m_UseCtRange = this->m_UseCtRange; AccessByItk_3(image, CalculateFirstOrderStatistics, mask, featureList, params); return featureList; } mitk::GIFFirstOrderStatistics::FeatureNameListType mitk::GIFFirstOrderStatistics::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("FirstOrder Minimum"); featureList.push_back("FirstOrder Maximum"); featureList.push_back("FirstOrder Mean"); featureList.push_back("FirstOrder Variance"); featureList.push_back("FirstOrder Sum"); featureList.push_back("FirstOrder Median"); featureList.push_back("FirstOrder Standard deviation"); featureList.push_back("FirstOrder No. of Voxel"); return featureList; -} \ No newline at end of file +} diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp index 5591176612..2dcf34b024 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp @@ -1,222 +1,237 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include // ITK #include #include // STL #include template void CalculateGrayLevelRunLengthFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFGrayLevelRunLength::FeatureListType & featureList, mitk::GIFGrayLevelRunLength::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::Statistics::EnhancedScalarImageToRunLengthFeaturesFilter FilterType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef typename FilterType::RunLengthFeaturesFilterType TextureFilterType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); typename FilterType::Pointer filter = FilterType::New(); typename FilterType::OffsetVector::Pointer newOffset = FilterType::OffsetVector::New(); auto oldOffsets = filter->GetOffsets(); auto oldOffsetsIterator = oldOffsets->Begin(); while (oldOffsetsIterator != oldOffsets->End()) { bool continueOuterLoop = false; typename FilterType::OffsetType offset = oldOffsetsIterator->Value(); for (unsigned int i = 0; i < VImageDimension; ++i) { if (params.m_Direction == i + 2 && offset[i] != 0) { continueOuterLoop = true; } } if (params.m_Direction == 1) { offset[0] = 0; offset[1] = 0; offset[2] = 1; newOffset->push_back(offset); break; } oldOffsetsIterator++; if (continueOuterLoop) continue; newOffset->push_back(offset); } filter->SetOffsets(newOffset); // All features are required typename FilterType::FeatureNameVectorPointer requestedFeatures = FilterType::FeatureNameVector::New(); requestedFeatures->push_back(TextureFilterType::ShortRunEmphasis); requestedFeatures->push_back(TextureFilterType::LongRunEmphasis); requestedFeatures->push_back(TextureFilterType::GreyLevelNonuniformity); requestedFeatures->push_back(TextureFilterType::RunLengthNonuniformity); requestedFeatures->push_back(TextureFilterType::LowGreyLevelRunEmphasis); requestedFeatures->push_back(TextureFilterType::HighGreyLevelRunEmphasis); requestedFeatures->push_back(TextureFilterType::ShortRunLowGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::ShortRunHighGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::LongRunLowGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::LongRunHighGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::RunPercentage); requestedFeatures->push_back(TextureFilterType::NumberOfRuns); + requestedFeatures->push_back(TextureFilterType::GreyLevelVariance); + requestedFeatures->push_back(TextureFilterType::RunLengthVariance); + requestedFeatures->push_back(TextureFilterType::RunEntropy); typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New(); minMaxComputer->SetImage(itkImage); minMaxComputer->Compute(); filter->SetInput(itkImage); filter->SetMaskImage(maskImage); filter->SetRequestedFeatures(requestedFeatures); int rangeOfPixels = params.m_Range; if (rangeOfPixels < 2) rangeOfPixels = 256; if (params.m_UseCtRange) { filter->SetPixelValueMinMax((TPixel)(-1024.5),(TPixel)(3096.5)); filter->SetNumberOfBinsPerAxis(3096.5+1024.5); } else { filter->SetPixelValueMinMax(minMaxComputer->GetMinimum(),minMaxComputer->GetMaximum()); filter->SetNumberOfBinsPerAxis(rangeOfPixels); } filter->SetDistanceValueMinMax(0,rangeOfPixels); filter->Update(); auto featureMeans = filter->GetFeatureMeans (); auto featureStd = filter->GetFeatureStandardDeviations(); std::ostringstream ss; ss << rangeOfPixels; std::string strRange = ss.str(); for (std::size_t i = 0; i < featureMeans->size(); ++i) { switch (i) { case TextureFilterType::ShortRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::LongRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::GreyLevelNonuniformity : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelNonuniformity Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelNonuniformity Std.",featureStd->ElementAt(i))); break; case TextureFilterType::RunLengthNonuniformity : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthNonuniformity Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthNonuniformity Std.",featureStd->ElementAt(i))); break; case TextureFilterType::LowGreyLevelRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LowGreyLevelRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LowGreyLevelRunEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::HighGreyLevelRunEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") HighGreyLevelRunEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") HighGreyLevelRunEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::ShortRunLowGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunLowGreyLevelEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::ShortRunHighGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") ShortRunHighGreyLevelEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::LongRunLowGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunLowGreyLevelEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::LongRunHighGreyLevelEmphasis : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") LongRunHighGreyLevelEmphasis Std.",featureStd->ElementAt(i))); break; case TextureFilterType::RunPercentage : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunPercentage Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunPercentage Std.",featureStd->ElementAt(i))); break; case TextureFilterType::NumberOfRuns : featureList.push_back(std::make_pair("RunLength. ("+ strRange+") NumberOfRuns Means",featureMeans->ElementAt(i))); featureList.push_back(std::make_pair("RunLength. ("+ strRange+") NumberOfRuns Std.",featureStd->ElementAt(i))); break; + case TextureFilterType::GreyLevelVariance : + featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelVariance Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("RunLength. ("+ strRange+") GreyLevelVariance Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::RunLengthVariance : + featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthVariance Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunLengthVariance Std.",featureStd->ElementAt(i))); + break; + case TextureFilterType::RunEntropy : + featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunEntropy Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("RunLength. ("+ strRange+") RunEntropy Std.",featureStd->ElementAt(i))); + break; default: break; } } } mitk::GIFGrayLevelRunLength::GIFGrayLevelRunLength(): m_Range(1.0), m_UseCtRange(false), m_Direction(0) { } mitk::GIFGrayLevelRunLength::FeatureListType mitk::GIFGrayLevelRunLength::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; ParameterStruct params; params.m_UseCtRange=m_UseCtRange; params.m_Range = m_Range; params.m_Direction = m_Direction; AccessByItk_3(image, CalculateGrayLevelRunLengthFeatures, mask, featureList,params); return featureList; } mitk::GIFGrayLevelRunLength::FeatureNameListType mitk::GIFGrayLevelRunLength::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("RunLength. ShortRunEmphasis Means"); featureList.push_back("RunLength. ShortRunEmphasis Std."); featureList.push_back("RunLength. LongRunEmphasis Means"); featureList.push_back("RunLength. LongRunEmphasis Std."); featureList.push_back("RunLength. GreyLevelNonuniformity Means"); featureList.push_back("RunLength. GreyLevelNonuniformity Std."); featureList.push_back("RunLength. RunLengthNonuniformity Means"); featureList.push_back("RunLength. RunLengthNonuniformity Std."); featureList.push_back("RunLength. LowGreyLevelRunEmphasis Means"); featureList.push_back("RunLength. LowGreyLevelRunEmphasis Std."); featureList.push_back("RunLength. HighGreyLevelRunEmphasis Means"); featureList.push_back("RunLength. HighGreyLevelRunEmphasis Std."); featureList.push_back("RunLength. ShortRunLowGreyLevelEmphasis Means"); featureList.push_back("RunLength. ShortRunLowGreyLevelEmphasis Std."); featureList.push_back("RunLength. ShortRunHighGreyLevelEmphasis Means"); featureList.push_back("RunLength. ShortRunHighGreyLevelEmphasis Std."); featureList.push_back("RunLength. LongRunHighGreyLevelEmphasis Means"); featureList.push_back("RunLength. LongRunHighGreyLevelEmphasis Std."); return featureList; -} \ No newline at end of file +} diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp index 5a2ce2cfa0..e9d0d9debe 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp @@ -1,173 +1,307 @@ /*=================================================================== 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::GIFVolumetricStatistics::FeatureListType & featureList) { 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()); std::vector borderPoints; while(!iterator.IsAtEnd()) { if (iterator.GetCenterPixel() == 0) { ++iterator; continue; } bool border = false; for (int i = 0; i < (int)(iterator.Size()); ++i) { if (iterator.GetPixel(i) == 0) { border = true; break; } } if (border) { auto centerIndex = iterator.GetIndex(); PointType centerPoint; mask->TransformIndexToPhysicalPoint(centerIndex, centerPoint ); borderPoints.push_back(centerPoint); } ++iterator; } 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)); } mitk::GIFVolumetricStatistics::GIFVolumetricStatistics() { } mitk::GIFVolumetricStatistics::FeatureListType mitk::GIFVolumetricStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask) { FeatureListType featureList; AccessByItk_2(image, CalculateVolumeStatistic, mask, featureList); AccessByItk_1(mask, CalculateLargestDiameter, 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 compactness1 = pixelVolume / ( std::sqrt(pi) * std::pow(meshSurf, 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 compactness3 = pixelVolume / ( std::sqrt(pi) * std::pow(meshSurf, 3.0/2.0)); double sphericity=std::pow(pi,1/3.0) *std::pow(6*pixelVolume, 2.0/3.0) / meshSurf; double surfaceToVolume = meshSurf / pixelVolume; double sphericalDisproportion = meshSurf / 4 / pi / std::pow(3.0 / 4.0 / pi * pixelVolume, 2.0 / 3.0); + double asphericity = std::pow(compactness2,(-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]; + + mitk::Point3D geom; + mitk::Point3D weighted; + unsigned int noOfPx = 0; + double totalGrayValue = 0.0; + + std::vector pointsForPCA; + + MITK_WARN << xx << " " << yy << " " << zz; + 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) + { + geom[0] += x; + geom[1] += y; + geom[2] += z; + + weighted[0] += pxImage*x; + weighted[1] += pxImage*y; + weighted[2] += pxImage*z; + + noOfPx++; + totalGrayValue += pxImage; + + cv::Point3d tmp; + tmp.x = x; + tmp.y = y; + tmp.z = z; + + pointsForPCA.push_back(tmp); + } + + } + } + } + + geom[0] = geom[0] / noOfPx; + geom[1] = geom[1] / noOfPx; + geom[2] = geom[2] / noOfPx; + + weighted[0] = weighted[0] / totalGrayValue; + weighted[1] = weighted[1] / totalGrayValue; + weighted[2] = weighted[2] / totalGrayValue; + + double shift = std::sqrt( + (geom[0] - weighted[0])*(geom[0] - weighted[0]) + + (geom[1] - weighted[1])*(geom[1] - weighted[1]) + + (geom[2] - weighted[2])*(geom[2] - weighted[2]) + ); + + //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; + } + + //Perform PCA analysis + cv::PCA pca_analysis(data_pts, cv::Mat(), CV_PCA_DATA_AS_ROW); + + //Store the eigenvalues + std::vector eigen_val(3); + for (int i = 0; i < 3; ++i) + { + eigen_val[i] = pca_analysis.eigenvalues.at(0, i); + } + + std::sort (eigen_val.begin(), eigen_val.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; + 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 Center of mass Shift",shift)); + 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 Flatnes",flatness)); + 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; }