Medical Imaging Interaction Toolkit  2018.4.99-87d68d9f
Medical Imaging Interaction Toolkit
itkCoocurenceMatrixFeatureFunctor.h
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 #ifndef itkCooccurenceMatrixFeatureFunctor_h
14 #define itkCooccurenceMatrixFeatureFunctor_h
15 
16 #include "itkConstNeighborhoodIterator.h"
17 
18 #include <itkHistogramToTextureFeaturesFilter.h>
19 #include <itkHistogram.h>
20 #include <bitset>
21 
22 /*
23  * To Do:
24  * - Enable normalization of GLCM
25  */
26 
27 namespace itk
28 {
29 
30 namespace Functor
31 {
32 
36 template< typename TNeighborhoodType, typename TPixelOutputType>
38 {
39 public:
40  // Functor defines
41 
42  typedef itk::Statistics::Histogram< double, itk::Statistics::DenseFrequencyContainer2 > HistogramType;
43  typedef typename TNeighborhoodType::OffsetType OffsetType;
44  typedef typename HistogramType::SizeType SizeType;
45  typedef typename HistogramType::MeasurementVectorType MeasurementVectorType;
46  typedef typename itk::Statistics::HistogramToTextureFeaturesFilter<HistogramType> HistoToFeatureFilter;
47 
48  static const unsigned int OutputCount = 8;
49  typedef vnl_vector_fixed<TPixelOutputType, OutputCount> OutputVectorType;
50  typedef TNeighborhoodType NeighborhoodType;
52  {
54  };
55 
56  static const char * GetFeatureName(unsigned int f )
57  {
58  static const char * const FeatureNames[] = {"ENERGY", "ENTROPY", "CORRELATION", "INERTIA", "CLUSTERSHADE", "CLUSTERPROMINENCE", "HARALICKCORRELATION", "INVERSEDIFFERENCEMOMENT"};
59  return FeatureNames[f];
60  }
61 
63  {
71  DIR_n1_x1_x1 = 128,
72  DIR_x1_x0_n1 = 256,
73  DIR_x1_x1_n1 = 512,
74  DIR_x0_x1_n1 = 1024,
75  DIR_n1_x1_n1 = 2048,
76  DIR_x0_x0_x1 = 4096,
78  };
79 
80  static std::map<OffsetDirection, OffsetType > CreateStdOffsets()
81  {
82  std::map<OffsetDirection,OffsetType> offsetMap;
83 
84  {
85  itk::Offset<3> off = {{1,0,0}};
86  offsetMap[DIR_x1_x0_x0] = off;
87  }
88  {
89  itk::Offset<3> off = {{1,1,0}};
90  offsetMap[DIR_x1_x1_x0] = off;
91  }
92  {
93  itk::Offset<3> off = {{0,1,0}};
94  offsetMap[DIR_x0_x1_x0] = off;
95  }
96  {
97  itk::Offset<3> off = {{-1,1,0}};
98  offsetMap[DIR_n1_x1_x0] = off;
99  }
100  {
101  itk::Offset<3> off = {{1,0,1}};
102  offsetMap[DIR_x1_x0_x1]= off;
103  }
104  {
105  itk::Offset<3> off = {{1,1,1}};
106  offsetMap[DIR_x1_x1_x1] = off;
107  }
108  {
109  itk::Offset<3> off = {{0,1,1}};
110  offsetMap[DIR_x0_x1_x1] = off;
111  }
112  {
113  itk::Offset<3> off = {{-1,1,1}};
114  offsetMap[DIR_n1_x1_x1] = off;
115  }
116  {
117  itk::Offset<3> off = {{1,0,-1}};
118  offsetMap[DIR_x1_x0_n1] = off;
119  }
120  {
121  itk::Offset<3> off = {{1,1,-1}};
122  offsetMap[DIR_x1_x1_n1] = off;
123  }
124  {
125  itk::Offset<3> off = {{0,1,-1}};
126  offsetMap[DIR_x0_x1_n1] = off;
127  }
128  {
129  itk::Offset<3> off = {{-1,1,-1}};
130  offsetMap[DIR_n1_x1_n1] = off;
131  }
132  {
133  itk::Offset<3> off = {{0,0,1}};
134  offsetMap[DIR_x0_x0_x1] = off;
135  }
136  return offsetMap;
137  }
138 
140  std::map<OffsetDirection, OffsetType> m_offsetMap;
141  unsigned int m_levels;
142 
143  void DirectionFlags(int flags)
144  {
145  DirectionFlags(static_cast<OffsetDirection>(flags));
146  }
147 
149  {
150  m_DirectionFlags = flags;
151  }
152 
153  void SetLevels(unsigned int lvl)
154  {
155  m_levels = lvl;
156  }
157 
159  {
160  m_offsetMap = CreateStdOffsets();
161  m_levels = 5;
162  m_DirectionFlags =
163  DIR_x1_x0_x0; /*| DIR_x1_x1_x0 | DIR_x0_x1_x0 |
164  DIR_n1_x1_x0 | DIR_x1_x0_x1 | DIR_x1_x1_x1 |
165  DIR_x0_x1_x1 | DIR_n1_x1_x1 | DIR_x1_x0_n1 |
166  DIR_x1_x1_n1 | DIR_x0_x1_n1 | DIR_n1_x1_n1 |
167  DIR_x0_x0_x1;*/
168 
169  std::cout << "NeighborhoodCooccurenceMatrix" << std::endl;
170  }
171 
172  inline OutputVectorType operator()(const TNeighborhoodType & it) const
173  {
176  for (unsigned int i = 0; i < it.Size(); ++i)
177  {
178  double value = it.GetPixel(i);
179  max = (value > max) ? value : max;
180  min = (value < min) ? value : min;
181  }
182  if ( min >= max)
183  {
184  OutputVectorType nullRes; nullRes.fill(0);
185  return nullRes;
186  }
187 
188  SizeType size;
189  MeasurementVectorType minBorder;
190  MeasurementVectorType maxBorder;
191  OffsetType g1, g2;
192 
193  size.SetSize(2); // 2D Historgram
194  size.Fill(m_levels); // 5 bins each dim
195 
196  minBorder.SetSize(2); // min range value
197  minBorder.Fill(min);
198 
199  maxBorder.SetSize(2); // max range value
200  maxBorder.Fill(max);
201 
202  MeasurementVectorType cooccur;
203  cooccur.SetSize(2);
204 
205  OutputVectorType output_vector;
206  output_vector.fill(0);
207 
208  double div_num_dirs = 0;
209 
210  //std::bitset<14> x(m_DirectionFlags);
211  //std::cout << "Direction flag " << x;
212 
213  for(typename std::map<OffsetDirection,OffsetType>::const_iterator dir_it = m_offsetMap.begin(),
214  end = m_offsetMap.end(); dir_it != end; dir_it ++){
215 
216  if(! (dir_it->first & m_DirectionFlags))
217  continue;
218 
219  div_num_dirs++;
220 
221  HistogramType::Pointer histogram = HistogramType::New();
222  histogram->SetMeasurementVectorSize(2);
223  histogram->Initialize(size, minBorder, maxBorder);
224 
225  for (unsigned int i = 0; i < it.Size(); ++i)
226  {
227  // grayvalue pair (g1,g2)
228  // g1 ~ g2 with ~ as relation (e.g. a offset calculation)
229  g1 = it.GetOffset(i);
230  g2 = g1 + dir_it->second;
231 
232  cooccur[0] = it.GetPixel(i);
233  cooccur[1] = it.GetImagePointer()->GetPixel(it.GetIndex(i)+dir_it->second);
234 
235  histogram->IncreaseFrequencyOfMeasurement(cooccur, 1);
236 
237  std::swap(cooccur[0],cooccur[1]);
238  histogram->IncreaseFrequencyOfMeasurement(cooccur, 1);
239  }
240 
241  //To do Normalize GLCM N_g Number of levels
242 
243  HistoToFeatureFilter::Pointer filter = HistoToFeatureFilter::New();
244  filter->SetInput(histogram);
245  filter->Update();
246 
247  output_vector[ENERGY] += filter->GetEnergy();
248  output_vector[ENTROPY] += filter->GetEntropy();
249  output_vector[CORRELATION] += filter->GetCorrelation();
250  output_vector[INERTIA] += filter->GetInertia();
251  output_vector[CLUSTERSHADE] += filter->GetClusterShade();
252  output_vector[CLUSTERPROMINENCE] += filter->GetClusterProminence();
253  output_vector[HARALICKCORRELATION] += filter->GetHaralickCorrelation();
254  output_vector[INVERSEDIFFERENCEMOMENT] += filter->GetInverseDifferenceMoment();
255  }
256 
257  //std::cout << "Number of directions " << div_num_dirs << std::endl;
258 // output_vector /= div_num_dirs;
259  return output_vector;
260  }
261 };
262 
263 }// end namespace functor
264 
265 } // end namespace itk
266 
267 #endif // itkNeighborhoodFunctors_h
OutputVectorType operator()(const TNeighborhoodType &it) const
itk::Statistics::HistogramToTextureFeaturesFilter< HistogramType > HistoToFeatureFilter
HistogramType::MeasurementVectorType MeasurementVectorType
vnl_vector_fixed< TPixelOutputType, OutputCount > OutputVectorType
static T max(T x, T y)
Definition: svm.cpp:56
static T min(T x, T y)
Definition: svm.cpp:53
Functor for texture feature calculation based on the Cooccurence matrix.
static void swap(T &x, T &y)
Definition: svm.cpp:58
static std::map< OffsetDirection, OffsetType > CreateStdOffsets()
itk::Statistics::Histogram< double, itk::Statistics::DenseFrequencyContainer2 > HistogramType