OTB  6.7.0
Orfeo Toolbox
otbGridResampleImageFilter.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 otbGridResampleImageFilter_hxx
22 #define otbGridResampleImageFilter_hxx
23 
25 
26 #include "otbStreamingTraits.h"
27 #include "otbImage.h"
28 
29 #include "itkNumericTraits.h"
30 #include "itkProgressReporter.h"
32 #include "itkContinuousIndex.h"
33 
34 namespace otb
35 {
36 
37 template <typename TInputImage, typename TOutputImage,
38  typename TInterpolatorPrecision>
41  : m_OutputStartIndex(),
42  m_OutputSize(),
43  m_OutputOrigin(),
44  m_OutputSpacing(),
45  m_EdgePaddingValue(),
46  m_CheckOutputBounds(true),
47  m_Interpolator(),
48  m_ReachableOutputRegion()
49 {
50  // Set linear interpolator as default
51  m_Interpolator = dynamic_cast<InterpolatorType *>(DefaultInterpolatorType::New().GetPointer());
52 
53  // Initialize EdgePaddingValue
54  m_EdgePaddingValue = itk::NumericTraits<OutputPixelType>::ZeroValue(m_EdgePaddingValue);
55 
56  // Initialize origin and spacing
57  m_OutputOrigin.Fill(0.);
58  m_OutputSpacing.Fill(1.);
59  m_OutputStartIndex.Fill(0);
60  m_OutputSize.Fill(0);
61 }
62 
63 
65 template <typename TInputImage, typename TOutputImage,
66  typename TInterpolatorPrecision>
67 void
70 {
71  this->SetOutputOrigin ( image->GetOrigin() );
72  this->SetOutputSpacing ( internal::GetSignedSpacing( image ) );
73  this->SetOutputStartIndex ( image->GetLargestPossibleRegion().GetIndex() );
74  this->SetOutputSize ( image->GetLargestPossibleRegion().GetSize() );
75 }
77 
78 template <typename TInputImage, typename TOutputImage,
79  typename TInterpolatorPrecision>
80 void
83 {
84  // call the superclass' implementation of this method
85  Superclass::GenerateOutputInformation();
86 
87  // get pointers to the input and output
88  OutputImageType *outputPtr = this->GetOutput();
89  if ( !outputPtr )
90  {
91  return;
92  }
93 
94  // Fill output image data
95  typename TOutputImage::RegionType outputLargestPossibleRegion;
96  outputLargestPossibleRegion.SetSize(m_OutputSize);
97  outputLargestPossibleRegion.SetIndex(m_OutputStartIndex);
98 
99  outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
100  outputPtr->SetSignedSpacing(m_OutputSpacing);
101  outputPtr->SetOrigin(m_OutputOrigin);
102 
103  // TODO: Report no data value here
104 }
105 
106 template <typename TInputImage, typename TOutputImage,
107  typename TInterpolatorPrecision>
108 void
111 {
112  // call the superclass's implementation of this method
113  Superclass::GenerateInputRequestedRegion();
114 
115  // Check for input
116  if ( !this->GetInput() )
117  {
118  return;
119  }
120 
121  // get pointers to the input and output
122  typename TInputImage::Pointer inputPtr =
123  const_cast< TInputImage * >( this->GetInput() );
124 
125  // check for output
126  OutputImageType *outputPtr = this->GetOutput();
127  if ( !outputPtr )
128  {
129  return;
130  }
131 
132  typename TOutputImage::RegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
133 
134  IndexType outULIndex,outLRIndex;
135 
137 
138  ContinuousIndexType inULCIndex,inLRCIndex;
139 
140  // Get output image requested region corners as Index
141  outULIndex = outputRequestedRegion.GetIndex();
142  outLRIndex = outULIndex+outputRequestedRegion.GetSize();
143  outLRIndex[0]-=1;
144  outLRIndex[1]-=1;
145 
146  // Transform to physiscal points
147  PointType outULPoint,outLRPoint;
148  outputPtr->TransformIndexToPhysicalPoint(outULIndex,outULPoint);
149  outputPtr->TransformIndexToPhysicalPoint(outLRIndex,outLRPoint);
150 
151  // Transform to input image Index
152  inputPtr->TransformPhysicalPointToContinuousIndex(outULPoint,inULCIndex);
153  inputPtr->TransformPhysicalPointToContinuousIndex(outLRPoint,inLRCIndex);
154 
155  SizeType inSize;
156  IndexType inULIndex,inLRIndex;
157 
158  // Reorder index properly and compute size
159  for(unsigned int dim = 0; dim < ImageDimension;++dim)
160  {
161  if(inULCIndex[dim] > inLRCIndex[dim])
162  {
163  // swap
164  typename ContinuousIndexType::ValueType tmp(inULCIndex[dim]);
165  inULCIndex[dim]=inLRCIndex[dim];
166  inLRCIndex[dim]=tmp;
167  }
168 
169  // Ensure correct rounding of coordinates
170 
171  inULIndex[dim] = std::floor(inULCIndex[dim]);
172  inLRIndex[dim] = std::ceil(inLRCIndex[dim]);
173 
174  inSize[dim] = static_cast<typename SizeType::SizeValueType>(inLRIndex[dim]-inULIndex[dim])+1;
175  }
176 
177  // Build the input requested region
178  typename TInputImage::RegionType inputRequestedRegion;
179  inputRequestedRegion.SetIndex(inULIndex);
180  inputRequestedRegion.SetSize(inSize);
181 
182  // Compute the padding due to the interpolator
183  unsigned int interpolatorRadius =
185  inputRequestedRegion.PadByRadius(interpolatorRadius);
186 
187  // crop the input requested region at the input's largest possible region
188  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
189  {
190  inputPtr->SetRequestedRegion(inputRequestedRegion);
191  }
192  else
193  {
194 
195  // store what we tried to request (prior to trying to crop)
196  inputPtr->SetRequestedRegion(inputRequestedRegion);
197 
198  // build an exception
199  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
200  e.SetLocation(ITK_LOCATION);
201  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
202  e.SetDataObject(inputPtr);
203  throw e;
204  }
205 
206 }
207 
208 
209 template <typename TInputImage, typename TOutputImage,
210  typename TInterpolatorPrecision>
211 void
214 {
215  if ( !m_Interpolator )
216  {
217  itkExceptionMacro(<< "Interpolator not set");
218  }
219 
220  // Connect input image to interpolator
221  m_Interpolator->SetInputImage( this->GetInput() );
222 
223 
224  unsigned int nComponents
226  m_EdgePaddingValue );
227 
228  if (nComponents == 0)
229  {
230 
231  // Build a default value of the correct number of components
232  OutputPixelComponentType zeroComponent
234 
235  nComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
236 
237  itk::NumericTraits<OutputPixelType>::SetLength(m_EdgePaddingValue, nComponents );
238  for (unsigned int n=0; n<nComponents; n++)
239  {
240  OutputPixelConvertType::SetNthComponent( n, m_EdgePaddingValue,
241  zeroComponent );
242  }
243  }
244 
245  // Compute ReachableOutputRegion
246  // InputImage buffered region corresponds to a region of the output
247  // image. Computing it beforehand allows saving IsInsideBuffer
248  // calls in the interpolation loop
249 
250  // Compute the padding due to the interpolator
251 
252 
253  IndexType inUL = this->GetInput()->GetBufferedRegion().GetIndex();
254  IndexType inLR = this->GetInput()->GetBufferedRegion().GetIndex() + this->GetInput()->GetBufferedRegion().GetSize();
255  inLR[0]-=1;
256  inLR[1]-=1;
257 
258  // We should take interpolation radius into account here, but this
259  // does not match the IsInsideBuffer method
260  // unsigned int interpolatorRadius =
261  // StreamingTraits<typename Superclass::InputImageType>::CalculateNeededRadiusForInterpolator(this->GetInterpolator());
262  // inUL[0]+=interpolatorRadius;
263  // inUL[1]+=interpolatorRadius;
264  // inLR[0]-=interpolatorRadius;
265  // inLR[1]-=interpolatorRadius;
266 
267  PointType inULp, inLRp;
268  this->GetInput()->TransformIndexToPhysicalPoint(inUL,inULp);
269  this->GetInput()->TransformIndexToPhysicalPoint(inLR,inLRp);
270 
271  inULp-=0.5*this->GetInput()->GetSignedSpacing();
272  inLRp+=0.5*this->GetInput()->GetSignedSpacing();
273 
276  this->GetOutput()->TransformPhysicalPointToContinuousIndex(inULp,outUL);
277  this->GetOutput()->TransformPhysicalPointToContinuousIndex(inLRp,outLR);
278 
279  IndexType outputIndex;
280  // This needs to take into account negative spacing
281  outputIndex[0] = std::ceil(std::min(outUL[0],outLR[0]));
282  outputIndex[1] = std::ceil(std::min(outUL[1],outLR[1]));
283 
284  SizeType outputSize;
285  outputSize[0] = std::floor(std::max(outUL[0],outLR[0])) - outputIndex[0] + 1;
286  outputSize[1] = std::floor(std::max(outUL[1],outLR[1])) - outputIndex[1] + 1;
287 
288  m_ReachableOutputRegion.SetIndex(outputIndex);
289  m_ReachableOutputRegion.SetSize(outputSize);
290 }
291 
292 template <typename TInputImage, typename TOutputImage,
293  typename TInterpolatorPrecision>
294 void
296 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
297 {
298  // Get the output pointers
299  OutputImageType *outputPtr = this->GetOutput();
300 
301  // Get this input pointers
302  const InputImageType *inputPtr = this->GetInput();
303 
304  // Min/max values of the output pixel type AND these values
305  // represented as the output type of the interpolator
308 
309  const InterpolatorComponentType minOutputValue = static_cast< InterpolatorComponentType >( minValue );
310  const InterpolatorComponentType maxOutputValue = static_cast< InterpolatorComponentType >( maxValue );
311 
312  // Iterator on the output region for current thread
313  OutputImageRegionType regionToCompute = outputRegionForThread;
314  bool cropSucceed = regionToCompute.Crop(m_ReachableOutputRegion);
315 
316  // Fill thread buffer
317  itk::ImageScanlineIterator<OutputImageType> outItFull(outputPtr,outputRegionForThread);
318  outItFull.GoToBegin();
319  while(!outItFull.IsAtEnd())
320  {
321  while(!outItFull.IsAtEndOfLine())
322  {
323  outItFull.Set(m_EdgePaddingValue);
324  ++outItFull;
325  }
326  outItFull.NextLine();
327  }
328 
329  if(!cropSucceed)
330  return;
331 
332  itk::ImageScanlineIterator<OutputImageType> outIt(outputPtr, regionToCompute);
333 
334  // Support for progress methods/callbacks
335  itk::ProgressReporter progress( this,
336  threadId,
337  regionToCompute.GetSize()[1]);
338 
339  // Temporary variables for loop
340  PointType outPoint;
341  ContinuousInputIndexType inCIndex;
342  InterpolatorOutputType interpolatorValue; //(this->GetOutput()->GetNumberOfComponentsPerPixel());
343  OutputPixelType outputValue; //(this->GetOutput()->GetNumberOfComponentsPerPixel());
344 
345  // TODO: assert outputPtr->GetSignedSpacing() != 0 here
346  assert(outputPtr->GetSignedSpacing()[0]!=0&&"Null spacing will cause division by zero.");
347  const double delta = outputPtr->GetSignedSpacing()[0]/inputPtr->GetSignedSpacing()[0];
348 
349  // Iterate through the output region
350  outIt.GoToBegin();
351 
352  while(!outIt.IsAtEnd())
353  {
354  // Map output index to input continuous index
355  outputPtr->TransformIndexToPhysicalPoint(outIt.GetIndex(),outPoint);
356  inputPtr->TransformPhysicalPointToContinuousIndex(outPoint,inCIndex);
357 
358  while(!outIt.IsAtEndOfLine())
359  {
360  // Interpolate
361  interpolatorValue = m_Interpolator->EvaluateAtContinuousIndex(inCIndex);
362 
363  // Cast and check bounds
364  this->CastPixelWithBoundsChecking(interpolatorValue,minOutputValue,maxOutputValue,outputValue);
365 
366  // Set output value
367  outIt.Set(outputValue);
368 
369  // move one pixel forward
370  ++outIt;
371 
372  // Update input position
373  inCIndex[0]+=delta;
374  }
375 
376  // Report progress
377  progress.CompletedPixel();
378 
379  // Move to next line
380  outIt.NextLine();
381  }
382 
383 }
384 
385 template <typename TInputImage, typename TOutputImage,
386  typename TInterpolatorPrecision>
387 void
390 {
391  // Disconnect input image from the interpolator
392  m_Interpolator->SetInputImage(nullptr);
393 }
394 
395 
396 template <typename TInputImage, typename TOutputImage,
397  typename TInterpolatorPrecision>
400 ::GetMTime(void) const
401 {
403 
404  if ( m_Interpolator )
405  {
406  if ( latestTime < m_Interpolator->GetMTime() )
407  {
408  latestTime = m_Interpolator->GetMTime();
409  }
410  }
411 
412  return latestTime;
413 }
414 
415 
416 template <typename TInputImage, typename TOutputImage,
417  typename TInterpolatorPrecision>
418 void
420 ::PrintSelf(std::ostream& os, itk::Indent indent) const
421 {
422  Superclass::PrintSelf(os, indent);
423 
424  os << indent << "EdgePaddingValue: "
425  << static_cast< typename itk::NumericTraits< OutputPixelType >::PrintType >
426  ( m_EdgePaddingValue )
427  << std::endl;
428  os << indent << "OutputStartIndex: " << m_OutputStartIndex << std::endl;
429  os << indent << "OutputSize: " << m_OutputSize << std::endl;
430  os << indent << "OutputOrigin: " << m_OutputOrigin << std::endl;
431  os << indent << "OutputSpacing: " << m_OutputSpacing << std::endl;
432  os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
433  os << indent << "CheckOutputBounds: " << ( m_CheckOutputBounds ? "On" : "Off" )
434  << std::endl;
435 }
436 
437 
438 
439 } // namespace otb
440 
441 
442 #endif
const IndexValueType * GetIndex() const
static unsigned int CalculateNeededRadiusForInterpolator(const InterpolationType *interpolator)
void SetDataObject(DataObject *dobj)
void SetOutputParametersFromImage(const ImageBaseType *image)
const IndexType & GetIndex() const
itk::SizeValueType SizeValueType
virtual void SetDescription(const std::string &s)
virtual void SetLocation(const std::string &s)
itk::ModifiedTimeType GetMTime(void) const override
virtual const RegionType & GetLargestPossibleRegion() const
unsigned long ModifiedTimeType
const IndexType GetIndex() const
static T ZeroValue()
void SetSize(const SizeValueType val[VDimension])
ImageType::SpacingType GetSignedSpacing(const ImageType *input)
Definition: otbImage.h:41
static ITK_CONSTEXPR_FUNC T max(const T &)
OutputPixelConvertType::ComponentType OutputPixelComponentType
TInputImage InputImageType
static void SetLength(T &m, const unsigned int s)
const SizeType & GetSize() const
static unsigned int GetNumberOfComponents()
virtual ModifiedTimeType GetMTime() const
unsigned int ThreadIdType
TOutputImage::PixelType OutputPixelType
InterpolatorConvertType::ComponentType InterpolatorComponentType
static ITK_CONSTEXPR_FUNC T NonpositiveMin()
void Set(const PixelType &value) const
InterpolatorType::OutputType InterpolatorOutputType
OutputImageType::RegionType OutputImageRegionType
bool IsAtEnd(void) const
TOutputImage OutputImageType
void PrintSelf(std::ostream &os, itk::Indent indent) const override
virtual const PointType & GetOrigin() const
void SetIndex(const IndexValueType val[VIndexDimension])
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) override