Medical Imaging Interaction Toolkit  2018.4.99-b585543d
Medical Imaging Interaction Toolkit
itkFirstOrderStatisticsFeatureFunctor.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 itkNeighborhoodFirstOrderStatistics_h
14 #define itkNeighborhoodFirstOrderStatistics_h
15 
16 #include "itkConstNeighborhoodIterator.h"
17 
18 namespace itk
19 {
20 
21 namespace Functor
22 {
23 
24 template< typename TNeighborhoodType, typename TPixelOutputType >
26 {
27  static const unsigned int OutputCount = 6;
28  typedef vnl_vector_fixed<TPixelOutputType, OutputCount> OutputVectorType;
29  typedef TNeighborhoodType NeighborhoodType;
31  {
33  };
34 
35  NeighborhoodFirstOrderStatistics(){std::cout << "NeighborhoodFirstOrderStatistics" << std::endl;}
36 
37  static const char * GetFeatureName(unsigned int f )
38  {
39  static const char * const FeatureNames[] = {"MEAN","VARIANCE","SKEWNESS", "KURTOSIS", "MIN", "MAX"};
40  return FeatureNames[f];
41  }
42 
43  inline OutputVectorType operator()(const NeighborhoodType & it) const
44  {
45  double mean = 0;
46  double sqmean = 0;
47  double min = 10000000;
48  double max = -10000000;
49 
50  for (unsigned int i = 0; i < it.Size(); ++i)
51  {
52  double value = it.GetPixel(i);
53  mean += value;
54  sqmean += (value * value);
55  max = (value > max) ? value : max;
56  min = (value < min) ? value : min;
57  }
58  mean /= it.Size();
59  sqmean /= it.Size();
60 
61  double variance = sqmean - mean * mean;
62 
63  double stdderivation = std::sqrt(variance);
64  double skewness = 0;
65  double kurtosis = 0;
66 
67  for (unsigned int i = 0; i < it.Size(); ++i)
68  {
69  double value = it.GetPixel(i);
70  double tmpskewness = (value - mean) / stdderivation;
71  skewness += tmpskewness * tmpskewness * tmpskewness;
72  kurtosis += (value - mean) * (value - mean) * (value - mean) * (value - mean);
73  }
74  skewness /= it.Size();
75  kurtosis /= it.Size();
76  kurtosis /= (variance * variance);
77 
78  if(skewness!=skewness){
79  skewness = 0;
80  }
81  if(kurtosis!=kurtosis){
82  kurtosis = 0;
83  }
84 
85  OutputVectorType output_vector;
86  output_vector[MEAN] = mean;
87  output_vector[VARIANCE] = variance;
88  output_vector[SKEWNESS] = skewness;
89  output_vector[KURTOSIS] = kurtosis;
90  output_vector[MIN] = min;
91  output_vector[MAX] = max;
92 
93  return output_vector;
94 
95  }
96 };
97 
98 }// end namespace functor
99 
100 } // end namespace itk
101 #endif
vnl_vector_fixed< TPixelOutputType, OutputCount > OutputVectorType
OutputVectorType operator()(const NeighborhoodType &it) const
static T max(T x, T y)
Definition: svm.cpp:56
static T min(T x, T y)
Definition: svm.cpp:53