BPFTmplDef.hxx

Go to the documentation of this file.
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 
Generated by  doxygen 1.6.3