ProteoWizard
interpolation.hpp
Go to the documentation of this file.
1//
2// $Id$
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
29namespace 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_;
43 LinearInterpolate(TReal epsilon = std::numeric_limits<TReal>::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_;
100 CubicInterpolate(TReal epsilon = std::numeric_limits<TReal>::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
152 TReal epsilon = std::numeric_limits<TReal>::epsilon()
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 >
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 >
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
const double epsilon
Definition DiffTest.cpp:41
static void interpolateCubicHermite(YInputIterator begY, YInputIterator endY, XInputIterator begX, XInputIterator endX, OutputIterator out, TFunctor &functor, int start_index=0)
Cubic or Hermite interpolation worker.
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.
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:40
CosineInterpolate Functor Linear interpolation results in discontinuities at each point.
TReal operator()(TReal y1, TReal y2, TReal mu)
operator
CubicInterpolate(TReal epsilon=std::numeric_limits< TReal >::epsilon())
HermiteInterpolate(TReal tension, TReal bias, TReal epsilon=std::numeric_limits< TReal >::epsilon())
TReal operator()(TReal y0, TReal y1, TReal y2, TReal y3, TReal mu)
operator
LinearInterpolate(TReal epsilon=std::numeric_limits< TReal >::epsilon())