Medical Imaging Interaction Toolkit  2018.4.99-07c45cb1
Medical Imaging Interaction Toolkit
mitkGIFGreyLevelDistanceZone.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 <mitkIOUtil.h>
20 #include <mitkITKImageImport.h>
21 
22 // ITK
23 #include <itkMinimumMaximumImageCalculator.h>
24 #include <itkImageRegionIteratorWithIndex.h>
25 #include <itkBinaryCrossStructuringElement.h>
26 #include <itkBinaryErodeImageFilter.h>
27 #include <itkAddImageFilter.h>
28 
29 namespace mitk{
30  struct GreyLevelDistanceZoneMatrixHolder
31  {
32  public:
33  GreyLevelDistanceZoneMatrixHolder(mitk::IntensityQuantifier::Pointer quantifier, int number, int maxSize);
34 
35  int IntensityToIndex(double intensity);
36 
37  int m_NumberOfBins;
38  int m_MaximumSize;
39  int m_NumerOfVoxels;
40  Eigen::MatrixXd m_Matrix;
42 
43  };
44 }
45 
46 static
48  std::string prefix,
50 
51 
52 
53 mitk::GreyLevelDistanceZoneMatrixHolder::GreyLevelDistanceZoneMatrixHolder(mitk::IntensityQuantifier::Pointer quantifier, int number, int maxSize) :
54  m_NumberOfBins(number),
55  m_MaximumSize(maxSize),
56  m_NumerOfVoxels(0),
57  m_Quantifier(quantifier)
58 {
59  m_Matrix.resize(number, maxSize);
60  m_Matrix.fill(0);
61 }
62 
63 int mitk::GreyLevelDistanceZoneMatrixHolder::IntensityToIndex(double intensity)
64 {
65  return m_Quantifier->IntensityToIndex(intensity);
66 }
67 
68 
69 template<typename TPixel, unsigned int VImageDimension>
70 int
71 CalculateGlSZMatrix(itk::Image<TPixel, VImageDimension>* itkImage,
72  itk::Image<unsigned short, VImageDimension>* mask,
73  itk::Image<unsigned short, VImageDimension>* distanceImage,
74  std::vector<itk::Offset<VImageDimension> > offsets,
75  bool estimateLargestRegion,
76  mitk::GreyLevelDistanceZoneMatrixHolder &holder)
77 {
78  typedef itk::Image<TPixel, VImageDimension> ImageType;
79  typedef itk::Image<unsigned short, VImageDimension> MaskImageType;
80  typedef typename ImageType::IndexType IndexType;
81 
82  typedef itk::ImageRegionIteratorWithIndex<ImageType> ConstIterType;
83  typedef itk::ImageRegionIteratorWithIndex<MaskImageType> ConstMaskIterType;
84 
85  auto region = mask->GetLargestPossibleRegion();
86  typename MaskImageType::RegionType newRegion;
87  newRegion.SetSize(region.GetSize());
88  newRegion.SetIndex(region.GetIndex());
89 
90  ConstIterType imageIter(itkImage, itkImage->GetLargestPossibleRegion());
91  ConstMaskIterType maskIter(mask, mask->GetLargestPossibleRegion());
92 
93  typename MaskImageType::Pointer visitedImage = MaskImageType::New();
94  visitedImage->SetRegions(newRegion);
95  visitedImage->Allocate();
96  visitedImage->FillBuffer(0);
97 
98  int largestRegion = 0;
99  holder.m_NumberOfBins = 0;
100 
101  while (!maskIter.IsAtEnd())
102  {
103  if (maskIter.Value() > 0 )
104  {
105  auto startIntensityIndex = holder.IntensityToIndex(imageIter.Value());
106  std::vector<IndexType> indices;
107  indices.push_back(maskIter.GetIndex());
108  unsigned int steps = 0;
109  int smallestDistance = 500;
110 
111  while (indices.size() > 0)
112  {
113  auto currentIndex = indices.back();
114  indices.pop_back();
115 
116  if (!region.IsInside(currentIndex))
117  {
118  continue;
119  }
120 
121  auto wasVisited = visitedImage->GetPixel(currentIndex);
122  auto newIntensityIndex = holder.IntensityToIndex(itkImage->GetPixel(currentIndex));
123  auto isInMask = mask->GetPixel(currentIndex);
124 
125  if ((isInMask > 0) &&
126  (newIntensityIndex == startIntensityIndex) &&
127  (wasVisited < 1))
128  {
129  ++(holder.m_NumerOfVoxels);
130  smallestDistance = (smallestDistance > distanceImage->GetPixel(currentIndex)) ? distanceImage->GetPixel(currentIndex) : smallestDistance;
131  ++steps;
132  visitedImage->SetPixel(currentIndex, 1);
133  for (auto offset : offsets)
134  {
135  auto newIndex = currentIndex + offset;
136  indices.push_back(newIndex);
137  newIndex = currentIndex - offset;
138  indices.push_back(newIndex);
139  }
140 
141  }
142  }
143  if (steps > 0)
144  {
145  largestRegion = std::max<int>(steps, largestRegion);
146  steps = std::min<unsigned int>(steps, holder.m_MaximumSize);
147  if (!estimateLargestRegion)
148  {
149  holder.m_Matrix(startIntensityIndex, smallestDistance-1) += 1;
150  }
151  }
152  }
153  ++imageIter;
154  ++maskIter;
155  }
156  return largestRegion;
157 }
158 
159 template <typename TPixel, unsigned int VDimension>
161  itk::Image<TPixel, VDimension> *sourceImage,
162  mitk::Image::Pointer &resultImage,
163  int &maxDistance)
164 {
165  typedef itk::Image<TPixel, VDimension> ImageType;
166  typedef unsigned short MaskType;
167  typedef itk::Image<MaskType, VDimension> MaskImageType;
168 
169  typename MaskImageType::Pointer distanceImage = MaskImageType::New();
170  distanceImage->SetRegions(sourceImage->GetLargestPossibleRegion());
171  distanceImage->SetOrigin(sourceImage->GetOrigin());
172  distanceImage->SetSpacing(sourceImage->GetSpacing());
173  distanceImage->SetDirection(sourceImage->GetDirection());
174  distanceImage->Allocate();
175  distanceImage->FillBuffer(std::numeric_limits<MaskType>::max()-1);
176 
177  typename ImageType::SizeType radius;
178  radius.Fill(1);
179  itk::NeighborhoodIterator<ImageType> neighbourIter(radius, sourceImage, sourceImage->GetLargestPossibleRegion());
180  itk::NeighborhoodIterator<MaskImageType> distanceIter(radius, distanceImage, distanceImage->GetLargestPossibleRegion());
181 
182  bool imageChanged = true;
183  while (imageChanged)
184  {
185  imageChanged = false;
186  maxDistance = 0;
187  neighbourIter.GoToBegin();
188  distanceIter.GoToBegin();
189  while (!neighbourIter.IsAtEnd())
190  {
191  MaskType oldDistance = distanceIter.GetCenterPixel();
192  maxDistance = std::max<MaskType>(maxDistance, oldDistance);
193  if (neighbourIter.GetCenterPixel() < 1)
194  {
195  if (oldDistance > 0)
196  {
197  distanceIter.SetCenterPixel(0);
198  imageChanged = true;
199  }
200  }
201  else if (oldDistance>0) {
202  MaskType minimumDistance = oldDistance;
203  for (unsigned int i = 0; i < distanceIter.Size(); ++i)
204  {
205  minimumDistance = std::min<MaskType>(minimumDistance, 1+distanceIter.GetPixel(i));
206  }
207  if (minimumDistance != oldDistance)
208  {
209  distanceIter.SetCenterPixel(minimumDistance);
210  imageChanged = true;
211  }
212  }
213 
214  ++neighbourIter;
215  ++distanceIter;
216  }
217  }
218 
219  mitk::CastToMitkImage(distanceImage, resultImage);
220 }
221 
222 void erode(mitk::Image::Pointer input, mitk::Image::Pointer &output, int &maxDistance)
223 {
224  AccessByItk_2(input, itkErode2, output, maxDistance);
225 }
226 
227 
228 void erodeAndAdd(mitk::Image::Pointer input, mitk::Image::Pointer& finalOutput, int &maxDistance)
229 {
230  maxDistance = 0;
231  erode(input, finalOutput, maxDistance);
232 }
233 
234 
235 void static CalculateFeatures(
236  mitk::GreyLevelDistanceZoneMatrixHolder &holder,
238  )
239 {
240  auto SgzMatrix = holder.m_Matrix;
241  auto pgzMatrix = holder.m_Matrix;
242  auto pgMatrix = holder.m_Matrix;
243  auto pzMatrix = holder.m_Matrix;
244 
245  double Ns = pgzMatrix.sum();
246  pgzMatrix /= Ns;
247  pgMatrix.rowwise().normalize();
248  pzMatrix.colwise().normalize();
249 
250  for (int i = 0; i < pgzMatrix.rows(); ++i)
251  for (int j = 0; j < pgzMatrix.cols(); ++j)
252  {
253  if (pgzMatrix(i, j) != pgzMatrix(i, j))
254  pgzMatrix(i, j) = 0;
255  if (pgMatrix(i, j) != pgMatrix(i, j))
256  pgMatrix(i, j) = 0;
257  if (pzMatrix(i, j) != pzMatrix(i, j))
258  pzMatrix(i, j) = 0;
259  }
260 
261  Eigen::VectorXd SgVector = SgzMatrix.rowwise().sum();
262  Eigen::VectorXd SzVector = SgzMatrix.colwise().sum();
263 
264  for (int j = 0; j < SzVector.size(); ++j)
265  {
266  results.SmallDistanceEmphasis += SzVector(j) / (j+1) / (j+1);
267  results.LargeDistanceEmphasis += SzVector(j) * (j + 1.0) * (j + 1.0);
268  results.ZoneDistanceNonUniformity += SzVector(j) * SzVector(j);
269  results.ZoneDistanceNoneUniformityNormalized += SzVector(j) * SzVector(j);
270  }
271  for (int i = 0; i < SgVector.size(); ++i)
272  {
273  results.LowGreyLevelEmphasis += SgVector(i) / (i + 1) / (i + 1);
274  results.HighGreyLevelEmphasis += SgVector(i) * (i + 1) * (i + 1);
275  results.GreyLevelNonUniformity += SgVector(i)*SgVector(i);
276  results.GreyLevelNonUniformityNormalized += SgVector(i)*SgVector(i);
277  }
278 
279  for (int i = 0; i < SgzMatrix.rows(); ++i)
280  {
281  for (int j = 0; j < SgzMatrix.cols(); ++j)
282  {
283  results.SmallDistanceLowGreyLevelEmphasis += SgzMatrix(i, j) / (i + 1) / (i + 1) / (j + 1) / (j + 1);
284  results.SmallDistanceHighGreyLevelEmphasis += SgzMatrix(i, j) * (i + 1) * (i + 1) / (j + 1) / (j + 1);
285  results.LargeDistanceLowGreyLevelEmphasis += SgzMatrix(i, j) / (i + 1) / (i + 1) * (j + 1.0) * (j + 1.0);
286  results.LargeDistanceHighGreyLevelEmphasis += SgzMatrix(i, j) * (i + 1) * (i + 1) * (j + 1.0) * (j + 1.0);
287  results.ZonePercentage += SgzMatrix(i, j);
288 
289  results.GreyLevelMean += (i + 1)*pgzMatrix(i, j);
290  results.ZoneDistanceMean += (j + 1)*pgzMatrix(i, j);
291  if (pgzMatrix(i, j) > 0)
292  results.ZoneDistanceEntropy -= pgzMatrix(i, j) * std::log(pgzMatrix(i, j)) / std::log(2);
293  }
294  }
295 
296  for (int i = 0; i < SgzMatrix.rows(); ++i)
297  {
298  for (int j = 0; j < SgzMatrix.cols(); ++j)
299  {
300  results.GreyLevelVariance += (i + 1 - results.GreyLevelMean)*(i + 1 - results.GreyLevelMean)*pgzMatrix(i, j);
301  results.ZoneDistanceVariance += (j + 1 - results.ZoneDistanceMean)*(j + 1 - results.ZoneDistanceMean)*pgzMatrix(i, j);
302  }
303  }
304 
305  results.SmallDistanceEmphasis /= Ns;
306  results.LargeDistanceEmphasis /= Ns;
307  results.LowGreyLevelEmphasis /= Ns;
308  results.HighGreyLevelEmphasis /= Ns;
309 
310  results.SmallDistanceLowGreyLevelEmphasis /= Ns;
312  results.LargeDistanceLowGreyLevelEmphasis /= Ns;
314  results.GreyLevelNonUniformity /= Ns;
315  results.GreyLevelNonUniformityNormalized /= Ns*Ns;
316  results.ZoneDistanceNonUniformity /= Ns;
317  results.ZoneDistanceNoneUniformityNormalized /= Ns*Ns;
318 
319  results.ZonePercentage = Ns / holder.m_NumerOfVoxels;// results.ZonePercentage;
320 }
321 
322 template<typename TPixel, unsigned int VImageDimension>
323 static void
325 {
326  typedef itk::Image<unsigned short, VImageDimension> MaskType;
327  typedef itk::Neighborhood<TPixel, VImageDimension > NeighborhoodType;
328  typedef itk::Offset<VImageDimension> OffsetType;
329 
331  int maximumDistance = 0;
332  mitk::Image::Pointer mitkDistanceImage = mitk::Image::New();
333  erodeAndAdd(config.distanceMask, mitkDistanceImage, maximumDistance);
334  typename MaskType::Pointer distanceImage = MaskType::New();
335  mitk::CastToItkImage(mitkDistanceImage, distanceImage);
336 
337  typename MaskType::Pointer maskImage = MaskType::New();
338  mitk::CastToItkImage(mask, maskImage);
339 
340  //Find possible directions
341  std::vector < itk::Offset<VImageDimension> > offsetVector;
342  NeighborhoodType hood;
343  hood.SetRadius(1);
344  unsigned int centerIndex = hood.GetCenterNeighborhoodIndex();
345  OffsetType offset;
346  for (unsigned int d = 0; d < centerIndex; d++)
347  {
348  offset = hood.GetOffset(d);
349  bool useOffset = true;
350  for (unsigned int i = 0; i < VImageDimension; ++i)
351  {
352  if ((config.direction == i + 2) && offset[i] != 0)
353  {
354  useOffset = false;
355  }
356  }
357  if (useOffset)
358  {
359  offsetVector.push_back(offset);
360  }
361  }
362  if (config.direction == 1)
363  {
364  offsetVector.clear();
365  offset[0] = 0;
366  offset[1] = 0;
367  offset[2] = 1;
368  offsetVector.push_back(offset);
369  }
370 
371  MITK_INFO << "Maximum Distance: " << maximumDistance;
372  std::vector<mitk::GreyLevelDistanceZoneFeatures> resultVector;
373  mitk::GreyLevelDistanceZoneMatrixHolder holderOverall(config.Quantifier, config.Bins, maximumDistance + 1);
375  CalculateGlSZMatrix<TPixel, VImageDimension>(itkImage, maskImage, distanceImage, offsetVector, false, holderOverall);
376  CalculateFeatures(holderOverall, overallFeature);
377 
378  MatrixFeaturesTo(overallFeature, config.prefix, featureList);
379 }
380 
381 
382 static
384  std::string prefix,
386 {
387  featureList.push_back(std::make_pair(prefix + "Small Distance Emphasis", features.SmallDistanceEmphasis));
388  featureList.push_back(std::make_pair(prefix + "Large Distance Emphasis", features.LargeDistanceEmphasis));
389  featureList.push_back(std::make_pair(prefix + "Low Grey Level Emphasis", features.LowGreyLevelEmphasis));
390  featureList.push_back(std::make_pair(prefix + "High Grey Level Emphasis", features.HighGreyLevelEmphasis));
391  featureList.push_back(std::make_pair(prefix + "Small Distance Low Grey Level Emphasis", features.SmallDistanceLowGreyLevelEmphasis));
392  featureList.push_back(std::make_pair(prefix + "Small Distance High Grey Level Emphasis", features.SmallDistanceHighGreyLevelEmphasis));
393  featureList.push_back(std::make_pair(prefix + "Large Distance Low Grey Level Emphasis", features.LargeDistanceLowGreyLevelEmphasis));
394  featureList.push_back(std::make_pair(prefix + "Large Distance High Grey Level Emphasis", features.LargeDistanceHighGreyLevelEmphasis));
395  featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity", features.GreyLevelNonUniformity));
396  featureList.push_back(std::make_pair(prefix + "Grey Level Non-Uniformity Normalized", features.GreyLevelNonUniformityNormalized));
397  featureList.push_back(std::make_pair(prefix + "Distance Size Non-Uniformity", features.ZoneDistanceNonUniformity));
398  featureList.push_back(std::make_pair(prefix + "Distance Size Non-Uniformity Normalized", features.ZoneDistanceNoneUniformityNormalized));
399  featureList.push_back(std::make_pair(prefix + "Zone Percentage", features.ZonePercentage));
400  featureList.push_back(std::make_pair(prefix + "Grey Level Mean", features.GreyLevelMean));
401  featureList.push_back(std::make_pair(prefix + "Grey Level Variance", features.GreyLevelVariance));
402  featureList.push_back(std::make_pair(prefix + "Zone Distance Mean", features.ZoneDistanceMean));
403  featureList.push_back(std::make_pair(prefix + "Zone Distance Variance", features.ZoneDistanceVariance));
404  featureList.push_back(std::make_pair(prefix + "Zone Distance Entropy", features.ZoneDistanceEntropy));
405  featureList.push_back(std::make_pair(prefix + "Grey Level Entropy", features.ZoneDistanceEntropy));
406 }
407 
409 {
410  SetShortName("gldz");
411  SetLongName("distance-zone");
412  SetFeatureClassName("Grey Level Distance Zone");
413 }
414 
416 {
417  InitializeQuantifier(image, mask);
418  FeatureListType featureList;
419 
421  config.direction = GetDirection();
422 
423  if (GetMorphMask().IsNull())
424  {
425  config.distanceMask = mask->Clone();
426  }
427  else
428  {
429  config.distanceMask = GetMorphMask();
430  }
431 
432  config.MinimumIntensity = GetQuantifier()->GetMinimum();
433  config.MaximumIntensity = GetQuantifier()->GetMaximum();
434  config.Bins = GetQuantifier()->GetBins();
435  config.prefix = FeatureDescriptionPrefix();
436  config.Quantifier = GetQuantifier();
437 
438  AccessByItk_3(image, CalculateGreyLevelDistanceZoneFeatures, mask, featureList, config);
439 
440  return featureList;
441 }
442 
444 {
445  FeatureNameListType featureList;
446  return featureList;
447 }
448 
449 
450 
451 
453 {
454  std::string name = GetOptionPrefix();
455 
456  parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Grey Level Distance Zone", "Calculates the size zone based features.", us::Any());
457  AddQuantifierArguments(parser);
458 }
459 
460 void
462 {
463  auto parsedArgs = GetParameter();
464  std::string name = GetOptionPrefix();
465 
466  if (parsedArgs.count(GetLongName()))
467  {
468  InitializeQuantifierFromParameters(feature, mask);
469 
470  MITK_INFO << "Start calculating Grey Level Distance Zone ....";
471  auto localResults = this->CalculateFeatures(feature, maskNoNAN);
472  featureList.insert(featureList.end(), localResults.begin(), localResults.end());
473  MITK_INFO << "Finished calculating Grey Level Distance Zone.";
474  }
475 
476 }
477 
479 {
480  return QuantifierParameterString();
481 }
482 
#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 erodeAndAdd(mitk::Image::Pointer input, mitk::Image::Pointer &finalOutput, int &maxDistance)
void erode(mitk::Image::Pointer input, mitk::Image::Pointer &output, int &maxDistance)
static void MatrixFeaturesTo(mitk::GreyLevelDistanceZoneFeatures features, std::string prefix, mitk::GIFGreyLevelDistanceZone::FeatureListType &featureList)
DataCollection - Class to facilitate loading/accessing structured data.
FeatureNameListType GetFeatureNames() override
Returns a list of the names of all features that are calculated from this class.
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)
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...
std::vector< std::pair< std::string, double > > FeatureListType
int CalculateGlSZMatrix(itk::Image< TPixel, VImageDimension > *itkImage, itk::Image< unsigned short, VImageDimension > *mask, itk::Image< unsigned short, VImageDimension > *distanceImage, std::vector< itk::Offset< VImageDimension > > offsets, bool estimateLargestRegion, mitk::GreyLevelDistanceZoneMatrixHolder &holder)
static Vector3D offset
std::string FeatureDescriptionPrefix()
Returns a string that encodes the feature class name.
virtual mitk::Image::Pointer GetMorphMask() const
void itkErode2(itk::Image< TPixel, VDimension > *sourceImage, mitk::Image::Pointer &resultImage, int &maxDistance)
virtual std::string GetLongName() const
Definition: usAny.h:163
itk::Image< unsigned char, 3 > MaskImageType
Definition: CLBrainMask.cpp:32
virtual int GetDirection() const
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
static Pointer New()
void AddQuantifierArguments(mitkCommandLineParser &parser)
void CastToMitkImage(const itk::SmartPointer< ItkOutputImageType > &itkimage, itk::SmartPointer< mitk::Image > &mitkoutputimage)
Cast an itk::Image (with a specific type) to an mitk::Image.
Definition: mitkImageCast.h:74
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
FeatureListType CalculateFeatures(const Image::Pointer &image, const Image::Pointer &feature) override
Calculates the Cooccurence-Matrix based features for this class.
#define AccessByItk_2(mitkImage, itkImageTypeFunction, arg1, arg2)
virtual IntensityQuantifier::Pointer GetQuantifier()
void InitializeQuantifierFromParameters(const Image::Pointer &feature, const Image::Pointer &mask, unsigned int defaultBins=256)
static void CalculateGreyLevelDistanceZoneFeatures(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer mask, mitk::GIFGreyLevelDistanceZone::FeatureListType &featureList, mitk::GIFGreyLevelDistanceZone::GIFGreyLevelDistanceZoneConfiguration config)
void AddArguments(mitkCommandLineParser &parser) override
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.