Point Cloud Library (PCL)  1.9.1
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 
30 namespace pcl
31 {
32  namespace poisson
33  {
34 
35  /////////////////
36  // BSplineData //
37  /////////////////
38  // Support[i]:
39  // Odd: i +/- 0.5 * ( 1 + Degree )
40  // i - 0.5 * ( 1 + Degree ) < 0
41  // <=> i < 0.5 * ( 1 + Degree )
42  // i + 0.5 * ( 1 + Degree ) > 0
43  // <=> i > - 0.5 * ( 1 + Degree )
44  // i + 0.5 * ( 1 + Degree ) > r
45  // <=> i > r - 0.5 * ( 1 + Degree )
46  // i - 0.5 * ( 1 + Degree ) < r
47  // <=> i < r + 0.5 * ( 1 + Degree )
48  // Even: i + 0.5 +/- 0.5 * ( 1 + Degree )
49  // i - 0.5 * Degree < 0
50  // <=> i < 0.5 * Degree
51  // i + 1 + 0.5 * Degree > 0
52  // <=> i > -1 - 0.5 * Degree
53  // i + 1 + 0.5 * Degree > r
54  // <=> i > r - 1 - 0.5 * Degree
55  // i - 0.5 * Degree < r
56  // <=> i < r + 0.5 * Degree
57  template< int Degree > inline bool LeftOverlap( unsigned int depth , int offset )
58  {
59  offset <<= 1;
60  if( Degree & 1 ) return (offset < 1+Degree) && (offset > -1-Degree );
61  else return (offset < Degree) && (offset > -2-Degree );
62  }
63  template< int Degree > inline bool RightOverlap( unsigned int depth , int offset )
64  {
65  offset <<= 1;
66  int r = 1<<(depth+1);
67  if( Degree & 1 ) return (offset > 2-1-Degree) && (offset < 2+1+Degree );
68  else return (offset > 2-2-Degree) && (offset < 2+ Degree );
69  }
70  template< int Degree > inline int ReflectLeft( unsigned int depth , int offset )
71  {
72  if( Degree&1 ) return -offset;
73  else return -1-offset;
74  }
75  template< int Degree > inline int ReflectRight( unsigned int depth , int offset )
76  {
77  int r = 1<<(depth+1);
78  if( Degree&1 ) return r -offset;
79  else return r-1-offset;
80  }
81 
82  template< int Degree , class Real >
84  {
85  vvDotTable = dvDotTable = ddDotTable = NULL;
86  valueTables = dValueTables = NULL;
87  baseFunctions = NULL;
88  baseBSplines = NULL;
89  functionCount = sampleCount = 0;
90  }
91 
92  template< int Degree , class Real >
94  {
95  if( functionCount )
96  {
97  if( vvDotTable ) delete[] vvDotTable;
98  if( dvDotTable ) delete[] dvDotTable;
99  if( ddDotTable ) delete[] ddDotTable;
100 
101  if( valueTables ) delete[] valueTables;
102  if( dValueTables ) delete[] dValueTables;
103 
104  if( baseFunctions ) delete[] baseFunctions;
105  if( baseBSplines ) delete[] baseBSplines;
106  }
107  vvDotTable = dvDotTable = ddDotTable = NULL;
108  valueTables = dValueTables=NULL;
109  baseFunctions = NULL;
110  baseBSplines = NULL;
111  functionCount = 0;
112  }
113 
114  template<int Degree,class Real>
115  void BSplineData<Degree,Real>::set( int maxDepth , bool useDotRatios , bool reflectBoundary )
116  {
117  this->useDotRatios = useDotRatios;
118  this->reflectBoundary = reflectBoundary;
119 
120  depth = maxDepth;
121  // [Warning] This assumes that the functions spacing is dual
122  functionCount = BinaryNode< double >::CumulativeCenterCount( depth );
124  baseFunctions = new PPolynomial<Degree>[functionCount];
125  baseBSplines = new BSplineComponents[functionCount];
126 
127  baseFunction = PPolynomial< Degree >::BSpline();
128  for( int i=0 ; i<=Degree ; i++ ) baseBSpline[i] = Polynomial< Degree >::BSplineComponent( i ).shift( double(-(Degree+1)/2) + i - 0.5 );
129  dBaseFunction = baseFunction.derivative();
130  StartingPolynomial< Degree > sPolys[Degree+3];
131 
132  for( int i=0 ; i<Degree+3 ; i++ )
133  {
134  sPolys[i].start = double(-(Degree+1)/2) + i - 1.5;
135  sPolys[i].p *= 0;
136  if( i<=Degree ) sPolys[i].p += baseBSpline[i ].shift( -1 );
137  if( i>=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1];
138  for( int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
139  }
140  leftBaseFunction.set( sPolys , Degree+3 );
141  for( int i=0 ; i<Degree+3 ; i++ )
142  {
143  sPolys[i].start = double(-(Degree+1)/2) + i - 0.5;
144  sPolys[i].p *= 0;
145  if( i<=Degree ) sPolys[i].p += baseBSpline[i ];
146  if( i>=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1].shift( 1 );
147  for( int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
148  }
149  rightBaseFunction.set( sPolys , Degree+3 );
150  dLeftBaseFunction = leftBaseFunction.derivative();
151  dRightBaseFunction = rightBaseFunction.derivative();
152  leftBSpline = rightBSpline = baseBSpline;
153  leftBSpline [1] += leftBSpline[2].shift( -1 ) , leftBSpline[0] *= 0;
154  rightBSpline[1] += rightBSpline[0].shift( 1 ) , rightBSpline[2] *= 0;
155  double c , w;
156  for( int i=0 ; i<functionCount ; i++ )
157  {
159  baseFunctions[i] = baseFunction.scale(w).shift(c);
160  baseBSplines[i] = baseBSpline.scale(w).shift(c);
161  if( reflectBoundary )
162  {
163  int d , off , r;
165  r = 1<<d;
166  if ( off==0 ) baseFunctions[i] = leftBaseFunction.scale(w).shift(c);
167  else if( off==r-1 ) baseFunctions[i] = rightBaseFunction.scale(w).shift(c);
168  if ( off==0 ) baseBSplines[i] = leftBSpline.scale(w).shift(c);
169  else if( off==r-1 ) baseBSplines[i] = rightBSpline.scale(w).shift(c);
170  }
171  }
172  }
173  template<int Degree,class Real>
175  {
176  clearDotTables( flags );
177  int size = ( functionCount*functionCount + functionCount )>>1;
178  int fullSize = functionCount*functionCount;
179  if( flags & VV_DOT_FLAG )
180  {
181  vvDotTable = new Real[size];
182  memset( vvDotTable , 0 , sizeof(Real)*size );
183  }
184  if( flags & DV_DOT_FLAG )
185  {
186  dvDotTable = new Real[fullSize];
187  memset( dvDotTable , 0 , sizeof(Real)*fullSize );
188  }
189  if( flags & DD_DOT_FLAG )
190  {
191  ddDotTable = new Real[size];
192  memset( ddDotTable , 0 , sizeof(Real)*size );
193  }
194  double vvIntegrals[Degree+1][Degree+1];
195  double vdIntegrals[Degree+1][Degree ];
196  double dvIntegrals[Degree ][Degree+1];
197  double ddIntegrals[Degree ][Degree ];
198  int vvSums[Degree+1][Degree+1];
199  int vdSums[Degree+1][Degree ];
200  int dvSums[Degree ][Degree+1];
201  int ddSums[Degree ][Degree ];
202  SetBSplineElementIntegrals< Degree , Degree >( vvIntegrals );
203  SetBSplineElementIntegrals< Degree , Degree-1 >( vdIntegrals );
204  SetBSplineElementIntegrals< Degree-1 , Degree >( dvIntegrals );
205  SetBSplineElementIntegrals< Degree-1 , Degree-1 >( ddIntegrals );
206 
207  for( int d1=0 ; d1<=depth ; d1++ )
208  for( int off1=0 ; off1<(1<<d1) ; off1++ )
209  {
210  int ii = BinaryNode< Real >::CenterIndex( d1 , off1 );
212  BSplineElements< Degree-1 > db1;
213  b1.differentiate( db1 );
214 
215  int start1 , end1;
216 
217  start1 = -1;
218  for( int i=0 ; i<int(b1.size()) ; i++ ) for( int j=0 ; j<=Degree ; j++ )
219  {
220  if( b1[i][j] && start1==-1 ) start1 = i;
221  if( b1[i][j] ) end1 = i+1;
222  }
223  for( int d2=d1 ; d2<=depth ; d2++ )
224  {
225  for( int off2=0 ; off2<(1<<d2) ; off2++ )
226  {
227  int start2 = off2-Degree;
228  int end2 = off2+Degree+1;
229  if( start2>=end1 || start1>=end2 ) continue;
230  start2 = std::max< int >( start1 , start2 );
231  end2 = std::min< int >( end1 , end2 );
232  if( d1==d2 && off2<off1 ) continue;
233  int jj = BinaryNode< Real >::CenterIndex( d2 , off2 );
235  BSplineElements< Degree-1 > db2;
236  b2.differentiate( db2 );
237 
238  int idx = SymmetricIndex( ii , jj );
239  int idx1 = Index( ii , jj ) , idx2 = Index( jj , ii );
240 
241  memset( vvSums , 0 , sizeof( int ) * ( Degree+1 ) * ( Degree+1 ) );
242  memset( vdSums , 0 , sizeof( int ) * ( Degree+1 ) * ( Degree ) );
243  memset( dvSums , 0 , sizeof( int ) * ( Degree ) * ( Degree+1 ) );
244  memset( ddSums , 0 , sizeof( int ) * ( Degree ) * ( Degree ) );
245  for( int i=start2 ; i<end2 ; i++ )
246  {
247  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) vvSums[j][k] += b1[i][j] * b2[i][k];
248  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) vdSums[j][k] += b1[i][j] * db2[i][k];
249  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) dvSums[j][k] += db1[i][j] * b2[i][k];
250  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) ddSums[j][k] += db1[i][j] * db2[i][k];
251  }
252  double vvDot = 0 , dvDot = 0 , vdDot = 0 , ddDot = 0;
253  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) vvDot += vvIntegrals[j][k] * vvSums[j][k];
254  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) vdDot += vdIntegrals[j][k] * vdSums[j][k];
255  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) dvDot += dvIntegrals[j][k] * dvSums[j][k];
256  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) ddDot += ddIntegrals[j][k] * ddSums[j][k];
257  vvDot /= (1<<d2);
258  ddDot *= (1<<d2);
259  vvDot /= ( b1.denominator * b2.denominator );
260  dvDot /= ( b1.denominator * b2.denominator );
261  vdDot /= ( b1.denominator * b2.denominator );
262  ddDot /= ( b1.denominator * b2.denominator );
263  if( fabs(vvDot)<1e-15 ) continue;
264  if( flags & VV_DOT_FLAG ) vvDotTable [idx] = Real( vvDot );
265  if( useDotRatios )
266  {
267  if( flags & DV_DOT_FLAG ) dvDotTable[idx1] = Real( dvDot / vvDot );
268  if( flags & DV_DOT_FLAG ) dvDotTable[idx2] = Real( vdDot / vvDot );
269  if( flags & DD_DOT_FLAG ) ddDotTable[idx ] = Real( ddDot / vvDot );
270  }
271  else
272  {
273  if( flags & DV_DOT_FLAG ) dvDotTable[idx1] = Real( dvDot );
274  if( flags & DV_DOT_FLAG ) dvDotTable[idx2] = Real( dvDot );
275  if( flags & DD_DOT_FLAG ) ddDotTable[idx ] = Real( ddDot );
276  }
277  }
279  b = b1;
280  b.upSample( b1 );
281  b1.differentiate( db1 );
282  start1 = -1;
283  for( int i=0 ; i<int(b1.size()) ; i++ ) for( int j=0 ; j<=Degree ; j++ )
284  {
285  if( b1[i][j] && start1==-1 ) start1 = i;
286  if( b1[i][j] ) end1 = i+1;
287  }
288  }
289  }
290  }
291  template<int Degree,class Real>
293  {
294  if( (flags & VV_DOT_FLAG) && vvDotTable ) delete[] vvDotTable , vvDotTable = NULL;
295  if( (flags & DV_DOT_FLAG) && dvDotTable ) delete[] dvDotTable , dvDotTable = NULL;
296  if( (flags & DD_DOT_FLAG) && ddDotTable ) delete[] ddDotTable , ddDotTable = NULL;
297  }
298  template< int Degree , class Real >
299  void BSplineData< Degree , Real >::setSampleSpan( int idx , int& start , int& end , double smooth ) const
300  {
301  int d , off , res;
302  BinaryNode< double >::DepthAndOffset( idx , d , off );
303  res = 1<<d;
304  double _start = ( off + 0.5 - 0.5*(Degree+1) ) / res - smooth;
305  double _end = ( off + 0.5 + 0.5*(Degree+1) ) / res + smooth;
306  // (start)/(sampleCount-1) >_start && (start-1)/(sampleCount-1)<=_start
307  // => start > _start * (sampleCount-1 ) && start <= _start*(sampleCount-1) + 1
308  // => _start * (sampleCount-1) + 1 >= start > _start * (sampleCount-1)
309  start = int( floor( _start * (sampleCount-1) + 1 ) );
310  if( start<0 ) start = 0;
311  // (end)/(sampleCount-1)<_end && (end+1)/(sampleCount-1)>=_end
312  // => end < _end * (sampleCount-1 ) && end >= _end*(sampleCount-1) - 1
313  // => _end * (sampleCount-1) > end >= _end * (sampleCount-1) - 1
314  end = int( ceil( _end * (sampleCount-1) - 1 ) );
315  if( end>=sampleCount ) end = sampleCount-1;
316  }
317  template<int Degree,class Real>
318  void BSplineData<Degree,Real>::setValueTables( int flags , double smooth )
319  {
320  clearValueTables();
321  if( flags & VALUE_FLAG ) valueTables = new Real[functionCount*sampleCount];
322  if( flags & D_VALUE_FLAG ) dValueTables = new Real[functionCount*sampleCount];
323  PPolynomial<Degree+1> function;
324  PPolynomial<Degree> dFunction;
325  for( int i=0 ; i<functionCount ; i++ )
326  {
327  if(smooth>0)
328  {
329  function = baseFunctions[i].MovingAverage(smooth);
330  dFunction = baseFunctions[i].derivative().MovingAverage(smooth);
331  }
332  else
333  {
334  function = baseFunctions[i];
335  dFunction = baseFunctions[i].derivative();
336  }
337  for( int j=0 ; j<sampleCount ; j++ )
338  {
339  double x=double(j)/(sampleCount-1);
340  if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=Real( function(x));}
341  if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=Real(dFunction(x));}
342  }
343  }
344  }
345  template<int Degree,class Real>
346  void BSplineData<Degree,Real>::setValueTables(int flags,double valueSmooth,double normalSmooth){
347  clearValueTables();
348  if(flags & VALUE_FLAG){ valueTables=new Real[functionCount*sampleCount];}
349  if(flags & D_VALUE_FLAG){dValueTables=new Real[functionCount*sampleCount];}
350  PPolynomial<Degree+1> function;
351  PPolynomial<Degree> dFunction;
352  for(int i=0;i<functionCount;i++){
353  if(valueSmooth>0) { function=baseFunctions[i].MovingAverage(valueSmooth);}
354  else { function=baseFunctions[i];}
355  if(normalSmooth>0) {dFunction=baseFunctions[i].derivative().MovingAverage(normalSmooth);}
356  else {dFunction=baseFunctions[i].derivative();}
357 
358  for(int j=0;j<sampleCount;j++){
359  double x=double(j)/(sampleCount-1);
360  if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=Real( function(x));}
361  if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=Real(dFunction(x));}
362  }
363  }
364  }
365 
366 
367  template<int Degree,class Real>
369  if( valueTables){delete[] valueTables;}
370  if(dValueTables){delete[] dValueTables;}
371  valueTables=dValueTables=NULL;
372  }
373 
374  template<int Degree,class Real>
375  inline int BSplineData<Degree,Real>::Index( int i1 , int i2 ) const { return i1*functionCount+i2; }
376  template<int Degree,class Real>
377  inline int BSplineData<Degree,Real>::SymmetricIndex( int i1 , int i2 )
378  {
379  if( i1>i2 ) return ((i1*i1+i1)>>1)+i2;
380  else return ((i2*i2+i2)>>1)+i1;
381  }
382  template<int Degree,class Real>
383  inline int BSplineData<Degree,Real>::SymmetricIndex( int i1 , int i2 , int& index )
384  {
385  if( i1<i2 )
386  {
387  index = ((i2*i2+i2)>>1)+i1;
388  return 1;
389  }
390  else
391  {
392  index = ((i1*i1+i1)>>1)+i2;
393  return 0;
394  }
395  }
396 
397 
398  ////////////////////////
399  // BSplineElementData //
400  ////////////////////////
401  template< int Degree >
402  BSplineElements< Degree >::BSplineElements( int res , int offset , int boundary )
403  {
404  denominator = 1;
405  this->resize( res , BSplineElementCoefficients<Degree>() );
406 
407  for( int i=0 ; i<=Degree ; i++ )
408  {
409  int idx = -_off + offset + i;
410  if( idx>=0 && idx<res ) (*this)[idx][i] = 1;
411  }
412  if( boundary!=0 )
413  {
414  _addLeft( offset-2*res , boundary ) , _addRight( offset+2*res , boundary );
415  if( Degree&1 ) _addLeft( offset-res , boundary ) , _addRight( offset+res , boundary );
416  else _addLeft( -offset-1 , boundary ) , _addRight( -offset-1+2*res , boundary );
417  }
418  }
419  template< int Degree >
420  void BSplineElements< Degree >::_addLeft( int offset , int boundary )
421  {
422  int res = int( this->size() );
423  bool set = false;
424  for( int i=0 ; i<=Degree ; i++ )
425  {
426  int idx = -_off + offset + i;
427  if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set = true;
428  }
429  if( set ) _addLeft( offset-2*res , boundary );
430  }
431  template< int Degree >
432  void BSplineElements< Degree >::_addRight( int offset , int boundary )
433  {
434  int res = int( this->size() );
435  bool set = false;
436  for( int i=0 ; i<=Degree ; i++ )
437  {
438  int idx = -_off + offset + i;
439  if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set = true;
440  }
441  if( set ) _addRight( offset+2*res , boundary );
442  }
443  template< int Degree >
445  {
446  fprintf( stderr , "[ERROR] B-spline up-sampling not supported for degree %d\n" , Degree );
447  exit( 0 );
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
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