Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkGIFNeighbouringGreyLevelDependenceFeatures.cpp
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 
14 
15 // MITK
16 #include <mitkITKImageImport.h>
17 #include <mitkImageCast.h>
18 #include <mitkImageAccessByItk.h>
19 
20 // ITK
22 #include <itkMinimumMaximumImageCalculator.h>
23 #include <itkNeighborhoodIterator.h>
24 #include <itkImageRegionConstIterator.h>
25 
26 // STL
27 #include <sstream>
28 
29 namespace mitk
30 {
31  struct NGLDMMatrixHolder
32  {
33  public:
34  NGLDMMatrixHolder(double min, double max, int number, int depenence);
35 
36  int IntensityToIndex(double intensity);
37  double IndexToMinIntensity(int index);
38  double IndexToMeanIntensity(int index);
39  double IndexToMaxIntensity(int index);
40 
41  double m_MinimumRange;
42  double m_MaximumRange;
43  double m_Stepsize;
44  int m_NumberOfDependences;
45  int m_NumberOfBins;
46  Eigen::MatrixXd m_Matrix;
47 
48  int m_NeighbourhoodSize;
49  unsigned long m_NumberOfNeighbourVoxels;
50  unsigned long m_NumberOfDependenceNeighbourVoxels;
51  unsigned long m_NumberOfNeighbourhoods;
52  unsigned long m_NumberOfCompleteNeighbourhoods;
53  };
54 
55  struct NGLDMMatrixFeatures
56  {
57  NGLDMMatrixFeatures() :
58  LowDependenceEmphasis(0),
59  HighDependenceEmphasis(0),
60  LowGreyLevelCountEmphasis(0),
61  HighGreyLevelCountEmphasis(0),
62  LowDependenceLowGreyLevelEmphasis(0),
63  LowDependenceHighGreyLevelEmphasis(0),
64  HighDependenceLowGreyLevelEmphasis(0),
65  HighDependenceHighGreyLevelEmphasis(0),
66  GreyLevelNonUniformity(0),
67  GreyLevelNonUniformityNormalised(0),
68  DependenceCountNonUniformity(0),
69  DependenceCountNonUniformityNormalised(0),
70  DependenceCountPercentage(0),
71  GreyLevelVariance(0),
72  DependenceCountVariance(0),
73  DependenceCountEntropy(0),
74  DependenceCountEnergy(0),
75  MeanGreyLevelCount(0),
76  MeanDependenceCount(0),
77  ExpectedNeighbourhoodSize(0),
78  AverageNeighbourhoodSize(0),
79  AverageIncompleteNeighbourhoodSize(0),
80  PercentageOfCompleteNeighbourhoods(0),
81  PercentageOfDependenceNeighbours(0)
82  {
83  }
84 
85  public:
86  double LowDependenceEmphasis;
87  double HighDependenceEmphasis;
88  double LowGreyLevelCountEmphasis;
89  double HighGreyLevelCountEmphasis;
90  double LowDependenceLowGreyLevelEmphasis;
91  double LowDependenceHighGreyLevelEmphasis;
92  double HighDependenceLowGreyLevelEmphasis;
93  double HighDependenceHighGreyLevelEmphasis;
94 
95  double GreyLevelNonUniformity;
96  double GreyLevelNonUniformityNormalised;
97  double DependenceCountNonUniformity;
98  double DependenceCountNonUniformityNormalised;
99 
100  double DependenceCountPercentage;
101  double GreyLevelVariance;
102  double DependenceCountVariance;
103  double DependenceCountEntropy;
104  double DependenceCountEnergy;
105  double MeanGreyLevelCount;
106  double MeanDependenceCount;
107 
108  double ExpectedNeighbourhoodSize;
109  double AverageNeighbourhoodSize;
110  double AverageIncompleteNeighbourhoodSize;
111  double PercentageOfCompleteNeighbourhoods;
112  double PercentageOfDependenceNeighbours;
113 
114  };
115 }
116 
117 static
118 void MatrixFeaturesTo(mitk::NGLDMMatrixFeatures features,
119  std::string prefix,
121 
122 
123 mitk::NGLDMMatrixHolder::NGLDMMatrixHolder(double min, double max, int number, int depenence) :
124  m_MinimumRange(min),
125  m_MaximumRange(max),
126  m_Stepsize(0),
127  m_NumberOfDependences(depenence),
128  m_NumberOfBins(number),
129  m_NeighbourhoodSize(1),
130  m_NumberOfNeighbourVoxels(0),
131  m_NumberOfDependenceNeighbourVoxels(0),
132  m_NumberOfNeighbourhoods(0),
133  m_NumberOfCompleteNeighbourhoods(0)
134 {
135  m_Matrix.resize(number, depenence);
136  m_Matrix.fill(0);
137  m_Stepsize = (max - min) / (number);
138 }
139 
140 int mitk::NGLDMMatrixHolder::IntensityToIndex(double intensity)
141 {
142  return std::floor((intensity - m_MinimumRange) / m_Stepsize);
143 }
144 
145 double mitk::NGLDMMatrixHolder::IndexToMinIntensity(int index)
146 {
147  return m_MinimumRange + index * m_Stepsize;
148 }
149 double mitk::NGLDMMatrixHolder::IndexToMeanIntensity(int index)
150 {
151  return m_MinimumRange + (index+0.5) * m_Stepsize;
152 }
153 double mitk::NGLDMMatrixHolder::IndexToMaxIntensity(int index)
154 {
155  return m_MinimumRange + (index + 1) * m_Stepsize;
156 }
157 
158 template<typename TPixel, unsigned int VImageDimension>
159 void
160 CalculateNGLDMMatrix(itk::Image<TPixel, VImageDimension>* itkImage,
161  itk::Image<unsigned short, VImageDimension>* mask,
162  int alpha,
163  int range,
164  unsigned int direction,
165  mitk::NGLDMMatrixHolder &holder)
166 {
167  typedef itk::Image<TPixel, VImageDimension> ImageType;
168  typedef itk::Image<unsigned short, VImageDimension> MaskImageType;
169  typedef itk::NeighborhoodIterator<ImageType> ShapeIterType;
170  typedef itk::NeighborhoodIterator<MaskImageType> ShapeMaskIterType;
171 
172  holder.m_NumberOfCompleteNeighbourhoods = 0;
173  holder.m_NumberOfNeighbourhoods = 0;
174  holder.m_NumberOfNeighbourVoxels = 0;
175  holder.m_NumberOfDependenceNeighbourVoxels = 0;
176 
177  itk::Size<VImageDimension> radius;
178  radius.Fill(range);
179 
180  if ((direction > 1) && (direction - 2 <VImageDimension))
181  {
182  radius[direction - 2] = 0;
183  }
184 
185  ShapeIterType imageIter(radius, itkImage, itkImage->GetLargestPossibleRegion());
186  ShapeMaskIterType maskIter(radius, mask, mask->GetLargestPossibleRegion());
187 
188  auto region = mask->GetLargestPossibleRegion();
189 
190  auto center = imageIter.Size() / 2;
191  auto iterSize = imageIter.Size();
192  holder.m_NeighbourhoodSize = iterSize-1;
193  while (!maskIter.IsAtEnd())
194  {
195  int sameValues = 0;
196  bool completeNeighbourhood = true;
197 
198  int i = holder.IntensityToIndex(imageIter.GetCenterPixel());
199 
200  if ((imageIter.GetCenterPixel() != imageIter.GetCenterPixel()) ||
201  (maskIter.GetCenterPixel() < 1))
202  {
203  ++imageIter;
204  ++maskIter;
205  continue;
206  }
207 
208  for (unsigned int position = 0; position < iterSize; ++position)
209  {
210  if (position == center)
211  {
212  continue;
213  }
214  if ( ! region.IsInside(maskIter.GetIndex(position)))
215  {
216  completeNeighbourhood = false;
217  continue;
218  }
219  bool isInBounds;
220  auto jIntensity = imageIter.GetPixel(position, isInBounds);
221  auto jMask = maskIter.GetPixel(position, isInBounds);
222  if (jMask < 1 || (jIntensity != jIntensity) || ( ! isInBounds))
223  {
224  completeNeighbourhood = false;
225  continue;
226  }
227 
228  int j = holder.IntensityToIndex(jIntensity);
229  holder.m_NumberOfNeighbourVoxels += 1;
230  if (std::abs(i - j) <= alpha)
231  {
232  holder.m_NumberOfDependenceNeighbourVoxels += 1;
233  ++sameValues;
234  }
235  }
236  holder.m_Matrix(i, sameValues) += 1;
237  holder.m_NumberOfNeighbourhoods += 1;
238  if (completeNeighbourhood)
239  {
240  holder.m_NumberOfCompleteNeighbourhoods += 1;
241  }
242 
243  ++imageIter;
244  ++maskIter;
245  }
246 
247 }
248 
250  mitk::NGLDMMatrixHolder &holder,
251  mitk::NGLDMMatrixFeatures & results
252  )
253 {
254  auto sijMatrix = holder.m_Matrix;
255  auto piMatrix = holder.m_Matrix;
256  auto pjMatrix = holder.m_Matrix;
257  // double Ng = holder.m_NumberOfBins;
258  // int NgSize = holder.m_NumberOfBins;
259  double Ns = sijMatrix.sum();
260  piMatrix.rowwise().normalize();
261  pjMatrix.colwise().normalize();
262 
263  for (int i = 0; i < holder.m_NumberOfBins; ++i)
264  {
265  double sj = 0;
266  for (int j = 0; j < holder.m_NumberOfDependences; ++j)
267  {
268  double iInt = i+1 ;// holder.IndexToMeanIntensity(i);
269  double sij = sijMatrix(i, j);
270  double k = j + 1;
271  double pij = sij / Ns;
272 
273  results.LowDependenceEmphasis += sij / k / k;
274  results.HighDependenceEmphasis += sij * k*k;
275  if (iInt != 0)
276  {
277  results.LowGreyLevelCountEmphasis += sij / iInt / iInt;
278  }
279  results.HighGreyLevelCountEmphasis += sij * iInt*iInt;
280  if (iInt != 0)
281  {
282  results.LowDependenceLowGreyLevelEmphasis += sij / k / k / iInt / iInt;
283  }
284  results.LowDependenceHighGreyLevelEmphasis += sij * iInt*iInt / k / k;
285  if (iInt != 0)
286  {
287  results.HighDependenceLowGreyLevelEmphasis += sij *k * k / iInt / iInt;
288  }
289  results.HighDependenceHighGreyLevelEmphasis += sij * k*k*iInt*iInt;
290 
291 
292  results.MeanGreyLevelCount += iInt * pij;
293  results.MeanDependenceCount += k * pij;
294  if (pij > 0)
295  {
296  results.DependenceCountEntropy -= pij * std::log(pij) / std::log(2);
297  }
298  results.DependenceCountEnergy += pij*pij;
299  sj += sij;
300  }
301  results.GreyLevelNonUniformity += sj*sj;
302  results.GreyLevelNonUniformityNormalised += sj*sj;
303  }
304 
305  for (int j = 0; j < holder.m_NumberOfDependences; ++j)
306  {
307  double si = 0;
308  for (int i = 0; i < holder.m_NumberOfBins; ++i)
309  {
310  double sij = sijMatrix(i, j);
311  si += sij;
312  }
313  results.DependenceCountNonUniformity += si*si;
314  results.DependenceCountNonUniformityNormalised += si*si;
315  }
316  for (int i = 0; i < holder.m_NumberOfBins; ++i)
317  {
318  for (int j = 0; j < holder.m_NumberOfDependences; ++j)
319  {
320  double iInt = i + 1;// holder.IndexToMeanIntensity(i);
321  double sij = sijMatrix(i, j);
322  double k = j + 1;
323  double pij = sij / Ns;
324 
325  results.GreyLevelVariance += (iInt - results.MeanGreyLevelCount)* (iInt - results.MeanGreyLevelCount) * pij;
326  results.DependenceCountVariance += (k - results.MeanDependenceCount)* (k - results.MeanDependenceCount) * pij;
327  }
328  }
329  results.LowDependenceEmphasis /= Ns;
330  results.HighDependenceEmphasis /= Ns;
331  results.LowGreyLevelCountEmphasis /= Ns;
332  results.HighGreyLevelCountEmphasis /= Ns;
333  results.LowDependenceLowGreyLevelEmphasis /= Ns;
334  results.LowDependenceHighGreyLevelEmphasis /= Ns;
335  results.HighDependenceLowGreyLevelEmphasis /= Ns;
336  results.HighDependenceHighGreyLevelEmphasis /= Ns;
337 
338  results.GreyLevelNonUniformity /= Ns;
339  results.GreyLevelNonUniformityNormalised /= (Ns*Ns);
340  results.DependenceCountNonUniformity /= Ns;
341  results.DependenceCountNonUniformityNormalised /= (Ns*Ns);
342 
343  results.DependenceCountPercentage = 1;
344 
345  results.ExpectedNeighbourhoodSize = holder.m_NeighbourhoodSize;
346  results.AverageNeighbourhoodSize = holder.m_NumberOfNeighbourVoxels / (1.0 * holder.m_NumberOfNeighbourhoods);
347  results.AverageIncompleteNeighbourhoodSize = (holder.m_NumberOfNeighbourVoxels - holder.m_NumberOfCompleteNeighbourhoods* holder.m_NeighbourhoodSize) / (1.0 * (holder.m_NumberOfNeighbourhoods - holder.m_NumberOfCompleteNeighbourhoods));
348  results.PercentageOfCompleteNeighbourhoods = (1.0*holder.m_NumberOfCompleteNeighbourhoods) / (1.0 * holder.m_NumberOfNeighbourhoods);
349  results.PercentageOfDependenceNeighbours = holder.m_NumberOfDependenceNeighbourVoxels / (1.0 * holder.m_NumberOfNeighbourVoxels);
350 }
351 
352 template<typename TPixel, unsigned int VImageDimension>
353 void
355 {
356  typedef itk::Image<unsigned short, VImageDimension> MaskType;
357 
358  double rangeMin = config.MinimumIntensity;
359  double rangeMax = config.MaximumIntensity;
360  int numberOfBins = config.Bins;
361 
362  typename MaskType::Pointer maskImage = MaskType::New();
363  mitk::CastToItkImage(mask, maskImage);
364 
365  std::vector<mitk::NGLDMMatrixFeatures> resultVector;
366  int numberofDependency = 37;
367  if (VImageDimension == 2)
368  numberofDependency = 37;
369 
370  mitk::NGLDMMatrixHolder holderOverall(rangeMin, rangeMax, numberOfBins, numberofDependency);
371  mitk::NGLDMMatrixFeatures overallFeature;
372  CalculateNGLDMMatrix<TPixel, VImageDimension>(itkImage, maskImage, config.alpha, config.range, config.direction, holderOverall);
373  LocalCalculateFeatures(holderOverall, overallFeature);
374 
375  MatrixFeaturesTo(overallFeature, config.FeatureEncoding, featureList);
376 }
377 
378 
379 static
380 void MatrixFeaturesTo(mitk::NGLDMMatrixFeatures features,
381  std::string prefix,
383 {
384  featureList.push_back(std::make_pair(prefix + "Low Dependence Emphasis", features.LowDependenceEmphasis));
385  featureList.push_back(std::make_pair(prefix + "High Dependence Emphasis", features.HighDependenceEmphasis));
386  featureList.push_back(std::make_pair(prefix + "Low Grey Level Count Emphasis", features.LowGreyLevelCountEmphasis));
387  featureList.push_back(std::make_pair(prefix + "High Grey Level Count Emphasis", features.HighGreyLevelCountEmphasis));
388  featureList.push_back(std::make_pair(prefix + "Low Dependence Low Grey Level Emphasis", features.LowDependenceLowGreyLevelEmphasis));
389  featureList.push_back(std::make_pair(prefix + "Low Dependence High Grey Level Emphasis", features.LowDependenceHighGreyLevelEmphasis));
390  featureList.push_back(std::make_pair(prefix + "High Dependence Low Grey Level Emphasis", features.HighDependenceLowGreyLevelEmphasis));
391  featureList.push_back(std::make_pair(prefix + "High Dependence High Grey Level Emphasis", features.HighDependenceHighGreyLevelEmphasis));
392 
393  featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity", features.GreyLevelNonUniformity));
394  featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity Normalised", features.GreyLevelNonUniformityNormalised));
395  featureList.push_back(std::make_pair(prefix + "Dependence Count Non-Uniformity", features.DependenceCountNonUniformity));
396  featureList.push_back(std::make_pair(prefix + "Dependence Count Non-Uniformity Normalised", features.DependenceCountNonUniformityNormalised));
397 
398  featureList.push_back(std::make_pair(prefix + "Dependence Count Percentage", features.DependenceCountPercentage));
399  featureList.push_back(std::make_pair(prefix + "Grey Level Mean", features.MeanGreyLevelCount));
400  featureList.push_back(std::make_pair(prefix + "Grey Level Variance", features.GreyLevelVariance));
401  featureList.push_back(std::make_pair(prefix + "Dependence Count Mean", features.MeanDependenceCount));
402  featureList.push_back(std::make_pair(prefix + "Dependence Count Variance", features.DependenceCountVariance));
403  featureList.push_back(std::make_pair(prefix + "Dependence Count Entropy", features.DependenceCountEntropy));
404  featureList.push_back(std::make_pair(prefix + "Dependence Count Energy", features.DependenceCountEnergy));
405 
406  featureList.push_back(std::make_pair(prefix + "Expected Neighbourhood Size", features.ExpectedNeighbourhoodSize));
407  featureList.push_back(std::make_pair(prefix + "Average Neighbourhood Size", features.AverageNeighbourhoodSize));
408  featureList.push_back(std::make_pair(prefix + "Average Incomplete Neighbourhood Size", features.AverageIncompleteNeighbourhoodSize));
409  featureList.push_back(std::make_pair(prefix + "Percentage of complete Neighbourhoods", features.PercentageOfCompleteNeighbourhoods));
410  featureList.push_back(std::make_pair(prefix + "Percentage of Dependence Neighbour Voxels", features.PercentageOfDependenceNeighbours));
411 
412 }
413 
415 m_Range(1.0)
416 {
417  SetShortName("ngld");
418  SetLongName("neighbouring-grey-level-dependence");
419  SetFeatureClassName("Neighbouring Grey Level Dependence");
420 }
421 
423 {
424  FeatureListType featureList;
425  InitializeQuantifier(image, mask);
426 
428  config.direction = GetDirection();
429  config.range = m_Range;
430  config.alpha = 0;
431 
432  config.MinimumIntensity = GetQuantifier()->GetMinimum();
433  config.MaximumIntensity = GetQuantifier()->GetMaximum();
434  config.Bins = GetQuantifier()->GetBins();
435 
437 
438  AccessByItk_3(image, CalculateCoocurenceFeatures, mask, featureList,config);
439 
440  return featureList;
441 }
442 
444 {
445  FeatureNameListType featureList;
446 
447  return featureList;
448 }
449 
450 
451 
453 {
454  std::string name = GetOptionPrefix();
455 
456  parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Calculate Neighbouring Grey Level Dependence Features", "Calculate Neighbouring grey level dependence based features", us::Any());
457  parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::String, "NGLD Range", "Define the range that is used (Semicolon-separated)", us::Any());
458 
459  AddQuantifierArguments(parser);
460 }
461 
462 void
464 {
465  auto parsedArgs = GetParameter();
466  std::string name = GetOptionPrefix();
467 
468  if (parsedArgs.count(GetLongName()))
469  {
470  std::vector<double> ranges;
471  if (parsedArgs.count(name + "::range"))
472  {
473  ranges = SplitDouble(parsedArgs[name + "::range"].ToString(), ';');
474  }
475  else
476  {
477  ranges.push_back(1);
478  }
479  for (double range : ranges)
480  {
481  InitializeQuantifierFromParameters(feature, maskNoNAN);
482  this->SetRange(range);
483  MITK_INFO << "Start calculating NGLD";
484  auto localResults = this->CalculateFeatures(feature, maskNoNAN);
485  featureList.insert(featureList.end(), localResults.begin(), localResults.end());
486  MITK_INFO << "Finished calculating NGLD";
487  }
488  }
489 }
490 
492 {
493  std::ostringstream ss;
494  ss << m_Range;
495  std::string strRange = ss.str();
496  return QuantifierParameterString() + "_Range-"+ss.str();
497 }
498 
void CalculateNGLDMMatrix(itk::Image< TPixel, VImageDimension > *itkImage, itk::Image< unsigned short, VImageDimension > *mask, int alpha, int range, unsigned int direction, mitk::NGLDMMatrixHolder &holder)
#define AccessByItk_3(mitkImage, itkImageTypeFunction, arg1, arg2, arg3)
float k(1.0)
std::string GetCurrentFeatureEncoding() override
Adds an additional Separator to the name of the feature, which encodes the used parameters.
#define MITK_INFO
Definition: mitkLogMacros.h:18
itk::Image< unsigned char, 3 > ImageType
virtual void SetLongName(std::string _arg)
DataCollection - Class to facilitate loading/accessing structured data.
void addArgument(const std::string &longarg, const std::string &shortarg, Type type, const std::string &argLabel, const std::string &argHelp=std::string(), const us::Any &defaultValue=us::Any(), bool optional=true, bool ignoreRest=false, bool deprecated=false, mitkCommandLineParser::Channel channel=mitkCommandLineParser::Channel::None)
virtual void SetShortName(std::string _arg)
std::vector< std::pair< std::string, double > > FeatureListType
std::string FeatureDescriptionPrefix()
Returns a string that encodes the feature class name.
void CalculateFeaturesUsingParameters(const Image::Pointer &feature, const Image::Pointer &mask, const Image::Pointer &maskNoNAN, FeatureListType &featureList) override
Calculates the feature of this abstact interface. Does not necessarily considers the parameter settin...
virtual std::string GetLongName() const
void LocalCalculateFeatures(mitk::NGLDMMatrixHolder &holder, mitk::NGLDMMatrixFeatures &results)
Definition: usAny.h:163
itk::Image< unsigned char, 3 > MaskImageType
Definition: CLBrainMask.cpp:32
virtual int GetDirection() const
static T max(T x, T y)
Definition: svm.cpp:56
mitk::Image::Pointer image
void AddQuantifierArguments(mitkCommandLineParser &parser)
static T min(T x, T y)
Definition: svm.cpp:53
std::vector< double > SplitDouble(std::string str, char delimiter)
virtual void SetFeatureClassName(std::string _arg)
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
static void MatrixFeaturesTo(mitk::NGLDMMatrixFeatures features, std::string prefix, mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType &featureList)
void CalculateCoocurenceFeatures(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer mask, mitk::GIFNeighbouringGreyLevelDependenceFeature::FeatureListType &featureList, mitk::GIFNeighbouringGreyLevelDependenceFeature::GIFNeighbouringGreyLevelDependenceFeatureConfiguration config)
mitk::Image::Pointer mask
FeatureNameListType GetFeatureNames() override
Returns a list of the names of all features that are calculated from this class.
virtual IntensityQuantifier::Pointer GetQuantifier()
void InitializeQuantifierFromParameters(const Image::Pointer &feature, const Image::Pointer &mask, unsigned int defaultBins=256)
FeatureListType CalculateFeatures(const Image::Pointer &image, const Image::Pointer &feature) override
Calculates the Cooccurence-Matrix based features for this class.
void InitializeQuantifier(const Image::Pointer &feature, const Image::Pointer &mask, unsigned int defaultBins=256)
virtual ParameterTypes GetParameter() const