23 #ifndef INTERPOLATIONUTILITIES_H 24 #define INTERPOLATIONUTILITIES_H 39 template<
typename TReal>
46 inline TReal operator()
58 return ( y1 * (1-mu) + y2 * mu ) ;
68 template<
typename TReal>
80 return(y1*( 1. -mu2 )+y2*mu2);
96 template<
typename TReal>
105 inline TReal operator()
124 TReal a0,a1,a2,a3,mu2;
126 a0 = y3 - y2 - y0 + y1;
130 return(a0*mu*mu2 + a1*mu2 + a2*mu + a3);
143 template<
typename TReal>
153 ): tension_(tension), bias_(bias), epsilon_(
epsilon)
168 else if(-(mu - 1.) < epsilon_)
178 m0 = (y1-y0)*(1+bias_)*(1-tension_)/2;
179 m0 += (y2-y1)*(1-bias_)*(1-tension_)/2;
180 m1 = (y2-y1)*(1+bias_)*(1-tension_)/2;
181 m1 += (y3-y2)*(1-bias_)*(1-tension_)/2;
182 a0 = 2*mu3 - 3*mu2 + 1;
183 a1 = mu3 - 2*mu2 + mu;
186 return( a0*y1 + a1*m0 + a2*m1 + a3*y2 );
194 typename YInputIterator,
195 typename XInputIterator,
196 typename OutputIterator,
210 typedef typename std::iterator_traits<OutputIterator>::value_type TReal;
212 size_t nrY = std::distance( begY , endY );
213 OutputIterator outI = out;
215 for(
unsigned int i = 0 ; begX != endX ; ++i , ++begX, ++outI )
217 double xd = *begX - start_index;
218 int index =
static_cast<int>(floor(xd));
220 TReal mu = xd -
static_cast<double>(index);
230 TReal y4 = *(begY+1);
231 *outI = functor(y1 , y2 , y3 , y4 , mu );
238 TReal y3 = *(begY+1);
239 TReal y4 = *(begY+2);
240 *outI = functor(y1 , y2 , y3 , y4 , mu );
242 else if( index > 0 && index < static_cast<int>(nrY - 2) )
244 YInputIterator begTmp = (begY + index - 1);
246 TReal y2 = *(begTmp + 1);
247 TReal y3 = *(begTmp + 2);
248 TReal y4 = *(begTmp + 3);
249 *outI = functor(y1 , y2 , y3 , y4 , mu );
251 else if(index == static_cast<int>(nrY-2) )
253 YInputIterator begTmp = (begY + index - 1);
255 TReal y2 = *(begTmp+1);
256 TReal y3 = *(begTmp+2);
258 *outI = functor(y1 , y2 , y3 , y4 ,mu);
260 else if(index == static_cast<int>(nrY-1) )
262 YInputIterator begTmp = (begY + index - 1);
264 TReal y2 = *(begTmp+1);
265 TReal y3 = *(begTmp+1);
266 TReal y4 = *(begTmp+1);
267 *outI = functor(y1 , y2 , y3 , y4 , mu );
278 typename YInputIterator,
279 typename XInputIterator,
280 typename OutputIterator,
289 OutputIterator out_p,
290 TFunctor & interpolator,
294 typedef typename std::iterator_traits<OutputIterator>::value_type TReal ;
295 size_t nrX = std::distance(x_p,endX);
296 size_t nrY = std::distance(y_p,endY);
299 for(
unsigned int i = 0 ; i < nrX; ++i, ++x_p, ++out_p)
301 xd = * x_p - start_index;
302 double indexd = floor(xd);
303 int index =
static_cast<int>( indexd );
304 assert(fabs (index - indexd) < 0.001);
310 }
else if( index < static_cast<int>(nrY-1) )
312 TReal mu = xd - indexd;
313 YInputIterator y1_p = (y_p + index);
315 TReal y2 = *(++y1_p);
316 *out_p = interpolator(y1,y2,mu);
320 *out_p = *(y_p + (nrY-1));
CubicInterpolate(TReal epsilon=std::numeric_limits< TReal >::epsilon())
CosineInterpolate Functor Linear interpolation results in discontinuities at each point...
HermiteInterpolate(TReal tension, TReal bias, TReal epsilon=std::numeric_limits< TReal >::epsilon())
static void interpolateCubicHermite(YInputIterator begY, YInputIterator endY, XInputIterator begX, XInputIterator endX, OutputIterator out, TFunctor &functor, int start_index=0)
Cubic or Hermite interpolation worker.
CubicInterpolate Functor.
LinearInterpolate Functor.
LinearInterpolate(TReal epsilon=std::numeric_limits< TReal >::epsilon())
const double PI(3.14159265358979323846264338327950288)
the ratio of the circumference of a circle to its diameter;
EQUISPACEINTERPOL Interpolation on a equidistantly spaced grid.
TReal operator()(TReal y1, TReal y2, TReal mu)
operator
static void interpolateLinearCosine(YInputIterator y_p, YInputIterator endY, XInputIterator x_p, XInputIterator endX, OutputIterator out_p, TFunctor &interpolator, int start_index=0)
Linear cubic interpolator worker.
TReal operator()(TReal y1, TReal y2, TReal mu)