Cubic or Hermite interpolation worker.
209 {
210 typedef typename std::iterator_traits<OutputIterator>::value_type TReal;
211
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
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
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) )
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) )
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) )
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 }
274 }