diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx index 767a315238..a60da12f62 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx @@ -1,330 +1,330 @@ /*=================================================================== 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 __itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter_hxx #define __itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter_hxx #include "itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.h" #include "itkNumericTraits.h" #include "vnl/vnl_math.h" namespace itk { namespace Statistics { //constructor template EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter() : 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 < 5; i++ ) { this->ProcessObject::SetNthOutput( i, this->MakeOutput( i ) ); } } template void EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter< THistogram> ::SetInput(const HistogramType *histogram ) { this->ProcessObject::SetNthInput( 0, const_cast( histogram ) ); } template const typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::HistogramType * EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter< THistogram> ::GetInput() const { if ( this->GetNumberOfInputs() < 1 ) { return ITK_NULLPTR; } return itkDynamicCastInDebugMode(this->ProcessObject::GetInput(0) ); } template typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::DataObjectPointer EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed( idx ) ) { return MeasurementObjectType::New().GetPointer(); } template void EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter< THistogram>:: GenerateData( void ) { const HistogramType * inputHistogram = this->GetInput(); //The nis this->m_TotalNumberOfValidVoxels = static_cast ( inputHistogram->GetTotalFrequency() ); MeasurementType coarseness = NumericTraits::ZeroValue(); MeasurementType contrast = NumericTraits::ZeroValue(); MeasurementType busyness = NumericTraits::ZeroValue(); MeasurementType complexity = NumericTraits::ZeroValue(); MeasurementType strength = NumericTraits::ZeroValue(); typedef typename HistogramType::ConstIterator HistogramIterator; long unsigned int Ng = inputHistogram->GetSize(0); std::cout << "No of bins: " << Ng << " total number of valid voxels: " << this->m_TotalNumberOfValidVoxels << std::endl; double sumPiSi = 0.0; double sumSi = 0.0; unsigned int Ngp = 0; unsigned int Nv = 0; //Calculate sum(pi*si). for ( HistogramIterator hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { IndexType index = hit.GetIndex(); MeasurementType ni = hit.GetFrequency(); double si = m_siMatrix[index[0]]; if ( ni == 0 ) { continue; } sumPiSi += ni*si; sumSi += si; Ngp++; Nv += ni; } - sumPiSi /= Nv; + sumPiSi /= (Nv > 0) ? Nv : 1; //Calculate the other features //pi = ni / Nv, so we can simply calculate with ni instead of pi double pipj = 0.0; double ipijpj = 0.0; double complexitySum = 0.0; double strengthSum = 0.0; for ( HistogramIterator hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { MeasurementVectorType measurement = hit.GetMeasurementVector(); IndexType index = hit.GetIndex(); //MITK_WARN << "current index " << index; MeasurementType ni = hit.GetFrequency(); double si = m_siMatrix[index[0]]; double i = measurement[0]; for ( HistogramIterator hit2 = inputHistogram->Begin(); hit2 != inputHistogram->End(); ++hit2 ) { MeasurementVectorType measurement2 = hit2.GetMeasurementVector(); IndexType index2 = hit2.GetIndex(); MeasurementType nj= hit2.GetFrequency(); double sj = m_siMatrix[index2[0]]; double j = measurement2[0]; pipj += ni*nj*(i-j)*(i-j); ipijpj += std::abs(i*ni - j*nj); if(ni != 0 && nj != 0) complexitySum += std::abs(i-j) * (ni*si + nj*sj)/(ni+nj); strengthSum += (ni+nj)*(i-j)*(i-j); } } //Calculate Coarseness coarseness = (sumPiSi == 0) ? 1e6 : 1.0 / sumPiSi; contrast = (Ngp <= 1 || Nv == 0) ? 0 : (1.0/(Ngp * (Ngp - 1))*pipj / Nv / Nv) * (sumSi / Nv / Nv); busyness = (Ngp <= 1 || ipijpj == 0) ? 0 : sumPiSi / ipijpj / Nv; complexity = (Nv == 0) ? 0: complexitySum / Nv / Nv; strength = (sumSi == 0 || Nv == 0) ? 0 : strengthSum / Nv / sumSi; //Calculate the other features. MeasurementObjectType* CoarsenessOutputObject = static_cast( this->ProcessObject::GetOutput( 0 ) ); CoarsenessOutputObject->Set( coarseness ); MeasurementObjectType* ContrastOutputObject = static_cast( this->ProcessObject::GetOutput( 1 ) ); ContrastOutputObject->Set( contrast ); MeasurementObjectType* BusynessOutputObject = static_cast( this->ProcessObject::GetOutput( 2 ) ); BusynessOutputObject->Set( busyness ); MeasurementObjectType* ComplexityOutputObject = static_cast( this->ProcessObject::GetOutput( 3 ) ); ComplexityOutputObject->Set( complexity ); MeasurementObjectType* StrengthOutputObject = static_cast( this->ProcessObject::GetOutput( 4 ) ); StrengthOutputObject->Set( strength ); } template const typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::MeasurementObjectType* EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetCoarsenessOutput() const { return itkDynamicCastInDebugMode(this->ProcessObject::GetOutput( 0 ) ); } template const typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::MeasurementObjectType* EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetContrastOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 1 ) ); } template const typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::MeasurementObjectType* EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetBusynessOutput() const { return itkDynamicCastInDebugMode(this->ProcessObject::GetOutput( 2 ) ); } template const typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::MeasurementObjectType* EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetComplexityOutput() const { return itkDynamicCastInDebugMode(this->ProcessObject::GetOutput( 3 ) ); } template const typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::MeasurementObjectType* EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetStrengthOutput() const { return itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 4 ) ); } template typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::MeasurementType EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetCoarseness() const { return this->GetCoarsenessOutput()->Get(); } template typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::MeasurementType EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetContrast() const { return this->GetContrastOutput()->Get(); } template typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::MeasurementType EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetBusyness() const { return this->GetBusynessOutput()->Get(); } template typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::MeasurementType EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetComplexity() const { return this->GetComplexityOutput()->Get(); } template typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter::MeasurementType EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetStrength() const { return this->GetStrengthOutput()->Get(); } template typename EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter< THistogram>::MeasurementType EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetFeature( NeighbourhoodGreyLevelDifferenceFeatureName feature ) { switch( feature ) { case Coarseness: return this->GetCoarseness(); case Contrast: return this->GetContrast(); case Busyness: return this->GetBusyness(); case Complexity: return this->GetComplexity(); case Strength: return this->GetStrength(); default: return 0; } } template< typename THistogram> void EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter< THistogram>:: PrintSelf(std::ostream& os, Indent indent) const { Superclass::PrintSelf( os,indent ); } } // end of namespace Statistics } // end of namespace itk #endif diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx index fa37529483..6957524b15 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToTextureFeaturesFilter.hxx @@ -1,760 +1,760 @@ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ #ifndef __itkEnhancedHistogramToTextureFeaturesFilter_hxx #define __itkEnhancedHistogramToTextureFeaturesFilter_hxx #include "itkEnhancedHistogramToTextureFeaturesFilter.h" #include "itkNumericTraits.h" #include "vnl/vnl_math.h" #include "itkMath.h" #define itkMacroGLCMFeature(name, id) \ template< typename THistogram > \ const \ typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * \ EnhancedHistogramToTextureFeaturesFilter< THistogram > \ ::Get##name##Output() const \ { \ return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(id) ); \ } \ \ template< typename THistogram > \ typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType \ EnhancedHistogramToTextureFeaturesFilter< THistogram > \ ::Get##name() const \ { \ return this->Get##name##Output()->Get(); \ } namespace itk { namespace Statistics { //constructor template< typename THistogram > EnhancedHistogramToTextureFeaturesFilter< THistogram >::EnhancedHistogramToTextureFeaturesFilter(void) { this->ProcessObject::SetNumberOfRequiredInputs(1); // allocate the data objects for the outputs which are // just decorators real types for ( int i = 0; i < 28; ++i ) { this->ProcessObject::SetNthOutput( i, this->MakeOutput(i) ); } } template< typename THistogram > void EnhancedHistogramToTextureFeaturesFilter< THistogram > ::SetInput(const HistogramType *histogram) { this->ProcessObject::SetNthInput( 0, const_cast< HistogramType * >( histogram ) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::HistogramType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetInput() const { return itkDynamicCastInDebugMode< const HistogramType * >( this->GetPrimaryInput() ); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::DataObjectPointer EnhancedHistogramToTextureFeaturesFilter< THistogram > ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed(idx) ) { return MeasurementObjectType::New().GetPointer(); } template< typename THistogram > void EnhancedHistogramToTextureFeaturesFilter< THistogram >::GenerateData(void) { typedef typename HistogramType::ConstIterator HistogramIterator; const HistogramType *inputHistogram = this->GetInput(); //Normalize the absolute frequencies and populate the relative frequency //container TotalRelativeFrequencyType totalFrequency = static_cast< TotalRelativeFrequencyType >( inputHistogram->GetTotalFrequency() ); m_RelativeFrequencyContainer.clear(); typename HistogramType::SizeValueType binsPerAxis = inputHistogram->GetSize(0); std::vector sumP; std::vector diffP; sumP.resize(2*binsPerAxis,0.0); diffP.resize(binsPerAxis,0.0); double numberOfPixels = 0; for ( HistogramIterator hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { AbsoluteFrequencyType frequency = hit.GetFrequency(); RelativeFrequencyType relativeFrequency = (totalFrequency > 0) ? frequency / totalFrequency : 0; m_RelativeFrequencyContainer.push_back(relativeFrequency); IndexType index = inputHistogram->GetIndex( hit.GetInstanceIdentifier() ); sumP[index[0] + index[1]] += relativeFrequency; diffP[std::abs(index[0] - index[1])] += relativeFrequency; //if (index[1] == 0) numberOfPixels += frequency; } // Now get the various means and variances. This is takes two passes // through the histogram. double pixelMean; double marginalMean; double marginalDevSquared; double pixelVariance; this->ComputeMeansAndVariances(pixelMean, marginalMean, marginalDevSquared, pixelVariance); // Finally compute the texture features. Another one pass. MeasurementType energy = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType entropy = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType correlation = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inverseDifferenceMoment = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inertia = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType clusterShade = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType clusterProminence = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType haralickCorrelation = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType autocorrelation = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType contrast = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType dissimilarity = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType maximumProbability = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inverseVariance = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType homogeneity1 = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType clusterTendency = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType variance = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType sumAverage = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType sumEntropy = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType sumVariance = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType diffAverage = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType diffEntropy = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType diffVariance = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inverseDifferenceMomentNormalized = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inverseDifferenceNormalized = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType inverseDifference = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType averageProbability = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType firstMeasureOfInformationCorrelation = NumericTraits< MeasurementType >::ZeroValue(); MeasurementType secondMeasureOfInformationCorrelation = NumericTraits< MeasurementType >::ZeroValue(); double pixelVarianceSquared = pixelVariance * pixelVariance; // Variance is only used in correlation. If variance is 0, then // (index[0] - pixelMean) * (index[1] - pixelMean) // should be zero as well. In this case, set the variance to 1. in // order to avoid NaN correlation. if( Math::FloatAlmostEqual( pixelVarianceSquared, 0.0, 4, 2*NumericTraits::epsilon() ) ) { pixelVarianceSquared = 1.; } const double log2 = std::log(2.0); typename RelativeFrequencyContainerType::const_iterator rFreqIterator = m_RelativeFrequencyContainer.begin(); - IndexType globalIndex; + IndexType globalIndex(2); for ( HistogramIterator hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { RelativeFrequencyType frequency = *rFreqIterator; ++rFreqIterator; if ( frequency == 0 ) { continue; // no use doing these calculations if we're just multiplying by // zero. } IndexType index = inputHistogram->GetIndex( hit.GetInstanceIdentifier() ); globalIndex = index; energy += frequency * frequency; entropy -= ( frequency > 0.0001 ) ? frequency *std::log(frequency) / log2:0; correlation += ( ( index[0] - pixelMean ) * ( index[1] - pixelMean ) * frequency ) / pixelVarianceSquared; inverseDifferenceMoment += frequency / ( 1.0 + ( index[0] - index[1] ) * ( index[0] - index[1] ) ); inertia += ( index[0] - index[1] ) * ( index[0] - index[1] ) * frequency; clusterShade += std::pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 3 ) * frequency; clusterProminence += std::pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 4 ) * frequency; haralickCorrelation += index[0] * index[1] * frequency; // New Features added for Aerts compatibility autocorrelation +=index[0] * index[1] * frequency; contrast += (index[0] - index[1]) * (index[0] - index[1]) * frequency; dissimilarity += std::abs(index[0] - index[1]) * frequency; maximumProbability = std::max(maximumProbability, frequency); averageProbability += frequency * index[0]; if (index[0] != index[1]) inverseVariance += frequency / ((index[0] - index[1])*(index[0] - index[1])); homogeneity1 +=frequency / ( 1.0 + std::abs( index[0] - index[1] )); clusterTendency += std::pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 2 ) * frequency; variance += std::pow( ( index[0] - pixelMean ), 2) * frequency; if (numberOfPixels > 0) { inverseDifferenceMomentNormalized += frequency / ( 1.0 + ( index[0] - index[1] ) * ( index[0] - index[1] ) / numberOfPixels / numberOfPixels); inverseDifferenceNormalized += frequency / ( 1.0 + std::abs( index[0] - index[1] ) / numberOfPixels ); } else { inverseDifferenceMomentNormalized = 0; inverseDifferenceNormalized = 0; } inverseDifference += frequency / ( 1.0 + std::abs( index[0] - index[1] ) ); } for (int i = 0; i < (int)sumP.size();++i) { double frequency = sumP[i]; sumAverage += i * frequency; sumEntropy -= ( frequency > 0.0001 ) ? frequency *std::log(frequency) / log2:0; } for (int i = 0; i < (int)sumP.size();++i) { double frequency = sumP[i]; sumVariance += (i-sumAverage)*(i-sumAverage) * frequency; } for (int i = 0; i < (int)diffP.size();++i) { double frequency = diffP[i]; diffAverage += i * frequency; diffEntropy -= ( frequency > 0.0001 ) ? frequency *std::log(frequency) / log2:0; } for (int i = 0; i < (int)diffP.size();++i) { double frequency = diffP[i]; sumVariance += (i-diffAverage)*(i-diffAverage) * frequency; } if (marginalDevSquared > 0) { haralickCorrelation = ( haralickCorrelation - marginalMean * marginalMean ) / marginalDevSquared; } else { haralickCorrelation =0; } //Calculate the margin probs std::vector pi_margins; std::vector pj_margins; //pi. for ( int i = 1; i <= inputHistogram->GetSize(0); i++ ) { double pi_tmp= 0.0; for( int j = 1; j <= inputHistogram->GetSize(1); j++ ) { globalIndex[0] = i; globalIndex[1] = j; pi_tmp += inputHistogram->GetFrequency(globalIndex); } pi_margins.push_back(pi_tmp); } //pj. for ( int j = 1; j <= inputHistogram->GetSize(1); j++ ) { double pj_tmp= 0.0; for( int i = 1; i <= inputHistogram->GetSize(0); i++ ) { globalIndex[0] = i; globalIndex[1] = j; pj_tmp += inputHistogram->GetFrequency(globalIndex); } pj_margins.push_back(pj_tmp); } //Calculate HX double hx = 0.0; for ( int i = 0; i < inputHistogram->GetSize(0); i++ ) { double pi_margin = pi_margins[i]; hx -= ( pi_margin > 0.0001 ) ? pi_margin *std::log(pi_margin) / log2:0; } //Calculate HXY1 double hxy1 = 0.0; for ( int i = 0; i < inputHistogram->GetSize(0); i++ ) { for ( int j = 0; j < inputHistogram->GetSize(1); j++ ) { globalIndex[0] = i; globalIndex[1] = j; double pi_margin = pi_margins[i]; double pj_margin = pj_margins[j]; double p_ij = inputHistogram->GetFrequency(globalIndex); hxy1 -= ( (pi_margin * pj_margin) > 0.0001 ) ? p_ij *std::log(pi_margin * pj_margin) / log2:0; } } //Calculate HXY2 double hxy2 = 0.0; for ( int i = 0; i < inputHistogram->GetSize(0); i++ ) { for ( int j = 0; j < inputHistogram->GetSize(1); j++ ) { globalIndex[0] = i; globalIndex[1] = j; double pi_margin = pi_margins[i]; double pj_margin = pj_margins[j]; hxy1 -= ( (pi_margin * pj_margin) > 0.0001 ) ? (pi_margin * pj_margin) *std::log(pi_margin * pj_margin) / log2:0; } } firstMeasureOfInformationCorrelation = (entropy - hxy1) / hx; secondMeasureOfInformationCorrelation = (entropy > hxy2) ? 0 : std::sqrt(1 - std::exp(-2*(hxy2 - entropy))); MeasurementObjectType *energyOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(0) ); energyOutputObject->Set(energy); MeasurementObjectType *entropyOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(1) ); entropyOutputObject->Set(entropy); MeasurementObjectType *correlationOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(2) ); correlationOutputObject->Set(correlation); MeasurementObjectType *inverseDifferenceMomentOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(3) ); inverseDifferenceMomentOutputObject->Set(inverseDifferenceMoment); MeasurementObjectType *inertiaOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(4) ); inertiaOutputObject->Set(inertia); MeasurementObjectType *clusterShadeOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(5) ); clusterShadeOutputObject->Set(clusterShade); MeasurementObjectType *clusterProminenceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(6) ); clusterProminenceOutputObject->Set(clusterProminence); MeasurementObjectType *haralickCorrelationOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(7) ); haralickCorrelationOutputObject->Set(haralickCorrelation); MeasurementObjectType *autocorrelationOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(8) ); autocorrelationOutputObject->Set(autocorrelation); MeasurementObjectType *contrastOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(9) ); contrastOutputObject->Set(contrast); MeasurementObjectType *dissimilarityOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(10) ); dissimilarityOutputObject->Set(dissimilarity); MeasurementObjectType *maximumProbabilityOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(11) ); maximumProbabilityOutputObject->Set(maximumProbability); MeasurementObjectType *inverseVarianceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(12) ); inverseVarianceOutputObject->Set(inverseVariance); MeasurementObjectType *homogeneity1OutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(13) ); homogeneity1OutputObject->Set(homogeneity1); MeasurementObjectType *clusterTendencyOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(14) ); clusterTendencyOutputObject->Set(clusterTendency); MeasurementObjectType *varianceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(15) ); varianceOutputObject->Set(variance); MeasurementObjectType *sumAverageOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(16) ); sumAverageOutputObject->Set(sumAverage); MeasurementObjectType *sumEntropyOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(17) ); sumEntropyOutputObject->Set(sumEntropy); MeasurementObjectType *sumVarianceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(18) ); sumVarianceOutputObject->Set(sumVariance); MeasurementObjectType *diffAverageOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(19) ); diffAverageOutputObject->Set(diffAverage); MeasurementObjectType *diffEntropyOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(20) ); diffEntropyOutputObject->Set(diffEntropy); MeasurementObjectType *diffVarianceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(21) ); diffVarianceOutputObject->Set(diffVariance); MeasurementObjectType *inverseDifferenceMomentNormalizedOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(22) ); inverseDifferenceMomentNormalizedOutputObject->Set(inverseDifferenceMomentNormalized); MeasurementObjectType *inverseDifferenceNormalizedOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(23) ); inverseDifferenceNormalizedOutputObject->Set(inverseDifferenceNormalized); MeasurementObjectType *inverseDifferenceOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(24) ); inverseDifferenceOutputObject->Set(inverseDifference); MeasurementObjectType *jointAverageOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(25) ); jointAverageOutputObject->Set(averageProbability); MeasurementObjectType *firstMeasureOfInformationCorrelationOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(26) ); firstMeasureOfInformationCorrelationOutputObject->Set(firstMeasureOfInformationCorrelation); MeasurementObjectType *secondMeasureOfInformationCorrelationOutputObject = static_cast< MeasurementObjectType * >( this->ProcessObject::GetOutput(27) ); secondMeasureOfInformationCorrelationOutputObject->Set(secondMeasureOfInformationCorrelation); } template< typename THistogram > void EnhancedHistogramToTextureFeaturesFilter< THistogram >::ComputeMeansAndVariances(double & pixelMean, double & marginalMean, double & marginalDevSquared, double & pixelVariance) { // This function takes two passes through the histogram and two passes through // an array of the same length as a histogram axis. This could probably be // cleverly compressed to one pass, but it's not clear that that's necessary. typedef typename HistogramType::ConstIterator HistogramIterator; const HistogramType *inputHistogram = this->GetInput(); // Initialize everything typename HistogramType::SizeValueType binsPerAxis = inputHistogram->GetSize(0); double *marginalSums = new double[binsPerAxis]; for ( double *ms_It = marginalSums; ms_It < marginalSums + binsPerAxis; ms_It++ ) { *ms_It = 0; } pixelMean = 0; typename RelativeFrequencyContainerType::const_iterator rFreqIterator = m_RelativeFrequencyContainer.begin(); // Ok, now do the first pass through the histogram to get the marginal sums // and compute the pixel mean HistogramIterator hit = inputHistogram->Begin(); while ( hit != inputHistogram->End() ) { RelativeFrequencyType frequency = *rFreqIterator; IndexType index = inputHistogram->GetIndex( hit.GetInstanceIdentifier() ); pixelMean += index[0] * frequency; marginalSums[index[0]] += frequency; ++hit; ++rFreqIterator; } /* Now get the mean and deviaton of the marginal sums. Compute incremental mean and SD, a la Knuth, "The Art of Computer Programming, Volume 2: Seminumerical Algorithms", section 4.2.2. Compute mean and standard deviation using the recurrence relation: M(1) = x(1), M(k) = M(k-1) + (x(k) - M(k-1) ) / k S(1) = 0, S(k) = S(k-1) + (x(k) - M(k-1)) * (x(k) - M(k)) for 2 <= k <= n, then sigma = std::sqrt(S(n) / n) (or divide by n-1 for sample SD instead of population SD). */ marginalMean = marginalSums[0]; marginalDevSquared = 0; for ( unsigned int arrayIndex = 1; arrayIndex < binsPerAxis; arrayIndex++ ) { int k = arrayIndex + 1; double M_k_minus_1 = marginalMean; double S_k_minus_1 = marginalDevSquared; double x_k = marginalSums[arrayIndex]; double M_k = M_k_minus_1 + ( x_k - M_k_minus_1 ) / k; double S_k = S_k_minus_1 + ( x_k - M_k_minus_1 ) * ( x_k - M_k ); marginalMean = M_k; marginalDevSquared = S_k; } marginalDevSquared = marginalDevSquared / binsPerAxis; rFreqIterator = m_RelativeFrequencyContainer.begin(); // OK, now compute the pixel variances. pixelVariance = 0; for ( hit = inputHistogram->Begin(); hit != inputHistogram->End(); ++hit ) { RelativeFrequencyType frequency = *rFreqIterator; IndexType index = inputHistogram->GetIndex( hit.GetInstanceIdentifier() ); pixelVariance += ( index[0] - pixelMean ) * ( index[0] - pixelMean ) * frequency; ++rFreqIterator; } delete[] marginalSums; } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetEnergyOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(0) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetEntropyOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(1) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetCorrelationOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(2) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetInverseDifferenceMomentOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(3) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetInertiaOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(4) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetClusterShadeOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(5) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetClusterProminenceOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(6) ); } template< typename THistogram > const typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementObjectType * EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetHaralickCorrelationOutput() const { return static_cast< const MeasurementObjectType * >( this->ProcessObject::GetOutput(7) ); } itkMacroGLCMFeature(Autocorrelation,8) itkMacroGLCMFeature(Contrast,9) itkMacroGLCMFeature(Dissimilarity,10) itkMacroGLCMFeature(MaximumProbability,11) itkMacroGLCMFeature(InverseVariance,12) itkMacroGLCMFeature(Homogeneity1,13) itkMacroGLCMFeature(ClusterTendency,14) itkMacroGLCMFeature(Variance,15) itkMacroGLCMFeature(SumAverage,16) itkMacroGLCMFeature(SumEntropy,17) itkMacroGLCMFeature(SumVariance,18) itkMacroGLCMFeature(DifferenceAverage,19) itkMacroGLCMFeature(DifferenceEntropy,20) itkMacroGLCMFeature(DifferenceVariance,21) itkMacroGLCMFeature(InverseDifferenceMomentNormalized,22) itkMacroGLCMFeature(InverseDifferenceNormalized,23) itkMacroGLCMFeature(InverseDifference,24) itkMacroGLCMFeature(JointAverage,25) itkMacroGLCMFeature(FirstMeasureOfInformationCorrelation,26) itkMacroGLCMFeature(SecondMeasureOfInformationCorrelation,27) template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetEnergy() const { return this->GetEnergyOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetEntropy() const { return this->GetEntropyOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetCorrelation() const { return this->GetCorrelationOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetInverseDifferenceMoment() const { return this->GetInverseDifferenceMomentOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetInertia() const { return this->GetInertiaOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetClusterShade() const { return this->GetClusterShadeOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetClusterProminence() const { return this->GetClusterProminenceOutput()->Get(); } template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetHaralickCorrelation() const { return this->GetHaralickCorrelationOutput()->Get(); } #define itkMacroGLCMFeatureSwitch(name) \ case name : \ return this->Get##name() template< typename THistogram > typename EnhancedHistogramToTextureFeaturesFilter< THistogram >::MeasurementType EnhancedHistogramToTextureFeaturesFilter< THistogram > ::GetFeature(TextureFeatureName feature) { switch ( feature ) { itkMacroGLCMFeatureSwitch(Energy); itkMacroGLCMFeatureSwitch(Entropy); itkMacroGLCMFeatureSwitch(Correlation); itkMacroGLCMFeatureSwitch(InverseDifferenceMoment); itkMacroGLCMFeatureSwitch(Inertia); itkMacroGLCMFeatureSwitch(ClusterShade); itkMacroGLCMFeatureSwitch(ClusterProminence); itkMacroGLCMFeatureSwitch(HaralickCorrelation); itkMacroGLCMFeatureSwitch(Autocorrelation); itkMacroGLCMFeatureSwitch(Contrast); itkMacroGLCMFeatureSwitch(Dissimilarity); itkMacroGLCMFeatureSwitch(MaximumProbability); itkMacroGLCMFeatureSwitch(InverseVariance); itkMacroGLCMFeatureSwitch(Homogeneity1); itkMacroGLCMFeatureSwitch(ClusterTendency); itkMacroGLCMFeatureSwitch(Variance); itkMacroGLCMFeatureSwitch(SumAverage); itkMacroGLCMFeatureSwitch(SumEntropy); itkMacroGLCMFeatureSwitch(SumVariance); itkMacroGLCMFeatureSwitch(DifferenceAverage); itkMacroGLCMFeatureSwitch(DifferenceEntropy); itkMacroGLCMFeatureSwitch(DifferenceVariance); itkMacroGLCMFeatureSwitch(InverseDifferenceMomentNormalized); itkMacroGLCMFeatureSwitch(InverseDifferenceNormalized); itkMacroGLCMFeatureSwitch(InverseDifference); itkMacroGLCMFeatureSwitch(JointAverage); itkMacroGLCMFeatureSwitch(FirstMeasureOfInformationCorrelation); itkMacroGLCMFeatureSwitch(SecondMeasureOfInformationCorrelation); default: return 0; } } #undef itkMacroGLCMFeatureSwitch template< typename THistogram > void EnhancedHistogramToTextureFeaturesFilter< THistogram > ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); } } // end of namespace Statistics } // end of namespace itk #endif diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx index 7e70611f8a..719a32fa85 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx @@ -1,414 +1,410 @@ /*=================================================================== 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 __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter_hxx #define __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter_hxx #include "itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.h" #include "itkNeighborhood.h" #include #include "vnl/vnl_math.h" namespace itk { namespace Statistics { template EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter() { this->SetNumberOfRequiredInputs( 1 ); this->SetNumberOfRequiredOutputs( 1 ); for( int i = 0; i < 2; ++i ) { this->ProcessObject::SetNthOutput( i, this->MakeOutput( i ) ); } this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator = NeighbourhoodGreyLevelDifferenceMatrixFilterType::New(); this->m_FeatureMeans = FeatureValueVector::New(); this->m_FeatureStandardDeviations = FeatureValueVector::New(); // Set the requested features to the default value: // {Energy, Entropy, InverseDifferenceMoment, Inertia, ClusterShade, // ClusterProminence} FeatureNameVectorPointer requestedFeatures = FeatureNameVector::New(); // can't directly set this->m_RequestedFeatures since it is const! requestedFeatures->push_back( NeighbourhoodGreyLevelDifferenceFeaturesFilterType::Coarseness ); requestedFeatures->push_back( NeighbourhoodGreyLevelDifferenceFeaturesFilterType::Contrast ); requestedFeatures->push_back( NeighbourhoodGreyLevelDifferenceFeaturesFilterType::Busyness ); requestedFeatures->push_back( NeighbourhoodGreyLevelDifferenceFeaturesFilterType::Complexity ); requestedFeatures->push_back( NeighbourhoodGreyLevelDifferenceFeaturesFilterType::Strength ); requestedFeatures->push_back( 20 ); this->SetRequestedFeatures( requestedFeatures ); // Set the offset directions to their defaults: half of all the possible // directions 1 pixel away. (The other half is included by symmetry.) // We use a neighborhood iterator to calculate the appropriate offsets. typedef Neighborhood NeighborhoodType; NeighborhoodType hood; hood.SetRadius( 1 ); // select all "previous" neighbors that are face+edge+vertex // connected to the current pixel. do not include the center pixel. unsigned int centerIndex = hood.GetCenterNeighborhoodIndex(); OffsetVectorPointer offsets = OffsetVector::New(); for( unsigned int d = 0; d < centerIndex; d++ ) { OffsetType offset = hood.GetOffset( d ); offsets->push_back( offset ); } this->SetOffsets( offsets ); this->m_FastCalculations = false; } template typename EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::DataObjectPointer EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed(idx) ) { return FeatureValueVectorDataObjectType::New().GetPointer(); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GenerateData(void) { if ( this->m_FastCalculations ) { this->FastCompute(); } else { this->FullCompute(); } } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::FullCompute() { int numFeatures = this->m_RequestedFeatures->size(); double *features = new double [numFeatures]; unsigned long numberOfVoxels = 0; ImageRegionConstIterator voxelCountIter(this->GetMaskImage(),this->GetMaskImage()->GetLargestPossibleRegion()); while ( ! voxelCountIter.IsAtEnd() ) { if (voxelCountIter.Get() > 0) ++numberOfVoxels; ++voxelCountIter; } // For each offset, calculate each feature typename OffsetVector::ConstIterator offsetIt; int offsetNum, featureNum; typedef typename NeighbourhoodGreyLevelDifferenceFeaturesFilterType::NeighbourhoodGreyLevelDifferenceFeatureName InternalNeighbourhoodGreyLevelDifferenceFeatureName; std::vector tVector; for( offsetIt = this->m_Offsets->Begin(), offsetNum = 0; offsetIt != this->m_Offsets->End(); offsetIt++, offsetNum++ ) { - MITK_WARN << "Adding offset " << offsetIt.Value(); + //MITK_WARN << "Adding offset " << offsetIt.Value(); tVector.push_back(offsetIt.Value()); } this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->AddOffsets( tVector); this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->Update(); typename NeighbourhoodGreyLevelDifferenceFeaturesFilterType::Pointer NeighbourhoodGreyLevelDifferenceMatrixCalculator = NeighbourhoodGreyLevelDifferenceFeaturesFilterType::New(); NeighbourhoodGreyLevelDifferenceMatrixCalculator->SetInput( this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->GetOutput() ); NeighbourhoodGreyLevelDifferenceMatrixCalculator->SetSiMatrix( this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->GetSiMatrix() ); NeighbourhoodGreyLevelDifferenceMatrixCalculator->SetNumberOfVoxels(numberOfVoxels); NeighbourhoodGreyLevelDifferenceMatrixCalculator->Update(); - double y = NeighbourhoodGreyLevelDifferenceMatrixCalculator-> - GetFeature(NeighbourhoodGreyLevelDifferenceFeaturesFilterType::Coarseness); - MITK_WARN << y; - typename FeatureNameVector::ConstIterator fnameIt; for( fnameIt = this->m_RequestedFeatures->Begin(), featureNum = 0; fnameIt != this->m_RequestedFeatures->End(); fnameIt++, featureNum++ ) { InternalNeighbourhoodGreyLevelDifferenceFeatureName tn = ( InternalNeighbourhoodGreyLevelDifferenceFeatureName )fnameIt.Value(); double xx = NeighbourhoodGreyLevelDifferenceMatrixCalculator->GetFeature(tn); features[featureNum] = xx; } // Now get the mean and deviaton of each feature across the offsets. this->m_FeatureMeans->clear(); this->m_FeatureStandardDeviations->clear(); double *tempFeatureMeans = new double[numFeatures]; double *tempFeatureDevs = new double[numFeatures]; /*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). */ // Set up the initial conditions (k = 1) for( featureNum = 0; featureNum < numFeatures; featureNum++ ) { this->m_FeatureMeans->push_back( features[featureNum] ); //tempFeatureMeans[featureNum] = features[0][featureNum]; //tempFeatureDevs[featureNum] = 0; } // Zone through the recurrence (k = 2 ... N) /*for( offsetNum = 1; offsetNum < numOffsets; offsetNum++ ) { int k = offsetNum + 1; for( featureNum = 0; featureNum < numFeatures; featureNum++ ) { double M_k_minus_1 = tempFeatureMeans[featureNum]; double S_k_minus_1 = tempFeatureDevs[featureNum]; double x_k = features[offsetNum][featureNum]; 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 ); tempFeatureMeans[featureNum] = M_k; tempFeatureDevs[featureNum] = S_k; } } for( featureNum = 0; featureNum < numFeatures; featureNum++ ) { tempFeatureDevs[featureNum] = std::sqrt( tempFeatureDevs[featureNum] / numOffsets ); this->m_FeatureMeans->push_back( tempFeatureMeans[featureNum] ); this->m_FeatureStandardDeviations->push_back( tempFeatureDevs[featureNum] ); } */ FeatureValueVectorDataObjectType *meanOutputObject = itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 0 ) ); meanOutputObject->Set( this->m_FeatureMeans ); FeatureValueVectorDataObjectType *standardDeviationOutputObject = itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 1 ) ); standardDeviationOutputObject->Set( this->m_FeatureStandardDeviations ); delete[] tempFeatureMeans; delete[] tempFeatureDevs; /*for( int i = 0; i < numOffsets; i++ ) { delete[] features[i]; }*/ delete[] features; } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::FastCompute() { // Compute the feature for the first offset typename OffsetVector::ConstIterator offsetIt = this->m_Offsets->Begin(); this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->SetOffset( offsetIt.Value() ); this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->Update(); typename NeighbourhoodGreyLevelDifferenceFeaturesFilterType::Pointer NeighbourhoodGreyLevelDifferenceMatrixCalculator = NeighbourhoodGreyLevelDifferenceFeaturesFilterType::New(); NeighbourhoodGreyLevelDifferenceMatrixCalculator->SetInput( this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->GetOutput() ); NeighbourhoodGreyLevelDifferenceMatrixCalculator->SetSiMatrix( this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->GetSiMatrix() ); NeighbourhoodGreyLevelDifferenceMatrixCalculator->Update(); typedef typename NeighbourhoodGreyLevelDifferenceFeaturesFilterType::NeighbourhoodGreyLevelDifferenceFeatureName InternalNeighbourhoodGreyLevelDifferenceFeatureName; this->m_FeatureMeans->clear(); this->m_FeatureStandardDeviations->clear(); typename FeatureNameVector::ConstIterator fnameIt; for( fnameIt = this->m_RequestedFeatures->Begin(); fnameIt != this->m_RequestedFeatures->End(); fnameIt++ ) { this->m_FeatureMeans->push_back( NeighbourhoodGreyLevelDifferenceMatrixCalculator->GetFeature( ( InternalNeighbourhoodGreyLevelDifferenceFeatureName )fnameIt.Value() ) ); this->m_FeatureStandardDeviations->push_back( 0.0 ); } FeatureValueVectorDataObjectType *meanOutputObject = itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 0 ) ); meanOutputObject->Set( this->m_FeatureMeans ); FeatureValueVectorDataObjectType *standardDeviationOutputObject = itkDynamicCastInDebugMode( this->ProcessObject::GetOutput( 1 ) ); standardDeviationOutputObject->Set( this->m_FeatureStandardDeviations ); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::SetInput( const ImageType *image ) { // Process object is not const-correct so the const_cast is required here this->ProcessObject::SetNthInput( 0, const_cast( image ) ); this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->SetInput( image ); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::SetNumberOfBinsPerAxis( unsigned int numberOfBins ) { itkDebugMacro( "setting NumberOfBinsPerAxis to " << numberOfBins ); this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->SetNumberOfBinsPerAxis( numberOfBins ); this->Modified(); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::SetPixelValueMinMax( PixelType min, PixelType max ) { itkDebugMacro( "setting Min to " << min << "and Max to " << max ); this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->SetPixelValueMinMax( min, max ); this->Modified(); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::SetDistanceValueMinMax( double min, double max ) { itkDebugMacro( "setting Min to " << min << "and Max to " << max ); this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->SetDistanceValueMinMax( min, max ); this->Modified(); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::SetMaskImage( const ImageType *image ) { // Process object is not const-correct so the const_cast is required here this->ProcessObject::SetNthInput( 1, const_cast< ImageType * >( image ) ); this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->SetMaskImage( image ); } template const TImage * EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetInput() const { if ( this->GetNumberOfInputs() < 1 ) { return ITK_NULLPTR; } return static_cast( this->ProcessObject::GetInput( 0 ) ); } template const typename EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::FeatureValueVectorDataObjectType * EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetFeatureMeansOutput() const { return itkDynamicCastInDebugMode (this->ProcessObject::GetOutput( 0 ) ); } template const typename EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::FeatureValueVectorDataObjectType * EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetFeatureStandardDeviationsOutput() const { return itkDynamicCastInDebugMode< const FeatureValueVectorDataObjectType * > ( this->ProcessObject::GetOutput( 1 ) ); } template const TImage * EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::GetMaskImage() const { if ( this->GetNumberOfInputs() < 2 ) { return ITK_NULLPTR; } return static_cast< const ImageType *>( this->ProcessObject::GetInput( 1 ) ); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::SetInsidePixelValue( PixelType insidePixelValue ) { itkDebugMacro( "setting InsidePixelValue to " << insidePixelValue ); this->m_NeighbourhoodGreyLevelDifferenceMatrixGenerator->SetInsidePixelValue( insidePixelValue ); this->Modified(); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter ::PrintSelf(std::ostream & os, Indent indent) const { Superclass::PrintSelf(os, indent); os << indent << "RequestedFeatures: " << this->GetRequestedFeatures() << std::endl; os << indent << "FeatureStandardDeviations: " << this->GetFeatureStandardDeviations() << std::endl; os << indent << "FastCalculations: " << this->GetFastCalculations() << std::endl; os << indent << "Offsets: " << this->GetOffsets() << std::endl; os << indent << "FeatureMeans: " << this->GetFeatureMeans() << std::endl; } } // end of namespace Statistics } // end of namespace itk #endif diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx index 3328660621..530611c05a 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx @@ -1,400 +1,400 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*========================================================================= * * Copyright Insight Software Consortium * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0.txt * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. * *=========================================================================*/ #ifndef __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter_hxx #define __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter_hxx #include "itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.h" #include "itkConstNeighborhoodIterator.h" #include "itkNeighborhood.h" #include "vnl/vnl_math.h" #include "itkMacro.h" #include "itkRescaleIntensityImageFilter.h" #include "itkMaskImageFilter.h" #include "itkLabelStatisticsImageFilter.h" #include "itkScalarConnectedComponentImageFilter.h" #include "itkRelabelComponentImageFilter.h" #include "itkCastImageFilter.h" #include namespace itk { namespace Statistics { template EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter() : m_NumberOfBinsPerAxis( itkGetStaticConstMacro( DefaultBinsPerAxis ) ), m_Min( NumericTraits::NonpositiveMin() ), m_Max( NumericTraits::max() ), m_MinDistance( NumericTraits::ZeroValue() ), m_MaxDistance( NumericTraits::max() ), m_InsidePixelValue( NumericTraits::OneValue() ) { this->SetNumberOfRequiredInputs( 1 ); this->SetNumberOfRequiredOutputs( 1 ); const unsigned int measurementVectorSize = 1; this->ProcessObject::SetNthOutput( 0, this->MakeOutput( 0 ) ); this->ProcessObject::SetNthOutput( 1, this->MakeOutput( 1 ) ); HistogramType *output = const_cast( this->GetOutput() ); output->SetMeasurementVectorSize( measurementVectorSize ); m_siMatrix = new double[m_NumberOfBinsPerAxis]; for(int i = 0; i < m_NumberOfBinsPerAxis; i++) { m_siMatrix[i] = 0; } this->m_LowerBound.SetSize( measurementVectorSize ); this->m_UpperBound.SetSize( measurementVectorSize ); this->m_LowerBound[0] = this->m_Min; this->m_LowerBound[1] = this->m_MinDistance; this->m_UpperBound[0] = this->m_Max; this->m_UpperBound[1] = this->m_MaxDistance; } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::SetOffset( const OffsetType offset ) { OffsetVectorPointer offsetVector = OffsetVector::New(); offsetVector->push_back( offset ); this->SetOffsets( offsetVector ); MITK_WARN << "We now have " << this->GetOffsets()->size() << " offsets in matrixgenerator"; } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::AddOffsets( const std::vector _offsets ) { OffsetVectorPointer offsetVector = OffsetVector::New(); typename OffsetVector::ConstIterator offsets; //MITK_WARN << "We have " << this->GetOffsets()->size() << " offsets!"; for( int i = 0; i < _offsets.size(); i++) { offsetVector->push_back(_offsets[i]); auto k = _offsets[i]; this->NormalizeOffsetDirection(k); offsetVector->push_back(k); } this->SetOffsets( offsetVector ); MITK_WARN << "We now have " << this->GetOffsets()->size() << " offsets in matrixgenerator"; } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::SetInput( const ImageType *image ) { // Process object is not const-correct so the const_cast is required here this->ProcessObject::SetNthInput( 0, const_cast( image ) ); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::SetMaskImage( const ImageType *image ) { // Process object is not const-correct so the const_cast is required here this->ProcessObject::SetNthInput( 1, const_cast( image ) ); } template const TImageType * EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetInput() const { if( this->GetNumberOfInputs() < 1 ) { return ITK_NULLPTR; } return static_cast( this->ProcessObject::GetInput( 0 ) ); } template const TImageType * EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetMaskImage() const { if( this->GetNumberOfInputs() < 2 ) { return ITK_NULLPTR; } return static_cast( this->ProcessObject::GetInput( 1 ) ); } template const typename EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter::HistogramType * EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetOutput() const { const HistogramType *output = static_cast( this->ProcessObject::GetOutput( 0 ) ); return output; } template double* EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetSiMatrix() const { return m_siMatrix; } template typename EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter::DataObjectPointer EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed( idx ) ) { return HistogramType::New().GetPointer(); } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GenerateData() { HistogramType *output = static_cast( this->ProcessObject::GetOutput( 0 ) ); const ImageType * inputImage = this->GetInput(); const ImageType * maskImage = this->GetMaskImage(); // First, create an appropriate histogram with the right number of bins // and mins and maxes correct for the image type. typename HistogramType::SizeType size( output->GetMeasurementVectorSize() ); size.Fill( this->m_NumberOfBinsPerAxis ); this->m_LowerBound[0] = this->m_Min; this->m_LowerBound[1] = this->m_MinDistance; this->m_UpperBound[0] = this->m_Max; this->m_UpperBound[1] = this->m_MaxDistance; output->Initialize( size, this->m_LowerBound, this->m_UpperBound ); MeasurementVectorType run( output->GetMeasurementVectorSize() ); typename HistogramType::IndexType hIndex; //Cast the image to a float image - with no respect to the incoming image //to prevent some non-templated itk issues typedef itk::Image FloatImageType; typedef itk::CastImageFilter CastFilterType; typename CastFilterType::Pointer caster = CastFilterType::New(); caster->SetInput(inputImage); caster->Update(); typename FloatImageType::Pointer floatImage = caster->GetOutput(); //Cast the mask to an unsigned short image - with no respect to the incomimg maskimage //to prevent some non-templated itk issues typedef unsigned short LabelPixelType; typedef itk::Image LabelImageType; typedef itk::CastImageFilter MaskCastFilterType; typename MaskCastFilterType::Pointer maskCaster = MaskCastFilterType::New(); maskCaster->SetInput(maskImage); maskCaster->Update(); //Set all values out of the mask to nans. typedef itk::MaskImageFilter< FloatImageType, LabelImageType, FloatImageType > MaskFilterType; typename MaskFilterType::Pointer maskFilter = MaskFilterType::New(); maskFilter->SetInput(floatImage); maskFilter->SetMaskImage(maskCaster->GetOutput()); maskFilter->SetOutsideValue( std::numeric_limits::quiet_NaN()); maskFilter->Update(); FloatImageType::Pointer floatImageMasked = maskFilter->GetOutput(); typedef ConstNeighborhoodIterator NeighborhoodIteratorType; typename NeighborhoodIteratorType::RadiusType radius; radius.Fill( 1 ); NeighborhoodIteratorType neighborIt( radius, inputImage, inputImage->GetRequestedRegion() ); for( neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt ) { const PixelType centerPixelIntensity = neighborIt.GetCenterPixel(); IndexType centerIndex = neighborIt.GetIndex(); FloatImageType::IndexType cIndex; cIndex[0] = centerIndex[0]; cIndex[1] = centerIndex[1]; cIndex[2] = centerIndex[2]; float centerFloatPixel = floatImageMasked->GetPixel(cIndex); int px = 0; PixelType sum = 0.0; bool canCalculate = true; typename OffsetVector::ConstIterator offsets; for( offsets = this->GetOffsets()->Begin(); offsets != this->GetOffsets()->End(); offsets++ ) { OffsetType offset = offsets.Value(); IndexType index; index = centerIndex + offset; if(!inputImage->GetRequestedRegion().IsInside(index)) { canCalculate = false; break; } PixelType offsetPixel = inputImage->GetPixel(index); FloatImageType::IndexType fIndex; fIndex[0] = index[0]; fIndex[1] = index[1]; fIndex[2] = index[2]; float floatPixel = floatImageMasked->GetPixel(fIndex); //We have a nan here if(floatPixel != floatPixel || centerFloatPixel!= centerFloatPixel) { canCalculate = false; break; } sum += offsetPixel; px++; } //If we have a nan in the neighbourhood, continue if(!canCalculate) continue; PixelType mean = sum / px; double si = std::abs(mean-centerPixelIntensity); run[0] = centerPixelIntensity; //Check for NaN and inf if(run[0] == run[0] && !std::isinf(std::abs(run[0]))) { output->GetIndex( run, hIndex ); output->IncreaseFrequencyOfIndex( hIndex, 1 ); m_siMatrix[hIndex[0]] += si; } - MITK_WARN << " -> In this round we added: " << centerIndex << " with value " << centerPixelIntensity << " and si = " << si; - MITK_WARN << " -> Values are now siMatrix["< Values are now niMatrix: " << output->GetFrequency(hIndex) << "/" << run; + //MITK_WARN << " -> In this round we added: " << centerIndex << " with value " << centerPixelIntensity << " and si = " << si; + //MITK_WARN << " -> Values are now siMatrix["< Values are now niMatrix: " << output->GetFrequency(hIndex) << "/" << run; } } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::SetPixelValueMinMax( PixelType min, PixelType max ) { if( this->m_Min != min || this->m_Max != max ) { itkDebugMacro( "setting Min to " << min << "and Max to " << max ); this->m_Min = min; this->m_Max = max; this->Modified(); } } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::SetDistanceValueMinMax( RealType min, RealType max ) { if( this->m_MinDistance != min || this->m_MaxDistance != max ) { itkDebugMacro( "setting MinDistance to " << min << "and MaxDistance to " << max ); this->m_MinDistance = min; this->m_MaxDistance = max; this->Modified(); } } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::PrintSelf( std::ostream& os, Indent indent ) const { Superclass::PrintSelf( os,indent ); os << indent << "Offsets: " << this->GetOffsets() << std::endl; os << indent << "Min: " << this->m_Min << std::endl; os << indent << "Max: " << this->m_Max << std::endl; os << indent << "Min distance: " << this->m_MinDistance << std::endl; os << indent << "Max distance: " << this->m_MaxDistance << std::endl; os << indent << "NumberOfBinsPerAxis: " << this->m_NumberOfBinsPerAxis << std::endl; os << indent << "InsidePixelValue: " << this->m_InsidePixelValue << std::endl; } template void EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::NormalizeOffsetDirection(OffsetType &offset) { //MITK_WARN <<" -> NGTDM old offset = " << offset; itkDebugMacro("old offset = " << offset << std::endl); int sign = 1; bool metLastNonZero = false; for (int i = offset.GetOffsetDimension()-1; i>=0; i--) { if (metLastNonZero) { offset[i] *= sign; } else if (offset[i] != 0) { sign = (offset[i] > 0 ) ? 1 : -1; metLastNonZero = true; offset[i] *= sign; } } //MITK_WARN << " ->NGTDM new offset = " << offset; itkDebugMacro("new offset = " << offset << std::endl); } } // end of namespace Statistics } // end of namespace itk #endif diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp index 3f1c7d6236..de573335b9 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp @@ -1,211 +1,209 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include // MITK #include #include #include // ITK #include #include // STL #include template void CalculateFirstOrderStatistics(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFFirstOrderStatistics::FeatureListType & featureList, mitk::GIFFirstOrderStatistics::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::LabelStatisticsImageFilter FilterType; typedef typename FilterType::HistogramType HistogramType; typedef typename HistogramType::IndexType HIndexType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typename MaskType::Pointer maskImage = MaskType::New(); mitk::CastToItkImage(mask, maskImage); 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; } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp index 8319aadb6b..d185e134cf 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp @@ -1,226 +1,226 @@ /*=================================================================== 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 CalculateGrayLevelSizeZoneFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFGrayLevelSizeZone::FeatureListType & featureList, mitk::GIFGrayLevelSizeZone::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::Statistics::EnhancedScalarImageToSizeZoneFeaturesFilter FilterType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef typename FilterType::SizeZoneFeaturesFilterType 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::SmallZoneEmphasis); requestedFeatures->push_back(TextureFilterType::LargeZoneEmphasis); requestedFeatures->push_back(TextureFilterType::GreyLevelNonuniformity); requestedFeatures->push_back(TextureFilterType::GreyLevelNonuniformityNormalized); requestedFeatures->push_back(TextureFilterType::SizeZoneNonuniformity); requestedFeatures->push_back(TextureFilterType::SizeZoneNonuniformityNormalized); requestedFeatures->push_back(TextureFilterType::LowGreyLevelZoneEmphasis); requestedFeatures->push_back(TextureFilterType::HighGreyLevelZoneEmphasis); requestedFeatures->push_back(TextureFilterType::SmallZoneLowGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::SmallZoneHighGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::LargeZoneLowGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::LargeZoneHighGreyLevelEmphasis); requestedFeatures->push_back(TextureFilterType::ZonePercentage); requestedFeatures->push_back(TextureFilterType::GreyLevelVariance); requestedFeatures->push_back(TextureFilterType::SizeZoneVariance); requestedFeatures->push_back(TextureFilterType::ZoneEntropy); 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::SmallZoneEmphasis : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") SmallZoneEmphasis Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SmallZoneEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::LargeZoneEmphasis : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") LargeZoneEmphasis Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") LargeZoneEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::GreyLevelNonuniformity : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") GreyLevelNonuniformity Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") GreyLevelNonuniformity Means",featureMeans->ElementAt(i))); break; case TextureFilterType::GreyLevelNonuniformityNormalized : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") GreyLevelNonuniformityNormalized Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") GreyLevelNonuniformityNormalized Means",featureMeans->ElementAt(i))); break; case TextureFilterType::SizeZoneNonuniformity : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") SizeZoneNonuniformity Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SizeZoneNonuniformity Means",featureMeans->ElementAt(i))); break; case TextureFilterType::SizeZoneNonuniformityNormalized : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") SizeZoneNonuniformityNormalized Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SizeZoneNonuniformityNormalized Means",featureMeans->ElementAt(i))); break; case TextureFilterType::LowGreyLevelZoneEmphasis : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") LowGreyLevelZoneEmphasis Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") LowGreyLevelZoneEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::HighGreyLevelZoneEmphasis : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") HighGreyLevelZoneEmphasis Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") HighGreyLevelZoneEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::SmallZoneLowGreyLevelEmphasis : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") SmallZoneLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SmallZoneLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::SmallZoneHighGreyLevelEmphasis : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") SmallZoneHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SmallZoneHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::LargeZoneLowGreyLevelEmphasis : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") LargeZoneLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") LargeZoneLowGreyLevelEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::LargeZoneHighGreyLevelEmphasis : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") LargeZoneHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") LargeZoneHighGreyLevelEmphasis Means",featureMeans->ElementAt(i))); break; case TextureFilterType::ZonePercentage : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") ZonePercentage Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") ZonePercentage Means",featureMeans->ElementAt(i))); break; case TextureFilterType::GreyLevelVariance : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") GreyLevelVariance Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") GreyLevelVariance Means",featureMeans->ElementAt(i))); break; case TextureFilterType::SizeZoneVariance : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") SizeZoneVariance Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") SizeZoneVariance Means",featureMeans->ElementAt(i))); break; case TextureFilterType::ZoneEntropy : - featureList.push_back(std::make_pair("SizeZone. ("+ strRange+") ZoneEntropy Means",featureMeans->ElementAt(i))); + featureList.push_back(std::make_pair("SizeZone ("+ strRange+") ZoneEntropy Means",featureMeans->ElementAt(i))); break; default: break; } } } mitk::GIFGrayLevelSizeZone::GIFGrayLevelSizeZone(): m_Range(1.0), m_UseCtRange(false), m_Direction(0) { } mitk::GIFGrayLevelSizeZone::FeatureListType mitk::GIFGrayLevelSizeZone::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, CalculateGrayLevelSizeZoneFeatures, mask, featureList,params); return featureList; } mitk::GIFGrayLevelSizeZone::FeatureNameListType mitk::GIFGrayLevelSizeZone::GetFeatureNames() { FeatureNameListType featureList; - featureList.push_back("SizeZone. SmallZoneEmphasis Means"); - featureList.push_back("SizeZone. SmallZoneEmphasis Std."); - featureList.push_back("SizeZone. LargeZoneEmphasis Means"); - featureList.push_back("SizeZone. LargeZoneEmphasis Std."); - featureList.push_back("SizeZone. GreyLevelNonuniformity Means"); - featureList.push_back("SizeZone. GreyLevelNonuniformity Std."); - featureList.push_back("SizeZone. SizeZoneNonuniformity Means"); - featureList.push_back("SizeZone. SizeZoneNonuniformity Std."); - featureList.push_back("SizeZone. LowGreyLevelZoneEmphasis Means"); - featureList.push_back("SizeZone. LowGreyLevelZoneEmphasis Std."); - featureList.push_back("SizeZone. HighGreyLevelZoneEmphasis Means"); - featureList.push_back("SizeZone. HighGreyLevelZoneEmphasis Std."); - featureList.push_back("SizeZone. SmallZoneLowGreyLevelEmphasis Means"); - featureList.push_back("SizeZone. SmallZoneLowGreyLevelEmphasis Std."); - featureList.push_back("SizeZone. SmallZoneHighGreyLevelEmphasis Means"); - featureList.push_back("SizeZone. SmallZoneHighGreyLevelEmphasis Std."); - featureList.push_back("SizeZone. LargeZoneHighGreyLevelEmphasis Means"); - featureList.push_back("SizeZone. LargeZoneHighGreyLevelEmphasis Std."); + featureList.push_back("SizeZone SmallZoneEmphasis Means"); + featureList.push_back("SizeZone SmallZoneEmphasis Std."); + featureList.push_back("SizeZone LargeZoneEmphasis Means"); + featureList.push_back("SizeZone LargeZoneEmphasis Std."); + featureList.push_back("SizeZone GreyLevelNonuniformity Means"); + featureList.push_back("SizeZone GreyLevelNonuniformity Std."); + featureList.push_back("SizeZone SizeZoneNonuniformity Means"); + featureList.push_back("SizeZone SizeZoneNonuniformity Std."); + featureList.push_back("SizeZone LowGreyLevelZoneEmphasis Means"); + featureList.push_back("SizeZone LowGreyLevelZoneEmphasis Std."); + featureList.push_back("SizeZone HighGreyLevelZoneEmphasis Means"); + featureList.push_back("SizeZone HighGreyLevelZoneEmphasis Std."); + featureList.push_back("SizeZone SmallZoneLowGreyLevelEmphasis Means"); + featureList.push_back("SizeZone SmallZoneLowGreyLevelEmphasis Std."); + featureList.push_back("SizeZone SmallZoneHighGreyLevelEmphasis Means"); + featureList.push_back("SizeZone SmallZoneHighGreyLevelEmphasis Std."); + featureList.push_back("SizeZone LargeZoneHighGreyLevelEmphasis Means"); + featureList.push_back("SizeZone LargeZoneHighGreyLevelEmphasis Std."); return featureList; } diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp index 20ba05e44c..2f7511c0af 100644 --- a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp @@ -1,175 +1,174 @@ /*=================================================================== 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 CalculateGrayLevelNeighbourhoodGreyLevelDifferenceFeatures(itk::Image* itkImage, mitk::Image::Pointer mask, mitk::GIFNeighbourhoodGreyLevelDifference::FeatureListType & featureList, mitk::GIFNeighbourhoodGreyLevelDifference::ParameterStruct params) { typedef itk::Image ImageType; typedef itk::Image MaskType; typedef itk::Statistics::EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter FilterType; typedef itk::MinimumMaximumImageCalculator MinMaxComputerType; typedef typename FilterType::NeighbourhoodGreyLevelDifferenceFeaturesFilterType 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::Coarseness); requestedFeatures->push_back(TextureFilterType::Contrast); requestedFeatures->push_back(TextureFilterType::Busyness); requestedFeatures->push_back(TextureFilterType::Complexity); requestedFeatures->push_back(TextureFilterType::Strength); 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) { - MITK_WARN << "Writing output " << i; switch (i) { case TextureFilterType::Coarseness : featureList.push_back(std::make_pair("NeighbourhoodGreyLevelDifference ("+ strRange+") Coarseness Means",featureMeans->ElementAt(i))); break; case TextureFilterType::Contrast : featureList.push_back(std::make_pair("NeighbourhoodGreyLevelDifference ("+ strRange+") Contrast Means",featureMeans->ElementAt(i))); break; case TextureFilterType::Busyness : featureList.push_back(std::make_pair("NeighbourhoodGreyLevelDifference ("+ strRange+") Busyness Means",featureMeans->ElementAt(i))); break; case TextureFilterType::Complexity : featureList.push_back(std::make_pair("NeighbourhoodGreyLevelDifference ("+ strRange+") Complexity Means",featureMeans->ElementAt(i))); break; case TextureFilterType::Strength : featureList.push_back(std::make_pair("NeighbourhoodGreyLevelDifference ("+ strRange+") Strength Means",featureMeans->ElementAt(i))); break; default: break; } } } mitk::GIFNeighbourhoodGreyLevelDifference::GIFNeighbourhoodGreyLevelDifference(): m_Range(1.0), m_UseCtRange(false), m_Direction(0) { } mitk::GIFNeighbourhoodGreyLevelDifference::FeatureListType mitk::GIFNeighbourhoodGreyLevelDifference::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, CalculateGrayLevelNeighbourhoodGreyLevelDifferenceFeatures, mask, featureList,params); return featureList; } mitk::GIFNeighbourhoodGreyLevelDifference::FeatureNameListType mitk::GIFNeighbourhoodGreyLevelDifference::GetFeatureNames() { FeatureNameListType featureList; featureList.push_back("NeighbourhoodGreyLevelDifference. Coarseness Means"); featureList.push_back("NeighbourhoodGreyLevelDifference. Coarseness Std."); featureList.push_back("NeighbourhoodGreyLevelDifference. Contrast Means"); featureList.push_back("NeighbourhoodGreyLevelDifference. Contrast Std."); featureList.push_back("NeighbourhoodGreyLevelDifference. Busyness Means"); featureList.push_back("NeighbourhoodGreyLevelDifference. Busyness Std."); featureList.push_back("NeighbourhoodGreyLevelDifference. Complexity Means"); featureList.push_back("NeighbourhoodGreyLevelDifference. Complexity Std."); featureList.push_back("NeighbourhoodGreyLevelDifference. Strength Means"); featureList.push_back("NeighbourhoodGreyLevelDifference. Strength Std."); return featureList; }