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