Medical Imaging Interaction Toolkit  2018.4.99-c0f884b2
Medical Imaging Interaction Toolkit
mitkGIFGreyLevelSizeZone.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 <itkImageRegionIteratorWithIndex.h>
22 
23 // STL
24 
25 namespace mitk
26 {
27  struct GreyLevelSizeZoneMatrixHolder
28  {
29  public:
30  GreyLevelSizeZoneMatrixHolder(double min, double max, int number, int maxSize);
31 
32  int IntensityToIndex(double intensity);
33  double IndexToMinIntensity(int index);
34  double IndexToMeanIntensity(int index);
35  double IndexToMaxIntensity(int index);
36 
37  double m_MinimumRange;
38  double m_MaximumRange;
39  double m_Stepsize;
40  int m_NumberOfBins;
41  int m_MaximumSize;
42  Eigen::MatrixXd m_Matrix;
43 
44  };
45 
46  struct GreyLevelSizeZoneFeatures
47  {
48  GreyLevelSizeZoneFeatures() :
49  SmallZoneEmphasis(0),
50  LargeZoneEmphasis(0),
51  LowGreyLevelEmphasis(0),
52  HighGreyLevelEmphasis(0),
53  SmallZoneLowGreyLevelEmphasis(0),
54  SmallZoneHighGreyLevelEmphasis(0),
55  LargeZoneLowGreyLevelEmphasis(0),
56  LargeZoneHighGreyLevelEmphasis(0),
57  GreyLevelNonUniformity(0),
58  GreyLevelNonUniformityNormalized(0),
59  ZoneSizeNonUniformity(0),
60  ZoneSizeNoneUniformityNormalized(0),
61  ZonePercentage(0),
62  GreyLevelMean(0),
63  GreyLevelVariance(0),
64  ZoneSizeMean(0),
65  ZoneSizeVariance(0),
66  ZoneSizeEntropy(0)
67  {
68  }
69 
70  public:
71  double SmallZoneEmphasis;
72  double LargeZoneEmphasis;
73  double LowGreyLevelEmphasis;
74  double HighGreyLevelEmphasis;
75  double SmallZoneLowGreyLevelEmphasis;
76  double SmallZoneHighGreyLevelEmphasis;
77  double LargeZoneLowGreyLevelEmphasis;
78  double LargeZoneHighGreyLevelEmphasis;
79  double GreyLevelNonUniformity;
80  double GreyLevelNonUniformityNormalized;
81  double ZoneSizeNonUniformity;
82  double ZoneSizeNoneUniformityNormalized;
83  double ZonePercentage;
84  double GreyLevelMean;
85  double GreyLevelVariance;
86  double ZoneSizeMean;
87  double ZoneSizeVariance;
88  double ZoneSizeEntropy;
89  };
90 }
91 
92 static
93 void MatrixFeaturesTo(mitk::GreyLevelSizeZoneFeatures features,
94  std::string prefix,
96 
97 
98 
99 mitk::GreyLevelSizeZoneMatrixHolder::GreyLevelSizeZoneMatrixHolder(double min, double max, int number, int maxSize) :
100  m_MinimumRange(min),
101  m_MaximumRange(max),
102  m_NumberOfBins(number),
103  m_MaximumSize(maxSize)
104 {
105  m_Matrix.resize(number, maxSize);
106  m_Matrix.fill(0);
107  m_Stepsize = (max - min) / (number);
108 }
109 
110 int mitk::GreyLevelSizeZoneMatrixHolder::IntensityToIndex(double intensity)
111 {
112  return std::floor((intensity - m_MinimumRange) / m_Stepsize);
113 }
114 
115 double mitk::GreyLevelSizeZoneMatrixHolder::IndexToMinIntensity(int index)
116 {
117  return m_MinimumRange + index * m_Stepsize;
118 }
119 double mitk::GreyLevelSizeZoneMatrixHolder::IndexToMeanIntensity(int index)
120 {
121  return m_MinimumRange + (index+0.5) * m_Stepsize;
122 }
123 double mitk::GreyLevelSizeZoneMatrixHolder::IndexToMaxIntensity(int index)
124 {
125  return m_MinimumRange + (index + 1) * m_Stepsize;
126 }
127 
128 template<typename TPixel, unsigned int VImageDimension>
129 static int
130 CalculateGlSZMatrix(itk::Image<TPixel, VImageDimension>* itkImage,
131  itk::Image<unsigned short, VImageDimension>* mask,
132  std::vector<itk::Offset<VImageDimension> > offsets,
133  bool estimateLargestRegion,
134  mitk::GreyLevelSizeZoneMatrixHolder &holder)
135 {
136  typedef itk::Image<TPixel, VImageDimension> ImageType;
137  typedef itk::Image<unsigned short, VImageDimension> MaskImageType;
138  typedef typename ImageType::IndexType IndexType;
139 
140  typedef itk::ImageRegionIteratorWithIndex<ImageType> ConstIterType;
141  typedef itk::ImageRegionIteratorWithIndex<MaskImageType> ConstMaskIterType;
142 
143  auto region = mask->GetLargestPossibleRegion();
144  typename MaskImageType::RegionType newRegion;
145  newRegion.SetSize(region.GetSize());
146  newRegion.SetIndex(region.GetIndex());
147 
148  ConstIterType imageIter(itkImage, itkImage->GetLargestPossibleRegion());
149  ConstMaskIterType maskIter(mask, mask->GetLargestPossibleRegion());
150 
151  typename MaskImageType::Pointer visitedImage = MaskImageType::New();
152  visitedImage->SetRegions(newRegion);
153  visitedImage->Allocate();
154  visitedImage->FillBuffer(0);
155 
156  int largestRegion = 0;
157 
158  while (!maskIter.IsAtEnd())
159  {
160  if (maskIter.Value() > 0 )
161  {
162  auto startIntensityIndex = holder.IntensityToIndex(imageIter.Value());
163  std::vector<IndexType> indices;
164  indices.push_back(maskIter.GetIndex());
165  unsigned int steps = 0;
166 
167  while (indices.size() > 0)
168  {
169  auto currentIndex = indices.back();
170  indices.pop_back();
171 
172  if (!region.IsInside(currentIndex))
173  {
174  continue;
175  }
176 
177  auto wasVisited = visitedImage->GetPixel(currentIndex);
178  auto newIntensityIndex = holder.IntensityToIndex(itkImage->GetPixel(currentIndex));
179  auto isInMask = mask->GetPixel(currentIndex);
180 
181  if ((isInMask > 0) &&
182  (newIntensityIndex == startIntensityIndex) &&
183  (wasVisited < 1))
184  {
185  ++steps;
186  visitedImage->SetPixel(currentIndex, 1);
187  for (auto offset : offsets)
188  {
189  auto newIndex = currentIndex + offset;
190  indices.push_back(newIndex);
191  newIndex = currentIndex - offset;
192  indices.push_back(newIndex);
193  }
194 
195  }
196  }
197  if (steps > 0)
198  {
199  largestRegion = std::max<int>(steps, largestRegion);
200  steps = std::min<unsigned int>(steps, holder.m_MaximumSize);
201  if (!estimateLargestRegion)
202  {
203  holder.m_Matrix(startIntensityIndex, steps - 1) += 1;
204  }
205  }
206  }
207  ++imageIter;
208  ++maskIter;
209  }
210  return largestRegion;
211 }
212 
213 static void CalculateFeatures(
214  mitk::GreyLevelSizeZoneMatrixHolder &holder,
215  mitk::GreyLevelSizeZoneFeatures & results
216  )
217 {
218  auto SgzMatrix = holder.m_Matrix;
219  auto pgzMatrix = holder.m_Matrix;
220  auto pgMatrix = holder.m_Matrix;
221  auto pzMatrix = holder.m_Matrix;
222 
223  double Ns = pgzMatrix.sum();
224  pgzMatrix /= Ns;
225  pgMatrix.rowwise().normalize();
226  pzMatrix.colwise().normalize();
227 
228  for (int i = 0; i < holder.m_NumberOfBins; ++i)
229  for (int j = 0; j < holder.m_NumberOfBins; ++j)
230  {
231  if (pgzMatrix(i, j) != pgzMatrix(i, j))
232  pgzMatrix(i, j) = 0;
233  if (pgMatrix(i, j) != pgMatrix(i, j))
234  pgMatrix(i, j) = 0;
235  if (pzMatrix(i, j) != pzMatrix(i, j))
236  pzMatrix(i, j) = 0;
237  }
238 
239  Eigen::VectorXd SgVector = SgzMatrix.rowwise().sum();
240  Eigen::VectorXd SzVector = SgzMatrix.colwise().sum();
241 
242  for (int j = 0; j < SzVector.size(); ++j)
243  {
244  results.SmallZoneEmphasis += SzVector(j) / (j + 1) / (j + 1);
245  results.LargeZoneEmphasis += SzVector(j) * (j + 1.0) * (j + 1.0);
246  results.ZoneSizeNonUniformity += SzVector(j) * SzVector(j);
247  results.ZoneSizeNoneUniformityNormalized += SzVector(j) * SzVector(j);
248  }
249  for (int i = 0; i < SgVector.size(); ++i)
250  {
251  results.LowGreyLevelEmphasis += SgVector(i) / (i + 1) / (i + 1);
252  results.HighGreyLevelEmphasis += SgVector(i) * (i + 1) * (i + 1);
253  results.GreyLevelNonUniformity += SgVector(i)*SgVector(i);
254  results.GreyLevelNonUniformityNormalized += SgVector(i)*SgVector(i);
255  }
256 
257  for (int i = 0; i < SgzMatrix.rows(); ++i)
258  {
259  for (int j = 0; j < SgzMatrix.cols(); ++j)
260  {
261  results.SmallZoneLowGreyLevelEmphasis += SgzMatrix(i, j) / (i + 1) / (i + 1) / (j + 1) / (j + 1);
262  results.SmallZoneHighGreyLevelEmphasis += SgzMatrix(i, j) * (i + 1) * (i + 1) / (j + 1) / (j + 1);
263  results.LargeZoneLowGreyLevelEmphasis += SgzMatrix(i, j) / (i + 1) / (i + 1) * (j + 1.0) * (j + 1.0);
264  results.LargeZoneHighGreyLevelEmphasis += SgzMatrix(i, j) * (i + 1) * (i + 1) * (j + 1.0) * (j + 1.0);
265  results.ZonePercentage += SgzMatrix(i, j)*(j + 1);
266 
267  results.GreyLevelMean += (i + 1)*pgzMatrix(i, j);
268  results.ZoneSizeMean += (j + 1)*pgzMatrix(i, j);
269  if (pgzMatrix(i, j) > 0)
270  results.ZoneSizeEntropy -= pgzMatrix(i, j) * std::log(pgzMatrix(i, j)) / std::log(2);
271  }
272  }
273 
274  for (int i = 0; i < SgzMatrix.rows(); ++i)
275  {
276  for (int j = 0; j < SgzMatrix.cols(); ++j)
277  {
278  results.GreyLevelVariance += (i + 1 - results.GreyLevelMean)*(i + 1 - results.GreyLevelMean)*pgzMatrix(i, j);
279  results.ZoneSizeVariance += (j + 1 - results.ZoneSizeMean)*(j + 1 - results.ZoneSizeMean)*pgzMatrix(i, j);
280  }
281  }
282 
283  results.SmallZoneEmphasis /= Ns;
284  results.LargeZoneEmphasis /= Ns;
285  results.LowGreyLevelEmphasis /= Ns;
286  results.HighGreyLevelEmphasis /= Ns;
287 
288  results.SmallZoneLowGreyLevelEmphasis /= Ns;
289  results.SmallZoneHighGreyLevelEmphasis /= Ns;
290  results.LargeZoneLowGreyLevelEmphasis /= Ns;
291  results.LargeZoneHighGreyLevelEmphasis /= Ns;
292  results.GreyLevelNonUniformity /= Ns;
293  results.GreyLevelNonUniformityNormalized /= Ns*Ns;
294  results.ZoneSizeNonUniformity /= Ns;
295  results.ZoneSizeNoneUniformityNormalized /= Ns*Ns;
296 
297  results.ZonePercentage = Ns / results.ZonePercentage;
298 }
299 
300 template<typename TPixel, unsigned int VImageDimension>
301 static void
303 {
304  typedef itk::Image<unsigned short, VImageDimension> MaskType;
305  typedef itk::Neighborhood<TPixel, VImageDimension > NeighborhoodType;
306  typedef itk::Offset<VImageDimension> OffsetType;
307 
309  double rangeMin = config.MinimumIntensity;
310  double rangeMax = config.MaximumIntensity;
311  int numberOfBins = config.Bins;
312 
313  typename MaskType::Pointer maskImage = MaskType::New();
314  mitk::CastToItkImage(mask, maskImage);
315 
316  //Find possible directions
317  std::vector < itk::Offset<VImageDimension> > offsetVector;
318  NeighborhoodType hood;
319  hood.SetRadius(1);
320  unsigned int centerIndex = hood.GetCenterNeighborhoodIndex();
321  OffsetType offset;
322  for (unsigned int d = 0; d < centerIndex; d++)
323  {
324  offset = hood.GetOffset(d);
325  bool useOffset = true;
326  for (unsigned int i = 0; i < VImageDimension; ++i)
327  {
328  if ((config.direction == i + 2) && offset[i] != 0)
329  {
330  useOffset = false;
331  }
332  }
333  if (useOffset)
334  {
335  offsetVector.push_back(offset);
336  MITK_INFO << offset;
337  }
338  }
339  if (config.direction == 1)
340  {
341  offsetVector.clear();
342  offset[0] = 0;
343  offset[1] = 0;
344  offset[2] = 1;
345  offsetVector.push_back(offset);
346  }
347 
348  std::vector<mitk::GreyLevelSizeZoneFeatures> resultVector;
349  mitk::GreyLevelSizeZoneMatrixHolder tmpHolder(rangeMin, rangeMax, numberOfBins, 3);
350  int largestRegion = CalculateGlSZMatrix<TPixel, VImageDimension>(itkImage, maskImage, offsetVector, true, tmpHolder);
351  mitk::GreyLevelSizeZoneMatrixHolder holderOverall(rangeMin, rangeMax, numberOfBins,largestRegion);
352  mitk::GreyLevelSizeZoneFeatures overallFeature;
353  CalculateGlSZMatrix<TPixel, VImageDimension>(itkImage, maskImage, offsetVector, false, holderOverall);
354  CalculateFeatures(holderOverall, overallFeature);
355 
356  MatrixFeaturesTo(overallFeature, config.prefix, featureList);
357 }
358 
359 
360 static
361 void MatrixFeaturesTo(mitk::GreyLevelSizeZoneFeatures features,
362  std::string prefix,
364 {
365  featureList.push_back(std::make_pair(prefix + "Small Zone Emphasis", features.SmallZoneEmphasis));
366  featureList.push_back(std::make_pair(prefix + "Large Zone Emphasis", features.LargeZoneEmphasis));
367  featureList.push_back(std::make_pair(prefix + "Low Grey Level Emphasis", features.LowGreyLevelEmphasis));
368  featureList.push_back(std::make_pair(prefix + "High Grey Level Emphasis", features.HighGreyLevelEmphasis));
369  featureList.push_back(std::make_pair(prefix + "Small Zone Low Grey Level Emphasis", features.SmallZoneLowGreyLevelEmphasis));
370  featureList.push_back(std::make_pair(prefix + "Small Zone High Grey Level Emphasis", features.SmallZoneHighGreyLevelEmphasis));
371  featureList.push_back(std::make_pair(prefix + "Large Zone Low Grey Level Emphasis", features.LargeZoneLowGreyLevelEmphasis));
372  featureList.push_back(std::make_pair(prefix + "Large Zone High Grey Level Emphasis", features.LargeZoneHighGreyLevelEmphasis));
373  featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity", features.GreyLevelNonUniformity));
374  featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity Normalized", features.GreyLevelNonUniformityNormalized));
375  featureList.push_back(std::make_pair(prefix + "Zone Size Non-Uniformity", features.ZoneSizeNonUniformity));
376  featureList.push_back(std::make_pair(prefix + "Zone Size Non-Uniformity Normalized", features.ZoneSizeNoneUniformityNormalized));
377  featureList.push_back(std::make_pair(prefix + "Zone Percentage", features.ZonePercentage));
378  featureList.push_back(std::make_pair(prefix + "Grey Level Mean", features.GreyLevelMean));
379  featureList.push_back(std::make_pair(prefix + "Grey Level Variance", features.GreyLevelVariance));
380  featureList.push_back(std::make_pair(prefix + "Zone Size Mean", features.ZoneSizeMean));
381  featureList.push_back(std::make_pair(prefix + "Zone Size Variance", features.ZoneSizeVariance));
382  featureList.push_back(std::make_pair(prefix + "Zone Size Entropy", features.ZoneSizeEntropy));
383 }
384 
386 {
387  SetShortName("glsz");
388  SetLongName("grey-level-sizezone");
389  SetFeatureClassName("Grey Level Size Zone");
390 }
391 
393 {
394  InitializeQuantifier(image, mask);
395 
396  FeatureListType featureList;
397 
399  config.direction = GetDirection();
400 
401  config.MinimumIntensity = GetQuantifier()->GetMinimum();
402  config.MaximumIntensity = GetQuantifier()->GetMaximum();
403  config.Bins = GetQuantifier()->GetBins();
404  config.prefix = FeatureDescriptionPrefix();
405 
406  AccessByItk_3(image, CalculateGreyLevelSizeZoneFeatures, mask, featureList, config);
407 
408  return featureList;
409 }
410 
412 {
413  FeatureNameListType featureList;
414  return featureList;
415 }
416 
417 
418 
419 
421 {
422  std::string name = GetOptionPrefix();
423 
424  parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Grey Level Size Zone", "Calculates the size zone based features.", us::Any());
425  AddQuantifierArguments(parser);
426 }
427 
428 void
430 {
431  auto parsedArgs = GetParameter();
432  std::string name = GetOptionPrefix();
433 
434  if (parsedArgs.count(GetLongName()))
435  {
436  InitializeQuantifierFromParameters(feature, maskNoNAN);
437 
438  MITK_INFO << "Start calculating Grey leve size zone ...";
439  auto localResults = this->CalculateFeatures(feature, maskNoNAN);
440  featureList.insert(featureList.end(), localResults.begin(), localResults.end());
441  MITK_INFO << "Finished calculating Grey level size zone ...";
442  }
443 
444 }
445 
447 {
448  return QuantifierParameterString();
449 }
450 
#define AccessByItk_3(mitkImage, itkImageTypeFunction, arg1, arg2, arg3)
FeatureListType CalculateFeatures(const Image::Pointer &image, const Image::Pointer &feature) override
Calculates the Cooccurence-Matrix based features for this class.
#define MITK_INFO
Definition: mitkLogMacros.h:18
itk::Image< unsigned char, 3 > ImageType
virtual void SetLongName(std::string _arg)
static void CalculateGreyLevelSizeZoneFeatures(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer mask, mitk::GIFGreyLevelSizeZone::FeatureListType &featureList, mitk::GIFGreyLevelSizeZone::GIFGreyLevelSizeZoneConfiguration config)
DataCollection - Class to facilitate loading/accessing structured data.
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.
virtual void SetShortName(std::string _arg)
std::vector< std::pair< std::string, double > > FeatureListType
static Vector3D offset
std::string FeatureDescriptionPrefix()
Returns a string that encodes the feature class name.
static void MatrixFeaturesTo(mitk::GreyLevelSizeZoneFeatures features, std::string prefix, mitk::GIFGreyLevelSizeZone::FeatureListType &featureList)
FeatureListType CalculateFeatures(const Image::Pointer &image, const Image::Pointer &feature) override
Calculates the Cooccurence-Matrix based features for this class.
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 AddArguments(mitkCommandLineParser &parser) override
virtual std::string GetLongName() const
Definition: usAny.h:163
itk::Image< unsigned char, 3 > MaskImageType
Definition: CLBrainMask.cpp:32
virtual int GetDirection() const
static T max(T x, T y)
Definition: svm.cpp:56
mitk::Image::Pointer image
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)
static int CalculateGlSZMatrix(itk::Image< TPixel, VImageDimension > *itkImage, itk::Image< unsigned short, VImageDimension > *mask, std::vector< itk::Offset< VImageDimension > > offsets, bool estimateLargestRegion, mitk::GreyLevelSizeZoneMatrixHolder &holder)
void InitializeQuantifier(const Image::Pointer &feature, const Image::Pointer &mask, unsigned int defaultBins=256)
virtual ParameterTypes GetParameter() const
std::string GetCurrentFeatureEncoding() override
Adds an additional Separator to the name of the feature, which encodes the used parameters.