diff --git a/Modules/Classification/CLUtilities/files.cmake b/Modules/Classification/CLUtilities/files.cmake index 22a91df6b0..7e3e0adefe 100644 --- a/Modules/Classification/CLUtilities/files.cmake +++ b/Modules/Classification/CLUtilities/files.cmake @@ -1,26 +1,27 @@ file(GLOB_RECURSE H_FILES RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}/include/*") set(CPP_FILES Algorithms/itkLabelSampler.cpp Algorithms/itkSmoothedClassProbabilites.cpp Algorithms/mitkRandomImageSampler.cpp Features/itkNeighborhoodFunctorImageFilter.cpp Features/itkLineHistogramBasedMassImageFilter.cpp GlobalImageFeatures/mitkGIFCooccurenceMatrix.cpp GlobalImageFeatures/mitkGIFGrayLevelRunLength.cpp GlobalImageFeatures/mitkGIFGrayLevelSizeZone.cpp GlobalImageFeatures/mitkGIFFirstOrderStatistics.cpp GlobalImageFeatures/mitkGIFVolumetricStatistics.cpp + GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp #GlobalImageFeatures/itkEnhancedScalarImageToRunLengthFeaturesFilter.hxx #GlobalImageFeatures/itkEnhancedScalarImageToRunLengthMatrixFilter.hxx #GlobalImageFeatures/itkEnhancedHistogramToRunLengthFeaturesFilter.hxx #GlobalImageFeatures/itkEnhancedHistogramToTextureFeaturesFilter.hxx #GlobalImageFeatures/itkEnhancedScalarImageToTextureFeaturesFilter.hxx mitkCLUtil.cpp ) set( TOOL_FILES ) diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.h b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.h new file mode 100644 index 0000000000..c3952de08c --- /dev/null +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.h @@ -0,0 +1,188 @@ +/*=================================================================== + +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_h +#define __itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter_h + +#include "itkHistogram.h" +#include "itkMacro.h" +#include "itkProcessObject.h" +#include "itkSimpleDataObjectDecorator.h" + +namespace itk { + namespace Statistics { + /** \class EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter + * \brief This class computes texture feature coefficients from a grey level + * Zone-length matrix. + * + * By default, Zone length features are computed for each spatial + * direction and then averaged afterward, so it is possible to access the + * standard deviations of the texture features. These values give a clue as + * to texture anisotropy. However, doing this is much more work, because it + * involved computing one for each offset given. To compute a single matrix + * using the first offset, call FastCalculationsOn(). If this is called, + * then the texture standard deviations will not be computed (and will be set + * to zero), but texture computation will be much faster. + * + * This class is templated over the input histogram type. + * + * Print references: + * M. M. Galloway. Texture analysis using gray level Zone lengths. Computer + * Graphics and Image Processing, 4:172-179, 1975. + * + * A. Chu, C. M. Sehgal, and J. F. Greenleaf. Use of gray value distribution of + * Zone lengths for texture analysis. Pattern Recognition Letters, 11:415-420, + * 1990. + * + * B. R. Dasarathy and E. B. Holder. Image characterizations based on joint + * gray-level Zone-length distributions. Pattern Recognition Letters, 12:490-502, + * 1991. + * + * IJ article: http://hdl.handle.net/1926/1374 + * + * \sa ScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter + * \sa ScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter + * \sa EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter + * + * \author: Nick Tustison + * \ingroup ITKStatistics + */ + + template< typename THistogram > + class EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter : public ProcessObject + { + public: + /** Standard typedefs */ + typedef EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter Self; + typedef ProcessObject Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Zone-time type information (and related methods). */ + itkTypeMacro( EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter, ProcessObject ); + + /** standard New() method support */ + itkNewMacro( Self ); + + typedef THistogram HistogramType; + typedef typename HistogramType::Pointer HistogramPointer; + typedef typename HistogramType::ConstPointer HistogramConstPointer; + typedef typename HistogramType::MeasurementType MeasurementType; + typedef typename HistogramType::MeasurementVectorType MeasurementVectorType; + typedef typename HistogramType::IndexType IndexType; + typedef typename HistogramType:: + TotalAbsoluteFrequencyType FrequencyType; + + /** Method to Set/Get the input Histogram */ + using Superclass::SetInput; + void SetInput (const HistogramType * histogram ); + const HistogramType * GetInput() const; + + /** Smart Pointer type to a DataObject. */ + typedef DataObject::Pointer DataObjectPointer; + + /** Type of DataObjects used for scalar outputs */ + typedef SimpleDataObjectDecorator MeasurementObjectType; + + /** Methods to return the short Zone emphasis. */ + MeasurementType GetCoarseness() const; + const MeasurementObjectType* GetCoarsenessOutput() const; + + /** Methods to return the long Zone emphasis. */ + MeasurementType GetContrast() const; + const MeasurementObjectType* GetContrastOutput() const; + + /** Methods to return the grey level nonuniformity. */ + MeasurementType GetBusyness() const; + const MeasurementObjectType* GetBusynessOutput() const; + + /** Methods to return the Zone length nonuniformity. */ + MeasurementType GetComplexity() const; + const MeasurementObjectType* GetComplexityOutput() const; + + /** Methods to return the low grey level Zone emphasis. */ + MeasurementType GetStrength() const; + const MeasurementObjectType* GetStrengthOutput() const; + + itkGetMacro( TotalNumberOfValidVoxels, unsigned long ); + + itkGetConstMacro(NumberOfVoxels, unsigned long); + itkSetMacro(NumberOfVoxels, unsigned long); + + void SetSiMatrix(double* matrix) + { + m_siMatrix = matrix; + } + + /** Zone-length feature types */ + typedef enum + { + Coarseness, + Contrast, + Busyness, + Complexity, + Strength + } NeighbourhoodGreyLevelDifferenceFeatureName; + + /** convenience method to access the Zone length values */ + MeasurementType GetFeature( NeighbourhoodGreyLevelDifferenceFeatureName name ); + + protected: + EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter(); + ~EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter() {}; + virtual void PrintSelf(std::ostream& os, Indent indent) const ITK_OVERRIDE; + + /** Make a DataObject to be used for output output. */ + typedef ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType; + using Superclass::MakeOutput; + virtual DataObjectPointer MakeOutput( DataObjectPointerArraySizeType ) ITK_OVERRIDE; + + virtual void GenerateData() ITK_OVERRIDE; + + private: + EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + unsigned long m_TotalNumberOfValidVoxels; + unsigned long m_NumberOfVoxels; + + double* m_siMatrix; + }; + } // end of namespace Statistics +} // end of namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx" +#endif + +#endif diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx new file mode 100644 index 0000000000..7ac99c20ca --- /dev/null +++ b/Modules/Classification/CLUtilities/include/itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx @@ -0,0 +1,339 @@ +/*=================================================================== + +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; + + + //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]; + + if(ni*nj != 0) + { + MITK_WARN << "===INFO==="; + MITK_WARN << "At index " << index << "/" << index2; + MITK_WARN << " -> We have a greyvalue i = " << i << ", j = " << j; + MITK_WARN << " -> si = " << si << ", sj = " << sj; + MITK_WARN << " -> ni = " << ni << ", nj = " << nj; + } + + 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/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.h b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.h new file mode 100644 index 0000000000..29c801ebd6 --- /dev/null +++ b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.h @@ -0,0 +1,245 @@ +/*=================================================================== + +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_h +#define __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter_h + +#include "itkDataObjectDecorator.h" + +#include "itkEnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter.h" +#include "itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.h" + +namespace itk +{ + namespace Statistics + { + /** \class EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter + * \brief This class computes run length descriptions from an image. + * + * By default, run length features are computed for each spatial + * direction and then averaged afterward, so it is possible to access the + * standard deviations of the texture features. These values give a clue as + * to texture anisotropy. However, doing this is much more work, because it + * involved computing one for each offset given. To compute a single + * matrix using the first offset, call FastCalculationsOn(). If this is called, + * then the texture standard deviations will not be computed (and will be set + * to zero), but texture computation will be much faster. + * + * This class is templated over the input image type. + * + * Template Parameters: + * The image type, and the type of histogram frequency container. If you are + * using a large number of bins per axis, a sparse frequency container may be + * advisable. The default is to use a dense frequency container. + * + * Inputs and parameters: + * -# An image + * -# A mask defining the region over which texture features will be + * calculated. (Optional) + * -# The pixel value that defines the "inside" of the mask. (Optional, defaults + * to 1 if a mask is set.) + * -# The set of features to be calculated. These features are defined + * in the HistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter class. + * -# The number of intensity bins. (Optional, defaults to 256.) + * -# The set of directions (offsets) to average across. (Optional, defaults to + * {(-1, 0), (-1, -1), (0, -1), (1, -1)} for 2D images and scales analogously + * for ND images.) + * -# The pixel intensity range over which the features will be calculated. + * (Optional, defaults to the full dynamic range of the pixel type.) + * -# The distance range over which the features will be calculated. + * (Optional, defaults to the full dynamic range of double type.) + * + * In general, the default parameter values should be sufficient. + * + * Outputs: + * (1) The average value of each feature. + * (2) The standard deviation in the values of each feature. + * + * Print references: + * M. M. Galloway. Texture analysis using gray level run lengths. Computer + * Graphics and Image Processing, 4:172-179, 1975. + * + * A. Chu, C. M. Sehgal, and J. F. Greenleaf. Use of gray value distribution of + * run lengths for texture analysis. Pattern Recognition Letters, 11:415-420, + * 1990. + * + * B. R. Dasarathy and E. B. Holder. Image characterizations based on joint + * gray-level run-length distributions. Pattern Recognition Letters, 12:490-502, + * 1991. + * + * IJ article: http://hdl.handle.net/1926/1374 + * + * \sa EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter + * \sa ScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter + * \sa HistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter + * + * \author: Nick Tustison + * \ingroup ITKStatistics + */ + + template< typename TImageType, + typename THistogramFrequencyContainer = DenseFrequencyContainer2 > + class EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter:public ProcessObject + { + public: + /** Standard typedefs */ + typedef EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter Self; + typedef ProcessObject Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + /** Run-time type information (and related methods). */ + itkTypeMacro(EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter, ProcessObject); + + /** standard New() method support */ + itkNewMacro(Self); + + typedef THistogramFrequencyContainer FrequencyContainerType; + typedef TImageType ImageType; + typedef typename ImageType::Pointer ImagePointer; + + typedef typename ImageType::PixelType PixelType; + typedef typename ImageType::OffsetType OffsetType; + typedef VectorContainer< unsigned char, OffsetType > OffsetVector; + typedef typename OffsetVector::Pointer OffsetVectorPointer; + typedef typename OffsetVector::ConstPointer OffsetVectorConstPointer; + + typedef EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter< + ImageType, FrequencyContainerType > NeighbourhoodGreyLevelDifferenceMatrixFilterType; + + typedef typename NeighbourhoodGreyLevelDifferenceMatrixFilterType::HistogramType + HistogramType; + + typedef EnhancedHistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter< HistogramType > + NeighbourhoodGreyLevelDifferenceFeaturesFilterType; + + typedef short NeighbourhoodGreyLevelDifferenceFeatureName; + typedef VectorContainer FeatureNameVector; + typedef typename FeatureNameVector::Pointer FeatureNameVectorPointer; + typedef typename FeatureNameVector::ConstPointer FeatureNameVectorConstPointer; + typedef VectorContainer< unsigned char, double > FeatureValueVector; + typedef typename FeatureValueVector::Pointer FeatureValueVectorPointer; + + /** Smart Pointer type to a DataObject. */ + typedef DataObject::Pointer DataObjectPointer; + + /** Type of DataObjects used for scalar outputs */ + typedef DataObjectDecorator< FeatureValueVector > + FeatureValueVectorDataObjectType; + + const FeatureValueVectorDataObjectType * GetFeatureMeansOutput() const; + + const FeatureValueVectorDataObjectType * GetFeatureStandardDeviationsOutput() + const; + + /** Connects the input image for which the features are going to be computed + */ + using Superclass::SetInput; + void SetInput(const ImageType *); + + const ImageType * GetInput() const; + + /** Return the feature means and deviations. */ + itkGetConstReferenceObjectMacro(FeatureMeans, FeatureValueVector); + itkGetConstReferenceObjectMacro(FeatureStandardDeviations, FeatureValueVector); + + /** Set the desired feature set. Optional, for default value see above. */ + itkSetConstObjectMacro(RequestedFeatures, FeatureNameVector); + itkGetConstObjectMacro(RequestedFeatures, FeatureNameVector); + + /** Set the offsets over which the co-occurrence pairs will be computed. + Optional; for default value see above. */ + itkSetConstObjectMacro(Offsets, OffsetVector); + itkGetConstObjectMacro(Offsets, OffsetVector); + + /** Set number of histogram bins along each axis. + Optional; for default value see above. */ + void SetNumberOfBinsPerAxis(unsigned int); + + /** Set the min and max (inclusive) pixel value that will be used for + feature calculations. Optional; for default value see above. */ + void SetPixelValueMinMax(PixelType min, PixelType max); + + /** Set the min and max (inclusive) pixel value that will be used for + feature calculations. Optional; for default value see above. */ + void SetDistanceValueMinMax( double min, double max ); + + /** Connects the mask image for which the histogram is going to be computed. + Optional; for default value see above. */ + void SetMaskImage(const ImageType *); + + const ImageType * GetMaskImage() const; + + /** Set the pixel value of the mask that should be considered "inside" the + object. Optional; for default value see above. */ + void SetInsidePixelValue(PixelType InsidePixelValue); + + itkGetConstMacro(FastCalculations, bool); + itkSetMacro(FastCalculations, bool); + itkBooleanMacro(FastCalculations); + + protected: + EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter(); + virtual ~EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter() {} + virtual void PrintSelf( std::ostream & os, Indent indent ) const ITK_OVERRIDE; + + void FastCompute(); + + void FullCompute(); + + /** This method causes the filter to generate its output. */ + virtual void GenerateData() ITK_OVERRIDE; + + /** Make a DataObject to be used for output output. */ + typedef ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType; + using Superclass::MakeOutput; + virtual DataObjectPointer MakeOutput(DataObjectPointerArraySizeType) ITK_OVERRIDE; + + private: + typename NeighbourhoodGreyLevelDifferenceMatrixFilterType::Pointer m_NeighbourhoodGreyLevelDifferenceMatrixGenerator; + + FeatureValueVectorPointer m_FeatureMeans; + FeatureValueVectorPointer m_FeatureStandardDeviations; + FeatureNameVectorConstPointer m_RequestedFeatures; + OffsetVectorConstPointer m_Offsets; + bool m_FastCalculations; + }; + } // end of namespace Statistics +} // end of namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx" +#endif + +#endif diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx new file mode 100644 index 0000000000..7e70611f8a --- /dev/null +++ b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter.hxx @@ -0,0 +1,414 @@ +/*=================================================================== + +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(); + 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/itkEnhancedScalarImageToSizeZoneMatrixFilter.h b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.h similarity index 91% copy from Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.h copy to Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.h index a654543193..697f13028f 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.h +++ b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.h @@ -1,292 +1,299 @@ /*=================================================================== 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 __itkEnhancedScalarImageToSizeZoneMatrixFilter_h -#define __itkEnhancedScalarImageToSizeZoneMatrixFilter_h +#ifndef __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter_h +#define __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter_h #include "itkImage.h" #include "itkHistogram.h" #include "itkNumericTraits.h" #include "itkVectorContainer.h" namespace itk { namespace Statistics { - /** \class EnhancedScalarImageToSizeZoneMatrixFilter + /** \class EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter * \brief This class computes a run length matrix (histogram) from * a given image and a mask image if provided. Run length matrces are * used for image texture description. * * This filters creates a grey-level run length matrix from a N-D scalar * image. This is another possible texture description. See the following * references. * M. M. Galloway. Texture analysis using gray level run lengths. Computer * Graphics and Image Processing, 4:172-179, 1975. * * A. Chu, C. M. Sehgal, and J. F. Greenleaf. Use of gray value distribution of * run lengths for texture analysis. Pattern Recognition Letters, 11:415-420, * 1990. * * B. R. Dasarathy and E. B. Holder. Image characterizations based on joint * gray-level run-length distributions. Pattern Recognition Letters, 12:490-502, * 1991. * * The basic idea is as follows: * Given an image and an offset (e.g. (1, -1) for a 2-d image), each element * in the joint histogram describes the frequency for a particular distance/ * intensity pair within a given image. This distance/intensity pair can be * described as follows: we start at a given voxel which has some intensity. * We then "jump" to neighboring pixels in increments provided by the offset(s) * as long as the pixel to which we are jumping is within the same intensity * bin as the original voxel. The distance component is given by the distance * from the original to the final voxel satisfying our jumping criteria. * * The offset (or offsets) along which the co-occurences are calculated can be * set by the user. Traditionally, only one offset is used per histogram, and * offset components in the range [-1, 1] are used. For rotation-invariant * features averages of features computed over several histograms with different * offsets are generally used, instead of computing features from one histogram * create with several offsets. Additionally, instead of using offsets of two or * more pixels in any direction, multi-resolution techniques (e.g. image * pyramids) are generally used to deal with texture at different spatial * resolutions. * * This class calculates a 2-d histogram of all the intensity/distance pairs in * the given image's requested region, for a given set of offsets. That is, if * a given offset falls outside of the requested region (or outside the mask) * at a particular point, that distance/intensity pair will not be added to * the matrix. * * The number of histogram bins on each axis can be set (defaults to 256). Also, * by default the histogram min and max corresponds to the largest and smallest * possible pixel value of that pixel type. To customize the histogram bounds * for a given image, the max and min pixel values that will be placed in the * histogram can be set manually. NB: The min and max are INCLUSIVE. * * Further, the type of histogram frequency container used is an optional * template parameter. By default, a dense container is used, but for images * with little texture or in cases where the user wants more histogram bins, * a sparse container can be used for the histogram instead. * * WARNING: This probably won't work for pixels of double or long-double type * unless you set the histogram min and max manually. This is because the largest * histogram bin by default has max value of the largest possible pixel value * plus 1. For double and long-double types, whose "RealType" as defined by the * NumericTraits class is the same, and thus cannot hold any larger values, * this would cause a float overflow. * * IJ article: http://hdl.handle.net/1926/1374 * - * \sa ScalarImageToSizeZoneFeaturesFilter - * \sa EnhancedScalarImageToSizeZoneMatrixFilter - * \sa HistogramToSizeZoneFeaturesFilter + * \sa ScalarImageToNeighbourhoodGreyLevelDifferenceFeaturesFilter + * \sa EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter + * \sa HistogramToNeighbourhoodGreyLevelDifferenceFeaturesFilter * * \author: Nick Tustison * \ingroup ITKStatistics */ template - class EnhancedScalarImageToSizeZoneMatrixFilter : public ProcessObject + class EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter : public ProcessObject { public: /** Standard typedefs */ - typedef EnhancedScalarImageToSizeZoneMatrixFilter Self; + typedef EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter Self; typedef ProcessObject Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; /** Run-time type information (and related methods). */ - itkTypeMacro( EnhancedScalarImageToSizeZoneMatrixFilter, ProcessObject ); + itkTypeMacro( EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter, ProcessObject ); /** standard New() method support */ itkNewMacro( Self ); typedef TImageType ImageType; typedef typename ImageType::Pointer ImagePointer; typedef typename ImageType::ConstPointer ImageConstPointer; typedef typename ImageType::PixelType PixelType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::RegionType RegionType; typedef typename ImageType::SizeType RadiusType; typedef typename ImageType::OffsetType OffsetType; typedef VectorContainer OffsetVector; typedef typename OffsetVector::Pointer OffsetVectorPointer; typedef typename ImageType::PointType PointType; typedef typename NumericTraits::RealType MeasurementType; typedef typename NumericTraits::RealType RealType; typedef Histogram HistogramType; typedef typename HistogramType::Pointer HistogramPointer; typedef typename HistogramType::ConstPointer HistogramConstPointer; typedef typename HistogramType::MeasurementVectorType MeasurementVectorType; /** ImageDimension constants */ itkStaticConstMacro( ImageDimension, unsigned int, TImageType::ImageDimension ); /** Specify the default number of bins per axis */ itkStaticConstMacro( DefaultBinsPerAxis, unsigned int, 256 ); /** * Set the offsets over which the intensity/distance pairs will be computed. * Invoking this function clears the previous offsets. * Note: for each individual offset in the OffsetVector, the rightmost non-zero * offset element must be positive. For example, in the offset list of a 2D image, * (1, 0) means the offset along x-axis. (1, 0) has to be set instead * of (-1, 0). This is required from the iterating order of pixel iterator. * */ itkSetObjectMacro( Offsets, OffsetVector ); /** * Set offset over which the intensity/distance pairs will be computed. * Invoking this function clears the previous offset(s). * Note: for each individual offset, the rightmost non-zero * offset element must be positive. For example, in the offset list of a 2D image, * (1, 0) means the offset along x-axis. (1, 0) has to be set instead * of (-1, 0). This is required from the iterating order of pixel iterator. * */ void SetOffset( const OffsetType offset ); + void AddOffsets( const std::vector offset ); + /** * Get the current offset(s). */ itkGetModifiableObjectMacro(Offsets, OffsetVector ); /** Set number of histogram bins along each axis */ itkSetMacro( NumberOfBinsPerAxis, unsigned int ); /** Get number of histogram bins along each axis */ itkGetConstMacro( NumberOfBinsPerAxis, unsigned int ); /** * Set the min and max (inclusive) pixel value that will be used in * generating the histogram. */ void SetPixelValueMinMax( PixelType min, PixelType max ); /** Get the min pixel value defining one dimension of the joint histogram. */ itkGetConstMacro( Min, PixelType ); /** Get the max pixel value defining one dimension of the joint histogram. */ itkGetConstMacro( Max, PixelType ); /** * Set the min and max (inclusive) pixel value that will be used in * generating the histogram. */ void SetDistanceValueMinMax( RealType min, RealType max ); /** * Get the min distance value defining one dimension of the joint histogram. */ itkGetConstMacro( MinDistance, RealType ); /** * Get the max distance value defining one dimension of the joint histogram. */ itkGetConstMacro( MaxDistance, RealType ); /** Method to set the input image */ using Superclass::SetInput; void SetInput( const ImageType *image ); /** Method to get the input image */ const ImageType * GetInput() const; /** Method to set the mask image */ void SetMaskImage( const ImageType *image ); /** Method to get the mask image */ const ImageType * GetMaskImage() const; /** method to get the Histogram */ const HistogramType * GetOutput() const; + /** method to get the Histogram */ + double* GetSiMatrix() const; + /** * Set the pixel value of the mask that should be considered "inside" the * object. Defaults to 1. */ itkSetMacro( InsidePixelValue, PixelType ); itkGetConstMacro( InsidePixelValue, PixelType ); protected: - EnhancedScalarImageToSizeZoneMatrixFilter(); - virtual ~EnhancedScalarImageToSizeZoneMatrixFilter() {}; + EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter(); + virtual ~EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter() {}; virtual void PrintSelf( std::ostream& os, Indent indent ) const ITK_OVERRIDE; /** Standard itk::ProcessObject subclass method. */ typedef DataObject::Pointer DataObjectPointer; typedef ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType; using Superclass::MakeOutput; virtual DataObjectPointer MakeOutput( DataObjectPointerArraySizeType idx ) ITK_OVERRIDE; /** This method causes the filter to generate its output. */ virtual void GenerateData() ITK_OVERRIDE; /** * Normalize the direction of the offset before it is applied. * The last non-zero dimension of the offest has to be positive in order * to match to scanning order of the iterator. Only the sign is changed. * For example, the input offset (-1, 0) will be normalized as * (1, 0). * */ void NormalizeOffsetDirection(OffsetType &offset); private: unsigned int m_NumberOfBinsPerAxis; PixelType m_Min; PixelType m_Max; RealType m_MinDistance; RealType m_MaxDistance; PixelType m_InsidePixelValue; MeasurementVectorType m_LowerBound; MeasurementVectorType m_UpperBound; OffsetVectorPointer m_Offsets; + + double * m_siMatrix; }; } // end of namespace Statistics } // end of namespace itk #ifndef ITK_MANUAL_INSTANTIATION -#include "itkEnhancedScalarImageToSizeZoneMatrixFilter.hxx" +#include "itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx" #endif #endif diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx similarity index 56% copy from Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.hxx copy to Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx index 8d234c374b..85b71ccdb5 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter.hxx @@ -1,431 +1,397 @@ /*=================================================================== 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 __itkEnhancedScalarImageToSizeZoneMatrixFilter_hxx -#define __itkEnhancedScalarImageToSizeZoneMatrixFilter_hxx +#ifndef __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter_hxx +#define __itkEnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter_hxx -#include "itkEnhancedScalarImageToSizeZoneMatrixFilter.h" +#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 -EnhancedScalarImageToSizeZoneMatrixFilter -::EnhancedScalarImageToSizeZoneMatrixFilter() : +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 = 2; + 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]; + m_siMatrix[0] = 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 -EnhancedScalarImageToSizeZoneMatrixFilter +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 -EnhancedScalarImageToSizeZoneMatrixFilter +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 -EnhancedScalarImageToSizeZoneMatrixFilter +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 * -EnhancedScalarImageToSizeZoneMatrixFilter +EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetInput() const { if( this->GetNumberOfInputs() < 1 ) { return ITK_NULLPTR; } return static_cast( this->ProcessObject::GetInput( 0 ) ); } template const TImageType * -EnhancedScalarImageToSizeZoneMatrixFilter +EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetMaskImage() const { if( this->GetNumberOfInputs() < 2 ) { return ITK_NULLPTR; } return static_cast( this->ProcessObject::GetInput( 1 ) ); } template -const typename EnhancedScalarImageToSizeZoneMatrixFilter::HistogramType * -EnhancedScalarImageToSizeZoneMatrixFilter +EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::GetOutput() const { const HistogramType *output = static_cast( this->ProcessObject::GetOutput( 0 ) ); return output; } + +template +double* EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter +::GetSiMatrix() const +{ + return m_siMatrix; +} + template -typename EnhancedScalarImageToSizeZoneMatrixFilter::DataObjectPointer -EnhancedScalarImageToSizeZoneMatrixFilter +EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed( idx ) ) { return HistogramType::New().GetPointer(); } template void -EnhancedScalarImageToSizeZoneMatrixFilter +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(); - MITK_WARN << "InputImage casted."; - //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(); - MITK_WARN << "MaskImage casted."; - - //Set all values out of the mask to (m_Min + m_Max) / 2. + //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((m_Max + m_Min) / 2); + maskFilter->SetOutsideValue( std::numeric_limits::quiet_NaN()); maskFilter->Update(); + FloatImageType::Pointer floatImageMasked = maskFilter->GetOutput(); - MITK_WARN << "InputImage masked."; + typedef ConstNeighborhoodIterator NeighborhoodIteratorType; + typename NeighborhoodIteratorType::RadiusType radius; + radius.Fill( 1 ); + NeighborhoodIteratorType neighborIt( radius, + inputImage, inputImage->GetRequestedRegion() ); - //Rescale intensity to match the size of the histogram - typedef itk::Image< unsigned int, 3 > OutputImageType; + 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; - typedef itk::RescaleIntensityImageFilter< FloatImageType,OutputImageType> RescalerType; - typename RescalerType::Pointer rescaler = RescalerType::New(); - //We use 0 for nans, all valid numbers will be 1 < x < size - rescaler->SetOutputMinimum( 1 ); - rescaler->SetOutputMaximum( size[0] ); - rescaler->SetInput(maskFilter->GetOutput()); - rescaler->Update(); + index = centerIndex + offset; - typename OutputImageType::Pointer rescaled = rescaler->GetOutput(); + if(!inputImage->GetRequestedRegion().IsInside(index)) + { + canCalculate = false; + break; + } - MITK_WARN << "Intensities rescaled."; + PixelType offsetPixel = inputImage->GetPixel(index); - //Write back the nans because they get lost during rescaling + FloatImageType::IndexType fIndex; + fIndex[0] = index[0]; + fIndex[1] = index[1]; + fIndex[2] = index[2]; - int xx = inputImage->GetLargestPossibleRegion().GetSize()[0]; - int yy = inputImage->GetLargestPossibleRegion().GetSize()[1]; - int zz = inputImage->GetLargestPossibleRegion().GetSize()[2]; + float floatPixel = floatImageMasked->GetPixel(fIndex); - for (int x = 0; x < xx; x++) - { - for (int y = 0; y < yy; y++) - { - for (int z = 0; z < zz; z++) + //We have a nan here + if(floatPixel != floatPixel || centerFloatPixel!= centerFloatPixel) { - FloatImageType::IndexType indexF; - indexF[0] = x; - indexF[1] = y; - indexF[2] = z; - - OutputImageType::IndexType indexO; - indexO[0] = x; - indexO[1] = y; - indexO[2] = z; - - //Is Pixel NaN? - if(floatImage->GetPixel(indexF) != floatImage->GetPixel(indexF)) - { - rescaled->SetPixel(indexO,-1); - } + canCalculate = false; + break; } + + sum += offsetPixel; + px++; } - } + //If we have a nan in the neighbourhood, continue + if(!canCalculate) + continue; - OutputImageType::IndexType indexO; - indexO[0] = 0; - indexO[1] = 2; - indexO[2] = 1; - MITK_WARN << "is -1: " << rescaled->GetPixel(indexO); - indexO[0] = 0; - indexO[1] = 0; - indexO[2] = 0; - MITK_WARN << "is 1: " << rescaled->GetPixel(indexO); - - PixelType distanceThreshold = 1; - - //Calculate the connected components - typedef itk::ScalarConnectedComponentImageFilter - ConnectedComponentImageFilterType; - - typename ConnectedComponentImageFilterType::Pointer connected = ConnectedComponentImageFilterType::New (); - connected->SetInput(rescaled); - connected->SetMaskImage(maskCaster->GetOutput()); - connected->SetDistanceThreshold(distanceThreshold); - connected->Update(); - - MITK_WARN << "Connected components calculated."; - - //Relabel the components - typedef itk::RelabelComponentImageFilter RelabelFilterType; - typename RelabelFilterType::Pointer relabel = RelabelFilterType::New(); - - typename RelabelFilterType::ObjectSizeType minSize = 1; - - relabel->SetInput(connected->GetOutput()); - relabel->SetMinimumObjectSize(minSize); - relabel->Update(); - - MITK_WARN << "Components relabeled."; - - //Get the stats of the componentes - typedef itk::LabelStatisticsImageFilter< FloatImageType, LabelImageType> LabelStatisticsImageFilterType; - typename LabelStatisticsImageFilterType::Pointer labelStatisticsImageFilter = - LabelStatisticsImageFilterType::New(); - labelStatisticsImageFilter->SetLabelInput( relabel->GetOutput() ); - labelStatisticsImageFilter->SetInput(floatImage); - labelStatisticsImageFilter->UseHistogramsOn(); // needed to compute median - labelStatisticsImageFilter->Update(); - - /* - std::cout << "Number of labels: " - << labelStatisticsImageFilter->GetNumberOfLabels() << std::endl; - std::cout << std::endl; - */ - typedef typename LabelStatisticsImageFilterType::ValidLabelValuesContainerType ValidLabelValuesType; - - for(typename ValidLabelValuesType::const_iterator vIt = labelStatisticsImageFilter->GetValidLabelValues().begin(); - vIt != labelStatisticsImageFilter->GetValidLabelValues().end(); - ++vIt) - { - if ( labelStatisticsImageFilter->HasLabel(*vIt) ) + 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]))) { - LabelPixelType labelValue = *vIt; -/* - MITK_INFO << "Label: " << *vIt; - - MITK_INFO << "\tmin: " - << labelStatisticsImageFilter->GetMinimum( labelValue ); - MITK_INFO << "\tmax: " - << labelStatisticsImageFilter->GetMaximum( labelValue ); - MITK_INFO << "\tmedian: " - << labelStatisticsImageFilter->GetMedian( labelValue ); - MITK_INFO << "\tmean: " - << labelStatisticsImageFilter->GetMean( labelValue ); - MITK_INFO << "\tsigma: " - << labelStatisticsImageFilter->GetSigma( labelValue ); - MITK_INFO << "\tvariance: " - << labelStatisticsImageFilter->GetVariance( labelValue ); - MITK_INFO << "\tsum: " - << labelStatisticsImageFilter->GetSum( labelValue ); - MITK_INFO << "\tcount: " - << labelStatisticsImageFilter->GetCount( labelValue ); - MITK_INFO << "\tregion: " - << labelStatisticsImageFilter->GetRegion( labelValue ); - -*/ - run[0] = (labelStatisticsImageFilter->GetMinimum( labelValue ) + labelStatisticsImageFilter->GetMinimum( labelValue )) /2; - run[1] = labelStatisticsImageFilter->GetCount( labelValue ); - - //Check for NaN and inf - if(run[0] == run[0] && !std::isinf(std::abs(run[0]))) - { - output->GetIndex( run, hIndex ); - output->IncreaseFrequencyOfIndex( hIndex, 1 ); - } + output->GetIndex( run, hIndex ); + output->IncreaseFrequencyOfIndex( hIndex, 1 ); + + m_siMatrix[hIndex[0]] += si; } + MITK_WARN << " -> In this round we added: " << centerIndex << " with value " << centerPixelIntensity << " and si = " << si; + MITK_WARN << " -> Values are now siMatrix["< Values are now niMatrix: " << output->GetFrequency(hIndex) << "/" << run; } } + template void -EnhancedScalarImageToSizeZoneMatrixFilter +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 -EnhancedScalarImageToSizeZoneMatrixFilter +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 -EnhancedScalarImageToSizeZoneMatrixFilter +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 -EnhancedScalarImageToSizeZoneMatrixFilter +EnhancedScalarImageToNeighbourhoodGreyLevelDifferenceMatrixFilter ::NormalizeOffsetDirection(OffsetType &offset) { - MITK_WARN <<" -> Size zone old offset = " << 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 << " -> size zone new offset = " << offset; + //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/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.h b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.h index a654543193..c0cf3f65f5 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.h +++ b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.h @@ -1,292 +1,283 @@ /*=================================================================== 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 __itkEnhancedScalarImageToSizeZoneMatrixFilter_h #define __itkEnhancedScalarImageToSizeZoneMatrixFilter_h #include "itkImage.h" #include "itkHistogram.h" #include "itkNumericTraits.h" #include "itkVectorContainer.h" namespace itk { namespace Statistics { /** \class EnhancedScalarImageToSizeZoneMatrixFilter * \brief This class computes a run length matrix (histogram) from * a given image and a mask image if provided. Run length matrces are * used for image texture description. * * This filters creates a grey-level run length matrix from a N-D scalar * image. This is another possible texture description. See the following * references. * M. M. Galloway. Texture analysis using gray level run lengths. Computer * Graphics and Image Processing, 4:172-179, 1975. * * A. Chu, C. M. Sehgal, and J. F. Greenleaf. Use of gray value distribution of * run lengths for texture analysis. Pattern Recognition Letters, 11:415-420, * 1990. * * B. R. Dasarathy and E. B. Holder. Image characterizations based on joint * gray-level run-length distributions. Pattern Recognition Letters, 12:490-502, * 1991. * * The basic idea is as follows: * Given an image and an offset (e.g. (1, -1) for a 2-d image), each element * in the joint histogram describes the frequency for a particular distance/ * intensity pair within a given image. This distance/intensity pair can be * described as follows: we start at a given voxel which has some intensity. * We then "jump" to neighboring pixels in increments provided by the offset(s) * as long as the pixel to which we are jumping is within the same intensity * bin as the original voxel. The distance component is given by the distance * from the original to the final voxel satisfying our jumping criteria. * * The offset (or offsets) along which the co-occurences are calculated can be * set by the user. Traditionally, only one offset is used per histogram, and * offset components in the range [-1, 1] are used. For rotation-invariant * features averages of features computed over several histograms with different * offsets are generally used, instead of computing features from one histogram * create with several offsets. Additionally, instead of using offsets of two or * more pixels in any direction, multi-resolution techniques (e.g. image * pyramids) are generally used to deal with texture at different spatial * resolutions. * * This class calculates a 2-d histogram of all the intensity/distance pairs in * the given image's requested region, for a given set of offsets. That is, if * a given offset falls outside of the requested region (or outside the mask) * at a particular point, that distance/intensity pair will not be added to * the matrix. * * The number of histogram bins on each axis can be set (defaults to 256). Also, * by default the histogram min and max corresponds to the largest and smallest * possible pixel value of that pixel type. To customize the histogram bounds * for a given image, the max and min pixel values that will be placed in the * histogram can be set manually. NB: The min and max are INCLUSIVE. * * Further, the type of histogram frequency container used is an optional * template parameter. By default, a dense container is used, but for images * with little texture or in cases where the user wants more histogram bins, * a sparse container can be used for the histogram instead. * * WARNING: This probably won't work for pixels of double or long-double type * unless you set the histogram min and max manually. This is because the largest * histogram bin by default has max value of the largest possible pixel value * plus 1. For double and long-double types, whose "RealType" as defined by the * NumericTraits class is the same, and thus cannot hold any larger values, * this would cause a float overflow. * * IJ article: http://hdl.handle.net/1926/1374 * * \sa ScalarImageToSizeZoneFeaturesFilter * \sa EnhancedScalarImageToSizeZoneMatrixFilter * \sa HistogramToSizeZoneFeaturesFilter * * \author: Nick Tustison * \ingroup ITKStatistics */ template class EnhancedScalarImageToSizeZoneMatrixFilter : public ProcessObject { public: /** Standard typedefs */ typedef EnhancedScalarImageToSizeZoneMatrixFilter Self; typedef ProcessObject Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; /** Run-time type information (and related methods). */ itkTypeMacro( EnhancedScalarImageToSizeZoneMatrixFilter, ProcessObject ); /** standard New() method support */ itkNewMacro( Self ); typedef TImageType ImageType; typedef typename ImageType::Pointer ImagePointer; typedef typename ImageType::ConstPointer ImageConstPointer; typedef typename ImageType::PixelType PixelType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::RegionType RegionType; typedef typename ImageType::SizeType RadiusType; typedef typename ImageType::OffsetType OffsetType; typedef VectorContainer OffsetVector; typedef typename OffsetVector::Pointer OffsetVectorPointer; typedef typename ImageType::PointType PointType; typedef typename NumericTraits::RealType MeasurementType; typedef typename NumericTraits::RealType RealType; typedef Histogram HistogramType; typedef typename HistogramType::Pointer HistogramPointer; typedef typename HistogramType::ConstPointer HistogramConstPointer; typedef typename HistogramType::MeasurementVectorType MeasurementVectorType; /** ImageDimension constants */ itkStaticConstMacro( ImageDimension, unsigned int, TImageType::ImageDimension ); /** Specify the default number of bins per axis */ itkStaticConstMacro( DefaultBinsPerAxis, unsigned int, 256 ); /** * Set the offsets over which the intensity/distance pairs will be computed. * Invoking this function clears the previous offsets. * Note: for each individual offset in the OffsetVector, the rightmost non-zero * offset element must be positive. For example, in the offset list of a 2D image, * (1, 0) means the offset along x-axis. (1, 0) has to be set instead * of (-1, 0). This is required from the iterating order of pixel iterator. * */ itkSetObjectMacro( Offsets, OffsetVector ); /** * Set offset over which the intensity/distance pairs will be computed. * Invoking this function clears the previous offset(s). * Note: for each individual offset, the rightmost non-zero * offset element must be positive. For example, in the offset list of a 2D image, * (1, 0) means the offset along x-axis. (1, 0) has to be set instead * of (-1, 0). This is required from the iterating order of pixel iterator. * */ void SetOffset( const OffsetType offset ); /** * Get the current offset(s). */ itkGetModifiableObjectMacro(Offsets, OffsetVector ); /** Set number of histogram bins along each axis */ itkSetMacro( NumberOfBinsPerAxis, unsigned int ); /** Get number of histogram bins along each axis */ itkGetConstMacro( NumberOfBinsPerAxis, unsigned int ); /** * Set the min and max (inclusive) pixel value that will be used in * generating the histogram. */ void SetPixelValueMinMax( PixelType min, PixelType max ); /** Get the min pixel value defining one dimension of the joint histogram. */ itkGetConstMacro( Min, PixelType ); /** Get the max pixel value defining one dimension of the joint histogram. */ itkGetConstMacro( Max, PixelType ); /** * Set the min and max (inclusive) pixel value that will be used in * generating the histogram. */ void SetDistanceValueMinMax( RealType min, RealType max ); /** * Get the min distance value defining one dimension of the joint histogram. */ itkGetConstMacro( MinDistance, RealType ); /** * Get the max distance value defining one dimension of the joint histogram. */ itkGetConstMacro( MaxDistance, RealType ); /** Method to set the input image */ using Superclass::SetInput; void SetInput( const ImageType *image ); /** Method to get the input image */ const ImageType * GetInput() const; /** Method to set the mask image */ void SetMaskImage( const ImageType *image ); /** Method to get the mask image */ const ImageType * GetMaskImage() const; /** method to get the Histogram */ const HistogramType * GetOutput() const; /** * Set the pixel value of the mask that should be considered "inside" the * object. Defaults to 1. */ itkSetMacro( InsidePixelValue, PixelType ); itkGetConstMacro( InsidePixelValue, PixelType ); protected: EnhancedScalarImageToSizeZoneMatrixFilter(); virtual ~EnhancedScalarImageToSizeZoneMatrixFilter() {}; virtual void PrintSelf( std::ostream& os, Indent indent ) const ITK_OVERRIDE; /** Standard itk::ProcessObject subclass method. */ typedef DataObject::Pointer DataObjectPointer; typedef ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType; using Superclass::MakeOutput; virtual DataObjectPointer MakeOutput( DataObjectPointerArraySizeType idx ) ITK_OVERRIDE; /** This method causes the filter to generate its output. */ virtual void GenerateData() ITK_OVERRIDE; - /** - * Normalize the direction of the offset before it is applied. - * The last non-zero dimension of the offest has to be positive in order - * to match to scanning order of the iterator. Only the sign is changed. - * For example, the input offset (-1, 0) will be normalized as - * (1, 0). - * */ - void NormalizeOffsetDirection(OffsetType &offset); - private: unsigned int m_NumberOfBinsPerAxis; PixelType m_Min; PixelType m_Max; RealType m_MinDistance; RealType m_MaxDistance; PixelType m_InsidePixelValue; MeasurementVectorType m_LowerBound; MeasurementVectorType m_UpperBound; OffsetVectorPointer m_Offsets; }; } // end of namespace Statistics } // end of namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkEnhancedScalarImageToSizeZoneMatrixFilter.hxx" #endif #endif diff --git a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.hxx b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.hxx index 8d234c374b..1dcd78b09f 100644 --- a/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.hxx +++ b/Modules/Classification/CLUtilities/include/itkEnhancedScalarImageToSizeZoneMatrixFilter.hxx @@ -1,431 +1,409 @@ /*=================================================================== 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 __itkEnhancedScalarImageToSizeZoneMatrixFilter_hxx #define __itkEnhancedScalarImageToSizeZoneMatrixFilter_hxx #include "itkEnhancedScalarImageToSizeZoneMatrixFilter.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 EnhancedScalarImageToSizeZoneMatrixFilter ::EnhancedScalarImageToSizeZoneMatrixFilter() : 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 = 2; this->ProcessObject::SetNthOutput( 0, this->MakeOutput( 0 ) ); HistogramType *output = const_cast( this->GetOutput() ); output->SetMeasurementVectorSize( measurementVectorSize ); 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 EnhancedScalarImageToSizeZoneMatrixFilter ::SetOffset( const OffsetType offset ) { OffsetVectorPointer offsetVector = OffsetVector::New(); offsetVector->push_back( offset ); this->SetOffsets( offsetVector ); } template void EnhancedScalarImageToSizeZoneMatrixFilter ::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 EnhancedScalarImageToSizeZoneMatrixFilter ::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 * EnhancedScalarImageToSizeZoneMatrixFilter ::GetInput() const { if( this->GetNumberOfInputs() < 1 ) { return ITK_NULLPTR; } return static_cast( this->ProcessObject::GetInput( 0 ) ); } template const TImageType * EnhancedScalarImageToSizeZoneMatrixFilter ::GetMaskImage() const { if( this->GetNumberOfInputs() < 2 ) { return ITK_NULLPTR; } return static_cast( this->ProcessObject::GetInput( 1 ) ); } template const typename EnhancedScalarImageToSizeZoneMatrixFilter::HistogramType * EnhancedScalarImageToSizeZoneMatrixFilter ::GetOutput() const { const HistogramType *output = static_cast( this->ProcessObject::GetOutput( 0 ) ); return output; } template typename EnhancedScalarImageToSizeZoneMatrixFilter::DataObjectPointer EnhancedScalarImageToSizeZoneMatrixFilter ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed( idx ) ) { return HistogramType::New().GetPointer(); } template void EnhancedScalarImageToSizeZoneMatrixFilter ::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(); - MITK_WARN << "InputImage casted."; + //MITK_WARN << "InputImage casted."; //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(); - MITK_WARN << "MaskImage casted."; + //MITK_WARN << "MaskImage casted."; //Set all values out of the mask to (m_Min + m_Max) / 2. typedef itk::MaskImageFilter< FloatImageType, LabelImageType, FloatImageType > MaskFilterType; typename MaskFilterType::Pointer maskFilter = MaskFilterType::New(); maskFilter->SetInput(floatImage); maskFilter->SetMaskImage(maskCaster->GetOutput()); maskFilter->SetOutsideValue((m_Max + m_Min) / 2); maskFilter->Update(); - MITK_WARN << "InputImage masked."; + //MITK_WARN << "InputImage masked."; //Rescale intensity to match the size of the histogram typedef itk::Image< unsigned int, 3 > OutputImageType; typedef itk::RescaleIntensityImageFilter< FloatImageType,OutputImageType> RescalerType; typename RescalerType::Pointer rescaler = RescalerType::New(); //We use 0 for nans, all valid numbers will be 1 < x < size rescaler->SetOutputMinimum( 1 ); rescaler->SetOutputMaximum( size[0] ); rescaler->SetInput(maskFilter->GetOutput()); rescaler->Update(); typename OutputImageType::Pointer rescaled = rescaler->GetOutput(); - MITK_WARN << "Intensities rescaled."; + //MITK_WARN << "Intensities rescaled."; //Write back the nans because they get lost during rescaling int xx = inputImage->GetLargestPossibleRegion().GetSize()[0]; int yy = inputImage->GetLargestPossibleRegion().GetSize()[1]; int zz = inputImage->GetLargestPossibleRegion().GetSize()[2]; for (int x = 0; x < xx; x++) { for (int y = 0; y < yy; y++) { for (int z = 0; z < zz; z++) { FloatImageType::IndexType indexF; indexF[0] = x; indexF[1] = y; indexF[2] = z; OutputImageType::IndexType indexO; indexO[0] = x; indexO[1] = y; indexO[2] = z; //Is Pixel NaN? if(floatImage->GetPixel(indexF) != floatImage->GetPixel(indexF)) { - rescaled->SetPixel(indexO,-1); + rescaled->SetPixel(indexO,0); } } } } + //All nans are now 0, the valid values are within [1,numberOfBins] + /* OutputImageType::IndexType indexO; indexO[0] = 0; indexO[1] = 2; indexO[2] = 1; - MITK_WARN << "is -1: " << rescaled->GetPixel(indexO); + MITK_WARN << "is 0: " << rescaled->GetPixel(indexO); indexO[0] = 0; indexO[1] = 0; indexO[2] = 0; MITK_WARN << "is 1: " << rescaled->GetPixel(indexO); + */ + + PixelType distanceThreshold = 1 - mitk::eps; - PixelType distanceThreshold = 1; //Calculate the connected components - typedef itk::ScalarConnectedComponentImageFilter + typedef itk::ScalarConnectedComponentImageFilter ConnectedComponentImageFilterType; typename ConnectedComponentImageFilterType::Pointer connected = ConnectedComponentImageFilterType::New (); connected->SetInput(rescaled); connected->SetMaskImage(maskCaster->GetOutput()); connected->SetDistanceThreshold(distanceThreshold); connected->Update(); + /* + indexO[0] = 0; + indexO[1] = 2; + indexO[2] = 1; + MITK_WARN << "is 0: " << (connected->GetOutput())->GetPixel(indexO); + indexO[0] = 0; + indexO[1] = 0; + indexO[2] = 0; + MITK_WARN << "is 1: " << (connected->GetOutput())->GetPixel(indexO); + MITK_WARN << "Connected components calculated."; + */ //Relabel the components - typedef itk::RelabelComponentImageFilter RelabelFilterType; + typedef itk::RelabelComponentImageFilter RelabelFilterType; typename RelabelFilterType::Pointer relabel = RelabelFilterType::New(); typename RelabelFilterType::ObjectSizeType minSize = 1; relabel->SetInput(connected->GetOutput()); relabel->SetMinimumObjectSize(minSize); relabel->Update(); - MITK_WARN << "Components relabeled."; + //MITK_WARN << "Components relabeled."; //Get the stats of the componentes - typedef itk::LabelStatisticsImageFilter< FloatImageType, LabelImageType> LabelStatisticsImageFilterType; + typedef itk::LabelStatisticsImageFilter< FloatImageType, OutputImageType> LabelStatisticsImageFilterType; typename LabelStatisticsImageFilterType::Pointer labelStatisticsImageFilter = LabelStatisticsImageFilterType::New(); labelStatisticsImageFilter->SetLabelInput( relabel->GetOutput() ); labelStatisticsImageFilter->SetInput(floatImage); labelStatisticsImageFilter->UseHistogramsOn(); // needed to compute median labelStatisticsImageFilter->Update(); /* std::cout << "Number of labels: " << labelStatisticsImageFilter->GetNumberOfLabels() << std::endl; std::cout << std::endl; */ typedef typename LabelStatisticsImageFilterType::ValidLabelValuesContainerType ValidLabelValuesType; for(typename ValidLabelValuesType::const_iterator vIt = labelStatisticsImageFilter->GetValidLabelValues().begin(); vIt != labelStatisticsImageFilter->GetValidLabelValues().end(); ++vIt) { if ( labelStatisticsImageFilter->HasLabel(*vIt) ) { LabelPixelType labelValue = *vIt; -/* - MITK_INFO << "Label: " << *vIt; - - MITK_INFO << "\tmin: " - << labelStatisticsImageFilter->GetMinimum( labelValue ); - MITK_INFO << "\tmax: " - << labelStatisticsImageFilter->GetMaximum( labelValue ); - MITK_INFO << "\tmedian: " - << labelStatisticsImageFilter->GetMedian( labelValue ); - MITK_INFO << "\tmean: " - << labelStatisticsImageFilter->GetMean( labelValue ); - MITK_INFO << "\tsigma: " - << labelStatisticsImageFilter->GetSigma( labelValue ); - MITK_INFO << "\tvariance: " - << labelStatisticsImageFilter->GetVariance( labelValue ); - MITK_INFO << "\tsum: " - << labelStatisticsImageFilter->GetSum( labelValue ); - MITK_INFO << "\tcount: " - << labelStatisticsImageFilter->GetCount( labelValue ); - MITK_INFO << "\tregion: " - << labelStatisticsImageFilter->GetRegion( labelValue ); - -*/ - run[0] = (labelStatisticsImageFilter->GetMinimum( labelValue ) + labelStatisticsImageFilter->GetMinimum( labelValue )) /2; + + run[0] = labelStatisticsImageFilter->GetMean( labelValue ); run[1] = labelStatisticsImageFilter->GetCount( labelValue ); //Check for NaN and inf if(run[0] == run[0] && !std::isinf(std::abs(run[0]))) { output->GetIndex( run, hIndex ); output->IncreaseFrequencyOfIndex( hIndex, 1 ); + + /* + MITK_INFO << "Adding a region:"; + MITK_INFO << "\tmin: " + << labelStatisticsImageFilter->GetMinimum( labelValue ); + MITK_INFO << "\tmax: " + << labelStatisticsImageFilter->GetMaximum( labelValue ); + MITK_INFO << "\tmean: " + << labelStatisticsImageFilter->GetMean( labelValue ); + MITK_INFO << "\tcount: " + << labelStatisticsImageFilter->GetCount( labelValue ); + */ } } } } template void EnhancedScalarImageToSizeZoneMatrixFilter ::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 EnhancedScalarImageToSizeZoneMatrixFilter ::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 EnhancedScalarImageToSizeZoneMatrixFilter ::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 -EnhancedScalarImageToSizeZoneMatrixFilter -::NormalizeOffsetDirection(OffsetType &offset) -{ - MITK_WARN <<" -> Size zone 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 << " -> size zone new offset = " << offset; - itkDebugMacro("new offset = " << offset << std::endl); -} } // end of namespace Statistics } // end of namespace itk #endif diff --git a/Modules/Classification/CLUtilities/include/mitkGIFNeighbourhoodGreyLevelDifference.h b/Modules/Classification/CLUtilities/include/mitkGIFNeighbourhoodGreyLevelDifference.h new file mode 100644 index 0000000000..d730ee5c30 --- /dev/null +++ b/Modules/Classification/CLUtilities/include/mitkGIFNeighbourhoodGreyLevelDifference.h @@ -0,0 +1,51 @@ +#ifndef mitkGIFNeighbourhoodGreyLevelDifference_h +#define mitkGIFNeighbourhoodGreyLevelDifference_h + +#include +#include +#include + +namespace mitk +{ + class MITKCLUTILITIES_EXPORT GIFNeighbourhoodGreyLevelDifference : public AbstractGlobalImageFeature + { + public: + mitkClassMacro(GIFNeighbourhoodGreyLevelDifference,AbstractGlobalImageFeature) + itkFactorylessNewMacro(Self) + itkCloneMacro(Self) + + GIFNeighbourhoodGreyLevelDifference(); + + /** + * \brief Calculates the Cooccurence-Matrix based features for this class. + */ + virtual FeatureListType CalculateFeatures(const Image::Pointer & image, const Image::Pointer &feature) override; + + /** + * \brief Returns a list of the names of all features that are calculated from this class + */ + virtual FeatureNameListType GetFeatureNames() override; + + itkGetConstMacro(Range,double); + itkSetMacro(Range, double); + + itkGetConstMacro(UseCtRange, bool); + itkSetMacro(UseCtRange, bool); + + itkGetConstMacro(Direction, unsigned int); + itkSetMacro(Direction, unsigned int); + + struct ParameterStruct + { + bool m_UseCtRange; + double m_Range; + unsigned int m_Direction; + }; + + private: + double m_Range; + bool m_UseCtRange; + unsigned int m_Direction; + }; +} +#endif //mitkGIFNeighbourhoodGreyLevelDifference_h diff --git a/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp new file mode 100644 index 0000000000..20ba05e44c --- /dev/null +++ b/Modules/Classification/CLUtilities/src/GlobalImageFeatures/mitkGIFNeighbourhoodGreyLevelDifference.cpp @@ -0,0 +1,175 @@ +/*=================================================================== + +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; +}