OTB  9.0.0
Orfeo Toolbox
otbHaralickTexturesImageFunction.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbHaralickTexturesImageFunction_hxx
22 #define otbHaralickTexturesImageFunction_hxx
23 
25 #include "itkImageRegionIteratorWithIndex.h"
26 #include "itkConstNeighborhoodIterator.h"
27 #include "itkNumericTraits.h"
28 #include <complex>
29 #include <cmath>
30 
31 namespace otb
32 {
33 
37 template <class TInputImage, class TCoordRep>
39  : m_NeighborhoodRadius(10), m_Offset(), m_NumberOfBinsPerAxis(8), m_InputImageMinimum(0), m_InputImageMaximum(255)
40 {
41 }
42 
43 template <class TInputImage, class TCoordRep>
44 void HaralickTexturesImageFunction<TInputImage, TCoordRep>::PrintSelf(std::ostream& os, itk::Indent indent) const
45 {
46  this->Superclass::PrintSelf(os, indent);
47  os << indent << " Neighborhood radius value : " << m_NeighborhoodRadius << std::endl;
48  os << indent << " Input image minimum value : " << m_InputImageMinimum << std::endl;
49  os << indent << " Input Image maximum value : " << m_InputImageMaximum << std::endl;
50  os << indent << " Number of bins per axis : " << m_NumberOfBinsPerAxis << std::endl;
51  os << indent << " Offset : " << m_Offset << std::endl;
52 }
53 
54 template <class TInputImage, class TCoordRep>
57 {
58  // Build textures vector
59  OutputType textures;
60 
61  // Initialize textures
62  textures.Fill(itk::NumericTraits<ScalarRealType>::Zero);
63 
64  // Check for input image
65  if (!this->GetInputImage())
66  {
67  return textures;
68  }
69 
70  // Check for out of buffer
71  if (!this->IsInsideBuffer(index))
72  {
73  return textures;
74  }
75 
76  const double log2 = std::log(2.0);
77  // Retrieve the input pointer
78  InputImagePointerType inputPtr = const_cast<InputImageType*>(this->GetInputImage());
79 
80  // Compute the region on which co-occurence will be estimated
81  typename InputRegionType::IndexType inputIndex;
82  typename InputRegionType::SizeType inputSize;
83 
84  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
85  {
86  inputIndex[dim] = index[dim] - m_NeighborhoodRadius;
87  inputSize[dim] = 2 * m_NeighborhoodRadius + 1;
88  }
89 
90  // Build the input region
91  InputRegionType inputRegion;
92  inputRegion.SetIndex(inputIndex);
93  inputRegion.SetSize(inputSize);
94  inputRegion.Crop(inputPtr->GetRequestedRegion());
95 
96  CooccurrenceIndexedListPointerType GLCIList = CooccurrenceIndexedListType::New();
97  GLCIList->Initialize(m_NumberOfBinsPerAxis, m_InputImageMinimum, m_InputImageMaximum);
98 
99  // Next, find the minimum radius that encloses all the offsets.
100  unsigned int minRadius = 0;
101  for (unsigned int i = 0; i < m_Offset.GetOffsetDimension(); i++)
102  {
103  unsigned int distance = std::abs(m_Offset[i]);
104  if (distance > minRadius)
105  {
106  minRadius = distance;
107  }
108  }
109  SizeType radius;
110  radius.Fill(minRadius);
111 
112 
113  typedef itk::ConstNeighborhoodIterator<InputImageType> NeighborhoodIteratorType;
114  NeighborhoodIteratorType neighborIt;
115  neighborIt = NeighborhoodIteratorType(radius, inputPtr, inputRegion);
116  for (neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt)
117  {
118  const InputPixelType centerPixelIntensity = neighborIt.GetCenterPixel();
119  bool pixelInBounds;
120  const InputPixelType pixelIntensity = neighborIt.GetPixel(m_Offset, pixelInBounds);
121  if (!pixelInBounds)
122  {
123  continue; // don't put a pixel in the value if it's out-of-bounds.
124  }
125  GLCIList->AddPixelPair(centerPixelIntensity, pixelIntensity);
126  }
127 
128 
129  double pixelMean = 0.;
130  double marginalMean;
131  double marginalDevSquared = 0.;
132  double pixelVariance = 0.;
133 
134  // Create and Initialize marginalSums
135  std::vector<double> marginalSums(m_NumberOfBinsPerAxis, 0);
136 
137  // get co-occurrence vector and totalfrequency
138  VectorType glcVector = GLCIList->GetVector();
139  double totalFrequency = static_cast<double>(GLCIList->GetTotalFrequency());
140 
141  // Normalize the co-occurrence indexed list and compute mean, marginalSum
142  typename VectorType::iterator it = glcVector.begin();
143  while (it != glcVector.end())
144  {
145  double frequency = (*it).second / totalFrequency;
146  CooccurrenceIndexType index2 = (*it).first;
147  pixelMean += index2[0] * frequency;
148  marginalSums[index2[0]] += frequency;
149  ++it;
150  }
151 
152  /* Now get the mean and deviaton of the marginal sums.
153  Compute incremental mean and SD, a la Knuth, "The Art of Computer
154  Programming, Volume 2: Seminumerical Algorithms", section 4.2.2.
155  Compute mean and standard deviation using the recurrence relation:
156  M(1) = x(1), M(k) = M(k-1) + (x(k) - M(k-1) ) / k
157  S(1) = 0, S(k) = S(k-1) + (x(k) - M(k-1)) * (x(k) - M(k))
158  for 2 <= k <= n, then
159  sigma = std::sqrt(S(n) / n) (or divide by n-1 for sample SD instead of
160  population SD).
161  */
162  std::vector<double>::const_iterator msIt = marginalSums.begin();
163  marginalMean = *msIt;
164  // Increment iterator to start with index 1
165  ++msIt;
166  for (int k = 2; msIt != marginalSums.end(); ++k, ++msIt)
167  {
168  double M_k_minus_1 = marginalMean;
169  double S_k_minus_1 = marginalDevSquared;
170  double x_k = *msIt;
171  double M_k = M_k_minus_1 + (x_k - M_k_minus_1) / k;
172  double S_k = S_k_minus_1 + (x_k - M_k_minus_1) * (x_k - M_k);
173  marginalMean = M_k;
174  marginalDevSquared = S_k;
175  }
176  marginalDevSquared = marginalDevSquared / m_NumberOfBinsPerAxis;
177 
178  VectorConstIteratorType constVectorIt;
179  constVectorIt = glcVector.begin();
180  while (constVectorIt != glcVector.end())
181  {
182  RelativeFrequencyType frequency = (*constVectorIt).second / totalFrequency;
183  CooccurrenceIndexType index2 = (*constVectorIt).first;
184  pixelVariance += (index2[0] - pixelMean) * (index2[0] - pixelMean) * frequency;
185  ++constVectorIt;
186  }
187 
188  double pixelVarianceSquared = pixelVariance * pixelVariance;
189  // Variance is only used in correlation. If variance is 0, then (index[0] - pixelMean) * (index[1] - pixelMean)
190  // should be zero as well. In this case, set the variance to 1. in order to
191  // avoid NaN correlation.
192  if (pixelVarianceSquared < GetPixelValueTolerance())
193  {
194  pixelVarianceSquared = 1.;
195  }
196 
197  // Compute textures
198  constVectorIt = glcVector.begin();
199  while (constVectorIt != glcVector.end())
200  {
201  CooccurrenceIndexType index2 = (*constVectorIt).first;
202  RelativeFrequencyType frequency = (*constVectorIt).second / totalFrequency;
203  textures[0] += frequency * frequency;
204  textures[1] -= (frequency > GetPixelValueTolerance()) ? frequency * std::log(frequency) / log2 : 0;
205  textures[2] += ((index2[0] - pixelMean) * (index2[1] - pixelMean) * frequency) / pixelVarianceSquared;
206  textures[3] += frequency / (1.0 + (index2[0] - index2[1]) * (index2[0] - index2[1]));
207  textures[4] += (index2[0] - index2[1]) * (index2[0] - index2[1]) * frequency;
208  textures[5] += std::pow((index2[0] - pixelMean) + (index2[1] - pixelMean), 3) * frequency;
209  textures[6] += std::pow((index2[0] - pixelMean) + (index2[1] - pixelMean), 4) * frequency;
210  textures[7] += index2[0] * index2[1] * frequency;
211  ++constVectorIt;
212  }
213  textures[7] = (fabs(marginalDevSquared) > 1E-8) ? (textures[7] - marginalMean * marginalMean) / marginalDevSquared : 0;
214 
215  // Return result
216  return textures;
217 }
218 
219 } // namespace otb
220 
221 #endif
otb::HaralickTexturesImageFunction::InputPixelType
InputImageType::PixelType InputPixelType
Definition: otbHaralickTexturesImageFunction.h:117
otb::HaralickTexturesImageFunction::InputImageType
TInputImage InputImageType
Definition: otbHaralickTexturesImageFunction.h:109
otb::HaralickTexturesImageFunction::CooccurrenceIndexType
CooccurrenceIndexedListType::IndexType CooccurrenceIndexType
Definition: otbHaralickTexturesImageFunction.h:129
otb::HaralickTexturesImageFunction::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbHaralickTexturesImageFunction.hxx:44
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::HaralickTexturesImageFunction::RelativeFrequencyType
CooccurrenceIndexedListType::RelativeFrequencyType RelativeFrequencyType
Definition: otbHaralickTexturesImageFunction.h:131
otb::HaralickTexturesImageFunction::OutputType
Superclass::OutputType OutputType
Definition: otbHaralickTexturesImageFunction.h:123
otb::HaralickTexturesImageFunction::EvaluateAtIndex
OutputType EvaluateAtIndex(const IndexType &index) const override
Definition: otbHaralickTexturesImageFunction.hxx:56
otb::HaralickTexturesImageFunction::VectorConstIteratorType
VectorType::const_iterator VectorConstIteratorType
Definition: otbHaralickTexturesImageFunction.h:135
otb::HaralickTexturesImageFunction::InputRegionType
InputImageType::RegionType InputRegionType
Definition: otbHaralickTexturesImageFunction.h:118
otb::HaralickTexturesImageFunction::HaralickTexturesImageFunction
HaralickTexturesImageFunction()
Definition: otbHaralickTexturesImageFunction.hxx:38
otb::HaralickTexturesImageFunction::VectorType
CooccurrenceIndexedListType::VectorType VectorType
Definition: otbHaralickTexturesImageFunction.h:132
otb::HaralickTexturesImageFunction::InputImagePointerType
InputImageType::Pointer InputImagePointerType
Definition: otbHaralickTexturesImageFunction.h:116
otb::HaralickTexturesImageFunction::SizeType
InputRegionType::SizeType SizeType
Definition: otbHaralickTexturesImageFunction.h:120
otb::HaralickTexturesImageFunction::CooccurrenceIndexedListPointerType
CooccurrenceIndexedListType::Pointer CooccurrenceIndexedListPointerType
Definition: otbHaralickTexturesImageFunction.h:127
otb::HaralickTexturesImageFunction::IndexType
Superclass::IndexType IndexType
Definition: otbHaralickTexturesImageFunction.h:113
otbHaralickTexturesImageFunction.h