OTB  5.11.0
Orfeo Toolbox
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 , m_SubsampleFactor()
43 , m_SubsampleOffset()
44 {
45  // There are 8 outputs corresponding to the 8 textures indices
46  this->SetNumberOfRequiredOutputs(8);
47 
48  // Create the 8 outputs
49  this->SetNthOutput(0, OutputImageType::New());
50  this->SetNthOutput(1, OutputImageType::New());
51  this->SetNthOutput(2, OutputImageType::New());
52  this->SetNthOutput(3, OutputImageType::New());
53  this->SetNthOutput(4, OutputImageType::New());
54  this->SetNthOutput(5, OutputImageType::New());
55  this->SetNthOutput(6, OutputImageType::New());
56  this->SetNthOutput(7, OutputImageType::New());
57 
58  this->m_SubsampleFactor.Fill(1);
59  this->m_SubsampleOffset.Fill(0);
60 }
61 
62 template <class TInputImage, class TOutputImage>
65 {}
66 
67 template <class TInputImage, class TOutputImage>
69 ::OutputImageType *
72 {
73  if (this->GetNumberOfOutputs() < 1)
74  {
75  return ITK_NULLPTR;
76  }
77  return static_cast<OutputImageType *>(this->GetOutput(0));
78 }
79 
80 template <class TInputImage, class TOutputImage>
82 ::OutputImageType *
85 {
86  if (this->GetNumberOfOutputs() < 2)
87  {
88  return ITK_NULLPTR;
89  }
90  return static_cast<OutputImageType *>(this->GetOutput(1));
91 }
92 
93 template <class TInputImage, class TOutputImage>
95 ::OutputImageType *
98 {
99  if (this->GetNumberOfOutputs() < 3)
100  {
101  return ITK_NULLPTR;
102  }
103  return static_cast<OutputImageType *>(this->GetOutput(2));
104 }
105 
106 template <class TInputImage, class TOutputImage>
108 ::OutputImageType *
111 {
112  if (this->GetNumberOfOutputs() < 4)
113  {
114  return ITK_NULLPTR;
115  }
116  return static_cast<OutputImageType *>(this->GetOutput(3));
117 }
118 
119 template <class TInputImage, class TOutputImage>
121 ::OutputImageType *
124 {
125  if (this->GetNumberOfOutputs() < 5)
126  {
127  return ITK_NULLPTR;
128  }
129  return static_cast<OutputImageType *>(this->GetOutput(4));
130 }
131 
132 template <class TInputImage, class TOutputImage>
134 ::OutputImageType *
137 {
138  if (this->GetNumberOfOutputs() < 6)
139  {
140  return ITK_NULLPTR;
141  }
142  return static_cast<OutputImageType *>(this->GetOutput(5));
143 }
144 
145 template <class TInputImage, class TOutputImage>
147 ::OutputImageType *
150 {
151  if (this->GetNumberOfOutputs() < 7)
152  {
153  return ITK_NULLPTR;
154  }
155  return static_cast<OutputImageType *>(this->GetOutput(6));
156 }
157 
158 template <class TInputImage, class TOutputImage>
160 ::OutputImageType *
163 {
164  if (this->GetNumberOfOutputs() < 8)
165  {
166  return ITK_NULLPTR;
167  }
168  return static_cast<OutputImageType *>(this->GetOutput(7));
169 }
170 
171 template <class TInputImage, class TOutputImage>
172 void
175 {
176  // First, call superclass implementation
177  Superclass::GenerateOutputInformation();
178 
179  // Compute output size, origin & spacing
180  InputRegionType inputRegion = this->GetInput()->GetLargestPossibleRegion();
181  OutputRegionType outputRegion;
182  outputRegion.SetIndex(0,0);
183  outputRegion.SetIndex(1,0);
184  outputRegion.SetSize(0, 1 + (inputRegion.GetSize(0) - 1 - m_SubsampleOffset[0]) / m_SubsampleFactor[0]);
185  outputRegion.SetSize(1, 1 + (inputRegion.GetSize(1) - 1 - m_SubsampleOffset[1]) / m_SubsampleFactor[1]);
186 
187  typename OutputImageType::SpacingType outSpacing = this->GetInput()->GetSpacing();
188  outSpacing[0] *= m_SubsampleFactor[0];
189  outSpacing[1] *= m_SubsampleFactor[1];
190 
191  typename OutputImageType::PointType outOrigin;
192  this->GetInput()->TransformIndexToPhysicalPoint(inputRegion.GetIndex()+m_SubsampleOffset,outOrigin);
193 
194  for (unsigned int i=0 ; i<this->GetNumberOfOutputs() ; i++)
195  {
196  OutputImagePointerType outputPtr = this->GetOutput(i);
197  outputPtr->SetLargestPossibleRegion(outputRegion);
198  outputPtr->SetOrigin(outOrigin);
199  outputPtr->SetSpacing(outSpacing);
200  }
201 }
202 
203 template <class TInputImage, class TOutputImage>
204 void
207 {
208  // First, call superclass implementation
209  Superclass::GenerateInputRequestedRegion();
210 
211  // Retrieve the input and output pointers
212  InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
213  OutputImagePointerType outputPtr = this->GetOutput();
214 
215  if (!inputPtr || !outputPtr)
216  {
217  return;
218  }
219 
220  // Retrieve the output requested region
221  // We use only the first output since requested regions for all outputs are enforced to be equal
222  // by the default GenerateOutputRequestedRegiont() implementation
223  OutputRegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
224 
225  typename OutputRegionType::IndexType outputIndex = outputRequestedRegion.GetIndex();
226  typename OutputRegionType::SizeType outputSize = outputRequestedRegion.GetSize();
227  typename InputRegionType::IndexType inputIndex;
228  typename InputRegionType::SizeType inputSize;
229  InputRegionType inputLargest = inputPtr->GetLargestPossibleRegion();
230 
231  // Convert index and size to full grid
232  outputIndex[0] = outputIndex[0] * m_SubsampleFactor[0] + m_SubsampleOffset[0] + inputLargest.GetIndex(0);
233  outputIndex[1] = outputIndex[1] * m_SubsampleFactor[1] + m_SubsampleOffset[1] + inputLargest.GetIndex(1);
234  outputSize[0] = 1 + (outputSize[0] - 1) * m_SubsampleFactor[0];
235  outputSize[1] = 1 + (outputSize[1] - 1) * m_SubsampleFactor[1];
236 
237  // First, apply offset
238  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
239  {
240  inputIndex[dim] = std::min(outputIndex[dim], outputIndex[dim] + m_Offset[dim]);
241  inputSize[dim] =
242  std::max(outputIndex[dim] + outputSize[dim], outputIndex[dim] + outputSize[dim] +
243  m_Offset[dim]) - inputIndex[dim];
244  }
245 
246  // Build the input requested region
247  InputRegionType inputRequestedRegion;
248  inputRequestedRegion.SetIndex(inputIndex);
249  inputRequestedRegion.SetSize(inputSize);
250 
251  // Apply the radius
252  inputRequestedRegion.PadByRadius(m_Radius);
253 
254  // Try to apply the requested region to the input image
255  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
256  {
257  inputPtr->SetRequestedRegion(inputRequestedRegion);
258  }
259  else
260  {
261  // Build an exception
262  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
263  e.SetLocation(ITK_LOCATION);
264  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
265  e.SetDataObject(inputPtr);
266  throw e;
267  }
268 }
269 
270 template <class TInputImage, class TOutputImage>
271 void
274 {
275  unsigned int minRadius = 0;
276  for ( unsigned int i = 0; i < m_Offset.GetOffsetDimension(); i++ )
277  {
278  unsigned int distance = vcl_abs(m_Offset[i]);
279  if ( distance > minRadius )
280  {
281  minRadius = distance;
282  }
283  }
284  m_NeighborhoodRadius.Fill(minRadius);
285 }
286 
287 template <class TInputImage, class TOutputImage>
288 void
290 ::ThreadedGenerateData(const OutputRegionType& outputRegionForThread, itk::ThreadIdType threadId)
291 {
292  // Retrieve the input and output pointers
293  InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
294  OutputImagePointerType energyPtr = this->GetEnergyOutput();
295  OutputImagePointerType entropyPtr = this->GetEntropyOutput();
296  OutputImagePointerType correlationPtr = this->GetCorrelationOutput();
297  OutputImagePointerType invDiffMomentPtr = this->GetInverseDifferenceMomentOutput();
298  OutputImagePointerType inertiaPtr = this->GetInertiaOutput();
299  OutputImagePointerType clusterShadePtr = this->GetClusterShadeOutput();
300  OutputImagePointerType clusterProminencePtr = this->GetClusterProminenceOutput();
301  OutputImagePointerType haralickCorPtr = this->GetHaralickCorrelationOutput();
302 
303  // Build output iterators
304  itk::ImageRegionIteratorWithIndex<OutputImageType> energyIt(energyPtr, outputRegionForThread);
305  itk::ImageRegionIterator<OutputImageType> entropyIt(entropyPtr, outputRegionForThread);
306  itk::ImageRegionIterator<OutputImageType> correlationIt(correlationPtr, outputRegionForThread);
307  itk::ImageRegionIterator<OutputImageType> invDiffMomentIt(invDiffMomentPtr, outputRegionForThread);
308  itk::ImageRegionIterator<OutputImageType> inertiaIt(inertiaPtr, outputRegionForThread);
309  itk::ImageRegionIterator<OutputImageType> clusterShadeIt(clusterShadePtr, outputRegionForThread);
310  itk::ImageRegionIterator<OutputImageType> clusterProminenceIt(clusterProminencePtr, outputRegionForThread);
311  itk::ImageRegionIterator<OutputImageType> haralickCorIt(haralickCorPtr, outputRegionForThread);
312 
313  // Go to begin
314  energyIt.GoToBegin();
315  entropyIt.GoToBegin();
316  correlationIt.GoToBegin();
317  invDiffMomentIt.GoToBegin();
318  inertiaIt.GoToBegin();
319  clusterShadeIt.GoToBegin();
320  clusterProminenceIt.GoToBegin();
321  haralickCorIt.GoToBegin();
322 
323  const double log2 = vcl_log(2.0);
324 
325  InputRegionType inputLargest = inputPtr->GetLargestPossibleRegion();
326 
327  // Set-up progress reporting
328  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
329 
330  // Iterate on outputs to compute textures
331  while (!energyIt.IsAtEnd()
332  && !entropyIt.IsAtEnd()
333  && !correlationIt.IsAtEnd()
334  && !invDiffMomentIt.IsAtEnd()
335  && !inertiaIt.IsAtEnd()
336  && !clusterShadeIt.IsAtEnd()
337  && !clusterProminenceIt.IsAtEnd()
338  && !haralickCorIt.IsAtEnd())
339  {
340  // Compute the region on which co-occurence will be estimated
341  typename InputRegionType::IndexType inputIndex;
342  typename InputRegionType::SizeType inputSize;
343 
344  // Convert index to full grid
345  typename OutputImageType::IndexType outIndex;
346 
347  // First, create an window for neighborhood iterator based on m_Radius
348  // For example, if xradius and yradius is 2. window size is 5x5 (2 * radius + 1).
349  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
350  {
351  outIndex[dim] = energyIt.GetIndex()[dim] * m_SubsampleFactor[dim] + m_SubsampleOffset[dim] + inputLargest.GetIndex(dim);
352  inputIndex[dim] = outIndex[dim] - m_Radius[dim];
353  inputSize[dim] = 2 * m_Radius[dim] + 1;
354  }
355 
356  // Build the input region
357  InputRegionType inputRegion;
358  inputRegion.SetIndex(inputIndex);
359  inputRegion.SetSize(inputSize);
360  inputRegion.Crop(inputPtr->GetRequestedRegion());
361 
362  CooccurrenceIndexedListPointerType GLCIList = CooccurrenceIndexedListType::New();
363  GLCIList->Initialize(m_NumberOfBinsPerAxis, m_InputImageMinimum, m_InputImageMaximum);
364 
365  typedef itk::ConstNeighborhoodIterator< InputImageType > NeighborhoodIteratorType;
366  NeighborhoodIteratorType neighborIt;
367  neighborIt = NeighborhoodIteratorType(m_NeighborhoodRadius, inputPtr, inputRegion);
368  for ( neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt )
369  {
370  const InputPixelType centerPixelIntensity = neighborIt.GetCenterPixel();
371  bool pixelInBounds;
372  const InputPixelType pixelIntensity = neighborIt.GetPixel(m_Offset, pixelInBounds);
373  if ( !pixelInBounds )
374  {
375  continue; // don't put a pixel in the co-occurrence list if the value is
376  // out of bounds
377  }
378  GLCIList->AddPixelPair(centerPixelIntensity, pixelIntensity);
379  }
380 
381  double pixelMean = 0.;
382  double marginalMean;
383  double marginalDevSquared = 0.;
384  double pixelVariance = 0.;
385 
386  //Create and Initialize marginalSums
387  std::vector<double> marginalSums(m_NumberOfBinsPerAxis, 0);
388 
389  //get co-occurrence vector and totalfrequency
390  VectorType glcVector = GLCIList->GetVector();
391  double totalFrequency = static_cast<double> (GLCIList->GetTotalFrequency());
392 
393  //Normalize the co-occurrence indexed list and compute mean, marginalSum
394  typename VectorType::iterator it = glcVector.begin();
395  while( it != glcVector.end())
396  {
397  double frequency = (*it).second / totalFrequency;
398  CooccurrenceIndexType index = (*it).first;
399  pixelMean += index[0] * frequency;
400  marginalSums[index[0]] += frequency;
401  ++it;
402  }
403 
404  /* Now get the mean and deviaton of the marginal sums.
405  Compute incremental mean and SD, a la Knuth, "The Art of Computer
406  Programming, Volume 2: Seminumerical Algorithms", section 4.2.2.
407  Compute mean and standard deviation using the recurrence relation:
408  M(1) = x(1), M(k) = M(k-1) + (x(k) - M(k-1) ) / k
409  S(1) = 0, S(k) = S(k-1) + (x(k) - M(k-1)) * (x(k) - M(k))
410  for 2 <= k <= n, then
411  sigma = vcl_sqrt(S(n) / n) (or divide by n-1 for sample SD instead of
412  population SD).
413  */
414  std::vector<double>::const_iterator msIt = marginalSums.begin();
415  marginalMean = *msIt;
416  //Increment iterator to start with index 1
417  ++msIt;
418  for(int k= 2; msIt != marginalSums.end(); ++k, ++msIt)
419  {
420  double M_k_minus_1 = marginalMean;
421  double S_k_minus_1 = marginalDevSquared;
422  double x_k = *msIt;
423  double M_k = M_k_minus_1 + ( x_k - M_k_minus_1 ) / k;
424  double S_k = S_k_minus_1 + ( x_k - M_k_minus_1 ) * ( x_k - M_k );
425  marginalMean = M_k;
426  marginalDevSquared = S_k;
427  }
428  marginalDevSquared = marginalDevSquared / m_NumberOfBinsPerAxis;
429 
430  VectorConstIteratorType constVectorIt;
431  constVectorIt = glcVector.begin();
432  while( constVectorIt != glcVector.end())
433  {
434  RelativeFrequencyType frequency = (*constVectorIt).second / totalFrequency;
435  CooccurrenceIndexType index = (*constVectorIt).first;
436  pixelVariance += ( index[0] - pixelMean ) * ( index[0] - pixelMean ) * frequency;
437  ++constVectorIt;
438  }
439 
440  double pixelVarianceSquared = pixelVariance * pixelVariance;
441  // Variance is only used in correlation. If variance is 0, then (index[0] - pixelMean) * (index[1] - pixelMean)
442  // should be zero as well. In this case, set the variance to 1. in order to
443  // avoid NaN correlation.
444  if(pixelVarianceSquared < GetPixelValueTolerance())
445  {
446  pixelVarianceSquared = 1.;
447  }
448 
449  //Initialize texture variables;
458 
459  //Compute textures
460  constVectorIt = glcVector.begin();
461  while( constVectorIt != glcVector.end())
462  {
463  CooccurrenceIndexType index = (*constVectorIt).first;
464  RelativeFrequencyType frequency = (*constVectorIt).second / totalFrequency;
465  energy += frequency * frequency;
466  entropy -= ( frequency > GetPixelValueTolerance() ) ? frequency *vcl_log(frequency) / log2 : 0;
467  correlation += ( ( index[0] - pixelMean ) * ( index[1] - pixelMean ) * frequency ) / pixelVarianceSquared;
468  inverseDifferenceMoment += frequency / ( 1.0 + ( index[0] - index[1] ) * ( index[0] - index[1] ) );
469  inertia += ( index[0] - index[1] ) * ( index[0] - index[1] ) * frequency;
470  clusterShade += vcl_pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 3 ) * frequency;
471  clusterProminence += vcl_pow( ( index[0] - pixelMean ) + ( index[1] - pixelMean ), 4 ) * frequency;
472  haralickCorrelation += index[0] * index[1] * frequency;
473  ++constVectorIt;
474  }
475 
476  haralickCorrelation = (fabs(marginalDevSquared) > 1E-8) ?
477  ( haralickCorrelation - marginalMean * marginalMean ) / marginalDevSquared : 0;
478 
479  // Fill outputs
480  energyIt.Set(energy);
481  entropyIt.Set(entropy);
482  correlationIt.Set(correlation);
483  invDiffMomentIt.Set(inverseDifferenceMoment);
484  inertiaIt.Set(inertia);
485  clusterShadeIt.Set(clusterShade);
486  clusterProminenceIt.Set(clusterProminence);
487  haralickCorIt.Set(haralickCorrelation);
488 
489  // Update progress
490  progress.CompletedPixel();
491 
492  // Increment iterators
493  ++energyIt;
494  ++entropyIt;
495  ++correlationIt;
496  ++invDiffMomentIt;
497  ++inertiaIt;
498  ++clusterShadeIt;
499  ++clusterProminenceIt;
500  ++haralickCorIt;
501  }
502 }
503 
504 } // End namespace otb
505 
506 #endif
itk::Size< Monteverdi_DIMENSION > SizeType
Definition: mvdTypes.h:144
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
void ThreadedGenerateData(const OutputRegionType &outputRegion, itk::ThreadIdType threadId) ITK_OVERRIDE
itk::Index< Monteverdi_DIMENSION > IndexType
Definition: mvdTypes.h:140
CooccurrenceIndexedListType::RelativeFrequencyType RelativeFrequencyType
This class compute 8 local Haralick textures features. The 8 output image channels are: Energy...
VectorType::const_iterator VectorConstIteratorType
unsigned int ThreadIdType
VectorImageType::SpacingType SpacingType
Definition: mvdTypes.h:188
bool IsAtEnd(void) const
TOutputImage OutputImageType
VectorImageType::PointType PointType
Definition: mvdTypes.h:196
const SizeValueType * GetSize() const