CLAM-Development
1.1
|
00001 /* 00002 * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG) 00003 * UNIVERSITAT POMPEU FABRA 00004 * 00005 * 00006 * This program is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program; if not, write to the Free Software 00018 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00019 * 00020 */ 00021 00022 00023 #ifndef _BPFTmplDef_ 00024 #define _BPFTmplDef_ 00025 00026 #include "GlobalEnums.hxx" 00027 00028 namespace CLAM 00029 { 00030 00031 00032 template <typename TX, typename TY> const TData BPFTmpl<TX,TY>::Infinity = TData(1e30); 00033 00037 template <class TX,class TY> 00038 BPFTmpl<TX,TY>::BPFTmpl() : 00039 mArray(0), 00040 mSearch(mArray), 00041 mClosestPoints(10), 00042 mc(2), 00043 md(2), 00044 mOrder(1), 00045 mLastIndex(-1), 00046 mLastX(0), 00047 mSplineTable(0), 00048 mIsSplineUpdated(false), 00049 mLeftDerivative((TData)Infinity), 00050 mRightDerivative((TData)Infinity) 00051 { 00052 mClosestPoints.SetSize(mClosestPoints.AllocatedSize()); 00053 mc.SetSize(2); 00054 md.SetSize(2); 00055 meInterpolation=EInterpolation::eLinear; 00056 } 00057 00062 template <class TX,class TY> 00063 BPFTmpl<TX,TY>::BPFTmpl(const EInterpolation& eInterpolation) : 00064 mArray(0), 00065 mSearch(mArray), 00066 mClosestPoints(10), 00067 mc(0), 00068 md(0), 00069 mLastIndex(-1), 00070 mLastX(0), 00071 mSplineTable(0), 00072 mIsSplineUpdated(false), 00073 mLeftDerivative(Infinity), 00074 mRightDerivative(Infinity) 00075 { 00076 mClosestPoints.SetSize(mClosestPoints.AllocatedSize()); 00077 SetIntpType(eInterpolation); 00078 } 00079 00080 00088 template <class TX,class TY> 00089 BPFTmpl<TX,TY>::BPFTmpl(TSize size) : 00090 mArray(size), 00091 mSearch(mArray), 00092 mClosestPoints(10), 00093 mc(0), 00094 md(0), 00095 mLastIndex(-1), 00096 mLastX(0), 00097 mSplineTable(0), 00098 mIsSplineUpdated(false), 00099 mLeftDerivative(Infinity), 00100 mRightDerivative(Infinity) 00101 { 00102 mClosestPoints.SetSize(mClosestPoints.AllocatedSize()); 00103 SetIntpType(EInterpolation::eLinear); 00104 SetSize(size); 00105 } 00106 00114 template <class TX,class TY> 00115 BPFTmpl<TX,TY>::BPFTmpl(TSize size,const EInterpolation& eInterpolation) : 00116 mArray(size), 00117 mSearch(mArray), 00118 mClosestPoints(10), 00119 mc(0), 00120 md(0), 00121 mLastIndex(-1), 00122 mLastX(0), 00123 mSplineTable(0), 00124 mIsSplineUpdated(false), 00125 mLeftDerivative(Infinity), 00126 mRightDerivative(Infinity) 00127 { 00128 mClosestPoints.SetSize(mClosestPoints.AllocatedSize()); 00129 SetIntpType(eInterpolation); 00130 SetSize(size); 00131 } 00132 00137 template <class TX,class TY> 00138 BPFTmpl<TX,TY>::BPFTmpl(const BPFTmpl<TX,TY>& orig) : 00139 meInterpolation(orig.meInterpolation), 00140 mArray(orig.mArray), 00141 mClosestPoints(orig.mClosestPoints.AllocatedSize()), 00142 mc(orig.mc.AllocatedSize()), 00143 md(orig.md.AllocatedSize()), 00144 mOrder(orig.mOrder), 00145 mLastIndex(-1), 00146 mLastX(0), 00147 mSplineTable(orig.mSplineTable), 00148 mIsSplineUpdated(orig.mIsSplineUpdated), 00149 mLeftDerivative(orig.mLeftDerivative), 00150 mRightDerivative(orig.mRightDerivative) 00151 { 00152 mClosestPoints.SetSize(orig.mClosestPoints.Size()); 00153 mc.SetSize(orig.mc.Size()); 00154 md.SetSize(orig.md.Size()); 00155 mSearch.Set(mArray); 00156 mnPoints=orig.mnPoints; 00157 00158 } 00159 00164 template <class TX,class TY> 00165 void BPFTmpl<TX,TY>::Init() 00166 { 00167 mArray.SetSize(AllocatedSize()); 00168 for(int i=0;i<AllocatedSize();i++) 00169 { 00170 mArray[i].SetX((TX)i); 00171 SetValue(i,0); 00172 } 00173 } 00174 00175 00183 template <class TX,class TY> 00184 TIndex BPFTmpl<TX,TY>::Insert(const PointTmpl<TX,TY> &point) 00185 { 00186 mIsSplineUpdated = false; 00187 if (mArray.Size() == 0 || point.GetX() > mArray[Size()-1].GetX()) 00188 { 00189 mArray.AddElem(point); 00190 return Size()-1; 00191 } 00192 00193 TIndex closestIndex = mSearch.Find(point); 00194 if (closestIndex == -1) 00195 { 00196 if (point.GetX() >= GetXValue(Size() - 1)) 00197 { 00198 mArray.AddElem(point); 00199 return Size()-1; 00200 } 00201 else 00202 { 00203 mArray.InsertElem(0, point); 00204 return 0; 00205 } 00206 } 00207 else 00208 { 00209 if (GetXValue(closestIndex) == point.GetX()) 00210 { 00211 SetValue(closestIndex, point.GetY()); 00212 return closestIndex; 00213 } 00214 else 00215 { 00216 mArray.InsertElem(closestIndex+1 ,point); 00217 return closestIndex+1; 00218 } 00219 } 00220 } 00221 00230 template <class TX,class TY> 00231 TIndex BPFTmpl<TX,TY>::Insert(const TX &x,const TX &y) 00232 { 00233 PointTmpl<TX,TY> tmpPoint(x,y); 00234 return Insert(tmpPoint); 00235 } 00236 00241 template <class TX,class TY> 00242 void BPFTmpl<TX,TY>::DeleteIndex(TIndex index){ 00243 mArray.DeleteElem(index); 00244 mIsSplineUpdated=false; 00245 } 00251 template <class TX,class TY> 00252 void BPFTmpl<TX,TY>::DeleteThroughXValue(const TX& x) 00253 { 00254 mArray.DeleteElem(GetPosition(x)); 00255 mIsSplineUpdated=false; 00256 } 00257 00263 template <class TX,class TY> 00264 void BPFTmpl<TX,TY>::DeleteBetweenIndex(TIndex leftIndex,TIndex rightIndex) 00265 { 00266 for (int i=leftIndex; i<=rightIndex; i++) 00267 { 00268 DeleteIndex(leftIndex); 00269 } 00270 mIsSplineUpdated=false; 00271 } 00277 template <class TX,class TY> 00278 void BPFTmpl<TX,TY>::DeleteBetweenXValues(const TX& leftX,const TX& rightX) 00279 { 00280 TIndex leftIndex=GetPosition(leftX); 00281 TIndex rightIndex=GetPosition(rightX); 00282 DeleteBetweenIndex(leftIndex,rightIndex); 00283 mIsSplineUpdated=false; 00284 } 00285 00291 template <class TX,class TY> 00292 void BPFTmpl<TX,TY>::SetIntpType(const EInterpolation& eInterpolation) 00293 { 00294 meInterpolation=eInterpolation; 00295 switch(eInterpolation) 00296 { 00297 case(EInterpolation::eLinear)://linear interpolation between two closest points 00298 { 00299 mOrder=1; 00300 break; 00301 } 00302 case(EInterpolation::ePolynomial2)://parabolic interpolation 00303 { 00304 mOrder=2; 00305 break; 00306 } 00307 case(EInterpolation::ePolynomial3)://3rd order polynomial interpolation 00308 { 00309 mOrder=3; 00310 break; 00311 } 00312 case(EInterpolation::ePolynomial4)://4th order polynomial interpolation 00313 { 00314 mOrder=4; 00315 break; 00316 } 00317 case(EInterpolation::ePolynomial5)://5th order polynomial interpolation 00318 { 00319 mOrder=5; 00320 break; 00321 } 00322 case(EInterpolation::ePolynomialn):/*nth order polynomial interpolation where n is number 00323 of points in the BPF-1. Must be less than 10*/ 00324 { 00325 CLAM_ASSERT(Size()<11,"BPF::SetIntpType:Cannot ser more than 10th order interpolation"); 00326 mOrder=Size()-1; 00327 break; 00328 } 00329 case(EInterpolation::eSpline ): // MRJ: Nice refactoring, but forgot about the Spline... 00330 { 00331 return; 00332 } 00333 default: 00334 { 00335 CLAM_ASSERT( false, "Unsupported interpolation method" ); 00336 } 00337 } 00338 const unsigned newSize = mOrder+1; 00339 mc.Resize(newSize); 00340 mc.SetSize(newSize); 00341 md.Resize(newSize); 00342 md.SetSize(newSize); 00343 } 00344 00345 00353 template <class TX,class TY> 00354 TY BPFTmpl<TX,TY>::GetValue(const TX& x,const EInterpolation& eInterpolation) const /*Gets value 00355 according to the interpolation type set*/ 00356 { 00357 // First we should check if the the x value belongs to the BPF itself 00358 // and there is no need to interpolate 00359 00360 if(x<mLastX) mLastIndex=0; 00361 00362 TIndex i=mSearch.Find( PointTmpl<TX,TY>(x,0.0), mLastIndex); 00363 if(i==-1) 00364 { 00365 // Outside of the BPF; get the first or last value 00366 if (GetXValue(0)>x) 00367 return GetValueFromIndex(0); 00368 else 00369 return GetValueFromIndex(Size()-1); 00370 } 00371 mLastIndex=i; 00372 mLastX=x; 00373 if(GetXValue(mLastIndex)==x) return GetValueFromIndex(mLastIndex); 00374 switch(eInterpolation) 00375 { 00376 case(EInterpolation::eStep)://returns previous point value 00377 { 00378 GetnClosest(mLastIndex); 00379 if(GetXValue(mClosestPoints[0])<=x) 00380 return GetValueFromIndex(mClosestPoints[0]); 00381 return GetValueFromIndex(mClosestPoints[0]+1); 00382 } 00383 case(EInterpolation::eRound)://return closest point value 00384 { 00385 GetnClosest(mLastIndex); 00386 return GetValueFromIndex(mClosestPoints[0]); 00387 } 00388 case(EInterpolation::eLinear)://linear interpolation between two closest points 00389 { 00390 TData error=0; 00391 if(GetXValue(mLastIndex)<=x) 00392 { 00393 mClosestPoints[0]=mLastIndex; 00394 mClosestPoints[1]=mLastIndex+1; 00395 } 00396 else 00397 { 00398 mClosestPoints[0]=mLastIndex-1; 00399 mClosestPoints[1]=mLastIndex; 00400 } 00401 return BPFPolInt(x,mClosestPoints,error); 00402 } 00403 case(EInterpolation::eSpline)://3rd order spline interpolation 00404 { 00405 CLAM_ASSERT(mIsSplineUpdated,"BPF::Spline table not updated"); 00406 return BPFSplineInt(x);//get actual value 00407 } 00408 case(EInterpolation::ePolynomial2)://parabolic interpolation 00409 case(EInterpolation::ePolynomial3)://3rd order polynomial interpolation 00410 case(EInterpolation::ePolynomial4)://4th order polynomial interpolation 00411 case(EInterpolation::ePolynomial5)://5th order polynomial interpolation 00412 { 00413 TData error=0; 00414 GetnClosest(mLastIndex); 00415 return BPFPolInt(x,mClosestPoints,error); 00416 } 00417 case(EInterpolation::ePolynomialn):/*nth order polynomial interpolation where n is number 00418 of points in the BPF-1*/ 00419 { 00420 Array<TIndex> indexArray(mArray.Size()); 00421 TData error=0; 00422 for(TIndex i=0; i<mArray.Size(); i++) 00423 { 00424 indexArray[i]=i; 00425 } 00426 return BPFPolInt(x,indexArray,error); 00427 } 00428 default: 00429 { 00430 CLAM_ASSERT(false, "Invalid BPF interpolation method."); 00431 return 0; 00432 } 00433 } 00434 } 00435 00445 template <class TX,class TY> 00446 void BPFTmpl<TX,TY>::GetnClosest(TIndex foundIndex) const 00447 { 00448 CLAM_ASSERT(Size() >= mOrder+1, "BPF size must be at least the interpolation order plus one"); 00449 int first = foundIndex-mOrder/2; 00450 // Fix the first index when it goes too 00451 if (first<0) first=0; 00452 // At least mOrder+1 fordward 00453 int safelimit = Size()-(mOrder+1); 00454 if (first>safelimit) first = safelimit; 00455 00456 for(int i=0; i<mOrder+1; i++) { 00457 mClosestPoints[i]=first+i; 00458 } 00459 00460 /* 00461 TIndex oP=GetPosition(x); 00462 int N = mArray.Size(); 00463 TIndex i, j; 00464 00465 TData elemToTake = (n-1.0)/2.0; 00466 00467 int toTakeToTheLeft = (int)elemToTake; 00468 int toTakeToTheRight = (int)ceil(elemToTake); 00469 00470 int dL = (oP) - toTakeToTheLeft; 00471 int dR = (N-oP) - toTakeToTheRight; 00472 00473 if (dL>= 0 && dR >= 0) 00474 { 00475 i = oP - toTakeToTheLeft; 00476 j = oP + toTakeToTheRight; 00477 00478 if (j == N) 00479 { 00480 j = N-1; 00481 i--; 00482 } 00483 } 00484 else if (dL<0 && dR >= 0) 00485 { 00486 i = 0; 00487 j = oP + toTakeToTheRight + (-1)*dL; 00488 } 00489 else if (dL>=0 && dR < 0) 00490 { 00491 j = N-1; 00492 i = oP - toTakeToTheLeft - (-1)*dR - 1 ; 00493 00494 } 00495 else if (dL<0 && dR < 0) 00496 { 00497 throw Err("Not enough Points to interpolate"); 00498 } 00499 00500 for (int k=i; k<=j; k++) 00501 { 00502 mClosestPoints[k-i] = k; 00503 }*/ 00504 00505 } 00506 00512 template <class TX,class TY> 00513 TIndex BPFTmpl<TX,TY>::GetPosition(const TX& x) const 00514 { 00515 PointTmpl<TX,TY> tmpPoint(x,0); 00516 return mSearch.Find(tmpPoint); 00517 } 00518 00522 template <class TX,class TY> 00523 BPFTmpl<TX,TY>& BPFTmpl<TX,TY>::operator=(const BPFTmpl<TX,TY> &originalBPF) 00524 { 00525 meInterpolation = originalBPF.meInterpolation; 00526 mArray = originalBPF.mArray; 00527 mSearch.Set(mArray); 00528 mSplineTable = originalBPF.mSplineTable; 00529 mIsSplineUpdated = originalBPF.mIsSplineUpdated; 00530 mnPoints=originalBPF.mnPoints; 00531 mOrder=originalBPF.mOrder; 00532 return *this; 00533 } 00534 00538 template <class TX,class TY> 00539 void BPFTmpl<TX,TY>::UpdateSplineTable() 00540 { 00541 if(!mIsSplineUpdated) 00542 CreateSplineTable(); // Create spline table if not updated 00543 mIsSplineUpdated=true; 00544 } 00545 00549 template <class TX,class TY> 00550 void BPFTmpl<TX,TY>::SetLeftDerivative(TData val) 00551 { 00552 mLeftDerivative = val; 00553 mIsSplineUpdated=false; 00554 } 00555 template <class TX,class TY> 00556 void BPFTmpl<TX,TY>::UnsetLeftDerivative(void) 00557 { 00558 if (mLeftDerivative < Infinity) { 00559 mLeftDerivative = Infinity; 00560 mIsSplineUpdated=false; 00561 } 00562 } 00566 template <class TX,class TY> 00567 void BPFTmpl<TX,TY>::SetRightDerivative(TData val) 00568 { 00569 mRightDerivative = val; 00570 mIsSplineUpdated=false; 00571 } 00572 template <class TX,class TY> 00573 void BPFTmpl<TX,TY>::UnsetRightDerivative(void) 00574 { 00575 if (mLeftDerivative < Infinity) { 00576 mRightDerivative = Infinity; 00577 mIsSplineUpdated=false; 00578 } 00579 } 00580 00585 template <class TX,class TY> 00586 TY BPFTmpl<TX,TY>::BPFPolInt 00587 (const TX& x,const Array<TIndex>& closestPointsIndex, TData &errorEstimate) 00588 const 00589 { 00590 int iClosest=0; 00591 TX dif=Abs(x-GetXValue(closestPointsIndex[0])); 00592 for(int i=0; i<=mOrder; i++) 00593 { 00594 TX dift=Abs(x-GetXValue(closestPointsIndex[i])); 00595 if(dift<dif) 00596 { 00597 iClosest=i; 00598 dif=dift; 00599 } 00600 md[i]=mc[i]=GetValueFromIndex(closestPointsIndex[i]); 00601 } 00602 00603 TY y=GetValueFromIndex(closestPointsIndex[iClosest]); 00604 for(int m=0; m<mOrder; m++) 00605 { 00606 for(int i=0; i<mOrder-m; i++) 00607 { 00608 TX ho=GetXValue(closestPointsIndex[i])-x; 00609 TX hp=GetXValue(closestPointsIndex[i+m+1])-x; 00610 TX w=mc[i+1]-md[i]; 00611 TX den=ho-hp; 00612 CLAM_ASSERT(den!=0, "Division by zero error interpolating BPF"); 00613 den=w/den; 00614 md[i]=hp*den; 00615 mc[i]=ho*den; 00616 } 00617 if ( 2*iClosest < mOrder-m ) 00618 { 00619 errorEstimate=mc[iClosest]; 00620 } 00621 else 00622 { 00623 errorEstimate=md[iClosest-1]; 00624 iClosest--; 00625 } 00626 y+=errorEstimate; 00627 } 00628 00629 return y; 00630 } 00635 template <class TX,class TY> 00636 void BPFTmpl<TX,TY>::CreateSplineTable() 00637 { 00638 int n=mArray.Size(); 00639 00640 CLAM_ASSERT(n > 2,"BPF size too small for spline"); 00641 00642 Array<TY> u(n); 00643 u.SetSize(n); 00644 mSplineTable.Resize(n); 00645 mSplineTable.SetSize(n); 00646 00647 if (mLeftDerivative >= Infinity) 00648 mSplineTable[0]=u[0]=0.0; // For a 'natural' spline 00649 else { 00650 mSplineTable[0] = -0.5; 00651 u[0]= (TData(3.0) / (GetXValue(1)-GetXValue(0)) ) * 00652 ( (GetValueFromIndex(1) - GetValueFromIndex(0)) / 00653 (GetXValue(1) - GetXValue(0)) 00654 - mLeftDerivative ); 00655 } 00656 00657 for(int i=2;i<=n-1;i++) 00658 { 00659 TY sig=(GetXValue(i-1)-GetXValue(i-2))/(GetXValue(i)-GetXValue(i-2)); 00660 TY p=sig*mSplineTable[i-2]+2; 00661 mSplineTable[i-1]=(sig-1)/p; 00662 u[i-1]=(GetValueFromIndex(i)-GetValueFromIndex(i-1))/(GetXValue(i)-GetXValue(i-1))- 00663 (GetValueFromIndex(i-1)-GetValueFromIndex(i-2))/(GetXValue(i-1)-GetXValue(i-2)); 00664 u[i-1]=(6*u[i-1]/(GetXValue(i)-GetXValue(i-2))-sig*u[i-2])/p; 00665 } 00666 00667 TY qn = 0.0; 00668 TY un = 0.0; 00669 if (mRightDerivative < Infinity) { 00670 qn = 0.5; 00671 un = (TData(3.0)/(GetXValue(n-1)-GetXValue(n-2))) * 00672 ( mRightDerivative - 00673 ( GetValueFromIndex(n-1) - GetValueFromIndex(n-2)) / 00674 ( GetXValue(n-1) - GetXValue(n-2))); 00675 } 00676 mSplineTable[n-1]=((un-qn*u[n-2])/(qn*mSplineTable[n-2]+1)); 00677 for(int k=n-1; k>=1; k--) 00678 { 00679 mSplineTable[k-1]=mSplineTable[k-1]*mSplineTable[k]+u[k-1]; 00680 } 00681 } 00686 template <class TX,class TY> 00687 TY BPFTmpl<TX,TY>::BPFSplineInt(const TX& x) const 00688 { 00689 int klo=1; 00690 int khi=mArray.Size(); 00691 while(khi-klo>1) 00692 { 00693 int k=(khi+klo)>>1; 00694 if(GetXValue(k-1)>x) khi=k; 00695 else klo=k; 00696 } 00697 TX h=GetXValue(khi-1)-GetXValue(klo-1); 00698 CLAM_ASSERT(h!=0.0, "Error interpolating Spline"); 00699 TX a=(GetXValue(khi-1)-x)/h; 00700 TX b=(x-GetXValue(klo-1))/h; 00701 return (a*GetValueFromIndex(klo-1)+b*GetValueFromIndex(khi-1)+((a*a*a-a)* 00702 mSplineTable[klo-1]+(b*b*b-b)*mSplineTable[khi-1])*(h*h)/6); 00703 } 00704 00705 template <class TX,class TY> 00706 void BPFTmpl<TX,TY>::StoreOn(Storage & storage) const 00707 { 00708 XMLComponentAdapter adapterInt(meInterpolation,"Interpolation",true); 00709 storage.Store(adapterInt); 00710 XMLComponentAdapter adapter(mArray, "Points", true); 00711 storage.Store(adapter); 00712 } 00713 template <class TX,class TY> 00714 void BPFTmpl<TX,TY>::LoadFrom(Storage & storage) 00715 { 00716 XMLComponentAdapter adapterInt(meInterpolation,"Interpolation",true); 00717 storage.Load(adapterInt); 00718 XMLComponentAdapter adapter(mArray, "Points", true); 00719 storage.Load(adapter); 00720 } 00721 00722 } // namespace CLAM 00723 00724 00725 #endif // _BPFTmplDef_ 00726