|
| 1 | +/*========================================================================= |
| 2 | + * |
| 3 | + * Copyright NumFOCUS |
| 4 | + * |
| 5 | + * Licensed under the Apache License, Version 2.0 (the "License"); |
| 6 | + * you may not use this file except in compliance with the License. |
| 7 | + * You may obtain a copy of the License at |
| 8 | + * |
| 9 | + * https://www.apache.org/licenses/LICENSE-2.0.txt |
| 10 | + * |
| 11 | + * Unless required by applicable law or agreed to in writing, software |
| 12 | + * distributed under the License is distributed on an "AS IS" BASIS, |
| 13 | + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 14 | + * See the License for the specific language governing permissions and |
| 15 | + * limitations under the License. |
| 16 | + * |
| 17 | + *=========================================================================*/ |
| 18 | +#ifndef itkCoocurrenceTextureFeaturesImageFilter_h |
| 19 | +#define itkCoocurrenceTextureFeaturesImageFilter_h |
| 20 | + |
| 21 | +#include "itkImageToImageFilter.h" |
| 22 | +#include "itkConstNeighborhoodIterator.h" |
| 23 | +#include "itkVectorContainer.h" |
| 24 | + |
| 25 | +namespace itk |
| 26 | +{ |
| 27 | +namespace Statistics |
| 28 | +{ |
| 29 | +/** \class CoocurrenceTextureFeaturesImageFilter |
| 30 | + * \brief This class computes texture descriptions for each voxel of |
| 31 | + * a given image and a mask image if provided. The output image can then be |
| 32 | + * displayed by using colormaps. |
| 33 | + * |
| 34 | + * This filter computes a N-D image where each voxel will contain |
| 35 | + * a vector of up to 8 scalars representing the texture features |
| 36 | + * (of the specified neighborhood) from a N-D scalar image. |
| 37 | + * The texture features are computed for each spatial |
| 38 | + * direction and averaged afterward. |
| 39 | + * |
| 40 | + * This filter use grey-level co-occurence matrix in order to compute description a la Haralick. (See |
| 41 | + * Haralick, R.M., K. Shanmugam and I. Dinstein. 1973. Textural Features for |
| 42 | + * Image Classification. IEEE Transactions on Systems, Man and Cybernetics. |
| 43 | + * SMC-3(6):610-620. See also Haralick, R.M. 1979. Statistical and Structural |
| 44 | + * Approaches to Texture. Proceedings of the IEEE, 67:786-804.) |
| 45 | + * |
| 46 | + * Template Parameters: |
| 47 | + * -# The input image type: a N dimensional image where the pixel type MUST be integer. |
| 48 | + * -# The output image type: a N dimensional image where the pixel type MUST be a vector of floating points or a |
| 49 | + *VectorImage. |
| 50 | + * |
| 51 | + * Inputs and parameters: |
| 52 | + * -# An image |
| 53 | + * -# A mask defining the region over which texture features will be |
| 54 | + * calculated. (Optional) |
| 55 | + * -# The pixel value that defines the "inside" of the mask. (Optional, defaults |
| 56 | + * to 1 if a mask is set.) |
| 57 | + * -# The number of intensity bins. (Optional, defaults to 256.) |
| 58 | + * -# The set of directions (offsets) to average across. (Optional, defaults to |
| 59 | + * {(-1, 0), (-1, -1), (0, -1), (1, -1)} for 2D images and scales analogously |
| 60 | + * for ND images.) |
| 61 | + * -# The pixel intensity range over which the features will be calculated. |
| 62 | + * (Optional, defaults to the full dynamic range of the pixel type.) |
| 63 | + * -# The size of the neighborhood radius. (Optional, defaults to 2.) |
| 64 | + * |
| 65 | + * Recommendations: |
| 66 | + * -# Input image: To improve the computation time, the useful data should take as much |
| 67 | + * space as possible in the input image. If the useful data is concentrated in one part of |
| 68 | + * the image a crop step should be considered prior to the usage of this filter. |
| 69 | + * -# Mask: Even if optional, the usage of a mask will greatly improve the computation time. |
| 70 | + * -# Number of intensity bins: The number of bins should be adapted to the type of results expected |
| 71 | + * and the intensity range. In addition a high number of bins will increase the |
| 72 | + * computation time. |
| 73 | + * -# Pixel intensity range: For better results the Pixel intensity should be adapted to the input image |
| 74 | + * intensity range. For example they could be the minimum and maximum intensity of the image, or 0 and |
| 75 | + * the maximum intensity (if the negative values are considered as noise). |
| 76 | + * |
| 77 | + * WARNING: This probably won't work for pixels of double or long-double type |
| 78 | + * unless you set the histogram min and max manually. This is because the largest |
| 79 | + * histogram bin by default has max value of the largest possible pixel value |
| 80 | + * plus 1. For double and long-double types, whose "RealType" as defined by the |
| 81 | + * NumericTraits class is the same, and thus cannot hold any larger values, |
| 82 | + * this would cause a float overflow. |
| 83 | + * |
| 84 | + * \sa HistogramToTextureFeaturesFilter |
| 85 | + * \sa ScalarImageToCooccurrenceMatrixFilte |
| 86 | + * \sa ScalarImageToTextureFeaturesFilter |
| 87 | + * |
| 88 | + * \author: Jean-Baptiste Vimort |
| 89 | + * \ingroup TextureFeatures |
| 90 | + **/ |
| 91 | + |
| 92 | +template <typename TInputImage, |
| 93 | + typename TOutputImage, |
| 94 | + typename TMaskImage = Image<unsigned char, TInputImage::ImageDimension>> |
| 95 | +class ITK_TEMPLATE_EXPORT CoocurrenceTextureFeaturesImageFilter : public ImageToImageFilter<TInputImage, TOutputImage> |
| 96 | +{ |
| 97 | +public: |
| 98 | + /** Standard type alias */ |
| 99 | + using Self = CoocurrenceTextureFeaturesImageFilter; |
| 100 | + using Superclass = ImageToImageFilter<TInputImage, TOutputImage>; |
| 101 | + using Pointer = SmartPointer<Self>; |
| 102 | + using ConstPointer = SmartPointer<const Self>; |
| 103 | + |
| 104 | + /** Run-time type information (and related methods). */ |
| 105 | + itkOverrideGetNameOfClassMacro(CoocurrenceTextureFeaturesImageFilter); |
| 106 | + |
| 107 | + /** standard New() method support */ |
| 108 | + itkNewMacro(Self); |
| 109 | + |
| 110 | + using InputImageType = TInputImage; |
| 111 | + using OutputImageType = TOutputImage; |
| 112 | + using MaskImageType = TMaskImage; |
| 113 | + |
| 114 | + using PixelType = typename InputImageType::PixelType; |
| 115 | + using MaskPixelType = typename MaskImageType::PixelType; |
| 116 | + using IndexType = typename InputImageType::IndexType; |
| 117 | + using PointType = typename InputImageType::PointType; |
| 118 | + |
| 119 | + using OffsetType = typename InputImageType::OffsetType; |
| 120 | + using OffsetVector = VectorContainer<unsigned char, OffsetType>; |
| 121 | + using OffsetVectorPointer = typename OffsetVector::Pointer; |
| 122 | + using OffsetVectorConstPointer = typename OffsetVector::ConstPointer; |
| 123 | + |
| 124 | + using InputRegionType = typename InputImageType::RegionType; |
| 125 | + using OutputRegionType = typename OutputImageType::RegionType; |
| 126 | + |
| 127 | + using NeighborhoodRadiusType = typename itk::ConstNeighborhoodIterator<InputImageType>::RadiusType; |
| 128 | + |
| 129 | + using MeasurementType = typename NumericTraits<PixelType>::RealType; |
| 130 | + using RealType = typename NumericTraits<PixelType>::RealType; |
| 131 | + |
| 132 | + /** Method to set/get the Neighborhood radius */ |
| 133 | + itkSetMacro(NeighborhoodRadius, NeighborhoodRadiusType); |
| 134 | + itkGetConstMacro(NeighborhoodRadius, NeighborhoodRadiusType); |
| 135 | + |
| 136 | + /** Method to set the mask image */ |
| 137 | + itkSetInputMacro(MaskImage, MaskImageType); |
| 138 | + |
| 139 | + /** Method to get the mask image */ |
| 140 | + itkGetInputMacro(MaskImage, MaskImageType); |
| 141 | + |
| 142 | + |
| 143 | + /** Specify the default number of bins per axis */ |
| 144 | + static constexpr unsigned int DefaultBinsPerAxis = 256; |
| 145 | + |
| 146 | + /** |
| 147 | + * Set the offsets over which the intensities pairs will be computed. |
| 148 | + * Invoking this function clears the previous offsets. |
| 149 | + * Note: for each individual offset in the OffsetVector, the rightmost non-zero |
| 150 | + * offset element must be positive. For example, in the offset list of a 2D image, |
| 151 | + * (1, 0) means the offset along x-axis. (1, 0) has to be set instead |
| 152 | + * of (-1, 0). This is required from the iterating order of pixel iterator. |
| 153 | + * |
| 154 | + */ |
| 155 | + itkSetObjectMacro(Offsets, OffsetVector); |
| 156 | + |
| 157 | + /** |
| 158 | + * Set offset over which the intensities pairs will be computed. |
| 159 | + * Invoking this function clears the previous offset(s). |
| 160 | + * Note: for each individual offset, the rightmost non-zero |
| 161 | + * offset element must be positive. For example, in the offset list of a 2D image, |
| 162 | + * (1, 0) means the offset along x-axis. (1, 0) has to be set instead |
| 163 | + * of (-1, 0). This is required from the iterating order of pixel iterator. |
| 164 | + * |
| 165 | + */ |
| 166 | + void |
| 167 | + SetOffset(const OffsetType offset); |
| 168 | + |
| 169 | + /** |
| 170 | + * Get the current offset(s). |
| 171 | + */ |
| 172 | + itkGetModifiableObjectMacro(Offsets, OffsetVector); |
| 173 | + |
| 174 | + /** Set number of histogram bins along each axis */ |
| 175 | + itkSetMacro(NumberOfBinsPerAxis, unsigned int); |
| 176 | + |
| 177 | + /** Get number of histogram bins along each axis */ |
| 178 | + itkGetConstMacro(NumberOfBinsPerAxis, unsigned int); |
| 179 | + |
| 180 | + /** Get the max pixel value defining one dimension of the joint histogram. */ |
| 181 | + itkGetConstMacro(HistogramMaximum, PixelType); |
| 182 | + itkSetMacro(HistogramMaximum, PixelType); |
| 183 | + |
| 184 | + /** Get the min pixel value defining one dimension of the joint histogram. */ |
| 185 | + itkGetConstMacro(HistogramMinimum, PixelType); |
| 186 | + itkSetMacro(HistogramMinimum, PixelType); |
| 187 | + |
| 188 | + /** |
| 189 | + * Set the pixel value of the mask that should be considered "inside" the |
| 190 | + * object. Defaults to 1. |
| 191 | + */ |
| 192 | + itkSetMacro(InsidePixelValue, MaskPixelType); |
| 193 | + itkGetConstMacro(InsidePixelValue, MaskPixelType); |
| 194 | + |
| 195 | + /** Set the calculator to normalize the histogram (divide all bins by the |
| 196 | + total frequency). Normalization is off by default. */ |
| 197 | + itkSetMacro(Normalize, bool); |
| 198 | + itkGetConstMacro(Normalize, bool); |
| 199 | + itkBooleanMacro(Normalize); |
| 200 | + |
| 201 | + using OutputPixelType = typename OutputImageType::PixelType; |
| 202 | + using OutputRealType = typename NumericTraits<OutputPixelType>::ScalarRealType; |
| 203 | + |
| 204 | +#ifdef ITK_USE_CONCEPT_CHECKING |
| 205 | + // Begin concept checking |
| 206 | + itkConceptMacro(OutputPixelTypeCheck, (Concept::IsFloatingPoint<OutputRealType>)); |
| 207 | + // End concept checking |
| 208 | +#endif |
| 209 | + |
| 210 | +protected: |
| 211 | + using HistogramIndexType = int; |
| 212 | + using DigitizedImageType = itk::Image<HistogramIndexType, TInputImage::ImageDimension>; |
| 213 | + using NeighborhoodIteratorType = typename itk::ConstNeighborhoodIterator<DigitizedImageType>; |
| 214 | + using NeighborIndexType = typename NeighborhoodIteratorType::NeighborIndexType; |
| 215 | + |
| 216 | + CoocurrenceTextureFeaturesImageFilter(); |
| 217 | + ~CoocurrenceTextureFeaturesImageFilter() override = default; |
| 218 | + |
| 219 | + bool |
| 220 | + IsInsideNeighborhood(const OffsetType & iteratedOffset); |
| 221 | + void |
| 222 | + ComputeFeatures(const vnl_matrix<unsigned int> & hist, |
| 223 | + const unsigned int totalNumberOfFreq, |
| 224 | + typename TOutputImage::PixelType & outputPixel); |
| 225 | + void |
| 226 | + ComputeMeansAndVariances(const vnl_matrix<unsigned int> & hist, |
| 227 | + const unsigned int totalNumberOfFreq, |
| 228 | + double & pixelMean, |
| 229 | + double & marginalMean, |
| 230 | + double & marginalDevSquared, |
| 231 | + double & pixelVariance); |
| 232 | + void |
| 233 | + PrintSelf(std::ostream & os, Indent indent) const override; |
| 234 | + |
| 235 | + /** This method causes the filter to generate its output. */ |
| 236 | + void |
| 237 | + BeforeThreadedGenerateData() override; |
| 238 | + void |
| 239 | + AfterThreadedGenerateData() override; |
| 240 | + void |
| 241 | + DynamicThreadedGenerateData(const OutputRegionType & outputRegionForThread) override; |
| 242 | + void |
| 243 | + GenerateOutputInformation() override; |
| 244 | + |
| 245 | +private: |
| 246 | + typename DigitizedImageType::Pointer m_DigitizedInputImage; |
| 247 | + |
| 248 | + NeighborhoodRadiusType m_NeighborhoodRadius; |
| 249 | + OffsetVectorPointer m_Offsets; |
| 250 | + unsigned int m_NumberOfBinsPerAxis; |
| 251 | + PixelType m_HistogramMinimum; |
| 252 | + PixelType m_HistogramMaximum; |
| 253 | + MaskPixelType m_InsidePixelValue; |
| 254 | + bool m_Normalize; |
| 255 | +}; |
| 256 | +} // end of namespace Statistics |
| 257 | +} // end of namespace itk |
| 258 | + |
| 259 | +#ifndef ITK_MANUAL_INSTANTIATION |
| 260 | +# include "itkCoocurrenceTextureFeaturesImageFilter.hxx" |
| 261 | +#endif |
| 262 | + |
| 263 | +#endif |
0 commit comments