Orfeo Toolbox  4.0
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 
23 #include "itkImageRegionIterator.h"
24 #include "itkProgressReporter.h"
25 
26 namespace otb
27 {
28 template <class TInputImage, class TOutputImage>
31  m_NumberOfBinsPerAxis(8),
32  m_InputImageMinimum(0),
33  m_InputImageMaximum(256)
34 {
35  // There are 1 output corresponding to the Pan Tex texture indice
36  this->SetNumberOfRequiredOutputs(1);
37 
38  //Fill the offset list for contrast computation
39  OffsetType off;
40  off[0] = 0;
41  off[1] = 1;
42  m_OffsetList.push_back(off); //(0, 1)
43  off[1] = 2;
44  m_OffsetList.push_back(off); //(0, 2)
45  off[0] = 1;
46  off[1] = -2;
47  m_OffsetList.push_back(off); //(1, -2)
48  off[1] = -1;
49  m_OffsetList.push_back(off); //(1, -1)
50  off[1] = 0;
51  m_OffsetList.push_back(off); //(1, 0)
52  off[1] = 1;
53  m_OffsetList.push_back(off); //(1, 1)
54  off[1] = 2;
55  m_OffsetList.push_back(off); //(1, 2)
56  off[0] = 2;
57  off[1] = -1;
58  m_OffsetList.push_back(off); //(2, -1)
59  off[1] = 0;
60  m_OffsetList.push_back(off); //(2, 0)
61  off[1] = 1;
62  m_OffsetList.push_back(off); //(2, 1)
63 }
64 
65 template <class TInputImage, class TOutputImage>
68 {}
69 
70 template <class TInputImage, class TOutputImage>
71 void
74 {
75  // First, call superclass implementation
76  Superclass::GenerateInputRequestedRegion();
77 
78  // Retrieve the input and output pointers
79  InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
80  OutputImagePointerType outputPtr = this->GetOutput();
81 
82  if (!inputPtr || !outputPtr)
83  {
84  return;
85  }
86 
87  // Retrieve the output requested region
88  // We use only the first output since requested regions for all outputs are enforced to be equal
89  // by the default GenerateOutputRequestedRegiont() implementation
90  OutputRegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
91 
92  // Build the input requested region
93  InputRegionType inputRequestedRegion = outputRequestedRegion;
94 
95  // Apply the radius
96  SizeType maxOffsetSize;
97  maxOffsetSize[0] = 2;
98  maxOffsetSize[1] = 2;
99  inputRequestedRegion.PadByRadius(m_Radius + maxOffsetSize);
100 
101  // Try to apply the requested region to the input image
102  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
103  {
104  inputPtr->SetRequestedRegion(inputRequestedRegion);
105  }
106  else
107  {
108  // Build an exception
109  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
110  e.SetLocation(ITK_LOCATION);
111  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
112  e.SetDataObject(inputPtr);
113  throw e;
114  }
115 }
116 
117 template <class TInputImage, class TOutputImage>
118 void
120 ::ThreadedGenerateData(const OutputRegionType& outputRegionForThread, itk::ThreadIdType threadId)
121 {
122  // Retrieve the input and output pointers
123  InputImagePointerType inputPtr = const_cast<InputImageType *> (this->GetInput());
124  OutputImagePointerType outputPtr = this->GetOutput();
125 
126  itk::ImageRegionIteratorWithIndex<OutputImageType> outputIt(outputPtr, outputRegionForThread);
127 
128  // Go to begin
129  outputIt.GoToBegin();
130 
131  // Set-up progress reporting
132  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
133 
134  // Iterate on outputs to compute textures
135  while (!outputIt.IsAtEnd())
136  {
137  // Initialise output value
138  double out = itk::NumericTraits<double>::max();
139 
140  // For each offset
141  typename OffsetListType::const_iterator offIt;
142  for (offIt = m_OffsetList.begin(); offIt != m_OffsetList.end(); ++offIt)
143  {
144  OffsetType currentOffset = *offIt;
145 
146  // Compute the region on which co-occurence will be estimated
147  typename InputRegionType::IndexType inputIndex, inputIndexWithTwiceOffset;
148  typename InputRegionType::SizeType inputSize, inputSizeWithTwiceOffset;
149 
150  // First, apply offset
151  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
152  {
153  inputIndex[dim] = std::min(
154  static_cast<int>(outputIt.GetIndex()[dim] - m_Radius[dim]),
155  static_cast<int>(outputIt.GetIndex()[dim] - m_Radius[dim] + currentOffset[dim])
156  );
157  inputSize[dim] = 2 * m_Radius[dim] + 1 + std::abs(currentOffset[dim]);
158 
159  inputIndexWithTwiceOffset[dim] = static_cast<int>(outputIt.GetIndex()[dim] - m_Radius[dim] - std::abs(currentOffset[dim]));
160  inputSizeWithTwiceOffset[dim] = inputSize[dim] + std::abs(currentOffset[dim]);
161  }
162 
163  // Build the input region
164  InputRegionType inputRegion;
165  inputRegion.SetIndex(inputIndex);
166  inputRegion.SetSize(inputSize);
167  inputRegion.Crop(inputPtr->GetRequestedRegion());
168 
169  InputRegionType inputRegionWithTwiceOffset;
170  inputRegionWithTwiceOffset.SetIndex(inputIndexWithTwiceOffset);
171  inputRegionWithTwiceOffset.SetSize(inputSizeWithTwiceOffset);
172  inputRegionWithTwiceOffset.Crop(inputPtr->GetRequestedRegion());
173 
174  /*********************************************************************************/
175  //Local copy of the input image around the processed pixel *outputIt
176  InputImagePointerType localInputImage = InputImageType::New();
177  localInputImage->SetRegions(inputRegionWithTwiceOffset);
178  localInputImage->Allocate();
179  typedef itk::ImageRegionIteratorWithIndex<InputImageType> ImageRegionIteratorType;
180  ImageRegionIteratorType itInputPtr(inputPtr, inputRegionWithTwiceOffset);
181  ImageRegionIteratorType itLocalInputImage(localInputImage, inputRegionWithTwiceOffset);
182  for (itInputPtr.GoToBegin(), itLocalInputImage.GoToBegin(); !itInputPtr.IsAtEnd(); ++itInputPtr, ++itLocalInputImage)
183  {
184  itLocalInputImage.Set(itInputPtr.Get());
185  }
186  /*********************************************************************************/
187 
188  InputImagePointerType maskImage = InputImageType::New();
189  maskImage->SetRegions(inputRegionWithTwiceOffset);
190  maskImage->Allocate();
191  maskImage->FillBuffer(0);
192 
193  ImageRegionIteratorType itMask(maskImage, inputRegion);
194  for (itMask.GoToBegin(); !itMask.IsAtEnd(); ++itMask)
195  {
196  itMask.Set(1);
197  }
198 
199 
200  // Build the co-occurence matrix generator
201  CoocurrenceMatrixGeneratorPointerType coOccurenceMatrixGenerator = CoocurrenceMatrixGeneratorType::New();
202  coOccurenceMatrixGenerator->SetInput(localInputImage);
203  coOccurenceMatrixGenerator->SetOffset(currentOffset);
204  coOccurenceMatrixGenerator->SetNumberOfBinsPerAxis(m_NumberOfBinsPerAxis);
205  coOccurenceMatrixGenerator->SetPixelValueMinMax(m_InputImageMinimum, m_InputImageMaximum);
206 
207  // Compute the co-occurence matrix
208  coOccurenceMatrixGenerator->SetMaskImage(maskImage);
209  coOccurenceMatrixGenerator->SetInsidePixelValue(1);
210  coOccurenceMatrixGenerator->Update();
211 
212  typename HistogramType::ConstPointer histo = coOccurenceMatrixGenerator->GetOutput();
213 
214  double inertia = 0;
215  typename HistogramType::TotalAbsoluteFrequencyType totalFrequency = histo->GetTotalFrequency();
216  typename HistogramType::ConstIterator itr = histo->Begin();
217  typename HistogramType::ConstIterator end = histo->End();
218  for(; itr != end; ++itr )
219  {
220  MeasurementType frequency = itr.GetFrequency();
221  if (frequency == 0)
222  {
223  continue; // no use doing these calculations if we're just multiplying by zero.
224  }
225  typename HistogramType::IndexType index = histo->GetIndex(itr.GetInstanceIdentifier());
226  inertia += (index[0] - index[1]) * (index[0] - index[1]) * frequency / totalFrequency;
227  }
228 
229  if (inertia < out)
230  {
231  out = inertia;
232  }
233  }
234 
235  outputIt.Set(out);
236  ++outputIt;
237  progress.CompletedPixel();
238  }
239 }
240 
241 } // End namespace otb
242 
243 #endif

Generated at Sat Mar 8 2014 16:17:40 for Orfeo Toolbox with doxygen 1.8.3.1