Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkGIFFirstOrderNumericStatistics.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 #include <itkImageRegionIterator.h>
24 
25 // STL
26 #include <sstream>
27 
28 struct FirstOrderNumericParameterStruct {
30  double MinimumIntensity;
31  double MaximumIntensity;
32  int Bins;
33  std::string prefix;
34 };
35 
36 template<typename TPixel, unsigned int VImageDimension>
37 void
38 CalculateFirstOrderStatistics(itk::Image<TPixel, VImageDimension>* itkImage, mitk::Image::Pointer mask, mitk::GIFFirstOrderNumericStatistics::FeatureListType & featureList, FirstOrderNumericParameterStruct params)
39 {
40  typedef itk::Image<TPixel, VImageDimension> ImageType;
41  typedef itk::Image<unsigned short, VImageDimension> MaskType;
42 
43  typename MaskType::Pointer maskImage = MaskType::New();
44  mitk::CastToItkImage(mask, maskImage);
45 
46  //
47  // Calculate the Volume of Voxel (Maximum to the 3th order)
48  //
49  double voxelVolume = 1;
50  for (unsigned int i = 0; i < std::min<unsigned int>(3, VImageDimension); ++i)
51  voxelVolume *= itkImage->GetSpacing()[i];
52 
53  //
54  // Calculate the Hypervolume of Voxel
55  //
56  double voxelSpace = 1;
57  for (unsigned int i = 0; i < VImageDimension; ++i)
58  voxelSpace *= itkImage->GetSpacing()[i];
59 
60  unsigned int numberOfBins = params.quantifier->GetBins();
61  std::vector<unsigned int> histogram;
62  histogram.resize(numberOfBins, 0);
63 
64 
65  double minimum = std::numeric_limits<double>::max();
66  double maximum = std::numeric_limits<double>::lowest();
67  double absoluteMinimum = std::numeric_limits<double>::max();
68  double absoluteMaximum = std::numeric_limits<double>::lowest();
69  double sum = 0;
70  double sumTwo= 0;
71  double sumThree = 0;
72  unsigned int numberOfVoxels = 0;
73 
74  itk::ImageRegionIterator<ImageType> imageIter(itkImage, itkImage->GetLargestPossibleRegion());
75  itk::ImageRegionIterator<MaskType> maskIter(maskImage, maskImage->GetLargestPossibleRegion());
76 
77  while (!imageIter.IsAtEnd())
78  {
79  double value = imageIter.Get();
80 
81  absoluteMinimum = std::min<double>(minimum, value);
82  absoluteMaximum = std::max<double>(maximum, value);
83  if (maskIter.Get() > 0)
84  {
85  minimum = std::min<double>(minimum, value);
86  maximum = std::max<double>(maximum, value);
87 
88  sum += value;
89  sumTwo += value * value;
90  sumThree += value * value*value;
91 
92  histogram[params.quantifier->IntensityToIndex(value)] += 1;
93 
94  ++numberOfVoxels;
95  }
96  ++maskIter;
97  ++imageIter;
98  }
99 
100  //
101  // Histogram based calculations
102  //
103  unsigned int passedValues = 0;
104  double doubleNoOfVoxes = numberOfVoxels;
105  double median = 0;
106  double lastIntensityWithValues = params.quantifier->IndexToMeanIntensity(0);
107  std::size_t modeIdx = 0;
108 
109  double entropy = 0;
110  double uniformity = 0;
111 
112  std::vector<double> percentiles;
113  percentiles.resize(20, 0);
114 
115  for (std::size_t idx = 0; idx < histogram.size(); ++idx)
116  {
117  unsigned int actualValues = histogram[idx];
118 
119  for (std::size_t percentileIdx = 0; percentileIdx < percentiles.size(); ++percentileIdx)
120  {
121  double threshold = doubleNoOfVoxes * (percentileIdx + 1) *1.0 / (percentiles.size());
122  if ((passedValues < threshold) & ((passedValues + actualValues) >= threshold))
123  {
124  // Lower Bound
125  if (passedValues == std::floor(threshold))
126  {
127  percentiles[percentileIdx] = 0.5*(lastIntensityWithValues + params.quantifier->IndexToMeanIntensity(idx));
128  }
129  else
130  {
131  percentiles[percentileIdx] = params.quantifier->IndexToMeanIntensity(idx);
132  }
133  }
134  }
135 
136  if ((passedValues < doubleNoOfVoxes * 0.5) & ((passedValues + actualValues) >= doubleNoOfVoxes * 0.5))
137  {
138  // Lower Bound
139  if (passedValues == std::floor(doubleNoOfVoxes * 0.5))
140  {
141  median = 0.5*(lastIntensityWithValues + params.quantifier->IndexToMeanIntensity(idx));
142  }
143  else
144  {
145  median = params.quantifier->IndexToMeanIntensity(idx);
146  }
147  }
148 
149  if (actualValues > histogram[modeIdx])
150  {
151  modeIdx = idx;
152  }
153 
154  if (actualValues > 0)
155  {
156  lastIntensityWithValues = params.quantifier->IndexToMeanIntensity(idx);
157 
158  double currentProbability = actualValues / (1.0 *numberOfVoxels);
159  uniformity += currentProbability * currentProbability;
160  entropy += currentProbability * std::log(currentProbability) / std::log(2);
161  }
162  passedValues += actualValues;
163  }
164  double p10 = percentiles[1];
165  //double p25idx = params.quantifier->IntensityToIndex(percentiles[4]);
166  //double p75idx = params.quantifier->IntensityToIndex(percentiles[14]);
167  double p25idx = percentiles[4];
168  double p75idx = percentiles[14];
169  double p90 = percentiles[17];
170 
171  double mean = sum / (numberOfVoxels);
172  double variance = sumTwo / (numberOfVoxels) - (mean*mean);
173  double energy = sumTwo;
174  double rootMeanSquare = std::sqrt(sumTwo / numberOfVoxels);
175 
176  double sumAbsoluteDistanceToMean = 0;
177  double sumAbsoluteDistanceToMedian = 0;
178  double sumRobust = 0;
179  double sumRobustSquare = 0;
180  double sumRobustAbsolulteDistanceToMean = 0;
181  double sumValueMinusMean = 0;
182  double sumValueMinusMeanThree = 0;
183  double sumValueMinusMeanFour = 0;
184  unsigned int numberOfRobustVoxel = 0;
185 
186 
187  maskIter.GoToBegin();
188  imageIter.GoToBegin();
189  while (!imageIter.IsAtEnd())
190  {
191  if (maskIter.Get() > 0)
192  {
193  double value = imageIter.Get();
194  double valueMinusMean = value - mean;
195 
196  sumAbsoluteDistanceToMean += std::abs<double>(valueMinusMean);
197  sumAbsoluteDistanceToMedian += std::abs<double>(value - median);
198  sumValueMinusMean += valueMinusMean;
199  sumValueMinusMeanThree += valueMinusMean * valueMinusMean * valueMinusMean;
200  sumValueMinusMeanFour += valueMinusMean * valueMinusMean * valueMinusMean * valueMinusMean;
201 
202  if ((p10 <= value) & (value <= p90))
203  {
204  sumRobust += value;
205  sumRobustSquare += value * value;
206  ++numberOfRobustVoxel;
207  }
208  }
209  ++maskIter;
210  ++imageIter;
211  }
212  double robustMean = sumRobust / numberOfRobustVoxel;
213  double robustVariance = sumRobustSquare / numberOfRobustVoxel - (robustMean * robustMean);
214 
215  maskIter.GoToBegin();
216  imageIter.GoToBegin();
217  while (!imageIter.IsAtEnd())
218  {
219  if (maskIter.Get() > 0)
220  {
221  double value = imageIter.Get();
222 
223  if ((p10 <= value) & (value <= p90))
224  {
225  sumRobustAbsolulteDistanceToMean += std::abs(value - robustMean);
226  }
227  }
228  ++maskIter;
229  ++imageIter;
230  }
231 
232 
233 
234  double meanAbsoluteDeviation = sumAbsoluteDistanceToMean / numberOfVoxels;
235  double medianAbsoluteDeviation = sumAbsoluteDistanceToMedian / numberOfVoxels;
236  double robustMeanAbsoluteDeviation = sumRobustAbsolulteDistanceToMean / numberOfRobustVoxel;
237  double skewness = sumValueMinusMeanThree / numberOfVoxels / variance / std::sqrt(variance);
238  double kurtosis = sumValueMinusMeanFour / numberOfVoxels / variance / variance;
239  double interquantileRange = p75idx - p25idx;
240  double coefficientOfVariation = std::sqrt(variance) / mean;
241  double quantileCoefficientOfDispersion = (p75idx - p25idx) / (p75idx + p25idx);
242  double coveredImageRange = (maximum - minimum)/ (absoluteMaximum - absoluteMinimum) ;
243 
244  featureList.push_back(std::make_pair(params.prefix + "Mean", mean));
245  featureList.push_back(std::make_pair(params.prefix + "Variance", variance));
246  featureList.push_back(std::make_pair(params.prefix + "Skewness", skewness));
247  featureList.push_back(std::make_pair(params.prefix + "Excess kurtosis", kurtosis-3));
248  featureList.push_back(std::make_pair(params.prefix + "Median", median));
249  featureList.push_back(std::make_pair(params.prefix + "Minimum", minimum));
250  featureList.push_back(std::make_pair(params.prefix + "05th Percentile", percentiles[0]));
251  featureList.push_back(std::make_pair(params.prefix + "10th Percentile", percentiles[1]));
252  featureList.push_back(std::make_pair(params.prefix + "15th Percentile", percentiles[2]));
253  featureList.push_back(std::make_pair(params.prefix + "20th Percentile", percentiles[3]));
254  featureList.push_back(std::make_pair(params.prefix + "25th Percentile", percentiles[4]));
255  featureList.push_back(std::make_pair(params.prefix + "30th Percentile", percentiles[5]));
256  featureList.push_back(std::make_pair(params.prefix + "35th Percentile", percentiles[6]));
257  featureList.push_back(std::make_pair(params.prefix + "40th Percentile", percentiles[7]));
258  featureList.push_back(std::make_pair(params.prefix + "45th Percentile", percentiles[8]));
259  featureList.push_back(std::make_pair(params.prefix + "50th Percentile", percentiles[9]));
260  featureList.push_back(std::make_pair(params.prefix + "55th Percentile", percentiles[10]));
261  featureList.push_back(std::make_pair(params.prefix + "60th Percentile", percentiles[11]));
262  featureList.push_back(std::make_pair(params.prefix + "65th Percentile", percentiles[12]));
263  featureList.push_back(std::make_pair(params.prefix + "70th Percentile", percentiles[13]));
264  featureList.push_back(std::make_pair(params.prefix + "75th Percentile", percentiles[14]));
265  featureList.push_back(std::make_pair(params.prefix + "80th Percentile", percentiles[15]));
266  featureList.push_back(std::make_pair(params.prefix + "85th Percentile", percentiles[16]));
267  featureList.push_back(std::make_pair(params.prefix + "90th Percentile", percentiles[17]));
268  featureList.push_back(std::make_pair(params.prefix + "95th Percentile", percentiles[18]));
269  featureList.push_back(std::make_pair(params.prefix + "Maximum", maximum));
270  featureList.push_back(std::make_pair(params.prefix + "Interquantile range", interquantileRange));
271  featureList.push_back(std::make_pair(params.prefix + "Range", maximum-minimum));
272  featureList.push_back(std::make_pair(params.prefix + "Mean absolute deviation", meanAbsoluteDeviation));
273  featureList.push_back(std::make_pair(params.prefix + "Robust mean absolute deviation", robustMeanAbsoluteDeviation));
274  featureList.push_back(std::make_pair(params.prefix + "Median absolute deviation", medianAbsoluteDeviation));
275  featureList.push_back(std::make_pair(params.prefix + "Coefficient of variation", coefficientOfVariation));
276  featureList.push_back(std::make_pair(params.prefix + "Quantile coefficient of dispersion", quantileCoefficientOfDispersion));
277  featureList.push_back(std::make_pair(params.prefix + "Energy", energy));
278  featureList.push_back(std::make_pair(params.prefix + "Root mean square", rootMeanSquare));
279 
280  featureList.push_back(std::make_pair(params.prefix + "Standard Deviation", std::sqrt(variance)));
281  featureList.push_back(std::make_pair(params.prefix + "Kurtosis", kurtosis));
282  featureList.push_back(std::make_pair(params.prefix + "Robust mean", robustMean));
283  featureList.push_back(std::make_pair(params.prefix + "Robust variance", robustVariance));
284  featureList.push_back(std::make_pair(params.prefix + "Covered image intensity range", coveredImageRange));
285  featureList.push_back(std::make_pair(params.prefix + "Mode index", modeIdx));
286  featureList.push_back(std::make_pair(params.prefix + "Mode value", params.quantifier->IndexToMeanIntensity(modeIdx)));
287  featureList.push_back(std::make_pair(params.prefix + "Mode probability", histogram[modeIdx] / (1.0*numberOfVoxels)));
288  featureList.push_back(std::make_pair(params.prefix + "Entropy", entropy));
289  featureList.push_back(std::make_pair(params.prefix + "Uniformtiy", uniformity));
290  featureList.push_back(std::make_pair(params.prefix + "Number of voxels", numberOfVoxels));
291  featureList.push_back(std::make_pair(params.prefix + "Sum of voxels", sum));
292  featureList.push_back(std::make_pair(params.prefix + "Voxel space", voxelSpace));
293  featureList.push_back(std::make_pair(params.prefix + "Voxel volume", voxelVolume));
294  featureList.push_back(std::make_pair(params.prefix + "Image Dimension", VImageDimension));
295 
296  return;
297 }
298 
300 {
301  SetShortName("fon");
302  SetLongName("first-order-numeric");
303  SetFeatureClassName("First Order Numeric");
304 }
305 
307 {
308  InitializeQuantifier(image, mask);
309  FeatureListType featureList;
310 
311  FirstOrderNumericParameterStruct params;
312 
313  params.quantifier = GetQuantifier();
314  params.MinimumIntensity = GetQuantifier()->GetMinimum();
315  params.MaximumIntensity = GetQuantifier()->GetMaximum();
316  params.Bins = GetQuantifier()->GetBins();
317  params.prefix = FeatureDescriptionPrefix();
318  AccessByItk_3(image, CalculateFirstOrderStatistics, mask, featureList, params);
319 
320  return featureList;
321 }
322 
324 {
325  FeatureNameListType featureList;
326  return featureList;
327 }
328 
329 
331 {
332  std::string name = GetOptionPrefix();
333 
334  parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use First Order Statistic (Numeric)", "calculates First Order Statistic (Numeric)", us::Any());
335  AddQuantifierArguments(parser);
336 }
337 
338 void
340 {
341  auto parsedArgs = GetParameter();
342  if (parsedArgs.count(GetLongName()))
343  {
344  InitializeQuantifierFromParameters(feature, maskNoNAN);
345  MITK_INFO << "Start calculating first order features ....";
346  auto localResults = this->CalculateFeatures(feature, maskNoNAN);
347  featureList.insert(featureList.end(), localResults.begin(), localResults.end());
348  MITK_INFO << "Finished calculating first order features....";
349  }
350 }
351 
353 {
354  return QuantifierParameterString();
355 }
std::string GetCurrentFeatureEncoding() override
Adds an additional Separator to the name of the feature, which encodes the used parameters.
#define AccessByItk_3(mitkImage, itkImageTypeFunction, arg1, arg2, arg3)
#define MITK_INFO
Definition: mitkLogMacros.h:18
itk::Image< unsigned char, 3 > ImageType
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...
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)
FeatureNameListType GetFeatureNames() override
Returns a list of the names of all features that are calculated from this class.
std::vector< std::pair< std::string, double > > FeatureListType
Definition: usAny.h:163
static T max(T x, T y)
Definition: svm.cpp:56
mitk::Image::Pointer image
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
void CalculateFirstOrderStatistics(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer mask, mitk::GIFFirstOrderNumericStatistics::FeatureListType &featureList, FirstOrderNumericParameterStruct params)
void CalculateFeatures(mitk::CoocurenceMatrixHolder &holder, mitk::CoocurenceMatrixFeatures &results)
void AddArguments(mitkCommandLineParser &parser) override
FeatureListType CalculateFeatures(const Image::Pointer &image, const Image::Pointer &feature) override
Calculates the First Order Features based on a binned version of the image.