21 #ifndef otbGenericInterpolateImageFunction_hxx
22 #define otbGenericInterpolateImageFunction_hxx
24 #include "vnl/vnl_math.h"
30 template <
class TInputImage,
class TFunction,
class TBoundaryCondition,
class TCoordRep>
35 m_OffsetTable =
nullptr;
36 m_WeightOffsetTable =
nullptr;
37 m_TablesHaveBeenGenerated =
false;
38 m_NormalizeWeight =
false;
43 template <
class TInputImage,
class TFunction,
class TBoundaryCondition,
class TCoordRep>
46 this->ResetOffsetTable();
50 template <
class TInputImage,
class TFunction,
class TBoundaryCondition,
class TCoordRep>
54 if (m_OffsetTable !=
nullptr)
56 delete[] m_OffsetTable;
57 m_OffsetTable =
nullptr;
62 if (m_WeightOffsetTable !=
nullptr)
64 for (
unsigned int i = 0; i < m_OffsetTableSize; ++i)
66 delete[] m_WeightOffsetTable[i];
68 delete[] m_WeightOffsetTable;
69 m_WeightOffsetTable =
nullptr;
73 template <
class TInputImage,
class TFunction,
class TBoundaryCondition,
class TCoordRep>
77 this->GetFunction().SetRadius(rad);
78 m_WindowSize = rad << 1;
82 template <
class TInputImage,
class TFunction,
class TBoundaryCondition,
class TCoordRep>
85 Superclass::Modified();
86 m_TablesHaveBeenGenerated =
false;
90 template <
class TInputImage,
class TFunction,
class TBoundaryCondition,
class TCoordRep>
94 m_OffsetTableSize = 1;
95 for (
unsigned dim = 0; dim < ImageDimension; ++dim)
97 m_OffsetTableSize *= m_WindowSize;
102 m_OffsetTable =
new unsigned int[m_OffsetTableSize];
105 m_WeightOffsetTable =
new unsigned int*[m_OffsetTableSize];
106 for (
unsigned int i = 0; i < m_OffsetTableSize; ++i)
108 m_WeightOffsetTable[i] =
new unsigned int[ImageDimension];
113 template <
class TInputImage,
class TFunction,
class TBoundaryCondition,
class TCoordRep>
118 radius.Fill(this->GetRadius());
119 if (this->GetInputImage() !=
nullptr)
124 unsigned int iOffset = 0;
125 int empty =
static_cast<int>(this->GetRadius());
128 for (
unsigned int iPos = 0; iPos < it.Size(); ++iPos)
131 typename IteratorType::OffsetType off = it.GetOffset(iPos);
135 for (
unsigned int dim = 0; dim < ImageDimension; ++dim)
137 if (off[dim] == -empty)
147 m_OffsetTable[iOffset] = iPos;
150 for (
unsigned int dim = 0; dim < ImageDimension; ++dim)
152 m_WeightOffsetTable[iOffset][dim] = off[dim] + this->GetRadius() - 1;
161 itkExceptionMacro(<<
"An input has to be set");
166 template <
class TInputImage,
class TFunction,
class TBoundaryCondition,
class TCoordRep>
170 this->ResetOffsetTable();
172 this->InitializeTables();
174 this->FillWeightOffsetTable();
175 m_TablesHaveBeenGenerated =
true;
180 template <
class TInputImage,
class TFunction,
class TBoundaryCondition,
class TCoordRep>
184 if (!m_TablesHaveBeenGenerated)
186 itkExceptionMacro(<<
"The Interpolation functor need to be explicitly intanciated with the method Initialize()");
191 double distance[ImageDimension];
195 for (
unsigned int dim = 0; dim < ImageDimension; ++dim)
200 if (index[dim] >= 0.0)
202 baseIndex[dim] = (long)index[dim];
206 long tIndex = (long)index[dim];
207 if (
double(tIndex) != index[dim])
211 baseIndex[dim] = tIndex;
213 distance[dim] = index[dim] - double(baseIndex[dim]);
218 radius.Fill(this->GetRadius());
220 nit.SetLocation(baseIndex);
222 const unsigned int twiceRadius =
static_cast<const unsigned int>(2 * this->GetRadius());
224 std::vector<std::vector<double>> xWeight;
225 xWeight.resize(ImageDimension);
226 for (
unsigned int cpt = 0; cpt < xWeight.size(); ++cpt)
228 xWeight[cpt].resize(twiceRadius);
231 for (
unsigned int dim = 0; dim < ImageDimension; ++dim)
234 double x = distance[dim] + this->GetRadius();
250 for (
unsigned int i = 0; i < m_WindowSize; ++i)
257 xWeight[dim][i] = m_Function(x);
261 if (m_NormalizeWeight ==
true)
263 for (
unsigned int dim = 0; dim < ImageDimension; ++dim)
267 for (
unsigned int i = 0; i < m_WindowSize; ++i)
269 sum += xWeight[dim][i];
274 for (
unsigned int i = 0; i < m_WindowSize; ++i)
276 xWeight[dim][i] = xWeight[dim][i] / sum;
285 itk::NumericTraits<RealType>::SetLength(xPixelValue, this->GetInputImage()->GetNumberOfComponentsPerPixel());
286 xPixelValue =
static_cast<RealType>(0.0);
288 for (
unsigned int j = 0; j < m_OffsetTableSize; ++j)
291 unsigned int off = m_OffsetTable[j];
298 for (
unsigned int dim = 0; dim < ImageDimension; ++dim)
300 xVal *= xWeight[dim][m_WeightOffsetTable[j][dim]];
311 template <
class TInputImage,
class TFunction,
class TBoundaryCondition,
class TCoordRep>
314 Superclass::PrintSelf(os, indent);