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