OTB  6.7.0
Orfeo Toolbox
otbScalarImageToAdvancedTexturesFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 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 otbScalarImageToAdvancedTexturesFilter_hxx
22 #define otbScalarImageToAdvancedTexturesFilter_hxx
23 
27 #include "itkImageRegionIterator.h"
28 #include "itkProgressReporter.h"
29 #include "itkNumericTraits.h"
30 #include <algorithm>
31 
32 namespace otb
33 {
34 template <class TInputImage, class TOutputImage>
37 : m_Radius()
38 , m_Offset()
39 , m_NeighborhoodRadius()
40 , m_NumberOfBinsPerAxis(8)
41 , m_InputImageMinimum(0)
42 , m_InputImageMaximum(255)
43 , m_SubsampleFactor()
44 , m_SubsampleOffset()
45 {
46  // There are 10 outputs corresponding to the 9 textures indices
47  this->SetNumberOfRequiredOutputs(10);
48 
49  // Create the 10 outputs
50  this->SetNthOutput(0, OutputImageType::New());
51  this->SetNthOutput(1, OutputImageType::New());
52  this->SetNthOutput(2, OutputImageType::New());
53  this->SetNthOutput(3, OutputImageType::New());
54  this->SetNthOutput(4, OutputImageType::New());
55  this->SetNthOutput(5, OutputImageType::New());
56  this->SetNthOutput(6, OutputImageType::New());
57  this->SetNthOutput(7, OutputImageType::New());
58  this->SetNthOutput(8, OutputImageType::New());
59  this->SetNthOutput(9, 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 template <class TInputImage, class TOutputImage>
71 ::OutputImageType *
74 {
75  if (this->GetNumberOfOutputs() < 1)
76  {
77  return nullptr;
78  }
79  return static_cast<OutputImageType *>(this->GetOutput(0));
80 }
81 
82 template <class TInputImage, class TOutputImage>
84 ::OutputImageType *
87 {
88  if (this->GetNumberOfOutputs() < 2)
89  {
90  return nullptr;
91  }
92  return static_cast<OutputImageType *>(this->GetOutput(1));
93 }
94 
95 template <class TInputImage, class TOutputImage>
97 ::OutputImageType *
100 {
101  if (this->GetNumberOfOutputs() < 3)
102  {
103  return nullptr;
104  }
105  return static_cast<OutputImageType *>(this->GetOutput(2));
106 }
107 
108 template <class TInputImage, class TOutputImage>
110 ::OutputImageType *
113 {
114  if (this->GetNumberOfOutputs() < 4)
115  {
116  return nullptr;
117  }
118  return static_cast<OutputImageType *>(this->GetOutput(3));
119 }
120 
121 template <class TInputImage, class TOutputImage>
123 ::OutputImageType *
126 {
127  if (this->GetNumberOfOutputs() < 5)
128  {
129  return nullptr;
130  }
131  return static_cast<OutputImageType *>(this->GetOutput(4));
132 }
133 
134 template <class TInputImage, class TOutputImage>
136 ::OutputImageType *
139 {
140  if (this->GetNumberOfOutputs() < 6)
141  {
142  return nullptr;
143  }
144  return static_cast<OutputImageType *>(this->GetOutput(5));
145 }
146 
147 template <class TInputImage, class TOutputImage>
149 ::OutputImageType *
152 {
153  if (this->GetNumberOfOutputs() < 7)
154  {
155  return nullptr;
156  }
157  return static_cast<OutputImageType *>(this->GetOutput(6));
158 }
159 
160 template <class TInputImage, class TOutputImage>
162 ::OutputImageType *
165 {
166  if (this->GetNumberOfOutputs() < 8)
167  {
168  return nullptr;
169  }
170  return static_cast<OutputImageType *>(this->GetOutput(7));
171 }
172 
173 template <class TInputImage, class TOutputImage>
175 ::OutputImageType *
178 {
179  if (this->GetNumberOfOutputs() < 9)
180  {
181  return nullptr;
182  }
183  return static_cast<OutputImageType *>(this->GetOutput(8));
184 }
185 
186 template <class TInputImage, class TOutputImage>
188 ::OutputImageType *
191 {
192  if (this->GetNumberOfOutputs() < 10)
193  {
194  return nullptr;
195  }
196  return static_cast<OutputImageType *>(this->GetOutput(9));
197 }
198 
199 template <class TInputImage, class TOutputImage>
200 void
203 {
204  // First, call superclass implementation
205  Superclass::GenerateOutputInformation();
206 
207  // Compute output size, origin & spacing
208  InputRegionType inputRegion = this->GetInput()->GetLargestPossibleRegion();
209  OutputRegionType outputRegion;
210  outputRegion.SetIndex(0,0);
211  outputRegion.SetIndex(1,0);
212  outputRegion.SetSize(0, 1 + (inputRegion.GetSize(0) - 1 - m_SubsampleOffset[0]) / m_SubsampleFactor[0]);
213  outputRegion.SetSize(1, 1 + (inputRegion.GetSize(1) - 1 - m_SubsampleOffset[1]) / m_SubsampleFactor[1]);
214 
215  typename OutputImageType::SpacingType outSpacing = this->GetInput()->GetSignedSpacing();
216  outSpacing[0] *= m_SubsampleFactor[0];
217  outSpacing[1] *= m_SubsampleFactor[1];
218 
219  typename OutputImageType::PointType outOrigin;
220  this->GetInput()->TransformIndexToPhysicalPoint(inputRegion.GetIndex()+m_SubsampleOffset,outOrigin);
221 
222  for (unsigned int i=0 ; i<this->GetNumberOfOutputs() ; i++)
223  {
224  OutputImagePointerType outputPtr = this->GetOutput(i);
225  outputPtr->SetLargestPossibleRegion(outputRegion);
226  outputPtr->SetOrigin(outOrigin);
227  outputPtr->SetSignedSpacing(outSpacing);
228  }
229 }
230 
231 template <class TInputImage, class TOutputImage>
232 void
235 {
236  // First, call superclass implementation
237  Superclass::GenerateInputRequestedRegion();
238 
239  // Retrieve the input and output pointers
240  InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
241  OutputImagePointerType outputPtr = this->GetOutput();
242 
243  if (!inputPtr || !outputPtr)
244  {
245  return;
246  }
247 
248  // Retrieve the output requested region
249  // We use only the first output since requested regions for all outputs are enforced to be equal
250  // by the default GenerateOutputRequestedRegiont() implementation
251  OutputRegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
252  InputRegionType inputLargest = inputPtr->GetLargestPossibleRegion();
253 
254  typename OutputRegionType::IndexType outputIndex = outputRequestedRegion.GetIndex();
255  typename OutputRegionType::SizeType outputSize = outputRequestedRegion.GetSize();
256  typename InputRegionType::IndexType inputIndex;
257  typename InputRegionType::SizeType inputSize;
258 
259  // Convert index and size to full grid
260  outputIndex[0] = outputIndex[0] * m_SubsampleFactor[0] + m_SubsampleOffset[0] + inputLargest.GetIndex(0);
261  outputIndex[1] = outputIndex[1] * m_SubsampleFactor[1] + m_SubsampleOffset[1] + inputLargest.GetIndex(1);
262  outputSize[0] = 1 + (outputSize[0] - 1) * m_SubsampleFactor[0];
263  outputSize[1] = 1 + (outputSize[1] - 1) * m_SubsampleFactor[1];
264 
265  // First, apply offset
266  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
267  {
268  inputIndex[dim] = std::min(outputIndex[dim], outputIndex[dim] + m_Offset[dim]);
269  inputSize[dim] =
270  std::max(outputIndex[dim] + outputSize[dim], outputIndex[dim] + outputSize[dim] +
271  m_Offset[dim]) - inputIndex[dim];
272  }
273 
274  // Build the input requested region
275  InputRegionType inputRequestedRegion;
276  inputRequestedRegion.SetIndex(inputIndex);
277  inputRequestedRegion.SetSize(inputSize);
278 
279  // Apply the radius
280  inputRequestedRegion.PadByRadius(m_Radius);
281 
282  // Try to apply the requested region to the input image
283  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
284  {
285  inputPtr->SetRequestedRegion(inputRequestedRegion);
286  }
287  else
288  {
289  // Build an exception
290  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
291  e.SetLocation(ITK_LOCATION);
292  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
293  e.SetDataObject(inputPtr);
294  throw e;
295  }
296 }
297 
298 template <class TInputImage, class TOutputImage>
299 void
302 {
303  unsigned int minRadius = 0;
304  for ( unsigned int i = 0; i < m_Offset.GetOffsetDimension(); i++ )
305  {
306  unsigned int distance = std::abs(m_Offset[i]);
307  if ( distance > minRadius )
308  {
309  minRadius = distance;
310  }
311  }
312  m_NeighborhoodRadius.Fill(minRadius);
313 }
314 
315 template <class TInputImage, class TOutputImage>
316 void
318 ::ThreadedGenerateData(const OutputRegionType& outputRegionForThread, itk::ThreadIdType threadId)
319 {
320  // Retrieve the input and output pointers
321  InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
322  OutputImagePointerType meanPtr = this->GetMeanOutput();
323  OutputImagePointerType variancePtr = this->GetVarianceOutput();
324  OutputImagePointerType dissimilarityPtr = this->GetDissimilarityOutput();
325  OutputImagePointerType sumAveragePtr = this->GetSumAverageOutput();
326  OutputImagePointerType sumVariancePtr = this->GetSumVarianceOutput();
327  OutputImagePointerType sumEntropytPtr = this->GetSumEntropyOutput();
328  OutputImagePointerType differenceEntropyPtr = this->GetDifferenceEntropyOutput();
329  OutputImagePointerType differenceVariancePtr = this->GetDifferenceVarianceOutput();
330  OutputImagePointerType ic1Ptr = this->GetIC1Output();
331  OutputImagePointerType ic2Ptr = this->GetIC2Output();
332 
333  // Build output iterators
334  itk::ImageRegionIteratorWithIndex<OutputImageType> varianceIt(variancePtr, outputRegionForThread);
335  itk::ImageRegionIterator<OutputImageType> meanIt(meanPtr, outputRegionForThread);
336  itk::ImageRegionIterator<OutputImageType> dissimilarityIt(dissimilarityPtr, outputRegionForThread);
337  itk::ImageRegionIterator<OutputImageType> sumAverageIt(sumAveragePtr, outputRegionForThread);
338  itk::ImageRegionIterator<OutputImageType> sumVarianceIt(sumVariancePtr, outputRegionForThread);
339  itk::ImageRegionIterator<OutputImageType> sumEntropytIt(sumEntropytPtr, outputRegionForThread);
340  itk::ImageRegionIterator<OutputImageType> differenceEntropyIt(differenceEntropyPtr, outputRegionForThread);
341  itk::ImageRegionIterator<OutputImageType> differenceVarianceIt(differenceVariancePtr, outputRegionForThread);
342  itk::ImageRegionIterator<OutputImageType> ic1It(ic1Ptr, outputRegionForThread);
343  itk::ImageRegionIterator<OutputImageType> ic2It(ic2Ptr, outputRegionForThread);
344 
345  // Go to begin
346  varianceIt.GoToBegin();
347  meanIt.GoToBegin();
348  dissimilarityIt.GoToBegin();
349  sumAverageIt.GoToBegin();
350  sumVarianceIt.GoToBegin();
351  sumEntropytIt.GoToBegin();
352  differenceEntropyIt.GoToBegin();
353  differenceVarianceIt.GoToBegin();
354  ic1It.GoToBegin();
355  ic2It.GoToBegin();
356 
357  const double log2 = std::log(2.0);
358  const unsigned int histSize = m_NumberOfBinsPerAxis;
359  const long unsigned int twiceHistSize = 2 * m_NumberOfBinsPerAxis;
360 
361  InputRegionType inputLargest = inputPtr->GetLargestPossibleRegion();
362 
363  // Set-up progress reporting
364  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
365 
366  // Iterate on outputs to compute textures
367  while (!varianceIt.IsAtEnd()
368  && !meanIt.IsAtEnd()
369  && !dissimilarityIt.IsAtEnd()
370  && !sumAverageIt.IsAtEnd()
371  && !sumVarianceIt.IsAtEnd()
372  && !sumEntropytIt.IsAtEnd()
373  && !differenceEntropyIt.IsAtEnd()
374  && !differenceVarianceIt.IsAtEnd()
375  && !ic1It.IsAtEnd()
376  && !ic2It.IsAtEnd())
377  {
378  // Compute the region on which co-occurence will be estimated
379  typename InputRegionType::IndexType inputIndex;
380  typename InputRegionType::SizeType inputSize;
381 
382  // Convert index to full grid
383  typename OutputImageType::IndexType outIndex;
384 
385  // First, create an window for neighborhood iterator based on m_Radius
386  // For example, if xradius and yradius is 2. window size is 5x5 (2 * radius + 1).
387  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
388  {
389  outIndex[dim] = varianceIt.GetIndex()[dim] * m_SubsampleFactor[dim]
390  + m_SubsampleOffset[dim] + inputLargest.GetIndex(dim);
391  inputIndex[dim] = outIndex[dim] - m_Radius[dim];
392  inputSize[dim] = 2 * m_Radius[dim] + 1;
393  }
394 
395  // Build the input region
396  InputRegionType inputRegion;
397  inputRegion.SetIndex(inputIndex);
398  inputRegion.SetSize(inputSize);
399  inputRegion.Crop(inputPtr->GetRequestedRegion());
400 
401  CooccurrenceIndexedListPointerType GLCIList = CooccurrenceIndexedListType::New();
402  GLCIList->Initialize(m_NumberOfBinsPerAxis, m_InputImageMinimum, m_InputImageMaximum);
403 
404  typedef itk::ConstNeighborhoodIterator< InputImageType > NeighborhoodIteratorType;
405  NeighborhoodIteratorType neighborIt;
406  neighborIt = NeighborhoodIteratorType(m_NeighborhoodRadius, inputPtr, inputRegion);
407  for ( neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt )
408  {
409  const InputPixelType centerPixelIntensity = neighborIt.GetCenterPixel();
410  bool pixelInBounds;
411  const InputPixelType pixelIntensity = neighborIt.GetPixel(m_Offset, pixelInBounds);
412  if ( !pixelInBounds )
413  {
414  continue; // don't put a pixel in the co-occurrence list if the value is
415  // out of bounds
416  }
417  GLCIList->AddPixelPair(centerPixelIntensity, pixelIntensity);
418  }
419 
430 
431  double Entropy = 0;
432 
433  typedef itk::Array<double> DoubleArrayType;
434  DoubleArrayType hx(histSize);
435  DoubleArrayType hy(histSize);
436  DoubleArrayType pdxy(twiceHistSize);
437 
438  for(long unsigned int i = 0; i < histSize; i++)
439  {
440  hx[i] = 0.0;
441  hy[i] = 0.0;
442  pdxy[i] = 0.0;
443  }
444  for(long unsigned int i = histSize; i < twiceHistSize; i++)
445  {
446  pdxy[i] = 0.0;
447  }
448 
449  /* hx.Fill(0.0); hy.Fill(0.0); pdxy.Fill(0.0); */
450  double hxy1 = 0;
451 
452  //get co-occurrence vector and totalfrequency
453  VectorType glcVector = GLCIList->GetVector();
454  double totalFrequency = static_cast<double> (GLCIList->GetTotalFrequency());
455 
456  VectorConstIteratorType constVectorIt;
457  //Normalize the GreyLevelCooccurrenceListType
458  //Compute Mean, Entropy (f12), hx, hy, pdxy
459  constVectorIt = glcVector.begin();
460  while( constVectorIt != glcVector.end())
461  {
462  CooccurrenceIndexType index = (*constVectorIt).first;
463  double frequency = (*constVectorIt).second / totalFrequency;
464  m_Mean += static_cast<double>(index[0]) * frequency;
465  Entropy -= (frequency > 0.0001) ? frequency * std::log(frequency) / log2 : 0.;
466  unsigned int i = index[1];
467  unsigned int j = index[0];
468  hx[j] += frequency;
469  hy[i] += frequency;
470 
471  if( i+j > histSize-1)
472  {
473  pdxy[i+j] += frequency;
474  }
475  if( i <= j )
476  {
477  pdxy[j-i] += frequency;
478  }
479  ++constVectorIt;
480  }
481 
482  //second pass over normalized co-occurrence list to find variance and pipj.
483  //pipj is needed to calculate f11
484  constVectorIt = glcVector.begin();
485  while( constVectorIt != glcVector.end())
486  {
487  double frequency = (*constVectorIt).second / totalFrequency;
488  CooccurrenceIndexType index = (*constVectorIt).first;
489  unsigned int i = index[1];
490  unsigned int j = index[0];
491  double index0 = static_cast<double>(index[0]);
492  m_Variance += ((index0 - m_Mean) * (index0 - m_Mean)) * frequency;
493  double pipj = hx[j] * hy[i];
494  hxy1 -= (pipj > 0.0001) ? frequency * std::log(pipj) : 0.;
495  ++constVectorIt;
496  }
497 
498  //iterate histSize to compute sumEntropy
499  double PSSquareCumul = 0;
500  for(long unsigned int k = histSize; k < twiceHistSize; k++)
501  {
502  m_SumAverage += k * pdxy[k];
503  m_SumEntropy -= (pdxy[k] > 0.0001) ? pdxy[k] * std::log(pdxy[k]) / log2 : 0;
504  PSSquareCumul += k * k * pdxy[k];
505  }
506  m_SumVariance = PSSquareCumul - m_SumAverage * m_SumAverage;
507 
508  double PDSquareCumul = 0;
509  double PDCumul = 0;
510  double hxCumul = 0;
511  double hyCumul = 0;
512 
513  for (long unsigned int i = 0; i < histSize; ++i)
514  {
515  double pdTmp = pdxy[i];
516  PDCumul += i * pdTmp;
517  m_DifferenceEntropy -= (pdTmp > 0.0001) ? pdTmp * std::log(pdTmp) / log2 : 0;
518  PDSquareCumul += i * i * pdTmp;
519 
520  //comput hxCumul and hyCumul
521  double marginalfreq = hx[i];
522  hxCumul += (marginalfreq > 0.0001) ? std::log (marginalfreq) * marginalfreq : 0;
523 
524  marginalfreq = hy[i];
525  hyCumul += (marginalfreq > 0.0001) ? std::log (marginalfreq) * marginalfreq : 0;
526  }
527  m_DifferenceVariance = PDSquareCumul - PDCumul * PDCumul;
528 
529  /* pipj computed below is totally different from earlier one which was used
530  * to compute hxy1. */
531  double hxy2 = 0;
532  for(unsigned int i = 0; i < histSize; ++i)
533  {
534  for(unsigned int j = 0; j < histSize; ++j)
535  {
536  double pipj = hx[j] * hy[i];
537  hxy2 -= (pipj > 0.0001) ? pipj * std::log(pipj) : 0.;
538  double frequency = GLCIList->GetFrequency(i,j, glcVector) / totalFrequency;
539  m_Dissimilarity+= ( static_cast<double>(j) - static_cast<double>(i) ) * (frequency * frequency);
540  }
541  }
542 
543  //Information measures of correlation 1 & 2
544  m_IC1 = (std::abs(std::max (hxCumul, hyCumul)) > 0.0001) ? (Entropy - hxy1) / (std::max (hxCumul, hyCumul)) : 0;
545  m_IC2 = 1 - std::exp (-2. * std::abs (hxy2 - Entropy));
546  m_IC2 = (m_IC2 >= 0) ? std::sqrt (m_IC2) : 0;
547 
548  // Fill outputs
549  meanIt.Set(m_Mean);
550  varianceIt.Set(m_Variance);
551  dissimilarityIt.Set(m_Dissimilarity);
552  sumAverageIt.Set(m_SumAverage);
553  sumVarianceIt.Set(m_SumVariance);
554  sumEntropytIt.Set(m_SumEntropy);
555  differenceEntropyIt.Set(m_DifferenceEntropy);
556  differenceVarianceIt.Set(m_DifferenceVariance);
557  ic1It.Set(m_IC1);
558  ic2It.Set(m_IC2);
559 
560  // Update progress
561  progress.CompletedPixel();
562 
563  // Increment iterators
564  ++varianceIt;
565  ++meanIt;
566  ++dissimilarityIt;
567  ++sumAverageIt;
568  ++sumVarianceIt;
569  ++sumEntropytIt;
570  ++differenceEntropyIt;
571  ++differenceVarianceIt;
572  ++ic1It;
573  ++ic2It;
574  }
575 
576 }
577 
578 } // End namespace otb
579 
580 #endif
itk::Size< Monteverdi_DIMENSION > SizeType
Definition: mvdTypes.h:137
void SetDataObject(DataObject *dobj)
CooccurrenceIndexedListType::PixelValueType PixelValueType
virtual void SetDescription(const std::string &s)
virtual void SetLocation(const std::string &s)
void Set(const PixelType &value) const
void Set(const PixelType &value) const
const IndexType & GetIndex() const
itk::Index< Monteverdi_DIMENSION > IndexType
Definition: mvdTypes.h:133
void ThreadedGenerateData(const OutputRegionType &outputRegion, itk::ThreadIdType threadId) override
unsigned int ThreadIdType
In this case, 10 advanced texture features will be processed. The 10 output image channels are: Mean...
VectorImageType::SpacingType SpacingType
Definition: mvdTypes.h:181
bool IsAtEnd(void) const
TOutputImage OutputImageType
VectorImageType::PointType PointType
Definition: mvdTypes.h:189
const SizeValueType * GetSize() const