Orfeo ToolBox  4.2
Orfeo ToolBox is not a black box
otbScalarImageToTexturesFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3 Program: ORFEO Toolbox
4 Language: C++
5 Date: $Date$
6 Version: $Revision$
7 
8 
9 Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10 See OTBCopyright.txt for details.
11 
12 
13 This software is distributed WITHOUT ANY WARRANTY; without even
14 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15 PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 
19 #ifndef __otbScalarImageToTexturesFilter_txx
20 #define __otbScalarImageToTexturesFilter_txx
21 
25 #include "itkImageRegionIterator.h"
26 #include "itkProgressReporter.h"
27 #include "itkNumericTraits.h"
28 #include <vector>
29 #include <cmath>
30 
31 namespace otb
32 {
33 template <class TInputImage, class TOutputImage>
36 : m_Radius()
37 , m_Offset()
38 , m_NeighborhoodRadius()
39 , m_NumberOfBinsPerAxis(8)
40 , m_InputImageMinimum(0)
41 , m_InputImageMaximum(255)
42 {
43  // There are 8 outputs corresponding to the 8 textures indices
44  this->SetNumberOfRequiredOutputs(8);
45 
46  // Create the 8 outputs
47  this->SetNthOutput(0, OutputImageType::New());
48  this->SetNthOutput(1, OutputImageType::New());
49  this->SetNthOutput(2, OutputImageType::New());
50  this->SetNthOutput(3, OutputImageType::New());
51  this->SetNthOutput(4, OutputImageType::New());
52  this->SetNthOutput(5, OutputImageType::New());
53  this->SetNthOutput(6, OutputImageType::New());
54  this->SetNthOutput(7, OutputImageType::New());
55 }
56 
57 template <class TInputImage, class TOutputImage>
60 {}
61 
62 template <class TInputImage, class TOutputImage>
64 ::OutputImageType *
67 {
68  if (this->GetNumberOfOutputs() < 1)
69  {
70  return 0;
71  }
72  return static_cast<OutputImageType *>(this->GetOutput(0));
73 }
74 
75 template <class TInputImage, class TOutputImage>
77 ::OutputImageType *
80 {
81  if (this->GetNumberOfOutputs() < 2)
82  {
83  return 0;
84  }
85  return static_cast<OutputImageType *>(this->GetOutput(1));
86 }
87 
88 template <class TInputImage, class TOutputImage>
90 ::OutputImageType *
93 {
94  if (this->GetNumberOfOutputs() < 3)
95  {
96  return 0;
97  }
98  return static_cast<OutputImageType *>(this->GetOutput(2));
99 }
100 
101 template <class TInputImage, class TOutputImage>
103 ::OutputImageType *
106 {
107  if (this->GetNumberOfOutputs() < 4)
108  {
109  return 0;
110  }
111  return static_cast<OutputImageType *>(this->GetOutput(3));
112 }
113 
114 template <class TInputImage, class TOutputImage>
116 ::OutputImageType *
119 {
120  if (this->GetNumberOfOutputs() < 5)
121  {
122  return 0;
123  }
124  return static_cast<OutputImageType *>(this->GetOutput(4));
125 }
126 
127 template <class TInputImage, class TOutputImage>
129 ::OutputImageType *
132 {
133  if (this->GetNumberOfOutputs() < 6)
134  {
135  return 0;
136  }
137  return static_cast<OutputImageType *>(this->GetOutput(5));
138 }
139 
140 template <class TInputImage, class TOutputImage>
142 ::OutputImageType *
145 {
146  if (this->GetNumberOfOutputs() < 7)
147  {
148  return 0;
149  }
150  return static_cast<OutputImageType *>(this->GetOutput(6));
151 }
152 
153 template <class TInputImage, class TOutputImage>
155 ::OutputImageType *
158 {
159  if (this->GetNumberOfOutputs() < 8)
160  {
161  return 0;
162  }
163  return static_cast<OutputImageType *>(this->GetOutput(7));
164 }
165 
166 template <class TInputImage, class TOutputImage>
167 void
170 {
171  // First, call superclass implementation
172  Superclass::GenerateInputRequestedRegion();
173 
174  // Retrieve the input and output pointers
175  InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
176  OutputImagePointerType outputPtr = this->GetOutput();
177 
178  if (!inputPtr || !outputPtr)
179  {
180  return;
181  }
182 
183  // Retrieve the output requested region
184  // We use only the first output since requested regions for all outputs are enforced to be equal
185  // by the default GenerateOutputRequestedRegiont() implementation
186  OutputRegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
187 
188  typename OutputRegionType::IndexType outputIndex = outputRequestedRegion.GetIndex();
189  typename OutputRegionType::SizeType outputSize = outputRequestedRegion.GetSize();
190  typename InputRegionType::IndexType inputIndex;
191  typename InputRegionType::SizeType inputSize;
192 
193  // First, apply offset
194  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
195  {
196  inputIndex[dim] = std::min(outputIndex[dim], outputIndex[dim] + m_Offset[dim]);
197  inputSize[dim] =
198  std::max(outputIndex[dim] + outputSize[dim], outputIndex[dim] + outputSize[dim] +
199  m_Offset[dim]) - inputIndex[dim];
200  }
201 
202  // Build the input requested region
203  InputRegionType inputRequestedRegion;
204  inputRequestedRegion.SetIndex(inputIndex);
205  inputRequestedRegion.SetSize(inputSize);
206 
207  // Apply the radius
208  inputRequestedRegion.PadByRadius(m_Radius);
209 
210  // Try to apply the requested region to the input image
211  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
212  {
213  inputPtr->SetRequestedRegion(inputRequestedRegion);
214  }
215  else
216  {
217  // Build an exception
218  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
219  e.SetLocation(ITK_LOCATION);
220  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
221  e.SetDataObject(inputPtr);
222  throw e;
223  }
224 }
225 
226 template <class TInputImage, class TOutputImage>
227 void
230 {
231  unsigned int minRadius = 0;
232  for ( unsigned int i = 0; i < m_Offset.GetOffsetDimension(); i++ )
233  {
234  unsigned int distance = vcl_abs(m_Offset[i]);
235  if ( distance > minRadius )
236  {
237  minRadius = distance;
238  }
239  }
240  m_NeighborhoodRadius.Fill(minRadius);
241 }
242 
243 template <class TInputImage, class TOutputImage>
244 void
246 ::ThreadedGenerateData(const OutputRegionType& outputRegionForThread, itk::ThreadIdType threadId)
247 {
248  // Retrieve the input and output pointers
249  InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
250  OutputImagePointerType energyPtr = this->GetEnergyOutput();
251  OutputImagePointerType entropyPtr = this->GetEntropyOutput();
252  OutputImagePointerType correlationPtr = this->GetCorrelationOutput();
253  OutputImagePointerType invDiffMomentPtr = this->GetInverseDifferenceMomentOutput();
254  OutputImagePointerType inertiaPtr = this->GetInertiaOutput();
255  OutputImagePointerType clusterShadePtr = this->GetClusterShadeOutput();
256  OutputImagePointerType clusterProminencePtr = this->GetClusterProminenceOutput();
257  OutputImagePointerType haralickCorPtr = this->GetHaralickCorrelationOutput();
258 
259  // Build output iterators
260  itk::ImageRegionIteratorWithIndex<OutputImageType> energyIt(energyPtr, outputRegionForThread);
261  itk::ImageRegionIterator<OutputImageType> entropyIt(entropyPtr, outputRegionForThread);
262  itk::ImageRegionIterator<OutputImageType> correlationIt(correlationPtr, outputRegionForThread);
263  itk::ImageRegionIterator<OutputImageType> invDiffMomentIt(invDiffMomentPtr, outputRegionForThread);
264  itk::ImageRegionIterator<OutputImageType> inertiaIt(inertiaPtr, outputRegionForThread);
265  itk::ImageRegionIterator<OutputImageType> clusterShadeIt(clusterShadePtr, outputRegionForThread);
266  itk::ImageRegionIterator<OutputImageType> clusterProminenceIt(clusterProminencePtr, outputRegionForThread);
267  itk::ImageRegionIterator<OutputImageType> haralickCorIt(haralickCorPtr, outputRegionForThread);
268 
269  // Go to begin
270  energyIt.GoToBegin();
271  entropyIt.GoToBegin();
272  correlationIt.GoToBegin();
273  invDiffMomentIt.GoToBegin();
274  inertiaIt.GoToBegin();
275  clusterShadeIt.GoToBegin();
276  clusterProminenceIt.GoToBegin();
277  haralickCorIt.GoToBegin();
278 
279  const double log2 = vcl_log(2.0);
280 
281  // Set-up progress reporting
282  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
283 
284  // Iterate on outputs to compute textures
285  while (!energyIt.IsAtEnd()
286  && !entropyIt.IsAtEnd()
287  && !correlationIt.IsAtEnd()
288  && !invDiffMomentIt.IsAtEnd()
289  && !inertiaIt.IsAtEnd()
290  && !clusterShadeIt.IsAtEnd()
291  && !clusterProminenceIt.IsAtEnd()
292  && !haralickCorIt.IsAtEnd())
293  {
294  // Compute the region on which co-occurence will be estimated
295  typename InputRegionType::IndexType inputIndex;
296  typename InputRegionType::SizeType inputSize;
297 
298  // First, create an window for neighborhood iterator based on m_Radius
299  // For example, if xradius and yradius is 2. window size is 5x5 (2 * radius + 1).
300  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
301  {
302  inputIndex[dim] = energyIt.GetIndex()[dim] - m_Radius[dim];
303  inputSize[dim] = 2 * m_Radius[dim] + 1;
304  }
305 
306  // Build the input region
307  InputRegionType inputRegion;
308  inputRegion.SetIndex(inputIndex);
309  inputRegion.SetSize(inputSize);
310  inputRegion.Crop(inputPtr->GetRequestedRegion());
311 
312  CooccurrenceIndexedListPointerType GLCIList = CooccurrenceIndexedListType::New();
313  GLCIList->Initialize(m_NumberOfBinsPerAxis, m_InputImageMinimum, m_InputImageMaximum);
314 
315  typedef itk::ConstNeighborhoodIterator< InputImageType > NeighborhoodIteratorType;
316  NeighborhoodIteratorType neighborIt;
317  neighborIt = NeighborhoodIteratorType(m_NeighborhoodRadius, inputPtr, inputRegion);
318  for ( neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt )
319  {
320  const InputPixelType centerPixelIntensity = neighborIt.GetCenterPixel();
321  bool pixelInBounds;
322  const InputPixelType pixelIntensity = neighborIt.GetPixel(m_Offset, pixelInBounds);
323  if ( !pixelInBounds )
324  {
325  continue; // don't put a pixel in the co-occurrence list if the value is
326  // out of bounds
327  }
328  GLCIList->AddPixelPair(centerPixelIntensity, pixelIntensity);
329  }
330 
331  double pixelMean = 0.;
332  double marginalMean;
333  double marginalDevSquared = 0.;
334  double pixelVariance = 0.;
335 
336  //Create and Intialize marginalSums
337  std::vector<double> marginalSums(m_NumberOfBinsPerAxis, 0);
338 
339  //get co-occurrence vector and totalfrequency
340  VectorType glcVector = GLCIList->GetVector();
341  double totalFrequency = static_cast<double> (GLCIList->GetTotalFrequency());
342 
343  //Normalize the co-occurrence indexed list and compute mean, marginalSum
344  typename VectorType::iterator it = glcVector.begin();
345  while( it != glcVector.end())
346  {
347  double frequency = (*it).second / totalFrequency;
348  CooccurrenceIndexType index = (*it).first;
349  pixelMean += index[0] * frequency;
350  marginalSums[index[0]] += frequency;
351  ++it;
352  }
353 
354  /* Now get the mean and deviaton of the marginal sums.
355  Compute incremental mean and SD, a la Knuth, "The Art of Computer
356  Programming, Volume 2: Seminumerical Algorithms", section 4.2.2.
357  Compute mean and standard deviation using the recurrence relation:
358  M(1) = x(1), M(k) = M(k-1) + (x(k) - M(k-1) ) / k
359  S(1) = 0, S(k) = S(k-1) + (x(k) - M(k-1)) * (x(k) - M(k))
360  for 2 <= k <= n, then
361  sigma = vcl_sqrt(S(n) / n) (or divide by n-1 for sample SD instead of
362  population SD).
363  */
364  std::vector<double>::const_iterator msIt = marginalSums.begin();
365  marginalMean = *msIt;
366  //Increment iterator to start with index 1
367  ++msIt;
368  for(int k= 2; msIt != marginalSums.end(); ++k, ++msIt)
369  {
370  double M_k_minus_1 = marginalMean;
371  double S_k_minus_1 = marginalDevSquared;
372  double x_k = *msIt;
373  double M_k = M_k_minus_1 + ( x_k - M_k_minus_1 ) / k;
374  double S_k = S_k_minus_1 + ( x_k - M_k_minus_1 ) * ( x_k - M_k );
375  marginalMean = M_k;
376  marginalDevSquared = S_k;
377  }
378  marginalDevSquared = marginalDevSquared / m_NumberOfBinsPerAxis;
379 
380  VectorConstIteratorType constVectorIt;
381  constVectorIt = glcVector.begin();
382  while( constVectorIt != glcVector.end())
383  {
384  RelativeFrequencyType frequency = (*constVectorIt).second / totalFrequency;
385  CooccurrenceIndexType index = (*constVectorIt).first;
386  pixelVariance += ( index[0] - pixelMean ) * ( index[0] - pixelMean ) * frequency;
387  ++constVectorIt;
388  }
389 
390  double pixelVarianceSquared = pixelVariance * pixelVariance;
391  // Variance is only used in correlation. If variance is 0, then (index[0] - pixelMean) * (index[1] - pixelMean)
392  // should be zero as well. In this case, set the variance to 1. in order to
393  // avoid NaN correlation.
394  if(pixelVarianceSquared < GetPixelValueTolerance())
395  {
396  pixelVarianceSquared = 1.;
397  }
398 
399  //Initalize texture variables;
408 
409  //Compute textures
410  constVectorIt = glcVector.begin();
411  while( constVectorIt != glcVector.end())
412  {
413  CooccurrenceIndexType index = (*constVectorIt).first;
414  RelativeFrequencyType frequency = (*constVectorIt).second / totalFrequency;
415  energy += frequency * frequency;
416  entropy -= ( frequency > GetPixelValueTolerance() ) ? frequency *vcl_log(frequency) / log2 : 0;
417  correlation += ( ( index[0] - pixelMean ) * ( index[1] - pixelMean ) * frequency ) / pixelVarianceSquared;
418  inverseDifferenceMoment += frequency / ( 1.0 + ( index[0] - index[1] ) * ( index[0] - index[1] ) );
419  inertia += ( index[0] - index[1] ) * ( index[0] - index[1] ) * frequency;
420  clusterShade += vcl_pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 3 ) * frequency;
421  clusterProminence += vcl_pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 4 ) * frequency;
422  haralickCorrelation += index[0] * index[1] * frequency;
423  ++constVectorIt;
424  }
425 
426  haralickCorrelation = (fabs(marginalDevSquared) > 1E-8) ?
427  ( haralickCorrelation - marginalMean * marginalMean ) / marginalDevSquared : 0;
428 
429  // Fill outputs
430  energyIt.Set(energy);
431  entropyIt.Set(entropy);
432  correlationIt.Set(correlation);
433  invDiffMomentIt.Set(inverseDifferenceMoment);
434  inertiaIt.Set(inertia);
435  clusterShadeIt.Set(clusterShade);
436  clusterProminenceIt.Set(clusterProminence);
437  haralickCorIt.Set(haralickCorrelation);
438 
439  // Update progress
440  progress.CompletedPixel();
441 
442  // Increment iterators
443  ++energyIt;
444  ++entropyIt;
445  ++correlationIt;
446  ++invDiffMomentIt;
447  ++inertiaIt;
448  ++clusterShadeIt;
449  ++clusterProminenceIt;
450  ++haralickCorIt;
451  }
452 }
453 
454 } // End namespace otb
455 
456 #endif
void SetDataObject(DataObject *dobj)
virtual void SetDescription(const std::string &s)
CooccurrenceIndexedListType::PixelValueType PixelValueType
CooccurrenceIndexedListType::VectorType VectorType
virtual void SetLocation(const std::string &s)
void Set(const PixelType &value) const
void Set(const PixelType &value) const
const IndexType & GetIndex() const
virtual void ThreadedGenerateData(const OutputRegionType &outputRegion, itk::ThreadIdType threadId)
CooccurrenceIndexedListType::RelativeFrequencyType RelativeFrequencyType
This class compute 8 local Haralick textures features. The 8 output image channels are: Energy...
VectorType::const_iterator VectorConstIteratorType
bool IsAtEnd(void) const
unsigned int ThreadIdType