CLAM-Development  1.4.0
BPFTmplDef.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG)
3  * UNIVERSITAT POMPEU FABRA
4  *
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 #ifndef _BPFTmplDef_
24 #define _BPFTmplDef_
25 
26 #include "GlobalEnums.hxx"
27 
28 namespace CLAM
29 {
30 
31 
32  template <typename TX, typename TY> const TData BPFTmpl<TX,TY>::Infinity = TData(1e30);
33 
37  template <class TX,class TY>
39  mArray(0),
40  mSearch(mArray),
41  mClosestPoints(10),
42  mc(2),
43  md(2),
44  mOrder(1),
45  mLastIndex(-1),
46  mLastX(0),
47  mSplineTable(0),
48  mIsSplineUpdated(false),
49  mLeftDerivative((TData)Infinity),
50  mRightDerivative((TData)Infinity)
51  {
53  mc.SetSize(2);
54  md.SetSize(2);
56  }
57 
62  template <class TX,class TY>
63  BPFTmpl<TX,TY>::BPFTmpl(const EInterpolation& eInterpolation) :
64  mArray(0),
65  mSearch(mArray),
66  mClosestPoints(10),
67  mc(0),
68  md(0),
69  mLastIndex(-1),
70  mLastX(0),
71  mSplineTable(0),
72  mIsSplineUpdated(false),
73  mLeftDerivative(Infinity),
74  mRightDerivative(Infinity)
75  {
77  SetIntpType(eInterpolation);
78  }
79 
80 
88  template <class TX,class TY>
90  mArray(size),
91  mSearch(mArray),
92  mClosestPoints(10),
93  mc(0),
94  md(0),
95  mLastIndex(-1),
96  mLastX(0),
97  mSplineTable(0),
98  mIsSplineUpdated(false),
99  mLeftDerivative(Infinity),
100  mRightDerivative(Infinity)
101  {
104  SetSize(size);
105  }
106 
114  template <class TX,class TY>
115  BPFTmpl<TX,TY>::BPFTmpl(TSize size,const EInterpolation& eInterpolation) :
116  mArray(size),
117  mSearch(mArray),
118  mClosestPoints(10),
119  mc(0),
120  md(0),
121  mLastIndex(-1),
122  mLastX(0),
123  mSplineTable(0),
124  mIsSplineUpdated(false),
125  mLeftDerivative(Infinity),
126  mRightDerivative(Infinity)
127  {
129  SetIntpType(eInterpolation);
130  SetSize(size);
131  }
132 
137  template <class TX,class TY>
139  meInterpolation(orig.meInterpolation),
140  mArray(orig.mArray),
141  mClosestPoints(orig.mClosestPoints.AllocatedSize()),
142  mc(orig.mc.AllocatedSize()),
143  md(orig.md.AllocatedSize()),
144  mOrder(orig.mOrder),
145  mLastIndex(-1),
146  mLastX(0),
147  mSplineTable(orig.mSplineTable),
148  mIsSplineUpdated(orig.mIsSplineUpdated),
149  mLeftDerivative(orig.mLeftDerivative),
150  mRightDerivative(orig.mRightDerivative)
151  {
153  mc.SetSize(orig.mc.Size());
154  md.SetSize(orig.md.Size());
155  mSearch.Set(mArray);
156  mnPoints=orig.mnPoints;
157 
158  }
159 
164  template <class TX,class TY>
166  {
167  mArray.SetSize(AllocatedSize());
168  for(int i=0;i<AllocatedSize();i++)
169  {
170  mArray[i].SetX((TX)i);
171  SetValue(i,0);
172  }
173  }
174 
175 
183  template <class TX,class TY>
185  {
186  mIsSplineUpdated = false;
187  if (mArray.Size() == 0 || point.GetX() > mArray[Size()-1].GetX())
188  {
189  mArray.AddElem(point);
190  return Size()-1;
191  }
192 
193  TIndex closestIndex = mSearch.Find(point);
194  if (closestIndex == -1)
195  {
196  if (point.GetX() >= GetXValue(Size() - 1))
197  {
198  mArray.AddElem(point);
199  return Size()-1;
200  }
201  else
202  {
203  mArray.InsertElem(0, point);
204  return 0;
205  }
206  }
207  else
208  {
209  if (GetXValue(closestIndex) == point.GetX())
210  {
211  SetValue(closestIndex, point.GetY());
212  return closestIndex;
213  }
214  else
215  {
216  mArray.InsertElem(closestIndex+1 ,point);
217  return closestIndex+1;
218  }
219  }
220  }
221 
230  template <class TX,class TY>
231  TIndex BPFTmpl<TX,TY>::Insert(const TX &x,const TX &y)
232  {
233  PointTmpl<TX,TY> tmpPoint(x,y);
234  return Insert(tmpPoint);
235  }
236 
241  template <class TX,class TY>
243  mArray.DeleteElem(index);
244  mIsSplineUpdated=false;
245  }
251  template <class TX,class TY>
253  {
254  mArray.DeleteElem(GetPosition(x));
255  mIsSplineUpdated=false;
256  }
257 
263  template <class TX,class TY>
265  {
266  for (int i=leftIndex; i<=rightIndex; i++)
267  {
268  DeleteIndex(leftIndex);
269  }
270  mIsSplineUpdated=false;
271  }
277  template <class TX,class TY>
278  void BPFTmpl<TX,TY>::DeleteBetweenXValues(const TX& leftX,const TX& rightX)
279  {
280  TIndex leftIndex=GetPosition(leftX);
281  TIndex rightIndex=GetPosition(rightX);
282  DeleteBetweenIndex(leftIndex,rightIndex);
283  mIsSplineUpdated=false;
284  }
285 
291  template <class TX,class TY>
292  void BPFTmpl<TX,TY>::SetIntpType(const EInterpolation& eInterpolation)
293  {
294  meInterpolation=eInterpolation;
295  switch(eInterpolation)
296  {
297  case(EInterpolation::eLinear)://linear interpolation between two closest points
298  {
299  mOrder=1;
300  break;
301  }
302  case(EInterpolation::ePolynomial2)://parabolic interpolation
303  {
304  mOrder=2;
305  break;
306  }
307  case(EInterpolation::ePolynomial3)://3rd order polynomial interpolation
308  {
309  mOrder=3;
310  break;
311  }
312  case(EInterpolation::ePolynomial4)://4th order polynomial interpolation
313  {
314  mOrder=4;
315  break;
316  }
317  case(EInterpolation::ePolynomial5)://5th order polynomial interpolation
318  {
319  mOrder=5;
320  break;
321  }
322  case(EInterpolation::ePolynomialn):/*nth order polynomial interpolation where n is number
323  of points in the BPF-1. Must be less than 10*/
324  {
325  CLAM_ASSERT(Size()<11,"BPF::SetIntpType:Cannot ser more than 10th order interpolation");
326  mOrder=Size()-1;
327  break;
328  }
329  case(EInterpolation::eSpline ): // MRJ: Nice refactoring, but forgot about the Spline...
330  {
331  return;
332  }
333  default:
334  {
335  CLAM_ASSERT( false, "Unsupported interpolation method" );
336  }
337  }
338  const unsigned newSize = mOrder+1;
339  mc.Resize(newSize);
340  mc.SetSize(newSize);
341  md.Resize(newSize);
342  md.SetSize(newSize);
343  }
344 
345 
353  template <class TX,class TY>
354  TY BPFTmpl<TX,TY>::GetValue(const TX& x,const EInterpolation& eInterpolation) const /*Gets value
355  according to the interpolation type set*/
356  {
357  // First we should check if the the x value belongs to the BPF itself
358  // and there is no need to interpolate
359 
360  if(x<mLastX) mLastIndex=0;
361 
362  TIndex i=mSearch.Find( PointTmpl<TX,TY>(x,0.0), mLastIndex);
363  if(i==-1)
364  {
365  // Outside of the BPF; get the first or last value
366  if (GetXValue(0)>x)
367  return GetValueFromIndex(0);
368  else
369  return GetValueFromIndex(Size()-1);
370  }
371  mLastIndex=i;
372  mLastX=x;
373  if(GetXValue(mLastIndex)==x) return GetValueFromIndex(mLastIndex);
374  switch(eInterpolation)
375  {
376  case(EInterpolation::eStep)://returns previous point value
377  {
378  GetnClosest(mLastIndex);
379  if(GetXValue(mClosestPoints[0])<=x)
380  return GetValueFromIndex(mClosestPoints[0]);
381  return GetValueFromIndex(mClosestPoints[0]+1);
382  }
383  case(EInterpolation::eRound)://return closest point value
384  {
385  GetnClosest(mLastIndex);
386  return GetValueFromIndex(mClosestPoints[0]);
387  }
388  case(EInterpolation::eLinear)://linear interpolation between two closest points
389  {
390  TData error=0;
391  if(GetXValue(mLastIndex)<=x)
392  {
393  mClosestPoints[0]=mLastIndex;
394  mClosestPoints[1]=mLastIndex+1;
395  }
396  else
397  {
398  mClosestPoints[0]=mLastIndex-1;
399  mClosestPoints[1]=mLastIndex;
400  }
401  return BPFPolInt(x,mClosestPoints,error);
402  }
403  case(EInterpolation::eSpline)://3rd order spline interpolation
404  {
405  CLAM_ASSERT(mIsSplineUpdated,"BPF::Spline table not updated");
406  return BPFSplineInt(x);//get actual value
407  }
408  case(EInterpolation::ePolynomial2)://parabolic interpolation
409  case(EInterpolation::ePolynomial3)://3rd order polynomial interpolation
410  case(EInterpolation::ePolynomial4)://4th order polynomial interpolation
411  case(EInterpolation::ePolynomial5)://5th order polynomial interpolation
412  {
413  TData error=0;
414  GetnClosest(mLastIndex);
415  return BPFPolInt(x,mClosestPoints,error);
416  }
417  case(EInterpolation::ePolynomialn):/*nth order polynomial interpolation where n is number
418  of points in the BPF-1*/
419  {
420  Array<TIndex> indexArray(mArray.Size());
421  TData error=0;
422  for(TIndex i=0; i<mArray.Size(); i++)
423  {
424  indexArray[i]=i;
425  }
426  return BPFPolInt(x,indexArray,error);
427  }
428  default:
429  {
430  CLAM_ASSERT(false, "Invalid BPF interpolation method.");
431  return 0;
432  }
433  }
434  }
435 
445  template <class TX,class TY>
446  void BPFTmpl<TX,TY>::GetnClosest(TIndex foundIndex) const
447  {
448  CLAM_ASSERT(Size() >= mOrder+1, "BPF size must be at least the interpolation order plus one");
449  int first = foundIndex-mOrder/2;
450  // Fix the first index when it goes too
451  if (first<0) first=0;
452  // At least mOrder+1 fordward
453  int safelimit = Size()-(mOrder+1);
454  if (first>safelimit) first = safelimit;
455 
456  for(int i=0; i<mOrder+1; i++) {
457  mClosestPoints[i]=first+i;
458  }
459 
460  /*
461  TIndex oP=GetPosition(x);
462  int N = mArray.Size();
463  TIndex i, j;
464 
465  TData elemToTake = (n-1.0)/2.0;
466 
467  int toTakeToTheLeft = (int)elemToTake;
468  int toTakeToTheRight = (int)ceil(elemToTake);
469 
470  int dL = (oP) - toTakeToTheLeft;
471  int dR = (N-oP) - toTakeToTheRight;
472 
473  if (dL>= 0 && dR >= 0)
474  {
475  i = oP - toTakeToTheLeft;
476  j = oP + toTakeToTheRight;
477 
478  if (j == N)
479  {
480  j = N-1;
481  i--;
482  }
483  }
484  else if (dL<0 && dR >= 0)
485  {
486  i = 0;
487  j = oP + toTakeToTheRight + (-1)*dL;
488  }
489  else if (dL>=0 && dR < 0)
490  {
491  j = N-1;
492  i = oP - toTakeToTheLeft - (-1)*dR - 1 ;
493 
494  }
495  else if (dL<0 && dR < 0)
496  {
497  throw Err("Not enough Points to interpolate");
498  }
499 
500  for (int k=i; k<=j; k++)
501  {
502  mClosestPoints[k-i] = k;
503  }*/
504 
505  }
506 
512  template <class TX,class TY>
514  {
515  PointTmpl<TX,TY> tmpPoint(x,0);
516  return mSearch.Find(tmpPoint);
517  }
518 
522  template <class TX,class TY>
524  {
525  meInterpolation = originalBPF.meInterpolation;
526  mArray = originalBPF.mArray;
527  mSearch.Set(mArray);
528  mSplineTable = originalBPF.mSplineTable;
529  mIsSplineUpdated = originalBPF.mIsSplineUpdated;
530  mnPoints=originalBPF.mnPoints;
531  mOrder=originalBPF.mOrder;
532  return *this;
533  }
534 
538  template <class TX,class TY>
540  {
541  if(!mIsSplineUpdated)
542  CreateSplineTable(); // Create spline table if not updated
543  mIsSplineUpdated=true;
544  }
545 
549  template <class TX,class TY>
551  {
552  mLeftDerivative = val;
553  mIsSplineUpdated=false;
554  }
555  template <class TX,class TY>
557  {
558  if (mLeftDerivative < Infinity) {
559  mLeftDerivative = Infinity;
560  mIsSplineUpdated=false;
561  }
562  }
566  template <class TX,class TY>
568  {
569  mRightDerivative = val;
570  mIsSplineUpdated=false;
571  }
572  template <class TX,class TY>
574  {
575  if (mLeftDerivative < Infinity) {
576  mRightDerivative = Infinity;
577  mIsSplineUpdated=false;
578  }
579  }
580 
585  template <class TX,class TY>
587  (const TX& x,const Array<TIndex>& closestPointsIndex, TData &errorEstimate)
588  const
589  {
590  int iClosest=0;
591  TX dif=Abs(x-GetXValue(closestPointsIndex[0]));
592  for(int i=0; i<=mOrder; i++)
593  {
594  TX dift=Abs(x-GetXValue(closestPointsIndex[i]));
595  if(dift<dif)
596  {
597  iClosest=i;
598  dif=dift;
599  }
600  md[i]=mc[i]=GetValueFromIndex(closestPointsIndex[i]);
601  }
602 
603  TY y=GetValueFromIndex(closestPointsIndex[iClosest]);
604  for(int m=0; m<mOrder; m++)
605  {
606  for(int i=0; i<mOrder-m; i++)
607  {
608  TX ho=GetXValue(closestPointsIndex[i])-x;
609  TX hp=GetXValue(closestPointsIndex[i+m+1])-x;
610  TX w=mc[i+1]-md[i];
611  TX den=ho-hp;
612  CLAM_ASSERT(den!=0, "Division by zero error interpolating BPF");
613  den=w/den;
614  md[i]=hp*den;
615  mc[i]=ho*den;
616  }
617  if ( 2*iClosest < mOrder-m )
618  {
619  errorEstimate=mc[iClosest];
620  }
621  else
622  {
623  errorEstimate=md[iClosest-1];
624  iClosest--;
625  }
626  y+=errorEstimate;
627  }
628 
629  return y;
630  }
635  template <class TX,class TY>
637  {
638  int n=mArray.Size();
639 
640  CLAM_ASSERT(n > 2,"BPF size too small for spline");
641 
642  Array<TY> u(n);
643  u.SetSize(n);
644  mSplineTable.Resize(n);
645  mSplineTable.SetSize(n);
646 
647  if (mLeftDerivative >= Infinity)
648  mSplineTable[0]=u[0]=0.0; // For a 'natural' spline
649  else {
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))
654  - mLeftDerivative );
655  }
656 
657  for(int i=2;i<=n-1;i++)
658  {
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;
665  }
666 
667  TY qn = 0.0;
668  TY un = 0.0;
669  if (mRightDerivative < Infinity) {
670  qn = 0.5;
671  un = (TData(3.0)/(GetXValue(n-1)-GetXValue(n-2))) *
672  ( mRightDerivative -
673  ( GetValueFromIndex(n-1) - GetValueFromIndex(n-2)) /
674  ( GetXValue(n-1) - GetXValue(n-2)));
675  }
676  mSplineTable[n-1]=((un-qn*u[n-2])/(qn*mSplineTable[n-2]+1));
677  for(int k=n-1; k>=1; k--)
678  {
679  mSplineTable[k-1]=mSplineTable[k-1]*mSplineTable[k]+u[k-1];
680  }
681  }
686  template <class TX,class TY>
687  TY BPFTmpl<TX,TY>::BPFSplineInt(const TX& x) const
688  {
689  int klo=1;
690  int khi=mArray.Size();
691  while(khi-klo>1)
692  {
693  int k=(khi+klo)>>1;
694  if(GetXValue(k-1)>x) khi=k;
695  else klo=k;
696  }
697  TX h=GetXValue(khi-1)-GetXValue(klo-1);
698  CLAM_ASSERT(h!=0.0, "Error interpolating Spline");
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);
703  }
704 
705  template <class TX,class TY>
706  void BPFTmpl<TX,TY>::StoreOn(Storage & storage) const
707  {
708  XMLComponentAdapter adapterInt(meInterpolation,"Interpolation",true);
709  storage.Store(adapterInt);
710  XMLComponentAdapter adapter(mArray, "Points", true);
711  storage.Store(adapter);
712  }
713  template <class TX,class TY>
715  {
716  XMLComponentAdapter adapterInt(meInterpolation,"Interpolation",true);
717  storage.Load(adapterInt);
718  XMLComponentAdapter adapter(mArray, "Points", true);
719  storage.Load(adapter);
720  }
721 
722 } // namespace CLAM
723 
724 
725 #endif // _BPFTmplDef_
726