Point Cloud Library (PCL)  1.9.1-dev
bspline_data.hpp
1 /*
2 Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
3 All rights reserved.
4 
5 Redistribution and use in source and binary forms, with or without modification,
6 are permitted provided that the following conditions are met:
7 
8 Redistributions of source code must retain the above copyright notice, this list of
9 conditions and the following disclaimer. Redistributions in binary form must reproduce
10 the above copyright notice, this list of conditions and the following disclaimer
11 in the documentation and/or other materials provided with the distribution.
12 
13 Neither the name of the Johns Hopkins University nor the names of its contributors
14 may be used to endorse or promote products derived from this software without specific
15 prior written permission.
16 
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
18 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
19 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
20 SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
21 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
22 TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
23 BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
25 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
26 DAMAGE.
27 */
28 
29 #include "poisson_exceptions.h"
30 
31 namespace pcl
32 {
33  namespace poisson
34  {
35 
36  /////////////////
37  // BSplineData //
38  /////////////////
39  // Support[i]:
40  // Odd: i +/- 0.5 * ( 1 + Degree )
41  // i - 0.5 * ( 1 + Degree ) < 0
42  // <=> i < 0.5 * ( 1 + Degree )
43  // i + 0.5 * ( 1 + Degree ) > 0
44  // <=> i > - 0.5 * ( 1 + Degree )
45  // i + 0.5 * ( 1 + Degree ) > r
46  // <=> i > r - 0.5 * ( 1 + Degree )
47  // i - 0.5 * ( 1 + Degree ) < r
48  // <=> i < r + 0.5 * ( 1 + Degree )
49  // Even: i + 0.5 +/- 0.5 * ( 1 + Degree )
50  // i - 0.5 * Degree < 0
51  // <=> i < 0.5 * Degree
52  // i + 1 + 0.5 * Degree > 0
53  // <=> i > -1 - 0.5 * Degree
54  // i + 1 + 0.5 * Degree > r
55  // <=> i > r - 1 - 0.5 * Degree
56  // i - 0.5 * Degree < r
57  // <=> i < r + 0.5 * Degree
58  template< int Degree > inline bool LeftOverlap( unsigned int depth , int offset )
59  {
60  offset <<= 1;
61  if( Degree & 1 ) return (offset < 1+Degree) && (offset > -1-Degree );
62  else return (offset < Degree) && (offset > -2-Degree );
63  }
64  template< int Degree > inline bool RightOverlap( unsigned int depth , int offset )
65  {
66  offset <<= 1;
67  int r = 1<<(depth+1);
68  if( Degree & 1 ) return (offset > 2-1-Degree) && (offset < 2+1+Degree );
69  else return (offset > 2-2-Degree) && (offset < 2+ Degree );
70  }
71  template< int Degree > inline int ReflectLeft( unsigned int depth , int offset )
72  {
73  if( Degree&1 ) return -offset;
74  else return -1-offset;
75  }
76  template< int Degree > inline int ReflectRight( unsigned int depth , int offset )
77  {
78  int r = 1<<(depth+1);
79  if( Degree&1 ) return r -offset;
80  else return r-1-offset;
81  }
82 
83  template< int Degree , class Real >
85  {
86  vvDotTable = dvDotTable = ddDotTable = NULL;
87  valueTables = dValueTables = NULL;
88  baseFunctions = NULL;
89  baseBSplines = NULL;
90  functionCount = sampleCount = 0;
91  }
92 
93  template< int Degree , class Real >
95  {
96  if( functionCount )
97  {
98  if( vvDotTable ) delete[] vvDotTable;
99  if( dvDotTable ) delete[] dvDotTable;
100  if( ddDotTable ) delete[] ddDotTable;
101 
102  if( valueTables ) delete[] valueTables;
103  if( dValueTables ) delete[] dValueTables;
104 
105  if( baseFunctions ) delete[] baseFunctions;
106  if( baseBSplines ) delete[] baseBSplines;
107  }
108  vvDotTable = dvDotTable = ddDotTable = NULL;
109  valueTables = dValueTables=NULL;
110  baseFunctions = NULL;
111  baseBSplines = NULL;
112  functionCount = 0;
113  }
114 
115  template<int Degree,class Real>
116  void BSplineData<Degree,Real>::set( int maxDepth , bool useDotRatios , bool reflectBoundary )
117  {
118  this->useDotRatios = useDotRatios;
119  this->reflectBoundary = reflectBoundary;
120 
121  depth = maxDepth;
122  // [Warning] This assumes that the functions spacing is dual
123  functionCount = BinaryNode< double >::CumulativeCenterCount( depth );
125  baseFunctions = new PPolynomial<Degree>[functionCount];
126  baseBSplines = new BSplineComponents[functionCount];
127 
128  baseFunction = PPolynomial< Degree >::BSpline();
129  for( int i=0 ; i<=Degree ; i++ ) baseBSpline[i] = Polynomial< Degree >::BSplineComponent( i ).shift( double(-(Degree+1)/2) + i - 0.5 );
130  dBaseFunction = baseFunction.derivative();
131  StartingPolynomial< Degree > sPolys[Degree+3];
132 
133  for( int i=0 ; i<Degree+3 ; i++ )
134  {
135  sPolys[i].start = double(-(Degree+1)/2) + i - 1.5;
136  sPolys[i].p *= 0;
137  if( i<=Degree ) sPolys[i].p += baseBSpline[i ].shift( -1 );
138  if( i>=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1];
139  for( int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
140  }
141  leftBaseFunction.set( sPolys , Degree+3 );
142  for( int i=0 ; i<Degree+3 ; i++ )
143  {
144  sPolys[i].start = double(-(Degree+1)/2) + i - 0.5;
145  sPolys[i].p *= 0;
146  if( i<=Degree ) sPolys[i].p += baseBSpline[i ];
147  if( i>=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1].shift( 1 );
148  for( int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
149  }
150  rightBaseFunction.set( sPolys , Degree+3 );
151  dLeftBaseFunction = leftBaseFunction.derivative();
152  dRightBaseFunction = rightBaseFunction.derivative();
153  leftBSpline = rightBSpline = baseBSpline;
154  leftBSpline [1] += leftBSpline[2].shift( -1 ) , leftBSpline[0] *= 0;
155  rightBSpline[1] += rightBSpline[0].shift( 1 ) , rightBSpline[2] *= 0;
156  double c , w;
157  for( int i=0 ; i<functionCount ; i++ )
158  {
160  baseFunctions[i] = baseFunction.scale(w).shift(c);
161  baseBSplines[i] = baseBSpline.scale(w).shift(c);
162  if( reflectBoundary )
163  {
164  int d , off , r;
166  r = 1<<d;
167  if ( off==0 ) baseFunctions[i] = leftBaseFunction.scale(w).shift(c);
168  else if( off==r-1 ) baseFunctions[i] = rightBaseFunction.scale(w).shift(c);
169  if ( off==0 ) baseBSplines[i] = leftBSpline.scale(w).shift(c);
170  else if( off==r-1 ) baseBSplines[i] = rightBSpline.scale(w).shift(c);
171  }
172  }
173  }
174  template<int Degree,class Real>
176  {
177  clearDotTables( flags );
178  int size = ( functionCount*functionCount + functionCount )>>1;
179  int fullSize = functionCount*functionCount;
180  if( flags & VV_DOT_FLAG )
181  {
182  vvDotTable = new Real[size];
183  memset( vvDotTable , 0 , sizeof(Real)*size );
184  }
185  if( flags & DV_DOT_FLAG )
186  {
187  dvDotTable = new Real[fullSize];
188  memset( dvDotTable , 0 , sizeof(Real)*fullSize );
189  }
190  if( flags & DD_DOT_FLAG )
191  {
192  ddDotTable = new Real[size];
193  memset( ddDotTable , 0 , sizeof(Real)*size );
194  }
195  double vvIntegrals[Degree+1][Degree+1];
196  double vdIntegrals[Degree+1][Degree ];
197  double dvIntegrals[Degree ][Degree+1];
198  double ddIntegrals[Degree ][Degree ];
199  int vvSums[Degree+1][Degree+1];
200  int vdSums[Degree+1][Degree ];
201  int dvSums[Degree ][Degree+1];
202  int ddSums[Degree ][Degree ];
203  SetBSplineElementIntegrals< Degree , Degree >( vvIntegrals );
204  SetBSplineElementIntegrals< Degree , Degree-1 >( vdIntegrals );
205  SetBSplineElementIntegrals< Degree-1 , Degree >( dvIntegrals );
206  SetBSplineElementIntegrals< Degree-1 , Degree-1 >( ddIntegrals );
207 
208  for( int d1=0 ; d1<=depth ; d1++ )
209  for( int off1=0 ; off1<(1<<d1) ; off1++ )
210  {
211  int ii = BinaryNode< Real >::CenterIndex( d1 , off1 );
213  BSplineElements< Degree-1 > db1;
214  b1.differentiate( db1 );
215 
216  int start1 , end1;
217 
218  start1 = -1;
219  for( int i=0 ; i<int(b1.size()) ; i++ ) for( int j=0 ; j<=Degree ; j++ )
220  {
221  if( b1[i][j] && start1==-1 ) start1 = i;
222  if( b1[i][j] ) end1 = i+1;
223  }
224  for( int d2=d1 ; d2<=depth ; d2++ )
225  {
226  for( int off2=0 ; off2<(1<<d2) ; off2++ )
227  {
228  int start2 = off2-Degree;
229  int end2 = off2+Degree+1;
230  if( start2>=end1 || start1>=end2 ) continue;
231  start2 = std::max< int >( start1 , start2 );
232  end2 = std::min< int >( end1 , end2 );
233  if( d1==d2 && off2<off1 ) continue;
234  int jj = BinaryNode< Real >::CenterIndex( d2 , off2 );
235  BSplineElements< Degree > b2( 1<<d2 , off2 , reflectBoundary ? BSplineElements<Degree>::NEUMANN : BSplineElements< Degree>::NONE );
236  BSplineElements< Degree-1 > db2;
237  b2.differentiate( db2 );
238 
239  int idx = SymmetricIndex( ii , jj );
240  int idx1 = Index( ii , jj ) , idx2 = Index( jj , ii );
241 
242  memset( vvSums , 0 , sizeof( int ) * ( Degree+1 ) * ( Degree+1 ) );
243  memset( vdSums , 0 , sizeof( int ) * ( Degree+1 ) * ( Degree ) );
244  memset( dvSums , 0 , sizeof( int ) * ( Degree ) * ( Degree+1 ) );
245  memset( ddSums , 0 , sizeof( int ) * ( Degree ) * ( Degree ) );
246  for( int i=start2 ; i<end2 ; i++ )
247  {
248  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) vvSums[j][k] += b1[i][j] * b2[i][k];
249  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) vdSums[j][k] += b1[i][j] * db2[i][k];
250  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) dvSums[j][k] += db1[i][j] * b2[i][k];
251  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) ddSums[j][k] += db1[i][j] * db2[i][k];
252  }
253  double vvDot = 0 , dvDot = 0 , vdDot = 0 , ddDot = 0;
254  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) vvDot += vvIntegrals[j][k] * vvSums[j][k];
255  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) vdDot += vdIntegrals[j][k] * vdSums[j][k];
256  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) dvDot += dvIntegrals[j][k] * dvSums[j][k];
257  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) ddDot += ddIntegrals[j][k] * ddSums[j][k];
258  vvDot /= (1<<d2);
259  ddDot *= (1<<d2);
260  vvDot /= ( b1.denominator * b2.denominator );
261  dvDot /= ( b1.denominator * b2.denominator );
262  vdDot /= ( b1.denominator * b2.denominator );
263  ddDot /= ( b1.denominator * b2.denominator );
264  if( fabs(vvDot)<1e-15 ) continue;
265  if( flags & VV_DOT_FLAG ) vvDotTable [idx] = Real( vvDot );
266  if( useDotRatios )
267  {
268  if( flags & DV_DOT_FLAG ) dvDotTable[idx1] = Real( dvDot / vvDot );
269  if( flags & DV_DOT_FLAG ) dvDotTable[idx2] = Real( vdDot / vvDot );
270  if( flags & DD_DOT_FLAG ) ddDotTable[idx ] = Real( ddDot / vvDot );
271  }
272  else
273  {
274  if( flags & DV_DOT_FLAG ) dvDotTable[idx1] = Real( dvDot );
275  if( flags & DV_DOT_FLAG ) dvDotTable[idx2] = Real( dvDot );
276  if( flags & DD_DOT_FLAG ) ddDotTable[idx ] = Real( ddDot );
277  }
278  }
280  b = b1;
281  b.upSample( b1 );
282  b1.differentiate( db1 );
283  start1 = -1;
284  for( int i=0 ; i<int(b1.size()) ; i++ ) for( int j=0 ; j<=Degree ; j++ )
285  {
286  if( b1[i][j] && start1==-1 ) start1 = i;
287  if( b1[i][j] ) end1 = i+1;
288  }
289  }
290  }
291  }
292  template<int Degree,class Real>
294  {
295  if( (flags & VV_DOT_FLAG) && vvDotTable ) delete[] vvDotTable , vvDotTable = NULL;
296  if( (flags & DV_DOT_FLAG) && dvDotTable ) delete[] dvDotTable , dvDotTable = NULL;
297  if( (flags & DD_DOT_FLAG) && ddDotTable ) delete[] ddDotTable , ddDotTable = NULL;
298  }
299  template< int Degree , class Real >
300  void BSplineData< Degree , Real >::setSampleSpan( int idx , int& start , int& end , double smooth ) const
301  {
302  int d , off , res;
303  BinaryNode< double >::DepthAndOffset( idx , d , off );
304  res = 1<<d;
305  double _start = ( off + 0.5 - 0.5*(Degree+1) ) / res - smooth;
306  double _end = ( off + 0.5 + 0.5*(Degree+1) ) / res + smooth;
307  // (start)/(sampleCount-1) >_start && (start-1)/(sampleCount-1)<=_start
308  // => start > _start * (sampleCount-1 ) && start <= _start*(sampleCount-1) + 1
309  // => _start * (sampleCount-1) + 1 >= start > _start * (sampleCount-1)
310  start = int( floor( _start * (sampleCount-1) + 1 ) );
311  if( start<0 ) start = 0;
312  // (end)/(sampleCount-1)<_end && (end+1)/(sampleCount-1)>=_end
313  // => end < _end * (sampleCount-1 ) && end >= _end*(sampleCount-1) - 1
314  // => _end * (sampleCount-1) > end >= _end * (sampleCount-1) - 1
315  end = int( ceil( _end * (sampleCount-1) - 1 ) );
316  if( end>=sampleCount ) end = sampleCount-1;
317  }
318  template<int Degree,class Real>
319  void BSplineData<Degree,Real>::setValueTables( int flags , double smooth )
320  {
321  clearValueTables();
322  if( flags & VALUE_FLAG ) valueTables = new Real[functionCount*sampleCount];
323  if( flags & D_VALUE_FLAG ) dValueTables = new Real[functionCount*sampleCount];
324  PPolynomial<Degree+1> function;
325  PPolynomial<Degree> dFunction;
326  for( int i=0 ; i<functionCount ; i++ )
327  {
328  if(smooth>0)
329  {
330  function = baseFunctions[i].MovingAverage(smooth);
331  dFunction = baseFunctions[i].derivative().MovingAverage(smooth);
332  }
333  else
334  {
335  function = baseFunctions[i];
336  dFunction = baseFunctions[i].derivative();
337  }
338  for( int j=0 ; j<sampleCount ; j++ )
339  {
340  double x=double(j)/(sampleCount-1);
341  if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=Real( function(x));}
342  if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=Real(dFunction(x));}
343  }
344  }
345  }
346  template<int Degree,class Real>
347  void BSplineData<Degree,Real>::setValueTables(int flags,double valueSmooth,double normalSmooth){
348  clearValueTables();
349  if(flags & VALUE_FLAG){ valueTables=new Real[functionCount*sampleCount];}
350  if(flags & D_VALUE_FLAG){dValueTables=new Real[functionCount*sampleCount];}
351  PPolynomial<Degree+1> function;
352  PPolynomial<Degree> dFunction;
353  for(int i=0;i<functionCount;i++){
354  if(valueSmooth>0) { function=baseFunctions[i].MovingAverage(valueSmooth);}
355  else { function=baseFunctions[i];}
356  if(normalSmooth>0) {dFunction=baseFunctions[i].derivative().MovingAverage(normalSmooth);}
357  else {dFunction=baseFunctions[i].derivative();}
358 
359  for(int j=0;j<sampleCount;j++){
360  double x=double(j)/(sampleCount-1);
361  if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=Real( function(x));}
362  if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=Real(dFunction(x));}
363  }
364  }
365  }
366 
367 
368  template<int Degree,class Real>
370  if( valueTables){delete[] valueTables;}
371  if(dValueTables){delete[] dValueTables;}
372  valueTables=dValueTables=NULL;
373  }
374 
375  template<int Degree,class Real>
376  inline int BSplineData<Degree,Real>::Index( int i1 , int i2 ) const { return i1*functionCount+i2; }
377  template<int Degree,class Real>
378  inline int BSplineData<Degree,Real>::SymmetricIndex( int i1 , int i2 )
379  {
380  if( i1>i2 ) return ((i1*i1+i1)>>1)+i2;
381  else return ((i2*i2+i2)>>1)+i1;
382  }
383  template<int Degree,class Real>
384  inline int BSplineData<Degree,Real>::SymmetricIndex( int i1 , int i2 , int& index )
385  {
386  if( i1<i2 )
387  {
388  index = ((i2*i2+i2)>>1)+i1;
389  return 1;
390  }
391  else
392  {
393  index = ((i1*i1+i1)>>1)+i2;
394  return 0;
395  }
396  }
397 
398 
399  ////////////////////////
400  // BSplineElementData //
401  ////////////////////////
402  template< int Degree >
403  BSplineElements< Degree >::BSplineElements( int res , int offset , int boundary )
404  {
405  denominator = 1;
406  this->resize( res , BSplineElementCoefficients<Degree>() );
407 
408  for( int i=0 ; i<=Degree ; i++ )
409  {
410  int idx = -_off + offset + i;
411  if( idx>=0 && idx<res ) (*this)[idx][i] = 1;
412  }
413  if( boundary!=0 )
414  {
415  _addLeft( offset-2*res , boundary ) , _addRight( offset+2*res , boundary );
416  if( Degree&1 ) _addLeft( offset-res , boundary ) , _addRight( offset+res , boundary );
417  else _addLeft( -offset-1 , boundary ) , _addRight( -offset-1+2*res , boundary );
418  }
419  }
420  template< int Degree >
421  void BSplineElements< Degree >::_addLeft( int offset , int boundary )
422  {
423  int res = int( this->size() );
424  bool set = false;
425  for( int i=0 ; i<=Degree ; i++ )
426  {
427  int idx = -_off + offset + i;
428  if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set = true;
429  }
430  if( set ) _addLeft( offset-2*res , boundary );
431  }
432  template< int Degree >
433  void BSplineElements< Degree >::_addRight( int offset , int boundary )
434  {
435  int res = int( this->size() );
436  bool set = false;
437  for( int i=0 ; i<=Degree ; i++ )
438  {
439  int idx = -_off + offset + i;
440  if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set = true;
441  }
442  if( set ) _addRight( offset+2*res , boundary );
443  }
444  template< int Degree >
446  {
447  POISSON_THROW_EXCEPTION (pcl::poisson::PoissonBadArgumentException, "B-spline up-sampling not supported for degree " << Degree);
448  }
449  template<>
451  {
452  high.resize( size()*2 );
453  high.assign( high.size() , BSplineElementCoefficients<1>() );
454  for( int i=0 ; i<int(size()) ; i++ )
455  {
456  high[2*i+0][0] += 1 * (*this)[i][0];
457  high[2*i+0][1] += 0 * (*this)[i][0];
458  high[2*i+1][0] += 2 * (*this)[i][0];
459  high[2*i+1][1] += 1 * (*this)[i][0];
460 
461  high[2*i+0][0] += 1 * (*this)[i][1];
462  high[2*i+0][1] += 2 * (*this)[i][1];
463  high[2*i+1][0] += 0 * (*this)[i][1];
464  high[2*i+1][1] += 1 * (*this)[i][1];
465  }
466  high.denominator = denominator * 2;
467  }
468  template<>
470  {
471  // /----\
472  // / \
473  // / \ = 1 /--\ +3 /--\ +3 /--\ +1 /--\
474  // / \ / \ / \ / \ / \
475  // |----------| |----------| |----------| |----------| |----------|
476 
477  high.resize( size()*2 );
478  high.assign( high.size() , BSplineElementCoefficients<2>() );
479  for( int i=0 ; i<int(size()) ; i++ )
480  {
481  high[2*i+0][0] += 1 * (*this)[i][0];
482  high[2*i+0][1] += 0 * (*this)[i][0];
483  high[2*i+0][2] += 0 * (*this)[i][0];
484  high[2*i+1][0] += 3 * (*this)[i][0];
485  high[2*i+1][1] += 1 * (*this)[i][0];
486  high[2*i+1][2] += 0 * (*this)[i][0];
487 
488  high[2*i+0][0] += 3 * (*this)[i][1];
489  high[2*i+0][1] += 3 * (*this)[i][1];
490  high[2*i+0][2] += 1 * (*this)[i][1];
491  high[2*i+1][0] += 1 * (*this)[i][1];
492  high[2*i+1][1] += 3 * (*this)[i][1];
493  high[2*i+1][2] += 3 * (*this)[i][1];
494 
495  high[2*i+0][0] += 0 * (*this)[i][2];
496  high[2*i+0][1] += 1 * (*this)[i][2];
497  high[2*i+0][2] += 3 * (*this)[i][2];
498  high[2*i+1][0] += 0 * (*this)[i][2];
499  high[2*i+1][1] += 0 * (*this)[i][2];
500  high[2*i+1][2] += 1 * (*this)[i][2];
501  }
502  high.denominator = denominator * 4;
503  }
504 
505  template< int Degree >
507  {
508  d.resize( this->size() );
509  d.assign( d.size() , BSplineElementCoefficients< Degree-1 >() );
510  for( int i=0 ; i<int(this->size()) ; i++ ) for( int j=0 ; j<=Degree ; j++ )
511  {
512  if( j-1>=0 ) d[i][j-1] -= (*this)[i][j];
513  if( j<Degree ) d[i][j ] += (*this)[i][j];
514  }
515  d.denominator = denominator;
516  }
517  // If we were really good, we would implement this integral table to store
518  // rational values to improve precision...
519  template< int Degree1 , int Degree2 >
520  void SetBSplineElementIntegrals( double integrals[Degree1+1][Degree2+1] )
521  {
522  for( int i=0 ; i<=Degree1 ; i++ )
523  {
525  for( int j=0 ; j<=Degree2 ; j++ )
526  {
528  integrals[i][j] = ( p1 * p2 ).integral( 0 , 1 );
529  }
530  }
531  }
532 
533 
534  }
535 }
PPolynomial< Degree-1 > derivative(void) const
virtual void clearDotTables(int flags)
bool LeftOverlap(unsigned int depth, int offset)
static PPolynomial BSpline(double radius=0.5)
StartingPolynomial shift(double t) const
Definition: ppolynomial.hpp:58
Polynomial< Degree > p
Definition: ppolynomial.h:46
This file defines compatibility wrappers for low level I/O functions.
Definition: convolution.h:45
static int SymmetricIndex(int i1, int i2)
PPolynomial< Degree+1 > MovingAverage(double radius)
static int CumulativeCenterCount(int maxDepth)
Definition: binary_node.h:44
An exception that is thrown when the arguments number or type is wrong/unhandled. ...
void differentiate(BSplineElements< Degree-1 > &d) const
static int CenterCount(int depth)
Definition: binary_node.h:42
void _addLeft(int offset, int boundary)
void upSample(BSplineElements &high) const
virtual void setValueTables(int flags, double smooth=0)
int ReflectLeft(unsigned int depth, int offset)
void _addRight(int offset, int boundary)
virtual void clearValueTables(void)
static void CenterAndWidth(int depth, int offset, Real &center, Real &width)
Definition: binary_node.h:53
int Index(int i1, int i2) const
virtual void setDotTables(int flags)
bool RightOverlap(unsigned int depth, int offset)
int ReflectRight(unsigned int depth, int offset)
void set(int maxDepth, bool useDotRatios=true, bool reflectBoundary=false)
static int CenterIndex(int depth, int offSet)
Definition: binary_node.h:47
void SetBSplineElementIntegrals(double integrals[Degree1+1][Degree2+1])
static int CornerCount(int depth)
Definition: binary_node.h:43
static Polynomial BSplineComponent(int i)
Definition: polynomial.hpp:307
void setSampleSpan(int idx, int &start, int &end, double smooth=0) const
static void DepthAndOffset(int idx, int &depth, int &offset)
Definition: binary_node.h:64