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