32 template <
typename TX,
typename TY>
const TData BPFTmpl<TX,TY>::Infinity =
TData(1e30);
37 template <
class TX,
class TY>
48 mIsSplineUpdated(false),
49 mLeftDerivative((
TData)Infinity),
50 mRightDerivative((
TData)Infinity)
62 template <
class TX,
class TY>
72 mIsSplineUpdated(false),
73 mLeftDerivative(Infinity),
74 mRightDerivative(Infinity)
88 template <
class TX,
class TY>
98 mIsSplineUpdated(false),
99 mLeftDerivative(Infinity),
100 mRightDerivative(Infinity)
114 template <
class TX,
class TY>
124 mIsSplineUpdated(false),
125 mLeftDerivative(Infinity),
126 mRightDerivative(Infinity)
137 template <
class TX,
class TY>
139 meInterpolation(orig.meInterpolation),
141 mClosestPoints(orig.mClosestPoints.AllocatedSize()),
142 mc(orig.mc.AllocatedSize()),
143 md(orig.md.AllocatedSize()),
147 mSplineTable(orig.mSplineTable),
148 mIsSplineUpdated(orig.mIsSplineUpdated),
149 mLeftDerivative(orig.mLeftDerivative),
150 mRightDerivative(orig.mRightDerivative)
164 template <
class TX,
class TY>
167 mArray.SetSize(AllocatedSize());
168 for(
int i=0;i<AllocatedSize();i++)
170 mArray[i].SetX((TX)i);
183 template <
class TX,
class TY>
186 mIsSplineUpdated =
false;
187 if (mArray.Size() == 0 || point.
GetX() > mArray[Size()-1].GetX())
189 mArray.AddElem(point);
193 TIndex closestIndex = mSearch.Find(point);
194 if (closestIndex == -1)
196 if (point.
GetX() >= GetXValue(Size() - 1))
198 mArray.AddElem(point);
203 mArray.InsertElem(0, point);
209 if (GetXValue(closestIndex) == point.
GetX())
211 SetValue(closestIndex, point.
GetY());
216 mArray.InsertElem(closestIndex+1 ,point);
217 return closestIndex+1;
230 template <
class TX,
class TY>
234 return Insert(tmpPoint);
241 template <
class TX,
class TY>
243 mArray.DeleteElem(index);
244 mIsSplineUpdated=
false;
251 template <
class TX,
class TY>
254 mArray.DeleteElem(GetPosition(x));
255 mIsSplineUpdated=
false;
263 template <
class TX,
class TY>
266 for (
int i=leftIndex; i<=rightIndex; i++)
268 DeleteIndex(leftIndex);
270 mIsSplineUpdated=
false;
277 template <
class TX,
class TY>
280 TIndex leftIndex=GetPosition(leftX);
281 TIndex rightIndex=GetPosition(rightX);
282 DeleteBetweenIndex(leftIndex,rightIndex);
283 mIsSplineUpdated=
false;
291 template <
class TX,
class TY>
294 meInterpolation=eInterpolation;
295 switch(eInterpolation)
325 CLAM_ASSERT(Size()<11,
"BPF::SetIntpType:Cannot ser more than 10th order interpolation");
335 CLAM_ASSERT(
false,
"Unsupported interpolation method" );
338 const unsigned newSize = mOrder+1;
353 template <
class TX,
class TY>
360 if(x<mLastX) mLastIndex=0;
367 return GetValueFromIndex(0);
369 return GetValueFromIndex(Size()-1);
373 if(GetXValue(mLastIndex)==x)
return GetValueFromIndex(mLastIndex);
374 switch(eInterpolation)
378 GetnClosest(mLastIndex);
379 if(GetXValue(mClosestPoints[0])<=x)
380 return GetValueFromIndex(mClosestPoints[0]);
381 return GetValueFromIndex(mClosestPoints[0]+1);
385 GetnClosest(mLastIndex);
386 return GetValueFromIndex(mClosestPoints[0]);
391 if(GetXValue(mLastIndex)<=x)
393 mClosestPoints[0]=mLastIndex;
394 mClosestPoints[1]=mLastIndex+1;
398 mClosestPoints[0]=mLastIndex-1;
399 mClosestPoints[1]=mLastIndex;
401 return BPFPolInt(x,mClosestPoints,error);
405 CLAM_ASSERT(mIsSplineUpdated,
"BPF::Spline table not updated");
406 return BPFSplineInt(x);
414 GetnClosest(mLastIndex);
415 return BPFPolInt(x,mClosestPoints,error);
422 for(
TIndex i=0; i<mArray.Size(); i++)
426 return BPFPolInt(x,indexArray,error);
430 CLAM_ASSERT(
false,
"Invalid BPF interpolation method.");
445 template <
class TX,
class TY>
448 CLAM_ASSERT(Size() >= mOrder+1,
"BPF size must be at least the interpolation order plus one");
449 int first = foundIndex-mOrder/2;
451 if (first<0) first=0;
453 int safelimit = Size()-(mOrder+1);
454 if (first>safelimit) first = safelimit;
456 for(
int i=0; i<mOrder+1; i++) {
457 mClosestPoints[i]=first+i;
512 template <
class TX,
class TY>
516 return mSearch.Find(tmpPoint);
522 template <
class TX,
class TY>
526 mArray = originalBPF.
mArray;
531 mOrder=originalBPF.
mOrder;
538 template <
class TX,
class TY>
541 if(!mIsSplineUpdated)
543 mIsSplineUpdated=
true;
549 template <
class TX,
class TY>
552 mLeftDerivative = val;
553 mIsSplineUpdated=
false;
555 template <
class TX,
class TY>
558 if (mLeftDerivative < Infinity) {
559 mLeftDerivative = Infinity;
560 mIsSplineUpdated=
false;
566 template <
class TX,
class TY>
569 mRightDerivative = val;
570 mIsSplineUpdated=
false;
572 template <
class TX,
class TY>
575 if (mLeftDerivative < Infinity) {
576 mRightDerivative = Infinity;
577 mIsSplineUpdated=
false;
585 template <
class TX,
class TY>
591 TX dif=
Abs(x-GetXValue(closestPointsIndex[0]));
592 for(
int i=0; i<=mOrder; i++)
594 TX dift=
Abs(x-GetXValue(closestPointsIndex[i]));
600 md[i]=mc[i]=GetValueFromIndex(closestPointsIndex[i]);
603 TY y=GetValueFromIndex(closestPointsIndex[iClosest]);
604 for(
int m=0; m<mOrder; m++)
606 for(
int i=0; i<mOrder-m; i++)
608 TX ho=GetXValue(closestPointsIndex[i])-x;
609 TX hp=GetXValue(closestPointsIndex[i+m+1])-x;
612 CLAM_ASSERT(den!=0,
"Division by zero error interpolating BPF");
617 if ( 2*iClosest < mOrder-m )
619 errorEstimate=mc[iClosest];
623 errorEstimate=md[iClosest-1];
635 template <
class TX,
class TY>
640 CLAM_ASSERT(n > 2,
"BPF size too small for spline");
644 mSplineTable.Resize(n);
645 mSplineTable.SetSize(n);
647 if (mLeftDerivative >= Infinity)
648 mSplineTable[0]=u[0]=0.0;
650 mSplineTable[0] = -0.5;
651 u[0]= (
TData(3.0) / (GetXValue(1)-GetXValue(0)) ) *
652 ( (GetValueFromIndex(1) - GetValueFromIndex(0)) /
653 (GetXValue(1) - GetXValue(0))
657 for(
int i=2;i<=n-1;i++)
659 TY sig=(GetXValue(i-1)-GetXValue(i-2))/(GetXValue(i)-GetXValue(i-2));
660 TY p=sig*mSplineTable[i-2]+2;
661 mSplineTable[i-1]=(sig-1)/p;
662 u[i-1]=(GetValueFromIndex(i)-GetValueFromIndex(i-1))/(GetXValue(i)-GetXValue(i-1))-
663 (GetValueFromIndex(i-1)-GetValueFromIndex(i-2))/(GetXValue(i-1)-GetXValue(i-2));
664 u[i-1]=(6*u[i-1]/(GetXValue(i)-GetXValue(i-2))-sig*u[i-2])/p;
669 if (mRightDerivative < Infinity) {
671 un = (
TData(3.0)/(GetXValue(n-1)-GetXValue(n-2))) *
673 ( GetValueFromIndex(n-1) - GetValueFromIndex(n-2)) /
674 ( GetXValue(n-1) - GetXValue(n-2)));
676 mSplineTable[n-1]=((un-qn*u[n-2])/(qn*mSplineTable[n-2]+1));
677 for(
int k=n-1; k>=1; k--)
679 mSplineTable[k-1]=mSplineTable[k-1]*mSplineTable[k]+u[k-1];
686 template <
class TX,
class TY>
690 int khi=mArray.Size();
694 if(GetXValue(k-1)>x) khi=k;
697 TX h=GetXValue(khi-1)-GetXValue(klo-1);
699 TX a=(GetXValue(khi-1)-x)/h;
700 TX b=(x-GetXValue(klo-1))/h;
701 return (a*GetValueFromIndex(klo-1)+b*GetValueFromIndex(khi-1)+((a*a*a-a)*
702 mSplineTable[klo-1]+(b*b*b-b)*mSplineTable[khi-1])*(h*h)/6);
705 template <
class TX,
class TY>
709 storage.
Store(adapterInt);
711 storage.
Store(adapter);
713 template <
class TX,
class TY>
717 storage.
Load(adapterInt);
719 storage.
Load(adapter);
725 #endif // _BPFTmplDef_