Orfeo Toolbox  4.0
otbScalarImageToTexturesFilter.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 __otbScalarImageToTexturesFilter_txx
19 #define __otbScalarImageToTexturesFilter_txx
20 
23 #include "itkImageRegionIterator.h"
24 #include "itkProgressReporter.h"
25 
26 namespace otb
27 {
28 template <class TInputImage, class TOutputImage>
31  m_Offset(),
32  m_NumberOfBinsPerAxis(8),
33  m_InputImageMinimum(0),
34  m_InputImageMaximum(256)
35 {
36  // There are 8 outputs corresponding to the 8 textures indices
37  this->SetNumberOfRequiredOutputs(8);
38 
39  // Create the 8 outputs
40  this->SetNthOutput(0, OutputImageType::New());
41  this->SetNthOutput(1, OutputImageType::New());
42  this->SetNthOutput(2, OutputImageType::New());
43  this->SetNthOutput(3, OutputImageType::New());
44  this->SetNthOutput(4, OutputImageType::New());
45  this->SetNthOutput(5, OutputImageType::New());
46  this->SetNthOutput(6, OutputImageType::New());
47  this->SetNthOutput(7, OutputImageType::New());
48 }
49 
50 template <class TInputImage, class TOutputImage>
53 {}
54 
55 template <class TInputImage, class TOutputImage>
57 ::OutputImageType *
60 {
61  if (this->GetNumberOfOutputs() < 1)
62  {
63  return 0;
64  }
65  return static_cast<OutputImageType *>(this->GetOutput(0));
66 }
67 
68 template <class TInputImage, class TOutputImage>
70 ::OutputImageType *
73 {
74  if (this->GetNumberOfOutputs() < 2)
75  {
76  return 0;
77  }
78  return static_cast<OutputImageType *>(this->GetOutput(1));
79 }
80 
81 template <class TInputImage, class TOutputImage>
83 ::OutputImageType *
86 {
87  if (this->GetNumberOfOutputs() < 3)
88  {
89  return 0;
90  }
91  return static_cast<OutputImageType *>(this->GetOutput(2));
92 }
93 
94 template <class TInputImage, class TOutputImage>
96 ::OutputImageType *
99 {
100  if (this->GetNumberOfOutputs() < 4)
101  {
102  return 0;
103  }
104  return static_cast<OutputImageType *>(this->GetOutput(3));
105 }
106 
107 template <class TInputImage, class TOutputImage>
109 ::OutputImageType *
112 {
113  if (this->GetNumberOfOutputs() < 5)
114  {
115  return 0;
116  }
117  return static_cast<OutputImageType *>(this->GetOutput(4));
118 }
119 
120 template <class TInputImage, class TOutputImage>
122 ::OutputImageType *
125 {
126  if (this->GetNumberOfOutputs() < 6)
127  {
128  return 0;
129  }
130  return static_cast<OutputImageType *>(this->GetOutput(5));
131 }
132 
133 template <class TInputImage, class TOutputImage>
135 ::OutputImageType *
138 {
139  if (this->GetNumberOfOutputs() < 7)
140  {
141  return 0;
142  }
143  return static_cast<OutputImageType *>(this->GetOutput(6));
144 }
145 
146 template <class TInputImage, class TOutputImage>
148 ::OutputImageType *
151 {
152  if (this->GetNumberOfOutputs() < 8)
153  {
154  return 0;
155  }
156  return static_cast<OutputImageType *>(this->GetOutput(7));
157 }
158 
159 template <class TInputImage, class TOutputImage>
160 void
163 {
164  // First, call superclass implementation
165  Superclass::GenerateInputRequestedRegion();
166 
167  // Retrieve the input and output pointers
168  InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
169  OutputImagePointerType outputPtr = this->GetOutput();
170 
171  if (!inputPtr || !outputPtr)
172  {
173  return;
174  }
175 
176  // Retrieve the output requested region
177  // We use only the first output since requested regions for all outputs are enforced to be equal
178  // by the default GenerateOutputRequestedRegiont() implementation
179  OutputRegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
180 
181  typename OutputRegionType::IndexType outputIndex = outputRequestedRegion.GetIndex();
182  typename OutputRegionType::SizeType outputSize = outputRequestedRegion.GetSize();
183  typename InputRegionType::IndexType inputIndex;
184  typename InputRegionType::SizeType inputSize;
185 
186  // First, apply offset
187  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
188  {
189  inputIndex[dim] = std::min(outputIndex[dim], outputIndex[dim] + m_Offset[dim]);
190  inputSize[dim] =
191  std::max(outputIndex[dim] + outputSize[dim], outputIndex[dim] + outputSize[dim] +
192  m_Offset[dim]) - inputIndex[dim];
193  }
194 
195  // Build the input requested region
196  InputRegionType inputRequestedRegion;
197  inputRequestedRegion.SetIndex(inputIndex);
198  inputRequestedRegion.SetSize(inputSize);
199 
200  // Apply the radius
201  inputRequestedRegion.PadByRadius(m_Radius);
202 
203  // Try to apply the requested region to the input image
204  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
205  {
206  inputPtr->SetRequestedRegion(inputRequestedRegion);
207  }
208  else
209  {
210  // Build an exception
211  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
212  e.SetLocation(ITK_LOCATION);
213  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
214  e.SetDataObject(inputPtr);
215  throw e;
216  }
217 }
218 
219 template <class TInputImage, class TOutputImage>
220 void
222 ::ThreadedGenerateData(const OutputRegionType& outputRegionForThread, itk::ThreadIdType threadId)
223 {
224  // Retrieve the input and output pointers
225  InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
226  OutputImagePointerType energyPtr = this->GetEnergyOutput();
227  OutputImagePointerType entropyPtr = this->GetEntropyOutput();
228  OutputImagePointerType correlationPtr = this->GetCorrelationOutput();
229  OutputImagePointerType invDiffMomentPtr = this->GetInverseDifferenceMomentOutput();
230  OutputImagePointerType inertiaPtr = this->GetInertiaOutput();
231  OutputImagePointerType clusterShadePtr = this->GetClusterShadeOutput();
232  OutputImagePointerType clusterProminencePtr = this->GetClusterProminenceOutput();
233  OutputImagePointerType haralickCorPtr = this->GetHaralickCorrelationOutput();
234 
235  // Build output iterators
236  itk::ImageRegionIteratorWithIndex<OutputImageType> energyIt(energyPtr, outputRegionForThread);
237  itk::ImageRegionIterator<OutputImageType> entropyIt(entropyPtr, outputRegionForThread);
238  itk::ImageRegionIterator<OutputImageType> correlationIt(correlationPtr, outputRegionForThread);
239  itk::ImageRegionIterator<OutputImageType> invDiffMomentIt(invDiffMomentPtr, outputRegionForThread);
240  itk::ImageRegionIterator<OutputImageType> inertiaIt(inertiaPtr, outputRegionForThread);
241  itk::ImageRegionIterator<OutputImageType> clusterShadeIt(clusterShadePtr, outputRegionForThread);
242  itk::ImageRegionIterator<OutputImageType> clusterProminenceIt(clusterProminencePtr, outputRegionForThread);
243  itk::ImageRegionIterator<OutputImageType> haralickCorIt(haralickCorPtr, outputRegionForThread);
244 
245  // Go to begin
246  energyIt.GoToBegin();
247  entropyIt.GoToBegin();
248  correlationIt.GoToBegin();
249  invDiffMomentIt.GoToBegin();
250  inertiaIt.GoToBegin();
251  clusterShadeIt.GoToBegin();
252  clusterProminenceIt.GoToBegin();
253  haralickCorIt.GoToBegin();
254 
255 
256  // Build the co-occurence matrix generator
257  CoocurrenceMatrixGeneratorPointerType coOccurenceMatrixGenerator = CoocurrenceMatrixGeneratorType::New();
258  coOccurenceMatrixGenerator->SetOffset(m_Offset);
259  coOccurenceMatrixGenerator->SetNumberOfBinsPerAxis(m_NumberOfBinsPerAxis);
260  coOccurenceMatrixGenerator->SetPixelValueMinMax(m_InputImageMinimum, m_InputImageMaximum);
261 
262  // Build the texture calculator
263  TextureCoefficientsCalculatorPointerType texturesCalculator = TextureCoefficientsCalculatorType::New();
264 
265  // Set-up progress reporting
266  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
267 
268  // Iterate on outputs to compute textures
269  while (!energyIt.IsAtEnd()
270  && !entropyIt.IsAtEnd()
271  && !correlationIt.IsAtEnd()
272  && !invDiffMomentIt.IsAtEnd()
273  && !inertiaIt.IsAtEnd()
274  && !clusterShadeIt.IsAtEnd()
275  && !clusterProminenceIt.IsAtEnd()
276  && !haralickCorIt.IsAtEnd())
277  {
278  // Compute the region on which co-occurence will be estimated
279  typename InputRegionType::IndexType inputIndex, inputIndexWithTwiceOffset;
280  typename InputRegionType::SizeType inputSize, inputSizeWithTwiceOffset;
281 
282  // First, apply offset
283  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
284  {
285  inputIndex[dim] = std::min(
286  static_cast<int>(energyIt.GetIndex()[dim] - m_Radius[dim]),
287  static_cast<int>(energyIt.GetIndex()[dim] - m_Radius[dim] + m_Offset[dim])
288  );
289  inputSize[dim] = 2 * m_Radius[dim] + 1 + std::abs(m_Offset[dim]);
290 
291  inputIndexWithTwiceOffset[dim] = static_cast<int>(energyIt.GetIndex()[dim] - m_Radius[dim] - std::abs(m_Offset[dim]));
292  inputSizeWithTwiceOffset[dim] = inputSize[dim] + std::abs(m_Offset[dim]);
293  }
294 
295  // Build the input region
296  InputRegionType inputRegion;
297  inputRegion.SetIndex(inputIndex);
298  inputRegion.SetSize(inputSize);
299  inputRegion.Crop(inputPtr->GetRequestedRegion());
300 
301  InputRegionType inputRegionWithTwiceOffset;
302  inputRegionWithTwiceOffset.SetIndex(inputIndexWithTwiceOffset);
303  inputRegionWithTwiceOffset.SetSize(inputSizeWithTwiceOffset);
304  inputRegionWithTwiceOffset.Crop(inputPtr->GetRequestedRegion());
305 
306  /*********************************************************************************/
307  //Local copy of the input image around the processed pixel *outputIt
308  InputImagePointerType localInputImage = InputImageType::New();
309  localInputImage->SetRegions(inputRegionWithTwiceOffset);
310  localInputImage->Allocate();
311  typedef itk::ImageRegionIteratorWithIndex<InputImageType> ImageRegionIteratorType;
312  ImageRegionIteratorType itInputPtr(inputPtr, inputRegionWithTwiceOffset);
313  ImageRegionIteratorType itLocalInputImage(localInputImage, inputRegionWithTwiceOffset);
314  for (itInputPtr.GoToBegin(), itLocalInputImage.GoToBegin(); !itInputPtr.IsAtEnd(); ++itInputPtr, ++itLocalInputImage)
315  {
316  itLocalInputImage.Set(itInputPtr.Get());
317  }
318  /*********************************************************************************/
319 
320  InputImagePointerType maskImage = InputImageType::New();
321  maskImage->SetRegions(inputRegionWithTwiceOffset);
322  maskImage->Allocate();
323  maskImage->FillBuffer(0);
324 
325  ImageRegionIteratorType itMask(maskImage, inputRegion);
326  for (itMask.GoToBegin(); !itMask.IsAtEnd(); ++itMask)
327  {
328  itMask.Set(1);
329  }
330 
331  // Compute the co-occurence matrix
332  coOccurenceMatrixGenerator->SetInput(localInputImage);
333  coOccurenceMatrixGenerator->SetMaskImage(maskImage);
334  coOccurenceMatrixGenerator->SetInsidePixelValue(1);
335  coOccurenceMatrixGenerator->Update();
336 
337  // Compute textures indices
338  texturesCalculator->SetInput(coOccurenceMatrixGenerator->GetOutput());
339  texturesCalculator->Update();
340 
341  // Fill outputs
342  energyIt.Set(texturesCalculator->GetEnergy());
343  entropyIt.Set(texturesCalculator->GetEntropy());
344  correlationIt.Set(texturesCalculator->GetCorrelation());
345  invDiffMomentIt.Set(texturesCalculator->GetInverseDifferenceMoment());
346  inertiaIt.Set(texturesCalculator->GetInertia());
347  clusterShadeIt.Set(texturesCalculator->GetClusterShade());
348  clusterProminenceIt.Set(texturesCalculator->GetClusterProminence());
349  haralickCorIt.Set(texturesCalculator->GetHaralickCorrelation());
350 
351  // Update progress
352  progress.CompletedPixel();
353 
354  // Increment iterators
355  ++energyIt;
356  ++entropyIt;
357  ++correlationIt;
358  ++invDiffMomentIt;
359  ++inertiaIt;
360  ++clusterShadeIt;
361  ++clusterProminenceIt;
362  ++haralickCorIt;
363  }
364 
365 }
366 
367 } // End namespace otb
368 
369 #endif

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