OTB  6.7.0
Orfeo Toolbox
otbScalarImageToHigherOrderTexturesFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 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 otbScalarImageToHigherOrderTexturesFilter_hxx
22 #define otbScalarImageToHigherOrderTexturesFilter_hxx
23 
26 #include "itkImageRegionIterator.h"
27 #include "itkProgressReporter.h"
28 
29 namespace otb
30 {
31 template <class TInputImage, class TOutputImage>
34  m_Radius(),
35  m_NumberOfBinsPerAxis(8),
36  m_InputImageMinimum(0),
37  m_InputImageMaximum(255),
38  m_FastCalculations(false),
39  m_SubsampleFactor(),
40  m_SubsampleOffset()
41 {
42  // There are 10 outputs corresponding to the 8 textures indices
43  this->SetNumberOfRequiredOutputs(10);
44 
45  // Create the 11 outputs
46  this->SetNthOutput(0, OutputImageType::New());
47  this->SetNthOutput(1, OutputImageType::New());
48  this->SetNthOutput(2, OutputImageType::New());
49  this->SetNthOutput(3, OutputImageType::New());
50  this->SetNthOutput(4, OutputImageType::New());
51  this->SetNthOutput(5, OutputImageType::New());
52  this->SetNthOutput(6, OutputImageType::New());
53  this->SetNthOutput(7, OutputImageType::New());
54  this->SetNthOutput(8, OutputImageType::New());
55  this->SetNthOutput(9, OutputImageType::New());
56 
57  m_Radius.Fill(10);
58 
59  // Set the offset directions to their defaults: half of all the possible
60  // directions 1 pixel away. (The other half is included by symmetry.)
61  // We use a neighborhood iterator to calculate the appropriate offsets.
62  typedef itk::Neighborhood<typename InputImageType::PixelType,
63  InputImageType::ImageDimension> NeighborhoodType;
64  NeighborhoodType hood;
65  hood.SetRadius( 1 );
66 
67  // select all "previous" neighbors that are face+edge+vertex
68  // connected to the current pixel. do not include the center pixel.
69  unsigned int centerIndex = hood.GetCenterNeighborhoodIndex();
70  OffsetVectorPointer offsets = OffsetVector::New();
71  for( unsigned int d = 0; d < centerIndex; d++ )
72  {
73  OffsetType offset = hood.GetOffset( d );
74  offsets->push_back( offset );
75  }
76  this->SetOffsets( offsets );
77 
78  this->m_SubsampleFactor.Fill(1);
79  this->m_SubsampleOffset.Fill(0);
80 }
81 
82 template <class TInputImage, class TOutputImage>
85 {}
86 
87 template <class TInputImage, class TOutputImage>
89 ::OutputImageType *
92 {
93  return this->GetOutput(0);
94 }
95 
96 template <class TInputImage, class TOutputImage>
98 ::OutputImageType *
101 {
102  return this->GetOutput(1);
103 }
104 
105 template <class TInputImage, class TOutputImage>
107 ::OutputImageType *
110 {
111  return this->GetOutput(2);
112 }
113 
114 template <class TInputImage, class TOutputImage>
116 ::OutputImageType *
119 {
120  return this->GetOutput(3);
121 }
122 
123 template <class TInputImage, class TOutputImage>
125 ::OutputImageType *
128 {
129  return this->GetOutput(4);
130 }
131 
132 template <class TInputImage, class TOutputImage>
134 ::OutputImageType *
137 {
138  return this->GetOutput(5);
139 }
140 
141 template <class TInputImage, class TOutputImage>
143 ::OutputImageType *
146 {
147  return this->GetOutput(6);
148 }
149 
150 template <class TInputImage, class TOutputImage>
152 ::OutputImageType *
155 {
156  return this->GetOutput(7);
157 }
158 
159 template <class TInputImage, class TOutputImage>
161 ::OutputImageType *
164 {
165  return this->GetOutput(8);
166 }
167 
168 template <class TInputImage, class TOutputImage>
170 ::OutputImageType *
173 {
174  return this->GetOutput(9);
175 }
176 
177 template <class TInputImage, class TOutputImage>
178 void
180 ::SetOffset( const OffsetType offset )
181 {
182  OffsetVectorPointer offsetVector = OffsetVector::New();
183  offsetVector->push_back( offset );
184  this->SetOffsets( offsetVector );
185 }
186 
187 template <class TInputImage, class TOutputImage>
188 void
191 {
192  // Compute output size, origin & spacing
193  const InputImageType * inputPtr = this->GetInput();
194  InputRegionType inputRegion = inputPtr->GetLargestPossibleRegion();
195  OutputRegionType outputRegion;
196  outputRegion.SetIndex(0,0);
197  outputRegion.SetIndex(1,0);
198  outputRegion.SetSize(0, 1 + (inputRegion.GetSize(0) - 1 - m_SubsampleOffset[0]) / m_SubsampleFactor[0]);
199  outputRegion.SetSize(1, 1 + (inputRegion.GetSize(1) - 1 - m_SubsampleOffset[1]) / m_SubsampleFactor[1]);
200 
201  typename OutputImageType::SpacingType outSpacing = inputPtr->GetSignedSpacing();
202  outSpacing[0] *= m_SubsampleFactor[0];
203  outSpacing[1] *= m_SubsampleFactor[1];
204 
205  typename OutputImageType::PointType outOrigin;
206  inputPtr->TransformIndexToPhysicalPoint(inputRegion.GetIndex()+m_SubsampleOffset,outOrigin);
207 
208  for (unsigned int i=0 ; i < this->GetNumberOfOutputs() ; i++)
209  {
210  OutputImagePointerType outputPtr = this->GetOutput(i);
211  outputPtr->CopyInformation(inputPtr);
212  outputPtr->SetLargestPossibleRegion(outputRegion);
213  outputPtr->SetOrigin(outOrigin);
214  outputPtr->SetSignedSpacing(outSpacing);
215  }
216 }
217 
218 
219 template <class TInputImage, class TOutputImage>
220 void
223 {
224  // First, call superclass implementation
225  Superclass::GenerateInputRequestedRegion();
226 
227  // Retrieve the input and output pointers
228  InputImagePointerType inputPtr = const_cast<InputImageType *>(this->GetInput());
229  OutputImagePointerType outputPtr = this->GetOutput();
230 
231  if (!inputPtr || !outputPtr)
232  {
233  return;
234  }
235 
236  // Retrieve the output requested region
237  // We use only the first output since requested regions for all outputs are enforced to be equal
238  // by the default GenerateOutputRequestedRegiont() implementation
239  OutputRegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
240  typename OutputRegionType::IndexType outputIndex = outputRequestedRegion.GetIndex();
241  typename OutputRegionType::SizeType outputSize = outputRequestedRegion.GetSize();
242  InputRegionType inputLargest = inputPtr->GetLargestPossibleRegion();
243 
244  // Convert index and size to full grid
245  outputIndex[0] = outputIndex[0] * m_SubsampleFactor[0] + m_SubsampleOffset[0] + inputLargest.GetIndex(0);
246  outputIndex[1] = outputIndex[1] * m_SubsampleFactor[1] + m_SubsampleOffset[1] + inputLargest.GetIndex(1);
247  outputSize[0] = 1 + (outputSize[0] - 1) * m_SubsampleFactor[0];
248  outputSize[1] = 1 + (outputSize[1] - 1) * m_SubsampleFactor[1];
249 
250  InputRegionType inputRequestedRegion(outputIndex,outputSize);
251 
252  // Apply the radius
253  inputRequestedRegion.PadByRadius(m_Radius);
254 
255  // Try to apply the requested region to the input image
256  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
257  {
258  inputPtr->SetRequestedRegion(inputRequestedRegion);
259  }
260  else
261  {
262  // Build an exception
263  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
264  e.SetLocation(ITK_LOCATION);
265  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
266  e.SetDataObject(inputPtr);
267  throw e;
268  }
269 }
270 
271 template <class TInputImage, class TOutputImage>
272 void
274 ::ThreadedGenerateData(const OutputRegionType& outputRegionForThread, itk::ThreadIdType threadId)
275 {
276  // Retrieve the input and output pointers
277  const InputImageType * inputPtr = this->GetInput();
278 
279  typedef typename itk::ImageRegionIterator<OutputImageType> IteratorType;
280  std::vector<IteratorType> outputImagesIterators;
281 
282  for (unsigned int i = 0; i < this->GetNumberOfOutputs(); ++i)
283  {
284  outputImagesIterators.push_back( IteratorType(this->GetOutput(i), outputRegionForThread) );
285  outputImagesIterators[i].GoToBegin();
286  }
287 
288  // Compute the max possible run length (in physical unit)
289  typename InputImageType::PointType topLeftPoint;
290  typename InputImageType::PointType bottomRightPoint;
291  inputPtr->TransformIndexToPhysicalPoint( outputImagesIterators[0].GetIndex() - m_Radius, topLeftPoint );
292  inputPtr->TransformIndexToPhysicalPoint( outputImagesIterators[0].GetIndex() + m_Radius, bottomRightPoint );
293  double maxDistance = topLeftPoint.EuclideanDistanceTo(bottomRightPoint);
294 
295  InputRegionType inputLargest = inputPtr->GetLargestPossibleRegion();
296 
297  // Set-up progress reporting
298  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
299 
300  // Iterate on outputs to compute textures
301  while ( !outputImagesIterators[0].IsAtEnd() )
302  {
303  // Compute the region on which run-length matrix will be estimated
304  typename InputRegionType::IndexType inputIndex;
305  typename InputRegionType::SizeType inputSize;
306 
307  // Convert index to full grid
308  typename OutputImageType::IndexType outIndex;
309 
310  for (unsigned int dim = 0; dim < InputImageType::ImageDimension; ++dim)
311  {
312  outIndex[dim] = outputImagesIterators[0].GetIndex()[dim] * m_SubsampleFactor[dim]
313  + m_SubsampleOffset[dim] + inputLargest.GetIndex(dim);
314  inputIndex[dim] = outIndex[dim] - m_Radius[dim];
315  inputSize[dim] = 2 * m_Radius[dim] + 1;
316  }
317 
318  // Build the input region
319  InputRegionType inputRegion;
320  inputRegion.SetIndex(inputIndex);
321  inputRegion.SetSize(inputSize);
322 
323  inputRegion.Crop(inputPtr->GetBufferedRegion());
324 
325  // Create a local image corresponding to the input region
326  InputImagePointerType localInputImage = InputImageType::New();
327  localInputImage->SetRegions(inputRegion);
328  localInputImage->Allocate();
329  typedef itk::ImageRegionIteratorWithIndex<InputImageType> ImageRegionIteratorType;
330  typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> ImageRegionConstIteratorType;
331  ImageRegionConstIteratorType itInputPtr(inputPtr, inputRegion);
332  ImageRegionIteratorType itLocalInputImage(localInputImage, inputRegion);
333  for (itInputPtr.GoToBegin(), itLocalInputImage.GoToBegin();
334  !itInputPtr.IsAtEnd();
335  ++itInputPtr, ++itLocalInputImage)
336  {
337  itLocalInputImage.Set(itInputPtr.Get());
338  }
339 
340  typename ScalarImageToRunLengthFeaturesFilterType::Pointer runLengthFeatureCalculator = ScalarImageToRunLengthFeaturesFilterType::New();
341  runLengthFeatureCalculator->SetInput(localInputImage);
342  runLengthFeatureCalculator->SetOffsets(m_Offsets);
343  runLengthFeatureCalculator->SetNumberOfBinsPerAxis(m_NumberOfBinsPerAxis);
344  runLengthFeatureCalculator->SetPixelValueMinMax(m_InputImageMinimum, m_InputImageMaximum);
345  runLengthFeatureCalculator->SetDistanceValueMinMax(0, maxDistance);
346 
347  runLengthFeatureCalculator->Update();
348 
350  featuresMeans = *(runLengthFeatureCalculator->GetFeatureMeans().GetPointer());
351 
352  // Fill output
353  for (unsigned int i = 0; i < this->GetNumberOfOutputs(); ++i)
354  {
355  // Fill output
356  outputImagesIterators[i].Set(featuresMeans[i]);
357  // Increment iterators
358  ++outputImagesIterators[i];
359  }
360  // Update progress
361  progress.CompletedPixel();
362  }
363 }
364 
365 } // End namespace otb
366 
367 #endif
itk::Size< Monteverdi_DIMENSION > SizeType
Definition: mvdTypes.h:137
void SetDataObject(DataObject *dobj)
This class compute 10 local higher order statistics textures coefficients based on the grey level run...
virtual void SetDescription(const std::string &s)
virtual void SetLocation(const std::string &s)
ObjectType * GetPointer() const
void ThreadedGenerateData(const OutputRegionType &outputRegion, itk::ThreadIdType threadId) override
itk::Index< Monteverdi_DIMENSION > IndexType
Definition: mvdTypes.h:133
unsigned int ThreadIdType
VectorImageType::SpacingType SpacingType
Definition: mvdTypes.h:181
RealType EuclideanDistanceTo(const Point< TCoordRepB, NPointDimension > &pa) const
void SetRadius(const SizeType &)
VectorImageType::PointType PointType
Definition: mvdTypes.h:189
const SizeValueType * GetSize() const