Medical Imaging Interaction Toolkit  2018.4.99-6aa36ba9
Medical Imaging Interaction Toolkit
mitkGIFNeighbourhoodGreyToneDifferenceFeatures.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 <mitkPixelTypeMultiplex.h>
21 
22 // ITK
23 #include <itkImageRegionIteratorWithIndex.h>
24 #include <itkNeighborhoodIterator.h>
25 // STL
26 #include <limits>
27 
28 struct GIFNeighbourhoodGreyToneDifferenceParameter
29 {
30  int Range = 1;
32  std::string prefix;
33 };
34 
35 template<typename TPixel, unsigned int VImageDimension>
36 static void
37 CalculateIntensityPeak(itk::Image<TPixel, VImageDimension>* itkImage, mitk::Image::Pointer mask, GIFNeighbourhoodGreyToneDifferenceParameter params, mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::FeatureListType & featureList)
38 {
39  typedef itk::Image<TPixel, VImageDimension> ImageType;
40  typedef itk::Image<unsigned short, VImageDimension> MaskType;
41 
42  typename MaskType::Pointer itkMask = MaskType::New();
43  mitk::CastToItkImage(mask, itkMask);
44 
45  typename ImageType::SizeType regionSize;
46  regionSize.Fill(params.Range);
47 
48  itk::NeighborhoodIterator<ImageType> iter(regionSize, itkImage, itkImage->GetLargestPossibleRegion());
49  itk::NeighborhoodIterator<MaskType> iterMask(regionSize, itkMask, itkMask->GetLargestPossibleRegion());
50 
51  std::vector<double> pVector;
52  std::vector<double> sVector;
53  pVector.resize(params.quantifier->GetBins(), 0);
54  sVector.resize(params.quantifier->GetBins(), 0);
55 
56  int count = 0;
57  while (!iter.IsAtEnd())
58  {
59  if (iterMask.GetCenterPixel() > 0)
60  {
61  int localCount = 0;
62  double localMean = 0;
63  unsigned int localIndex = params.quantifier->IntensityToIndex(iter.GetCenterPixel());
64  for (itk::SizeValueType i = 0; i < iter.Size(); ++i)
65  {
66  if (i == (iter.Size() / 2))
67  continue;
68  if (iterMask.GetPixel(i) > 0)
69  {
70  ++localCount;
71  localMean += params.quantifier->IntensityToIndex(iter.GetPixel(i)) + 1;
72  }
73  }
74  if (localCount > 0)
75  {
76  localMean /= localCount;
77  }
78  localMean = std::abs<double>(localIndex + 1 - localMean);
79 
80  pVector[localIndex] += 1;
81  sVector[localIndex] += localMean;
82  ++count;
83 
84 
85  }
86  ++iterMask;
87  ++iter;
88  }
89 
90  unsigned int Ngp = 0;
91  for (unsigned int i = 0; i < params.quantifier->GetBins(); ++i)
92  {
93  if (pVector[i] > 0.1)
94  {
95  ++Ngp;
96  }
97  pVector[i] /= count;
98  }
99 
100  double sumS = 0;
101  double sumStimesP = 0;
102 
103  double contrastA = 0;
104  double busynessA = 0;
105  double complexity = 0;
106  double strengthA = 0;
107  for (unsigned int i = 0; i < params.quantifier->GetBins(); ++i)
108  {
109  sumS += sVector[i];
110  sumStimesP += pVector[i] * sVector[i];
111  for (unsigned int j = 0; j < params.quantifier->GetBins(); ++j)
112  {
113  double iMinusj = 1.0*i - 1.0*j;
114  contrastA += pVector[i] * pVector[j] * iMinusj*iMinusj;
115  if ((pVector[i] > 0) && (pVector[j] > 0))
116  {
117  busynessA += std::abs<double>((i + 1.0)*pVector[i] - (j + 1.0)*pVector[j]);
118  complexity += std::abs<double>(iMinusj)*(pVector[i] * sVector[i] + pVector[j] * sVector[j]) / (pVector[i] + pVector[j]);
119  strengthA += (pVector[i] + pVector[j])*iMinusj*iMinusj;
120  }
121  }
122  }
123  double coarsness = 1.0 / std::min<double>(sumStimesP, 1000000);
124  double contrast = 0;
125  double busyness = 0;
126  if (Ngp > 1)
127  {
128  contrast = contrastA / Ngp / (Ngp - 1) / count * sumS;
129  busyness = sumStimesP / busynessA;
130  }
131  complexity /= count;
132  double strength = 0;
133  if (sumS > 0)
134  {
135  strength = strengthA / sumS;
136  }
137 
138  std::string prefix = params.prefix;
139  featureList.push_back(std::make_pair(prefix + "Coarsness", coarsness));
140  featureList.push_back(std::make_pair(prefix + "Contrast", contrast));
141  featureList.push_back(std::make_pair(prefix + "Busyness", busyness));
142  featureList.push_back(std::make_pair(prefix + "Complexity", complexity));
143  featureList.push_back(std::make_pair(prefix + "Strength", strength));
144 }
145 
146 
148 m_Range(1)
149 {
150  SetLongName("neighbourhood-grey-tone-difference");
151  SetShortName("ngtd");
152  SetFeatureClassName("Neighbourhood Grey Tone Difference");
153 }
154 
155 mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::FeatureListType mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask)
156 {
157  FeatureListType featureList;
158 
160 
161  GIFNeighbourhoodGreyToneDifferenceParameter params;
162  params.Range = GetRange();
163  params.quantifier = GetQuantifier();
164  params.prefix = FeatureDescriptionPrefix();
165 
166  AccessByItk_3(image, CalculateIntensityPeak, mask, params, featureList);
167  return featureList;
168 }
169 
171 {
172  FeatureNameListType featureList;
173  return featureList;
174 }
175 
176 
178 {
179  AddQuantifierArguments(parser);
180  std::string name = GetOptionPrefix();
181 
182  parser.addArgument(GetLongName(), name, mitkCommandLineParser::Bool, "Use Neighbourhood Grey Tone Difference", "calculates Neighborhood Grey Tone based features", us::Any());
183  parser.addArgument(name + "::range", name + "::range", mitkCommandLineParser::Int, "Range for the local intensity", "Give the range that should be used for the local intensity in mm", us::Any());
184 
185 }
186 
187 void
189 {
190  InitializeQuantifierFromParameters(feature, mask);
191  std::string name = GetOptionPrefix();
192 
193  auto parsedArgs = GetParameter();
194  if (parsedArgs.count(GetLongName()))
195  {
196  MITK_INFO << "Start calculating Neighbourhood Grey Tone Difference features ....";
197  if (parsedArgs.count(name + "::range"))
198  {
199  int range = us::any_cast<int>(parsedArgs[name + "::range"]);
200  this->SetRange(range);
201  }
202  auto localResults = this->CalculateFeatures(feature, mask);
203  featureList.insert(featureList.end(), localResults.begin(), localResults.end());
204  MITK_INFO << "Finished calculating Neighbourhood Grey Tone Difference features....";
205  }
206 }
207 
209 {
210  std::ostringstream ss;
211  ss << m_Range;
212  std::string strRange = ss.str();
213  return QuantifierParameterString() + "_Range-" + ss.str();
214 }
215 
216 
217 
#define AccessByItk_3(mitkImage, itkImageTypeFunction, arg1, arg2, arg3)
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...
#define MITK_INFO
Definition: mitkLogMacros.h:18
itk::Image< unsigned char, 3 > ImageType
virtual void SetLongName(std::string _arg)
ValueType * any_cast(Any *operand)
Definition: usAny.h:377
FeatureListType CalculateFeatures(const Image::Pointer &image, const Image::Pointer &feature) override
Calculates the Cooccurence-Matrix based features for 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)
std::string FeatureDescriptionPrefix()
Returns a string that encodes the feature class name.
virtual std::string GetLongName() const
Definition: usAny.h:163
FeatureNameListType GetFeatureNames() override
Returns a list of the names of all features that are calculated from this class.
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)
mitkClassMacro(AbstractGlobalImageFeature, BaseData) typedef std typedef std::vector< std::string > FeatureNameListType
static void CalculateIntensityPeak(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer mask, GIFNeighbourhoodGreyToneDifferenceParameter params, mitk::GIFNeighbourhoodGreyToneDifferenceFeatures::FeatureListType &featureList)
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)
virtual ParameterTypes GetParameter() const