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

Generated at Sat Mar 8 2014 15:57:03 for Orfeo Toolbox with doxygen 1.8.3.1