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