OTB  6.7.0
Orfeo Toolbox
otbGenericInterpolateImageFunction.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 otbGenericInterpolateImageFunction_hxx
22 #define otbGenericInterpolateImageFunction_hxx
24 #include "vnl/vnl_math.h"
25 
26 namespace otb
27 {
28 
30 template<class TInputImage, class TFunction, class TBoundaryCondition, class TCoordRep>
33 {
34  m_WindowSize = 1;
35  this->SetRadius(1);
36  m_OffsetTable = nullptr;
37  m_WeightOffsetTable = nullptr;
38  m_TablesHaveBeenGenerated = false;
39  m_NormalizeWeight = false;
40 }
42 
44 template<class TInputImage, class TFunction, class TBoundaryCondition, class TCoordRep>
47 {
48  this->ResetOffsetTable();
49 }
50 
52 template<class TInputImage, class TFunction, class TBoundaryCondition, class TCoordRep>
53 void
56 {
57  // Clear the offset table
58  if (m_OffsetTable != nullptr)
59  {
60  delete[] m_OffsetTable;
61  m_OffsetTable = nullptr;
62  }
64 
65  // Clear the weights tales
66  if (m_WeightOffsetTable != nullptr)
67  {
68  for (unsigned int i = 0; i < m_OffsetTableSize; ++i)
69  {
70  delete[] m_WeightOffsetTable[i];
71  }
72  delete[] m_WeightOffsetTable;
73  m_WeightOffsetTable = nullptr;
74  }
75 }
76 
77 template<class TInputImage, class TFunction, class TBoundaryCondition, class TCoordRep>
78 void
80 ::SetRadius(unsigned int rad)
81 {
82  //m_Radius = rad;
83  this->GetFunction().SetRadius(rad);
84  m_WindowSize = rad << 1;
85  this->Modified();
86 }
87 
88 template<class TInputImage, class TFunction, class TBoundaryCondition, class TCoordRep>
89 void
91 ::Modified() const
92 {
93  Superclass::Modified();
94  m_TablesHaveBeenGenerated = false;
95 
96 }
97 
99 template<class TInputImage, class TFunction, class TBoundaryCondition, class TCoordRep>
100 void
103 {
104  // Compute the offset table size
105  m_OffsetTableSize = 1;
106  for (unsigned dim = 0; dim < ImageDimension; ++dim)
107  {
108  m_OffsetTableSize *= m_WindowSize;
109  }
111 
112  // Allocate the offset table
113  m_OffsetTable = new unsigned int[m_OffsetTableSize];
114 
115  // Allocate the weights tables
116  m_WeightOffsetTable = new unsigned int *[m_OffsetTableSize];
117  for (unsigned int i = 0; i < m_OffsetTableSize; ++i)
118  {
119  m_WeightOffsetTable[i] = new unsigned int[ImageDimension];
120  }
121 }
122 
124 template<class TInputImage, class TFunction, class TBoundaryCondition, class TCoordRep>
125 void
128 {
129  // Initialize the neighborhood
130  SizeType radius;
131  radius.Fill(this->GetRadius());
132  if (this->GetInputImage() != nullptr)
133  {
134  IteratorType it = IteratorType(radius, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
135  // Compute the offset tables (we ignore all the zero indices
136  // in the neighborhood)
137  unsigned int iOffset = 0;
138  int empty = static_cast<int>(this->GetRadius());
140 
141  for (unsigned int iPos = 0; iPos < it.Size(); ++iPos)
142  {
143  // Get the offset (index)
144  typename IteratorType::OffsetType off = it.GetOffset(iPos);
145 
146  // Check if the offset has zero weights
147  bool nonzero = true;
148  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
149  {
150  if (off[dim] == -empty)
151  {
152  nonzero = false;
153  break;
154  }
155  }
156  // Only use offsets with non-zero indices
157  if (nonzero)
158  {
159  // Set the offset index
160  m_OffsetTable[iOffset] = iPos;
161 
162  // Set the weight table indices
163  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
164  {
165  m_WeightOffsetTable[iOffset][dim] = off[dim] + this->GetRadius() - 1;
166  }
167  // Increment the index
168  iOffset++;
169  }
170  }
171  }
172  else
173  {
174  itkExceptionMacro(<< "An input has to be set");
175  }
176 }
177 
179 template<class TInputImage, class TFunction, class TBoundaryCondition, class TCoordRep>
180 void
183 {
184  // Delete existing tables
185  this->ResetOffsetTable();
186  // Tables initialization
187  this->InitializeTables();
188  // fill the weight table
189  this->FillWeightOffsetTable();
190  m_TablesHaveBeenGenerated = true;
191 }
193 
195 template<class TInputImage, class TFunction, class TBoundaryCondition, class TCoordRep>
199 {
200  if (!m_TablesHaveBeenGenerated)
201  {
202  itkExceptionMacro(<< "The Interpolation functor need to be explicitly intanciated with the method Initialize()");
203  }
204 
205  //unsigned int dim;
206  IndexType baseIndex;
207  double distance[ImageDimension];
208 
209  // Compute the integer index based on the continuous one by
210  // 'flooring' the index
211  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
212  {
213  // The following "if" block is equivalent to the following line without
214  // having to call floor.
215  // baseIndex[dim] = (long) std::floor(index[dim] );
216  if (index[dim] >= 0.0)
217  {
218  baseIndex[dim] = (long) index[dim];
219  }
220  else
221  {
222  long tIndex = (long) index[dim];
223  if (double(tIndex) != index[dim])
224  {
225  tIndex--;
226  }
227  baseIndex[dim] = tIndex;
228  }
229  distance[dim] = index[dim] - double(baseIndex[dim]);
230  }
231 
232  // Position the neighborhood at the index of interest
233  SizeType radius;
234  radius.Fill(this->GetRadius());
235  IteratorType nit = IteratorType(radius, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
236  nit.SetLocation(baseIndex);
237 
238  const unsigned int twiceRadius = static_cast<const unsigned int>(2 * this->GetRadius());
239  /* double xWeight[ImageDimension][ twiceRadius]; */
240  std::vector<std::vector<double> > xWeight;
241  xWeight.resize(ImageDimension);
242  for (unsigned int cpt = 0; cpt < xWeight.size(); ++cpt)
243  {
244  xWeight[cpt].resize(twiceRadius);
245  }
246 
247  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
248  {
249  // x is the offset, hence the parameter of the kernel
250  double x = distance[dim] + this->GetRadius();
251 
252  // If distance is zero, i.e. the index falls precisely on the
253  // pixel boundary, the weights form a delta function.
254  /*
255  if(distance[dim] == 0.0)
256  {
257  for( unsigned int i = 0; i < m_WindowSize; ++i)
258  {
259  xWeight[dim][i] = static_cast<int>(i) == (static_cast<int>(this->GetRadius()) - 1) ? 1. : 0.;
260  }
261  }
262  else
263  {
264  */
265  // i is the relative offset in dimension dim.
266  for (unsigned int i = 0; i < m_WindowSize; ++i)
267  {
268  // Increment the offset, taking it through the range
269  // (dist + rad - 1, ..., dist - rad), i.e. all x
270  // such that std::abs(x) <= rad
271  x -= 1.0;
272  // Compute the weight for this m
273  xWeight[dim][i] = m_Function(x);
274  }
275  //}
276  }
277  if (m_NormalizeWeight == true)
278  {
279  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
280  {
281  double sum = 0.;
282  // Compute the weights sum
283  for (unsigned int i = 0; i < m_WindowSize; ++i)
284  {
285  sum += xWeight[dim][i];
286  }
287  if (sum != 1.)
288  {
289  // Normalize the weights
290  for (unsigned int i = 0; i < m_WindowSize; ++i)
291  {
292  xWeight[dim][i] = xWeight[dim][i] / sum;
293  }
294  }
295  }
296  }
297 
298  // Iterate over the neighborhood, taking the correct set
299  // of weights in each dimension
300  RealType xPixelValue;
301  itk::NumericTraits<RealType>::SetLength(xPixelValue, this->GetInputImage()->GetNumberOfComponentsPerPixel());
302  xPixelValue=static_cast<RealType>(0.0);
303 
304  for (unsigned int j = 0; j < m_OffsetTableSize; ++j)
305  {
306  // Get the offset for this neighbor
307  unsigned int off = m_OffsetTable[j];
308 
309  // Get the intensity value at the pixel
310  RealType xVal = nit.GetPixel(off);
311 
312  // Multiply the intensity by each of the weights. Gotta hope
313  // that the compiler will unwrap this loop and pipeline this!
314  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
315  {
316  xVal *= xWeight[dim][m_WeightOffsetTable[j][dim]];
317  }
318 
319  // Increment the pixel value
320  xPixelValue += xVal;
321  }
322 
323  // Return the interpolated value
324  return static_cast<OutputType>(xPixelValue);
325 }
326 
327 template<class TInputImage, class TFunction, class TBoundaryCondition, class TCoordRep>
328 void
330 ::PrintSelf(std::ostream& os, itk::Indent indent) const
331 {
332  Superclass::PrintSelf(os, indent);
333 }
334 
335 } //namespace otb
336 
337 #endif
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Superclass::IndexType IndexType
Superclass::OutputType OutputType
NeighborIndexType Size() const
OffsetType GetOffset(NeighborIndexType i) const
void SetLocation(const IndexType &position)
virtual PixelType GetPixel(NeighborIndexType i, bool &IsInBounds) const
static void SetLength(T &m, const unsigned int s)
OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &index) const override
NumericTraits< typename TImageType::PixelType >::RealType RealType