Medical Imaging Interaction Toolkit  2018.4.99-936b789b
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  mitk::FeatureID id;
34 };
35 
36 template<typename TPixel, unsigned int VImageDimension>
37 void
38 CalculateFirstOrderStatistics(const itk::Image<TPixel, VImageDimension>* itkImage, const mitk::Image* 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::ImageRegionConstIterator<ImageType> imageIter(itkImage, itkImage->GetLargestPossibleRegion());
75  itk::ImageRegionConstIterator<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(mitk::CreateFeatureID(params.id, "Mean"),mean));
245  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Variance"),variance));
246  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Skewness"),skewness));
247  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Excess kurtosis"),kurtosis-3));
248  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Median"),median));
249  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Minimum"),minimum));
250  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "05th Percentile"),percentiles[0]));
251  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "10th Percentile"),percentiles[1]));
252  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "15th Percentile"),percentiles[2]));
253  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "20th Percentile"),percentiles[3]));
254  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "25th Percentile"),percentiles[4]));
255  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "30th Percentile"),percentiles[5]));
256  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "35th Percentile"),percentiles[6]));
257  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "40th Percentile"),percentiles[7]));
258  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "45th Percentile"),percentiles[8]));
259  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "50th Percentile"),percentiles[9]));
260  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "55th Percentile"),percentiles[10]));
261  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "60th Percentile"),percentiles[11]));
262  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "65th Percentile"),percentiles[12]));
263  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "70th Percentile"), percentiles[13]));
264  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "75th Percentile"), percentiles[14]));
265  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "80th Percentile"), percentiles[15]));
266  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "85th Percentile"), percentiles[16]));
267  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "90th Percentile"), percentiles[17]));
268  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "95th Percentile"), percentiles[18]));
269  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Maximum"), maximum));
270  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Interquantile range"), interquantileRange));
271  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Range"), maximum-minimum));
272  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Mean absolute deviation"), meanAbsoluteDeviation));
273  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Robust mean absolute deviation"), robustMeanAbsoluteDeviation));
274  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Median absolute deviation"), medianAbsoluteDeviation));
275  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Coefficient of variation"), coefficientOfVariation));
276  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Quantile coefficient of dispersion"), quantileCoefficientOfDispersion));
277  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Energy"), energy));
278  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Root mean square"), rootMeanSquare));
279 
280  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Standard Deviation"), std::sqrt(variance)));
281  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Kurtosis"), kurtosis));
282  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Robust mean"), robustMean));
283  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Robust variance"), robustVariance));
284  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Covered image intensity range"), coveredImageRange));
285  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Mode index"), modeIdx));
286  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Mode value"), params.quantifier->IndexToMeanIntensity(modeIdx)));
287  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Mode probability"), histogram[modeIdx] / (1.0*numberOfVoxels)));
288  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Entropy"), entropy));
289  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Uniformtiy"), uniformity));
290  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Number of voxels"), numberOfVoxels));
291  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Sum of voxels"), sum));
292  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Voxel space"), voxelSpace));
293  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "Voxel volume"), voxelVolume));
294  featureList.push_back(std::make_pair(mitk::CreateFeatureID(params.id, "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  this->AddQuantifierArguments(parser);
309 
310  std::string name = this->GetOptionPrefix();
311  parser.addArgument(this->GetLongName(), name, mitkCommandLineParser::Bool, "Use first order statistics (numericaly) as feature", "calculates first order statistics based features", us::Any());
312 }
313 
315 {
316  FeatureListType featureList;
317 
318  this->InitializeQuantifier(image, mask);
319 
320  MITK_INFO << "Start calculating first order features ....";
321 
322  FirstOrderNumericParameterStruct params;
323  params.quantifier = GetQuantifier();
324  params.MinimumIntensity = GetQuantifier()->GetMinimum();
325  params.MaximumIntensity = GetQuantifier()->GetMaximum();
326  params.Bins = GetQuantifier()->GetBins();
327  params.id = this->CreateTemplateFeatureID();
328  AccessByItk_3(image, CalculateFirstOrderStatistics, mask, featureList, params);
329 
330  MITK_INFO << "Finished calculating first order features....";
331 
332  return featureList;
333 }
334 
336 {
337  return Superclass::CalculateFeatures(image, maskNoNAN);
338 }
339 
#define AccessByItk_3(mitkImage, itkImageTypeFunction, arg1, arg2, arg3)
FeatureListType CalculateFeatures(const Image *image, const Image *mask, const Image *maskNoNAN) override
#define MITK_INFO
Definition: mitkLogMacros.h:18
itk::Image< unsigned char, 3 > ImageType
MITKCLCORE_EXPORT FeatureID CreateFeatureID(FeatureID templateID, std::string name)
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 CalculateFirstOrderStatistics(const itk::Image< TPixel, VImageDimension > *itkImage, const mitk::Image *mask, mitk::GIFFirstOrderNumericStatistics::FeatureListType &featureList, FirstOrderNumericParameterStruct params)
Image class for storing images.
Definition: mitkImage.h:72
Definition: usAny.h:163
std::vector< std::pair< FeatureID, double > > FeatureListType
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 AddArguments(mitkCommandLineParser &parser) const override
FeatureListType DoCalculateFeatures(const Image *image, const Image *mask) override
void CalculateFeatures(mitk::CoocurenceMatrixHolder &holder, mitk::CoocurenceMatrixFeatures &results)