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