OTB  9.0.0
Orfeo Toolbox
otbScalarImageToPanTexTextureFilter.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 otbScalarImageToPanTexTextureFilter_hxx
22 #define otbScalarImageToPanTexTextureFilter_hxx
23 
25 #include "itkImageRegionIteratorWithIndex.h"
26 #include "itkConstNeighborhoodIterator.h"
27 #include "itkImageRegionIterator.h"
28 #include "itkProgressReporter.h"
29 #include "itkNumericTraits.h"
30 
31 namespace otb
32 {
33 template <class TInputImage, class TOutputImage>
35  : m_Radius(), m_NumberOfBinsPerAxis(8), m_InputImageMinimum(0), m_InputImageMaximum(255)
36 {
37  // There are 1 output corresponding to the Pan Tex texture indice
38  this->SetNumberOfRequiredOutputs(1);
39 
40  // Ten offsets are selected for contrast computation (2 pixels displacement grid, and given that the
41  // co-occurance matrix is symmetric
42  m_OffsetList = { {0, 1}, {0, 2}, {1, -2}, {1, -1}, {1, 0}, {1, 1}, {1, 2}, {2, -1}, {2, 0}, {2, 1} };
43 }
44 
45 template <class TInputImage, class TOutputImage>
47 {
48  // First, call superclass implementation
49  Superclass::GenerateInputRequestedRegion();
50 
51  // Retrieve the input and output pointers
52  InputImagePointerType inputPtr = const_cast<InputImageType*>(this->GetInput());
53  OutputImagePointerType outputPtr = this->GetOutput();
54 
55  if (!inputPtr || !outputPtr)
56  {
57  return;
58  }
59 
60  // Retrieve the output requested region
61  // We use only the first output since requested regions for all outputs are enforced to be equal
62  // by the default GenerateOutputRequestedRegiont() implementation
63  OutputRegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
64 
65  // Build the input requested region
66  InputRegionType inputRequestedRegion = outputRequestedRegion;
67 
68  // Apply the radius
69  SizeType maxOffsetSize = {2, 2};
70  inputRequestedRegion.PadByRadius(m_Radius + maxOffsetSize);
71 
72  // Try to apply the requested region to the input image
73  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
74  {
75  inputPtr->SetRequestedRegion(inputRequestedRegion);
76  }
77  else
78  {
79  // Build an exception
80  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
81  e.SetLocation(ITK_LOCATION);
82  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
83  e.SetDataObject(inputPtr);
84  throw e;
85  }
86 }
87 
88 template <class TInputImage, class TOutputImage>
90  itk::ThreadIdType threadId)
91 {
92  // Retrieve the input and output pointers
93  InputImagePointerType inputPtr = const_cast<InputImageType*>(this->GetInput());
94  OutputImagePointerType outputPtr = this->GetOutput();
95 
96  itk::ImageRegionIteratorWithIndex<OutputImageType> outputIt(outputPtr, outputRegionForThread);
97  outputIt.GoToBegin();
98 
99  // Set-up progress reporting
100  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
101 
102  // Iterate on outputs to compute textures
103  while (!outputIt.IsAtEnd())
104  {
105  // Initialise output value
106  double out = itk::NumericTraits<double>::max();
107 
108  // For each offset
109  typename OffsetListType::const_iterator offIt;
110  for (offIt = m_OffsetList.begin(); offIt != m_OffsetList.end(); ++offIt)
111  {
112  OffsetType currentOffset = *offIt;
113 
114  // Compute the region on which co-occurence will be estimated
115  typename InputRegionType::IndexType inputIndex;
116  typename InputRegionType::SizeType inputSize;
117 
118  // First, create an window for neighborhood iterator based on m_Radius
119  // For example, if xradius and yradius is 2. window size is 5x5 (2 *
120  // radius + 1).
121  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
122  {
123  inputIndex[dim] = outputIt.GetIndex()[dim] - m_Radius[dim];
124  inputSize[dim] = 2 * m_Radius[dim] + 1;
125  }
126 
127  // Build the input region
128  InputRegionType inputRegion = {inputIndex, inputSize};
129  inputRegion.Crop(inputPtr->GetRequestedRegion());
130 
131  SizeType neighborhoodRadius;
133  unsigned int minRadius = 0;
134  for (unsigned int i = 0; i < currentOffset.GetOffsetDimension(); i++)
135  {
136  unsigned int distance = std::abs(currentOffset[i]);
137  if (distance > minRadius)
138  {
139  minRadius = distance;
140  }
141  }
142  neighborhoodRadius.Fill(minRadius);
144 
145  CooccurrenceIndexedListPointerType GLCIList = CooccurrenceIndexedListType::New();
146  GLCIList->Initialize(m_NumberOfBinsPerAxis, m_InputImageMinimum, m_InputImageMaximum);
147 
148  typedef itk::ConstNeighborhoodIterator<InputImageType> NeighborhoodIteratorType;
149  NeighborhoodIteratorType neighborIt = NeighborhoodIteratorType(neighborhoodRadius, inputPtr, inputRegion);
150 
151  for (neighborIt.GoToBegin(); !neighborIt.IsAtEnd(); ++neighborIt)
152  {
153  const InputPixelType centerPixelIntensity = neighborIt.GetCenterPixel();
154  bool pixelInBounds;
155  const InputPixelType pixelIntensity = neighborIt.GetPixel(currentOffset, pixelInBounds);
156  if (!pixelInBounds)
157  {
158  continue; // don't put a pixel in the histogram if it's out-of-bounds.
159  }
160  GLCIList->AddPixelPair(centerPixelIntensity, pixelIntensity);
161  }
162 
163  VectorConstIteratorType constVectorIt;
164  VectorType glcVector = GLCIList->GetVector();
165  double totalFrequency = static_cast<double>(GLCIList->GetTotalFrequency());
166 
167  // Compute inertia aka contrast
168  double inertia = 0;
169  constVectorIt = glcVector.begin();
170  while (constVectorIt != glcVector.end())
171  {
172  CooccurrenceIndexType index = constVectorIt->first;
173  RelativeFrequencyType frequency = constVectorIt->second / totalFrequency;
174  inertia += (index[0] - index[1]) * (index[0] - index[1]) * frequency;
175  ++constVectorIt;
176  }
177 
178  if (inertia < out)
179  {
180  out = inertia;
181  }
182  }
183 
184  outputIt.Set(out);
185  ++outputIt;
186  progress.CompletedPixel();
187  }
188 }
189 
190 } // End namespace otb
191 
192 #endif
otb::ScalarImageToPanTexTextureFilter::InputRegionType
InputImageType::RegionType InputRegionType
Definition: otbScalarImageToPanTexTextureFilter.h:69
otb::ScalarImageToPanTexTextureFilter::InputImagePointerType
InputImageType::Pointer InputImagePointerType
Definition: otbScalarImageToPanTexTextureFilter.h:67
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::ScalarImageToPanTexTextureFilter::VectorType
CooccurrenceIndexedListType::VectorType VectorType
Definition: otbScalarImageToPanTexTextureFilter.h:83
otb::ScalarImageToPanTexTextureFilter::VectorConstIteratorType
VectorType::const_iterator VectorConstIteratorType
Definition: otbScalarImageToPanTexTextureFilter.h:86
otbScalarImageToPanTexTextureFilter.h
otb::ScalarImageToPanTexTextureFilter::OutputRegionType
OutputImageType::RegionType OutputRegionType
Definition: otbScalarImageToPanTexTextureFilter.h:75
otb::ScalarImageToPanTexTextureFilter::OffsetType
InputImageType::OffsetType OffsetType
Definition: otbScalarImageToPanTexTextureFilter.h:71
otb::ScalarImageToPanTexTextureFilter::RelativeFrequencyType
CooccurrenceIndexedListType::RelativeFrequencyType RelativeFrequencyType
Definition: otbScalarImageToPanTexTextureFilter.h:82
otb::ScalarImageToPanTexTextureFilter::CooccurrenceIndexedListPointerType
CooccurrenceIndexedListType::Pointer CooccurrenceIndexedListPointerType
Definition: otbScalarImageToPanTexTextureFilter.h:78
otb::ScalarImageToPanTexTextureFilter::ScalarImageToPanTexTextureFilter
ScalarImageToPanTexTextureFilter()
Definition: otbScalarImageToPanTexTextureFilter.hxx:34
otb::ScalarImageToPanTexTextureFilter::m_OffsetList
OffsetListType m_OffsetList
Definition: otbScalarImageToPanTexTextureFilter.h:137
otb::ScalarImageToPanTexTextureFilter::SizeType
InputRegionType::SizeType SizeType
Definition: otbScalarImageToPanTexTextureFilter.h:70
otb::ScalarImageToPanTexTextureFilter::GenerateInputRequestedRegion
void GenerateInputRequestedRegion() override
Definition: otbScalarImageToPanTexTextureFilter.hxx:46
otb::ScalarImageToPanTexTextureFilter::ThreadedGenerateData
void ThreadedGenerateData(const OutputRegionType &outputRegion, itk::ThreadIdType threadId) override
Definition: otbScalarImageToPanTexTextureFilter.hxx:89
otb::ScalarImageToPanTexTextureFilter::InputImageType
TInpuImage InputImageType
Definition: otbScalarImageToPanTexTextureFilter.h:63
otb::ScalarImageToPanTexTextureFilter::OutputImagePointerType
OutputImageType::Pointer OutputImagePointerType
Definition: otbScalarImageToPanTexTextureFilter.h:74
otb::ScalarImageToPanTexTextureFilter::CooccurrenceIndexType
CooccurrenceIndexedListType::IndexType CooccurrenceIndexType
Definition: otbScalarImageToPanTexTextureFilter.h:80
otb::ScalarImageToPanTexTextureFilter::InputPixelType
InputImageType::PixelType InputPixelType
Definition: otbScalarImageToPanTexTextureFilter.h:68