Medical Imaging Interaction Toolkit  2018.4.99-87d68d9f
Medical Imaging Interaction Toolkit
mitkGIFCurvatureStatistic.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 #include <mitkPixelTypeMultiplex.h>
21 
22 // ITK
23 #include <itkLabelStatisticsImageFilter.h>
24 #include <itkNeighborhoodIterator.h>
25 
26 // VTK
27 #include <vtkSmartPointer.h>
28 #include <vtkImageMarchingCubes.h>
29 #include <vtkCurvatures.h>
30 #include <vtkPointData.h>
31 
32 // STL
33 #include <sstream>
34 #include <limits>
35 
36 static void calculateLocalStatistic(vtkDataArray* scalars, std::string name, std::string featureDescriptionPrefix, mitk::GIFCurvatureStatistic::FeatureListType & featureList)
37 {
38  int size = scalars->GetNumberOfTuples();
39  double minimum = std::numeric_limits<double>::max();
40  double maximum = std::numeric_limits<double>::lowest();
41  double mean1 = 0;
42  double mean2 = 0;
43  double mean3 = 0;
44  double mean4 = 0;
45  double mean1p = 0;
46  double mean2p = 0;
47  double mean3p = 0;
48  double mean4p = 0;
49  double mean1n = 0;
50  double mean2n = 0;
51  double mean3n = 0;
52  double mean4n = 0;
53  int countPositive = 0;
54  int countNegative = 0;
55 
56  for (int i = 0; i < size; ++i)
57  {
58  double actualValue = scalars->GetComponent(i, 0);
59  minimum = std::min<double>(minimum, scalars->GetComponent(i, 0));
60  maximum = std::max<double>(maximum, scalars->GetComponent(i, 0));
61  mean1 += actualValue;
62  mean2 += actualValue*actualValue;
63  mean3 += actualValue*actualValue*actualValue;
64  mean4 += actualValue*actualValue*actualValue*actualValue;
65  if (actualValue > 0)
66  {
67  mean1p += actualValue;
68  mean2p += actualValue*actualValue;
69  mean3p += actualValue*actualValue*actualValue;
70  mean4p += actualValue*actualValue*actualValue*actualValue;
71  countPositive++;
72  }
73  if (actualValue < 0)
74  {
75  mean1n += actualValue;
76  mean2n += actualValue*actualValue;
77  mean3n += actualValue*actualValue*actualValue;
78  mean4n += actualValue*actualValue*actualValue*actualValue;
79  countNegative++;
80  }
81  }
82  double mean = mean1 / size;
83  double stddev = std::sqrt(mean2 / size - mean*mean);
84  double skewness = ((mean3 / size) - 3 * mean*stddev*stddev - mean*mean*mean) / (stddev*stddev*stddev);
85 
86  double meanP = mean1p / countPositive;
87  double stddevP = std::sqrt(mean2p / countPositive - meanP*meanP);
88  double skewnessP = ((mean3p / countPositive) - 3 * meanP*stddevP*stddevP - meanP*meanP*meanP) / (stddevP*stddevP*stddevP);
89 
90  double meanN = mean1n / countNegative;
91  double stddevN = std::sqrt(mean2n / countNegative - meanN*meanN);
92  double skewnessN = ((mean3n / countNegative) - 3 * meanN*stddevN*stddevN - meanN*meanN*meanN) / (stddevN*stddevN*stddevN);
93 
94  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Minimum " + name + " Curvature", minimum));
95  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Maximum " + name + " Curvature", maximum));
96  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Mean " + name + " Curvature", mean));
97  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Standard Deviation " + name + " Curvature", stddev));
98  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Skewness " + name + " Curvature", skewness));
99  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Mean Positive " + name + " Curvature", meanP));
100  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Standard Deviation Positive " + name + " Curvature", stddevP));
101  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Skewness Positive " + name + " Curvature", skewnessP));
102  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Mean Negative " + name + " Curvature", meanN));
103  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Standard Deviation Negative " + name + " Curvature", stddevN));
104  featureList.push_back(std::make_pair(featureDescriptionPrefix + "Skewness Negative " + name + " Curvature", skewnessN));
105 
106 }
107 
109 {
110  SetLongName("curvature");
111  SetShortName("cur");
112  SetFeatureClassName("Curvature Feature");
113 }
114 
116 {
117  FeatureListType featureList;
118  if (image->GetDimension() < 3)
119  {
120  return featureList;
121  }
122 
123  vtkSmartPointer<vtkImageMarchingCubes> mesher = vtkSmartPointer<vtkImageMarchingCubes>::New();
124  vtkSmartPointer<vtkCurvatures> curvator = vtkSmartPointer<vtkCurvatures>::New();
125  mesher->SetInputData(mask->GetVtkImageData());
126  mesher->SetValue(0, 0.5);
127  curvator->SetInputConnection(mesher->GetOutputPort());
128  curvator->SetCurvatureTypeToMean();
129  curvator->Update();
130  vtkDataArray* scalars = curvator->GetOutput()->GetPointData()->GetScalars();
131  calculateLocalStatistic(scalars, "Mean", FeatureDescriptionPrefix(), featureList);
132 
133  curvator->SetCurvatureTypeToGaussian();
134  curvator->Update();
135  scalars = curvator->GetOutput()->GetPointData()->GetScalars();
136  calculateLocalStatistic(scalars, "Gaussian", FeatureDescriptionPrefix(), featureList);
137 
138  curvator->SetCurvatureTypeToMinimum();
139  curvator->Update();
140  scalars = curvator->GetOutput()->GetPointData()->GetScalars();
141  calculateLocalStatistic(scalars, "Minimum", FeatureDescriptionPrefix(), featureList);
142 
143  curvator->SetCurvatureTypeToMaximum();
144  curvator->Update();
145  scalars = curvator->GetOutput()->GetPointData()->GetScalars();
146  calculateLocalStatistic(scalars, "Maximum", FeatureDescriptionPrefix(), featureList);
147 
148  return featureList;
149 }
150 
152 {
153  FeatureNameListType featureList;
154  return featureList;
155 }
156 
157 
159 {
160  std::string name = GetOptionPrefix();
161 
162  parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Curvature of Surface as feature", "calculates shape curvature based features", us::Any());
163 }
164 
165 void
167 {
168  auto parsedArgs = GetParameter();
169  if (parsedArgs.count(GetLongName()))
170  {
171  MITK_INFO << "Start calculating volumetric features ....";
172  auto localResults = this->CalculateFeatures(feature, mask);
173  featureList.insert(featureList.end(), localResults.begin(), localResults.end());
174  MITK_INFO << "Finished calculating volumetric features....";
175  }
176 }
177 
#define MITK_INFO
Definition: mitkLogMacros.h:18
virtual void SetLongName(std::string _arg)
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.
FeatureListType CalculateFeatures(const Image::Pointer &image, const Image::Pointer &feature) override
Calculates the Cooccurence-Matrix based features for this class.
virtual std::string GetLongName() const
static void calculateLocalStatistic(vtkDataArray *scalars, std::string name, std::string featureDescriptionPrefix, mitk::GIFCurvatureStatistic::FeatureListType &featureList)
Definition: usAny.h:163
FeatureNameListType GetFeatureNames() override
Returns a list of the names of all features that are calculated from this class.
static T max(T x, T y)
Definition: svm.cpp:56
mitk::Image::Pointer image
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 void SetFeatureClassName(std::string _arg)
mitk::Image::Pointer mask
void AddArguments(mitkCommandLineParser &parser) override
virtual ParameterTypes GetParameter() const