Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkGIFFirstOrderStatistics.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
21 #include <itkLabelStatisticsImageFilter.h>
22 #include <itkMinimumMaximumImageCalculator.h>
23 
24 // STL
25 #include <sstream>
26 
27 template<typename TPixel, unsigned int VImageDimension>
28 void
30 {
31  typedef itk::Image<TPixel, VImageDimension> ImageType;
32  typedef itk::Image<unsigned short, VImageDimension> MaskType;
33  typedef itk::LabelStatisticsImageFilter<ImageType, MaskType> FilterType;
34  typedef typename FilterType::HistogramType HistogramType;
35  typedef typename HistogramType::IndexType HIndexType;
36  typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxComputerType;
37 
38  typename MaskType::Pointer maskImage = MaskType::New();
39  mitk::CastToItkImage(mask, maskImage);
40 
41  double voxelVolume = 1;
42  for (unsigned int i = 0; i < std::min<unsigned int>(3, VImageDimension); ++i)
43  voxelVolume *= itkImage->GetSpacing()[i];
44  double voxelSpace = 1;
45  for (unsigned int i = 0; i < VImageDimension; ++i)
46  voxelSpace *= itkImage->GetSpacing()[i];
47 
48  typename MinMaxComputerType::Pointer minMaxComputer = MinMaxComputerType::New();
49  minMaxComputer->SetImage(itkImage);
50  minMaxComputer->Compute();
51  double imageRange = minMaxComputer->GetMaximum() - minMaxComputer->GetMinimum();
52 
53  typename FilterType::Pointer labelStatisticsImageFilter = FilterType::New();
54  labelStatisticsImageFilter->SetInput( itkImage );
55  labelStatisticsImageFilter->SetLabelInput(maskImage);
56  labelStatisticsImageFilter->SetUseHistograms(true);
57 
58  double min = params.MinimumIntensity;
59  double max = params.MaximumIntensity;
60 
61  labelStatisticsImageFilter->SetHistogramParameters(params.Bins, min,max);
62  labelStatisticsImageFilter->Update();
63 
64  // --------------- Range --------------------
65  double range = labelStatisticsImageFilter->GetMaximum(1) - labelStatisticsImageFilter->GetMinimum(1);
66  // --------------- Uniformity, Entropy --------------------
67  double count = labelStatisticsImageFilter->GetCount(1);
68  //double std_dev = labelStatisticsImageFilter->GetSigma(1);
69  double mean = labelStatisticsImageFilter->GetMean(1);
70  double median = labelStatisticsImageFilter->GetMedian(1);
71  auto histogram = labelStatisticsImageFilter->GetHistogram(1);
72  bool histogramIsCalculated = histogram;
73 
74  HIndexType index;
75  index.SetSize(1);
76 
77  double uniformity = 0;
78  double entropy = 0;
79  double squared_sum = 0;
80  double kurtosis = 0;
81  double mean_absolut_deviation = 0;
82  double median_absolut_deviation = 0;
83  double skewness = 0;
84  double sum_prob = 0;
85  double binWidth = 0;
86  double p05th = 0, p10th = 0, p15th = 0, p20th = 0, p25th = 0, p30th = 0, p35th = 0, p40th = 0, p45th = 0, p50th = 0;
87  double p55th = 0, p60th = 0, p65th = 0, p70th = 0, p75th = 0, p80th = 0, p85th = 0, p90th = 0, p95th = 0;
88 
89  double voxelValue = 0;
90 
91  if (histogramIsCalculated)
92  {
93  binWidth = histogram->GetBinMax(0, 0) - histogram->GetBinMin(0, 0);
94  p05th = histogram->Quantile(0, 0.05);
95  p10th = histogram->Quantile(0, 0.10);
96  p15th = histogram->Quantile(0, 0.15);
97  p20th = histogram->Quantile(0, 0.20);
98  p25th = histogram->Quantile(0, 0.25);
99  p30th = histogram->Quantile(0, 0.30);
100  p35th = histogram->Quantile(0, 0.35);
101  p40th = histogram->Quantile(0, 0.40);
102  p45th = histogram->Quantile(0, 0.45);
103  p50th = histogram->Quantile(0, 0.50);
104  p55th = histogram->Quantile(0, 0.55);
105  p60th = histogram->Quantile(0, 0.60);
106  p65th = histogram->Quantile(0, 0.65);
107  p70th = histogram->Quantile(0, 0.70);
108  p75th = histogram->Quantile(0, 0.75);
109  p80th = histogram->Quantile(0, 0.80);
110  p85th = histogram->Quantile(0, 0.85);
111  p90th = histogram->Quantile(0, 0.90);
112  p95th = histogram->Quantile(0, 0.95);
113  }
114  double Log2=log(2);
115  double mode_bin;
116  double mode_value = 0;
117  double variance = 0;
118  if (histogramIsCalculated)
119  {
120  for (int i = 0; i < (int)(histogram->GetSize(0)); ++i)
121  {
122  index[0] = i;
123  double prob = histogram->GetFrequency(index);
124 
125 
126  if (prob < 0.00000001)
127  continue;
128 
129  voxelValue = histogram->GetBinMin(0, i) + binWidth * 0.5;
130 
131  if (prob > mode_value)
132  {
133  mode_value = prob;
134  mode_bin = voxelValue;
135  }
136 
137  sum_prob += prob;
138  squared_sum += prob * voxelValue*voxelValue;
139 
140  prob /= count;
141  mean_absolut_deviation += prob* std::abs(voxelValue - mean);
142  median_absolut_deviation += prob* std::abs(voxelValue - median);
143  variance += prob * (voxelValue - mean) * (voxelValue - mean);
144 
145  kurtosis += prob* (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean);
146  skewness += prob* (voxelValue - mean) * (voxelValue - mean) * (voxelValue - mean);
147 
148  uniformity += prob*prob;
149  if (prob > 0)
150  {
151  entropy += prob * std::log(prob) / Log2;
152  }
153  }
154  }
155  entropy = -entropy;
156 
157  double uncorrected_std_dev = std::sqrt(variance);
158  double rms = std::sqrt(squared_sum / count);
159  kurtosis = kurtosis / (variance * variance);
160  skewness = skewness / (variance * uncorrected_std_dev);
161  double coveredGrayValueRange = range / imageRange;
162  double coefficient_of_variation = (mean == 0) ? 0 : std::sqrt(variance) / mean;
163  double quantile_coefficient_of_dispersion = (p75th - p25th) / (p75th + p25th);
164 
165  //Calculate the robust mean absolute deviation
166  //First, set all frequencies to 0 that are <10th or >90th percentile
167  double meanRobust = 0.0;
168  double robustMeanAbsoluteDeviation = 0.0;
169  if (histogramIsCalculated)
170  {
171  for (int i = 0; i < (int)(histogram->GetSize(0)); ++i)
172  {
173  index[0] = i;
174  if (histogram->GetBinMax(0, i) < p10th)
175  {
176  histogram->SetFrequencyOfIndex(index, 0);
177  }
178  else if (histogram->GetBinMin(0, i) > p90th)
179  {
180  histogram->SetFrequencyOfIndex(index, 0);
181  }
182  }
183 
184  //Calculate the mean
185  for (int i = 0; i < (int)(histogram->GetSize(0)); ++i)
186  {
187  index[0] = i;
188  meanRobust += histogram->GetFrequency(index) * 0.5 * (histogram->GetBinMin(0, i) + histogram->GetBinMax(0, i));
189  }
190  meanRobust = meanRobust / histogram->GetTotalFrequency();
191  for (int i = 0; i < (int)(histogram->GetSize(0)); ++i)
192  {
193  index[0] = i;
194  robustMeanAbsoluteDeviation += std::abs(histogram->GetFrequency(index) *
195  ((0.5 * (histogram->GetBinMin(0, i) + histogram->GetBinMax(0, i)))
196  - meanRobust
197  ));
198  }
199  robustMeanAbsoluteDeviation = robustMeanAbsoluteDeviation / histogram->GetTotalFrequency();
200  }
201 
202  featureList.push_back(std::make_pair(params.prefix + "Mean", labelStatisticsImageFilter->GetMean(1)));
203  featureList.push_back(std::make_pair(params.prefix + "Unbiased Variance", labelStatisticsImageFilter->GetVariance(1))); //Siehe Definition von Unbiased Variance estimation. (Wird nicht durch n sondern durch n-1 normalisiert)
204  featureList.push_back(std::make_pair(params.prefix + "Biased Variance", variance));
205  featureList.push_back(std::make_pair(params.prefix + "Skewness", skewness));
206  featureList.push_back(std::make_pair(params.prefix + "Kurtosis", kurtosis));
207  featureList.push_back(std::make_pair(params.prefix + "Median", labelStatisticsImageFilter->GetMedian(1)));
208  featureList.push_back(std::make_pair(params.prefix + "Minimum", labelStatisticsImageFilter->GetMinimum(1)));
209  featureList.push_back(std::make_pair(params.prefix + "Maximum", labelStatisticsImageFilter->GetMaximum(1)));
210  featureList.push_back(std::make_pair(params.prefix + "Range", range));
211  featureList.push_back(std::make_pair(params.prefix + "Mean Absolute Deviation", mean_absolut_deviation));
212  featureList.push_back(std::make_pair(params.prefix + "Robust Mean Absolute Deviation", robustMeanAbsoluteDeviation));
213  featureList.push_back(std::make_pair(params.prefix + "Median Absolute Deviation", median_absolut_deviation));
214  featureList.push_back(std::make_pair(params.prefix + "Coefficient Of Variation", coefficient_of_variation));
215  featureList.push_back(std::make_pair(params.prefix + "Quantile Coefficient Of Dispersion", quantile_coefficient_of_dispersion));
216  featureList.push_back(std::make_pair(params.prefix + "Energy", squared_sum));
217  featureList.push_back(std::make_pair(params.prefix + "Root Mean Square", rms));
218 
219  typename HistogramType::MeasurementVectorType mv(1);
220  mv[0] = 0;
221  typename HistogramType::IndexType resultingIndex;
222  histogram->GetIndex(mv, resultingIndex);
223  featureList.push_back(std::make_pair(params.prefix + "Robust Mean", meanRobust));
224  featureList.push_back(std::make_pair(params.prefix + "Uniformity", uniformity));
225  featureList.push_back(std::make_pair(params.prefix + "Entropy", entropy));
226  featureList.push_back(std::make_pair(params.prefix + "Excess Kurtosis", kurtosis - 3));
227  featureList.push_back(std::make_pair(params.prefix + "Covered Image Intensity Range", coveredGrayValueRange));
228  featureList.push_back(std::make_pair(params.prefix + "Sum", labelStatisticsImageFilter->GetSum(1)));
229  featureList.push_back(std::make_pair(params.prefix + "Mode", mode_bin));
230  featureList.push_back(std::make_pair(params.prefix + "Mode Probability", mode_value));
231  featureList.push_back(std::make_pair(params.prefix + "Unbiased Standard deviation", labelStatisticsImageFilter->GetSigma(1)));
232  featureList.push_back(std::make_pair(params.prefix + "Biased Standard deviation", sqrt(variance)));
233  featureList.push_back(std::make_pair(params.prefix + "Number Of Voxels", labelStatisticsImageFilter->GetCount(1)));
234 
235  featureList.push_back(std::make_pair(params.prefix + "05th Percentile", p05th));
236  featureList.push_back(std::make_pair(params.prefix + "10th Percentile", p10th));
237  featureList.push_back(std::make_pair(params.prefix + "15th Percentile", p15th));
238  featureList.push_back(std::make_pair(params.prefix + "20th Percentile", p20th));
239  featureList.push_back(std::make_pair(params.prefix + "25th Percentile", p25th));
240  featureList.push_back(std::make_pair(params.prefix + "30th Percentile", p30th));
241  featureList.push_back(std::make_pair(params.prefix + "35th Percentile", p35th));
242  featureList.push_back(std::make_pair(params.prefix + "40th Percentile", p40th));
243  featureList.push_back(std::make_pair(params.prefix + "45th Percentile", p45th));
244  featureList.push_back(std::make_pair(params.prefix + "50th Percentile", p50th));
245  featureList.push_back(std::make_pair(params.prefix + "55th Percentile", p55th));
246  featureList.push_back(std::make_pair(params.prefix + "60th Percentile", p60th));
247  featureList.push_back(std::make_pair(params.prefix + "65th Percentile", p65th));
248  featureList.push_back(std::make_pair(params.prefix + "70th Percentile", p70th));
249  featureList.push_back(std::make_pair(params.prefix + "75th Percentile", p75th));
250  featureList.push_back(std::make_pair(params.prefix + "80th Percentile", p80th));
251  featureList.push_back(std::make_pair(params.prefix + "85th Percentile", p85th));
252  featureList.push_back(std::make_pair(params.prefix + "90th Percentile", p90th));
253  featureList.push_back(std::make_pair(params.prefix + "95th Percentile", p95th));
254  featureList.push_back(std::make_pair(params.prefix + "Interquartile Range", (p75th - p25th)));
255  featureList.push_back(std::make_pair(params.prefix + "Image Dimension", VImageDimension));
256  featureList.push_back(std::make_pair(params.prefix + "Voxel Space", voxelSpace));
257  featureList.push_back(std::make_pair(params.prefix + "Voxel Volume", voxelVolume));
258 }
259 
261 {
262  SetShortName("fo");
263  SetLongName("first-order");
264  SetFeatureClassName("First Order");
265 }
266 
268 {
269  InitializeQuantifier(image, mask);
270  FeatureListType featureList;
271 
272  ParameterStruct params;
273 
274  params.MinimumIntensity = GetQuantifier()->GetMinimum();
275  params.MaximumIntensity = GetQuantifier()->GetMaximum();
276  params.Bins = GetQuantifier()->GetBins();
277  params.prefix = FeatureDescriptionPrefix();
278  AccessByItk_3(image, CalculateFirstOrderStatistics, mask, featureList, params);
279 
280  return featureList;
281 }
282 
284 {
285  FeatureNameListType featureList;
286  featureList.push_back("First Order::Minimum");
287  featureList.push_back("First Order::Maximum");
288  featureList.push_back("First Order::Mean");
289  featureList.push_back("First Order::Variance");
290  featureList.push_back("First Order::Sum");
291  featureList.push_back("First Order::Median");
292  featureList.push_back("First Order::Standard deviation");
293  featureList.push_back("First Order::No. of Voxel");
294  return featureList;
295 }
296 
297 
299 {
300  std::string name = GetOptionPrefix();
301 
302  parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Volume-Statistic", "calculates volume based features", us::Any());
303  AddQuantifierArguments(parser);
304 }
305 
306 void
308 {
309  auto parsedArgs = GetParameter();
310  if (parsedArgs.count(GetLongName()))
311  {
312  InitializeQuantifierFromParameters(feature, maskNoNAN);
313  MITK_INFO << "Start calculating first order features ....";
314  auto localResults = this->CalculateFeatures(feature, maskNoNAN);
315  featureList.insert(featureList.end(), localResults.begin(), localResults.end());
316  MITK_INFO << "Finished calculating first order features....";
317  }
318 }
319 
321 {
322  return QuantifierParameterString();
323 }
#define AccessByItk_3(mitkImage, itkImageTypeFunction, arg1, arg2, arg3)
#define MITK_INFO
Definition: mitkLogMacros.h:18
itk::Image< unsigned char, 3 > ImageType
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)
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 SetShortName(std::string _arg)
std::vector< std::pair< std::string, double > > FeatureListType
FeatureNameListType GetFeatureNames() override
Returns a list of the names of all features that are calculated from this class.
std::string FeatureDescriptionPrefix()
Returns a string that encodes the feature class name.
virtual std::string GetLongName() const
void AddArguments(mitkCommandLineParser &parser) override
void CalculateFirstOrderStatistics(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer mask, mitk::GIFFirstOrderStatistics::FeatureListType &featureList, mitk::GIFFirstOrderStatistics::ParameterStruct params)
Definition: usAny.h:163
FeatureListType CalculateFeatures(const Image::Pointer &image, const Image::Pointer &feature) override
Calculates the First Order Features based on a binned version of the image.
static T max(T x, T y)
Definition: svm.cpp:56
mitk::Image::Pointer image
std::string GetCurrentFeatureEncoding() override
Adds an additional Separator to the name of the feature, which encodes the used parameters.
void AddQuantifierArguments(mitkCommandLineParser &parser)
static T min(T x, T y)
Definition: svm.cpp:53
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.
mitk::Image::Pointer mask
virtual IntensityQuantifier::Pointer GetQuantifier()
void InitializeQuantifierFromParameters(const Image::Pointer &feature, const Image::Pointer &mask, unsigned int defaultBins=256)
void InitializeQuantifier(const Image::Pointer &feature, const Image::Pointer &mask, unsigned int defaultBins=256)
virtual ParameterTypes GetParameter() const