OTB  6.7.0
Orfeo Toolbox
otbBCOInterpolateImageFunction.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 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 otbBCOInterpolateImageFunction_hxx
22 #define otbBCOInterpolateImageFunction_hxx
23 
25 
26 #include "itkNumericTraits.h"
27 
28 namespace otb
29 {
30 
31 template <class TInputImage, class TCoordRep>
33 ::PrintSelf(std::ostream& os, itk::Indent indent) const
34 {
35  Superclass::PrintSelf(os, indent);
36  os << indent << "Radius: " << m_Radius << std::endl;
37  os << indent << "Alpha: " << m_Alpha << std::endl;
38 }
39 
40 template <class TInputImage, class TCoordRep>
42 ::SetRadius(unsigned int radius)
43 {
44  if (radius < 2)
45  {
46  itkExceptionMacro(<< "Radius must be strictly greater than 1");
47  }
48  else
49  {
50  m_Radius = radius;
51  m_WinSize = 2*m_Radius+1;
52  }
53 }
54 
55 template <class TInputImage, class TCoordRep>
57 ::GetRadius() const
58 {
59  return m_Radius;
60 }
61 
62 template <class TInputImage, class TCoordRep>
64 ::SetAlpha(double alpha)
65 {
66  m_Alpha = alpha;
67 }
68 
69 template <class TInputImage, class TCoordRep>
71 ::GetAlpha() const
72 {
73  return m_Alpha;
74 }
75 
76 template<class TInputImage, class TCoordRep>
78 ::CoefContainerType
80 ::EvaluateCoef( const ContinuousIndexValueType & indexValue ) const
81 {
82  // Init BCO coefficient container
83 
84  CoefContainerType BCOCoef(m_WinSize, 0.);
85  double offset, dist, position, step;
86 
87  offset = indexValue - itk::Math::Floor<IndexValueType>(indexValue+0.5);
88 
89  // Compute BCO coefficients
90  step = 4./static_cast<double>(2*m_Radius);
91  position = - static_cast<double>(m_Radius) * step;
92 
93  double sum = 0.0;
94 
95  for ( unsigned int i = 0; i < m_WinSize; ++i)
96  {
97 
98  // Compute the BCO coefficients according to alpha.
99  dist = std::abs(position - offset*step);
100 
101  if( dist <= 2. )
102  {
103  if (dist <= 1.)
104  {
105  BCOCoef[i] = (m_Alpha + 2.)*std::abs(dist * dist * dist)
106  - (m_Alpha + 3.)*dist*dist + 1;
107  }
108  else
109  {
110  BCOCoef[i] = m_Alpha*std::abs(dist * dist * dist) - 5
111  *m_Alpha*dist*dist + 8*m_Alpha*std::abs(dist) - 4*m_Alpha;
112  }
113  }
114  else
115  {
116  BCOCoef[i] = 0;
117  }
118 
119  sum += BCOCoef[i];
120  position += step;
121  }
122 
123  for ( unsigned int i = 0; i < m_WinSize; ++i)
124  BCOCoef[i] = BCOCoef[i] / sum;
125 
126  return BCOCoef;
127 }
128 
129 template <class TInputImage, class TCoordRep>
131 ::PrintSelf(std::ostream& os, itk::Indent indent) const
132 {
133  Superclass::PrintSelf(os, indent);
134 }
135 
136 template <class TInputImage, class TCoordRep>
138 ::OutputType
141 {
142 
143  unsigned int dim;
144 
145  IndexType baseIndex;
146  IndexType neighIndex;
147 
149 
150  CoefContainerType BCOCoefX = this->EvaluateCoef(index[0]);
151  CoefContainerType BCOCoefY = this->EvaluateCoef(index[1]);
152 
153  // Compute base index = closet index
154  for( dim = 0; dim < ImageDimension; dim++ )
155  {
156  baseIndex[dim] = itk::Math::Floor< IndexValueType >( index[dim]+0.5 );
157  }
158 
159  for(unsigned int i = 0; i < this->m_WinSize; ++i )
160  {
161  RealType lineRes = 0.;
162  for(unsigned int j = 0; j < this->m_WinSize; ++j )
163  {
164  // get neighbor index
165  neighIndex[0] = baseIndex[0] + i - this->m_Radius;
166  neighIndex[1] = baseIndex[1] + j - this->m_Radius;
167 
168  if( neighIndex[0] > this->m_EndIndex[0] )
169  {
170  neighIndex[0] = this->m_EndIndex[0];
171  }
172  if( neighIndex[0] < this->m_StartIndex[0] )
173  {
174  neighIndex[0] = this->m_StartIndex[0];
175  }
176  if( neighIndex[1] > this->m_EndIndex[1] )
177  {
178  neighIndex[1] = this->m_EndIndex[1];
179  }
180  if( neighIndex[1] < this->m_StartIndex[1] )
181  {
182  neighIndex[1] = this->m_StartIndex[1];
183  }
184  lineRes += static_cast<RealType>( this->GetInputImage()->GetPixel( neighIndex ) ) * BCOCoefY[j];
185  }
186  value += lineRes*BCOCoefX[i];
187  }
188 
189 
190  return ( static_cast<OutputType>( value ) );
191 }
192 
193 template < typename TPixel, unsigned int VImageDimension, class TCoordRep >
195 ::PrintSelf(std::ostream& os, itk::Indent indent) const
196 {
197  Superclass::PrintSelf(os, indent);
198 }
199 
200 template < typename TPixel, unsigned int VImageDimension, class TCoordRep >
202 ::OutputType
204 ::EvaluateAtContinuousIndex( const ContinuousIndexType & index ) const
205 {
206  typedef typename itk::NumericTraits<InputPixelType>::ScalarRealType ScalarRealType;
207 
208 
209  unsigned int dim;
210  unsigned int componentNumber = this->GetInputImage()->GetNumberOfComponentsPerPixel();
211 
212  IndexType baseIndex;
213  IndexType neighIndex;
214 
215 
216  std::vector<ScalarRealType> lineRes(componentNumber);
217  OutputType output(componentNumber);
219 
220  CoefContainerType BCOCoefX = this->EvaluateCoef(index[0]);
221  CoefContainerType BCOCoefY = this->EvaluateCoef(index[1]);
222 
223  //Compute base index = closet index
224  for( dim = 0; dim < ImageDimension; dim++ )
225  {
226  baseIndex[dim] = itk::Math::Floor< IndexValueType >( index[dim]+0.5 );
227  }
228 
229  for(unsigned int i = 0; i < this->m_WinSize; ++i )
230  {
231  std::fill(lineRes.begin(), lineRes.end(), itk::NumericTraits<ScalarRealType>::Zero);
232  for(unsigned int j = 0; j < this->m_WinSize; ++j )
233  {
234  // get neighbor index
235  neighIndex[0] = baseIndex[0] + i - this->m_Radius;
236  neighIndex[1] = baseIndex[1] + j - this->m_Radius;
237  if( neighIndex[0] > this->m_EndIndex[0] )
238  {
239  neighIndex[0] = this->m_EndIndex[0];
240  }
241  if( neighIndex[0] < this->m_StartIndex[0] )
242  {
243  neighIndex[0] = this->m_StartIndex[0];
244  }
245  if( neighIndex[1] > this->m_EndIndex[1] )
246  {
247  neighIndex[1] = this->m_EndIndex[1];
248  }
249  if( neighIndex[1] < this->m_StartIndex[1] )
250  {
251  neighIndex[1] = this->m_StartIndex[1];
252  }
253 
254  const InputPixelType & pixel = this->GetInputImage()->GetPixel( neighIndex );
255  for( unsigned int k = 0; k<componentNumber; ++k)
256  {
257  lineRes[k] += pixel.GetElement(k) * BCOCoefY[j];
258  }
259  }
260  for( unsigned int k = 0; k<componentNumber; ++k)
261  {
262  output[k] += lineRes[k]*BCOCoefX[i];
263  }
264  }
265 
266  return ( output );
267 }
268 
269 } //namespace otb
270 
271 #endif
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Interpolate an image at specified positions using bicubic interpolation.
Superclass::IndexType IndexType
void PrintSelf(std::ostream &os, itk::Indent indent) const override
RealType ScalarRealType
OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &index) const override
virtual CoefContainerType EvaluateCoef(const ContinuousIndexValueType &indexValue) const
NumericTraits< typename TImageType::PixelType >::RealType RealType