ProteoWizard
interpolation.hpp
Go to the documentation of this file.
1 //
2 // $Id: interpolation.hpp 5313 2013-12-17 18:06:54Z chambm $
3 //
4 //
5 // Original author: Witold Wolski <wewolski@gmail.com>
6 //
7 // Copyright : ETH Zurich
8 //
9 // Licensed under the Apache License, Version 2.0 (the "License");
10 // you may not use this file except in compliance with the License.
11 // You may obtain a copy of the License at
12 //
13 // http://www.apache.org/licenses/LICENSE-2.0
14 //
15 // Unless required by applicable law or agreed to in writing, software
16 // distributed under the License is distributed on an "AS IS" BASIS,
17 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18 // See the License for the specific language governing permissions and
19 // limitations under the License.
20 //
21 
22 
23 #ifndef INTERPOLATIONUTILITIES_H
24 #define INTERPOLATIONUTILITIES_H
25 
27 #include <cmath>
28 
29 namespace ralab
30 {
31  namespace base
32  {
33  namespace base
34  {
35 
36  namespace utilities
37  {
38  /// LinearInterpolate Functor
39  template<typename TReal>
41  {
42  TReal epsilon_;
44  :epsilon_(epsilon){}
45 
46  inline TReal operator()
47  (
48  TReal y1, //!< y1
49  TReal y2, //!< y2
50  TReal mu //!< location parameter 0,1
51  )
52  {
53  if(mu < epsilon_)
54  return y1;
55  else if(-(mu - 1.) < epsilon_)
56  return y2;
57  else
58  return ( y1 * (1-mu) + y2 * mu ) ;
59  }
60  };//
61 
62  /// CosineInterpolate Functor
63  /// Linear interpolation results in discontinuities at each point.
64  /// Often a smoother interpolating function is desirable, perhaps the simplest is cosine interpolation.
65  /// A suitable orientated piece of a cosine function serves
66  /// to provide a smooth transition between adjacent segments.
67 
68  template<typename TReal>
70  {
71  /*!\brief operator */
72  inline TReal operator()(
73  TReal y1,//!< y1
74  TReal y2,//!< y2
75  TReal mu//!< location parameter in [0.,1.]
76  )
77  {
78  TReal mu2;
79  mu2 = (1. - cos(mu* ralab::constants::PI))/2;
80  return(y1*( 1. -mu2 )+y2*mu2);
81  }
82  };
83 
84 
85  /// CubicInterpolate Functor
86 
87  /// Cubic interpolation is the simplest method that offers true continuity between the segments.
88  /// As such it requires more than just the two endpoints of the segment but also the two points on either side of them.
89  /// So the function requires 4 points in all labeled y0, y1, y2, and y3, in the code below.
90  /// mu still behaves the same way for interpolating between the segment y1 to y2.
91  /// This doe s raise issues for how to interpolate between the first and last segments.
92  /// In the examples here I just haven't bothered.
93  /// A common solution is the dream up two extra points at the start and end of the sequence,
94  /// the new points are created so that they have a slope equal to the slope of the start or end segment.
95 
96  template<typename TReal>
98  {
99  TReal epsilon_;
101  ):epsilon_(epsilon)
102  {}
103 
104  /*!\brief operator */
105  inline TReal operator()
106  (
107  TReal y0,//!< y0
108  TReal y1,//!< y1
109  TReal y2,//!< y2
110  TReal y3,//!< y3
111  double mu//!< location parameter in [0.,1.]
112  )
113  {
114  if(mu < epsilon_)
115  {
116  return y1;
117  }
118  else if(-(mu - 1.) < epsilon_)
119  {
120  return y2;
121  }
122  else
123  {
124  TReal a0,a1,a2,a3,mu2;
125  mu2 = mu*mu;
126  a0 = y3 - y2 - y0 + y1;
127  a1 = y0 - y1 - a0;
128  a2 = y2 - y0;
129  a3 = y1;
130  return(a0*mu*mu2 + a1*mu2 + a2*mu + a3);
131  }
132  }
133  };
134 
135  /// HermiteInterpolation.
136  /// Hermite interpolation like cubic requires 4 points so that it can achieve a higher degree of continuity.
137  /// In addition it has nice tension and biasing controls.
138  /// Tension can be used to tighten up the curvature at the known points.
139  /// The bias is used to twist the curve about the known points.
140  /// The examples shown here have the default tension and bias values of 0,
141  /// it will be left as an exercise for the reader to explore different tension and bias values.
142 
143  template<typename TReal>
145  {
146  TReal tension_;
147  TReal bias_;
148  TReal epsilon_;
150  TReal tension,//!< 1 is high, 0 normal, -1 is low
151  TReal bias,//!< 0 is even, positive is towards first segment, negative towards the other
153  ): tension_(tension), bias_(bias), epsilon_(epsilon)
154  {}
155  /*!\brief operator */
156  inline TReal operator ()(
157  TReal y0,//!< y0
158  TReal y1,//!< y1
159  TReal y2,//!< y2
160  TReal y3,//!< y3
161  TReal mu //!< location
162  )
163  {
164  if(mu < epsilon_)
165  {
166  return y1;
167  }
168  else if(-(mu - 1.) < epsilon_)
169  {
170  return y2;
171  }
172  else
173  {
174  TReal m0,m1,mu2,mu3;
175  TReal a0,a1,a2,a3;
176  mu2 = mu * mu;
177  mu3 = mu2 * mu;
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;
184  a2 = mu3 - mu2;
185  a3 = -2*mu3 + 3*mu2;
186  return( a0*y1 + a1*m0 + a2*m1 + a3*y2 );
187  }
188  }
189  }; // HermiteInterpolation
190 
191 
192  /// Cubic or Hermite interpolation worker
193  template<
194  typename YInputIterator,
195  typename XInputIterator,
196  typename OutputIterator,
197  typename TFunctor
198  >
199  static void interpolateCubicHermite
200  (
201  YInputIterator begY,
202  YInputIterator endY,
203  XInputIterator begX,
204  XInputIterator endX,
205  OutputIterator out, //!< interpolated values, same length as x.
206  TFunctor & functor, //!< either CubicInterpolate or HermiteInterpolate
207  int start_index = 0
208  )
209  {
210  typedef typename std::iterator_traits<OutputIterator>::value_type TReal;
211  //size_t nrX = std::distance( begX , endX );
212  size_t nrY = std::distance( begY , endY );
213  OutputIterator outI = out;
214 
215  for( unsigned int i = 0 ; begX != endX ; ++i , ++begX, ++outI )
216  {
217  double xd = *begX - start_index;
218  int index = static_cast<int>(floor(xd));
219  //interpolate
220  TReal mu = xd - static_cast<double>(index);
221  if(index < -1)
222  {
223  *outI = *begY;
224  }
225  else if(index == -1)
226  {
227  TReal y1 = *begY;
228  TReal y2 = *begY;
229  TReal y3 = *begY;
230  TReal y4 = *(begY+1);
231  *outI = functor(y1 , y2 , y3 , y4 , mu );
232  }
233  //extrapolate
234  else if(index == 0 )
235  {
236  TReal y1 = 0;
237  TReal y2 = *begY;
238  TReal y3 = *(begY+1);
239  TReal y4 = *(begY+2);
240  *outI = functor(y1 , y2 , y3 , y4 , mu );
241  }
242  else if( index > 0 && index < static_cast<int>(nrY - 2) )//the normal case
243  {
244  YInputIterator begTmp = (begY + index - 1);
245  TReal y1 = *begTmp;
246  TReal y2 = *(begTmp + 1);
247  TReal y3 = *(begTmp + 2);
248  TReal y4 = *(begTmp + 3);
249  *outI = functor(y1 , y2 , y3 , y4 , mu );
250  }
251  else if(index == static_cast<int>(nrY-2) ) //you are getting out of range
252  {
253  YInputIterator begTmp = (begY + index - 1);
254  TReal y1 = *begTmp;
255  TReal y2 = *(begTmp+1);
256  TReal y3 = *(begTmp+2);
257  TReal y4 = 0 ;
258  *outI = functor(y1 , y2 , y3 , y4 ,mu);
259  }
260  else if(index == static_cast<int>(nrY-1) ) //you are even farther out...
261  {
262  YInputIterator begTmp = (begY + index - 1);
263  TReal y1 = *begTmp;
264  TReal y2 = *(begTmp+1);
265  TReal y3 = *(begTmp+1);
266  TReal y4 = *(begTmp+1);
267  *outI = functor(y1 , y2 , y3 , y4 , mu );
268  }
269  else
270  {
271  *outI = *(endY-1);
272  }
273  }//end for
274  }//end interpolate_cubic
275 
276  /// Linear cubic interpolator worker
277  template <
278  typename YInputIterator,
279  typename XInputIterator,
280  typename OutputIterator,
281  typename TFunctor
282  >
283  static void interpolateLinearCosine
284  (
285  YInputIterator y_p, //!< y values equidistantly spaced. spacing is [0,1,2, .... ,len(y)]
286  YInputIterator endY,
287  XInputIterator x_p, //!< points to interpolate at
288  XInputIterator endX,
289  OutputIterator out_p, //!< interpolated values, same length as x.
290  TFunctor & interpolator, //!< interpolation functor, either: CosineInterpolate, LinearInterpolate.
291  int start_index = 0 //!< if y values are placed on a grid with start_index != 0
292  )
293  {
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);
297  TReal xd;
298 
299  for(unsigned int i = 0 ; i < nrX; ++i, ++x_p, ++out_p)
300  {
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);
305 
306  //interpolate
307  if(index < 0 )
308  {
309  *out_p = *y_p;
310  }else if( index < static_cast<int>(nrY-1) )
311  {
312  TReal mu = xd - indexd;
313  YInputIterator y1_p = (y_p + index);
314  TReal y1 = *y1_p;
315  TReal y2 = *(++y1_p);
316  *out_p = interpolator(y1,y2,mu);
317  }
318  else
319  {
320  *out_p = *(y_p + (nrY-1));
321  }
322  }//end for
323  }//end interpolate cubic
324  }//end utilities
325  }//end base
326  }//base
327 }//ralab
328 
329 #endif
CubicInterpolate(TReal epsilon=std::numeric_limits< TReal >::epsilon())
CosineInterpolate Functor Linear interpolation results in discontinuities at each point...
const double epsilon
Definition: DiffTest.cpp:41
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.
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.
Definition: base.hpp:39
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)