Developer Documentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
BSplineData.inl
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 
30 // BSplineData //
32 // Support[i]:
33 // Odd: i +/- 0.5 * ( 1 + Degree )
34 // i - 0.5 * ( 1 + Degree ) < 0
35 // <=> i < 0.5 * ( 1 + Degree )
36 // i + 0.5 * ( 1 + Degree ) > 0
37 // <=> i > - 0.5 * ( 1 + Degree )
38 // i + 0.5 * ( 1 + Degree ) > r
39 // <=> i > r - 0.5 * ( 1 + Degree )
40 // i - 0.5 * ( 1 + Degree ) < r
41 // <=> i < r + 0.5 * ( 1 + Degree )
42 // Even: i + 0.5 +/- 0.5 * ( 1 + Degree )
43 // i - 0.5 * Degree < 0
44 // <=> i < 0.5 * Degree
45 // i + 1 + 0.5 * Degree > 0
46 // <=> i > -1 - 0.5 * Degree
47 // i + 1 + 0.5 * Degree > r
48 // <=> i > r - 1 - 0.5 * Degree
49 // i - 0.5 * Degree < r
50 // <=> i < r + 0.5 * Degree
51 template< int Degree > inline bool LeftOverlap( unsigned int depth , int offset )
52 {
53  offset <<= 1;
54  if( Degree & 1 ) return (offset < 1+Degree) && (offset > -1-Degree );
55  else return (offset < Degree) && (offset > -2-Degree );
56 }
57 template< int Degree > inline bool RightOverlap( unsigned int depth , int offset )
58 {
59  offset <<= 1;
60 // int r = 1<<(depth+1);
61  if( Degree & 1 ) return (offset > 2-1-Degree) && (offset < 2+1+Degree );
62  else return (offset > 2-2-Degree) && (offset < 2+ Degree );
63 }
64 template< int Degree > inline int ReflectLeft( unsigned int depth , int offset )
65 {
66  if( Degree&1 ) return -offset;
67  else return -1-offset;
68 }
69 template< int Degree > inline int ReflectRight( unsigned int depth , int offset )
70 {
71  int r = 1<<(depth+1);
72  if( Degree&1 ) return r -offset;
73  else return r-1-offset;
74 }
75 
76 template<int Degree,class Real>
77 BSplineData<Degree,Real>::BSplineData( void ) :
78  useDotRatios(false),
79  boundaryType(0),
80  depth(0),
81  functionCount(0),
82  sampleCount(0),
83  baseFunctions(NULL),
84  baseBSplines(NULL)
85 {
86  vvDotTable = dvDotTable = ddDotTable = NullPointer< Real >();
87  valueTables = dValueTables = NullPointer< Real >();
88 }
89 
90 template<int Degree,class Real>
91 BSplineData< Degree , Real >::~BSplineData(void)
92 {
93  if( functionCount )
94  {
95  if( vvDotTable ) DeletePointer( vvDotTable );
96  if( dvDotTable ) DeletePointer( dvDotTable );
97  if( ddDotTable ) DeletePointer( ddDotTable );
98  if( valueTables ) DeletePointer( valueTables );
99  if( dValueTables ) DeletePointer( dValueTables );
100  }
101  functionCount = 0;
102 }
103 
104 template<int Degree,class Real>
105 void BSplineData<Degree,Real>::set( int maxDepth , bool useDotRatios , int boundaryType )
106 {
107  this->useDotRatios = useDotRatios;
108  this->boundaryType = boundaryType;
109 
110  depth = maxDepth;
111  // [Warning] This assumes that the functions spacing is dual
112  functionCount = BinaryNode< double >::CumulativeCenterCount( depth );
113  sampleCount = BinaryNode< double >::CenterCount( depth ) + BinaryNode< double >::CornerCount( depth );
114  baseFunctions = NewPointer< PPolynomial< Degree > >( functionCount );
115  baseBSplines = NewPointer< BSplineComponents >( functionCount );
116 
117  baseFunction = PPolynomial< Degree >::BSpline();
118  for( int i=0 ; i<=Degree ; i++ ) baseBSpline[i] = Polynomial< Degree >::BSplineComponent( i ).shift( double(-(Degree+1)/2) + i - 0.5 );
119  dBaseFunction = baseFunction.derivative();
120  StartingPolynomial< Degree > sPolys[Degree+4];
121 
122  for( int i=0 ; i<Degree+3 ; i++ )
123  {
124  sPolys[i].start = double(-(Degree+1)/2) + i - 1.5;
125  sPolys[i].p *= 0;
126  if( i<=Degree ) sPolys[i].p += baseBSpline[i ].shift( -1 ) * boundaryType;
127  if( i>=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1];
128  for( int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
129  }
130  leftBaseFunction.set( sPolys , Degree+3 );
131  for( int i=0 ; i<Degree+3 ; i++ )
132  {
133  sPolys[i].start = double(-(Degree+1)/2) + i - 0.5;
134  sPolys[i].p *= 0;
135  if( i<=Degree ) sPolys[i].p += baseBSpline[i ];
136  if( i>=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1].shift( 1 ) * boundaryType;
137  for( int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
138  }
139  rightBaseFunction.set( sPolys , Degree+3 );
140  for( int i=0 ; i<Degree+4 ; i++ )
141  {
142  sPolys[i].start = double(-(Degree+1)/2) + i - 1.5;
143  sPolys[i].p *= 0;
144  if( i<=Degree ) sPolys[i].p += baseBSpline[i ].shift( -1 ) * boundaryType; // The left-shifted B-spline
145  if( i>=1 && i<=Degree+1 ) sPolys[i].p += baseBSpline[i-1]; // The centered B-Spline
146  if( i>=2 && i<=Degree+2 ) sPolys[i].p += baseBSpline[i-2].shift( 1 ) * boundaryType; // The right-shifted B-spline
147  for( int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
148  }
149  leftRightBaseFunction.set( sPolys , Degree+4 );
150 
151  dLeftBaseFunction = leftBaseFunction.derivative();
152  dRightBaseFunction = rightBaseFunction.derivative();
153  dLeftRightBaseFunction = leftRightBaseFunction.derivative();
154  leftRightBSpline = leftBSpline = rightBSpline = baseBSpline;
155  leftBSpline [1] += leftBSpline[2].shift( -1 ) , leftBSpline[0] *= 0;
156  rightBSpline[1] += rightBSpline[0].shift( 1 ) , rightBSpline[2] *= 0;
157  leftRightBSpline[1] += leftRightBSpline[2].shift( -1 ) + leftRightBSpline[0].shift( 1 ) , leftRightBSpline[0] *= 0 , leftRightBSpline[2] *= 0 ;
158 
159  double c , w;
160  for( int i=0 ; i<functionCount ; i++ )
161  {
162  BinaryNode< double >::CenterAndWidth( i , c , w );
163  baseFunctions[i] = baseFunction.scale(w).shift(c);
164  baseBSplines[i] = baseBSpline.scale(w).shift(c);
165  if( boundaryType )
166  {
167  int d , off , r;
168  BinaryNode< double >::DepthAndOffset( i , d , off );
169  r = 1<<d;
170  if ( off==0 && off==r-1 ) baseFunctions[i] = leftRightBaseFunction.scale(w).shift(c);
171  else if( off==0 ) baseFunctions[i] = leftBaseFunction.scale(w).shift(c);
172  else if( off==r-1 ) baseFunctions[i] = rightBaseFunction.scale(w).shift(c);
173  if ( off==0 && off==r-1 ) baseBSplines [i] = leftRightBSpline.scale(w).shift(c);
174  else if( off==0 ) baseBSplines [i] = leftBSpline.scale(w).shift(c);
175  else if( off==r-1 ) baseBSplines [i] = rightBSpline.scale(w).shift(c);
176  }
177  }
178 }
179 template<int Degree,class Real>
180 void BSplineData<Degree,Real>::setDotTables( int flags , bool inset )
181 {
182  clearDotTables( flags );
183  int size = ( functionCount*functionCount + functionCount )>>1;
184  int fullSize = functionCount*functionCount;
185  if( flags & VV_DOT_FLAG )
186  {
187  vvDotTable = NewPointer< Real >( size );
188  memset( vvDotTable , 0 , sizeof(Real)*size );
189  }
190  if( flags & DV_DOT_FLAG )
191  {
192  dvDotTable = NewPointer< Real >( fullSize );
193  memset( dvDotTable , 0 , sizeof(Real)*fullSize );
194  }
195  if( flags & DD_DOT_FLAG )
196  {
197  ddDotTable = NewPointer< Real >( size );
198  memset( ddDotTable , 0 , sizeof(Real)*size );
199  }
200  double vvIntegrals[Degree+1][Degree+1];
201  double vdIntegrals[Degree+1][Degree ];
202  double dvIntegrals[Degree ][Degree+1];
203  double ddIntegrals[Degree ][Degree ];
204  int vvSums[Degree+1][Degree+1];
205  int vdSums[Degree+1][Degree ];
206  int dvSums[Degree ][Degree+1];
207  int ddSums[Degree ][Degree ];
208  SetBSplineElementIntegrals< Degree , Degree >( vvIntegrals );
209  SetBSplineElementIntegrals< Degree , Degree-1 >( vdIntegrals );
210  SetBSplineElementIntegrals< Degree-1 , Degree >( dvIntegrals );
211  SetBSplineElementIntegrals< Degree-1 , Degree-1 >( ddIntegrals );
212 
213  for( int d1=0 ; d1<=depth ; d1++ )
214  for( int off1=0 ; off1<(1<<d1) ; off1++ )
215  {
216  int ii = BinaryNode< Real >::CenterIndex( d1 , off1 );
217  BSplineElements< Degree > b1( 1<<d1 , off1 , boundaryType , inset ? ( 1<<(d1-2) ) : 0 );
218  BSplineElements< Degree-1 > db1;
219  b1.differentiate( db1 );
220 
221  int start1 , end1;
222 
223  start1 = -1 , end1 = -1;
224  for( int i=0 ; i<int(b1.size()) ; i++ ) for( int j=0 ; j<=Degree ; j++ )
225  {
226  if( b1[i][j] && start1==-1 ) start1 = i;
227  if( b1[i][j] ) end1 = i+1;
228  }
229  if( start1==end1 ) continue;
230  for( int d2=d1 ; d2<=depth ; d2++ )
231  {
232  for( int off2=0 ; off2<(1<<d2) ; off2++ )
233  {
234  int start2 = off2-Degree;
235  int end2 = off2+Degree+1;
236  if( start2>=end1 || start1>=end2 ) continue;
237  start2 = std::max< int >( start1 , start2 );
238  end2 = std::min< int >( end1 , end2 );
239  if( d1==d2 && off2<off1 ) continue;
240  int jj = BinaryNode< Real >::CenterIndex( d2 , off2 );
241  BSplineElements< Degree > b2( 1<<d2 , off2 , boundaryType , inset ? ( 1<<(d2-2) ) : 0 );
242  BSplineElements< Degree-1 > db2;
243  b2.differentiate( db2 );
244 
245  int idx = SymmetricIndex( ii , jj );
246  int idx1 = Index( ii , jj ) , idx2 = Index( jj , ii );
247 
248  memset( vvSums , 0 , sizeof( int ) * ( Degree+1 ) * ( Degree+1 ) );
249  memset( vdSums , 0 , sizeof( int ) * ( Degree+1 ) * ( Degree ) );
250  memset( dvSums , 0 , sizeof( int ) * ( Degree ) * ( Degree+1 ) );
251  memset( ddSums , 0 , sizeof( int ) * ( Degree ) * ( Degree ) );
252  for( int i=start2 ; i<end2 ; i++ )
253  {
254  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) vvSums[j][k] += b1[i][j] * b2[i][k];
255  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) vdSums[j][k] += b1[i][j] * db2[i][k];
256  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) dvSums[j][k] += db1[i][j] * b2[i][k];
257  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) ddSums[j][k] += db1[i][j] * db2[i][k];
258  }
259  double vvDot = 0 , dvDot = 0 , vdDot = 0 , ddDot = 0;
260  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) vvDot += vvIntegrals[j][k] * vvSums[j][k];
261  for( int j=0 ; j<=Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) vdDot += vdIntegrals[j][k] * vdSums[j][k];
262  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k<=Degree ; k++ ) dvDot += dvIntegrals[j][k] * dvSums[j][k];
263  for( int j=0 ; j< Degree ; j++ ) for( int k=0 ; k< Degree ; k++ ) ddDot += ddIntegrals[j][k] * ddSums[j][k];
264  vvDot /= (1<<d2);
265  ddDot *= (1<<d2);
266  vvDot /= ( b1.denominator * b2.denominator );
267  dvDot /= ( b1.denominator * b2.denominator );
268  vdDot /= ( b1.denominator * b2.denominator );
269  ddDot /= ( b1.denominator * b2.denominator );
270  if( fabs(vvDot)<1e-15 ) continue;
271  if( flags & VV_DOT_FLAG ) vvDotTable [idx] = Real( vvDot );
272  if( useDotRatios )
273  {
274  if( flags & DV_DOT_FLAG ) dvDotTable[idx1] = Real( dvDot / vvDot );
275  if( flags & DV_DOT_FLAG ) dvDotTable[idx2] = Real( vdDot / vvDot );
276  if( flags & DD_DOT_FLAG ) ddDotTable[idx ] = Real( ddDot / vvDot );
277  }
278  else
279  {
280  if( flags & DV_DOT_FLAG ) dvDotTable[idx1] = Real( dvDot );
281  if( flags & DV_DOT_FLAG ) dvDotTable[idx2] = Real( dvDot );
282  if( flags & DD_DOT_FLAG ) ddDotTable[idx ] = Real( ddDot );
283  }
284  }
285  BSplineElements< Degree > b;
286  b = b1;
287  b.upSample( b1 );
288  b1.differentiate( db1 );
289  start1 = -1;
290  for( int i=0 ; i<int(b1.size()) ; i++ ) for( int j=0 ; j<=Degree ; j++ )
291  {
292  if( b1[i][j] && start1==-1 ) start1 = i;
293  if( b1[i][j] ) end1 = i+1;
294  }
295  }
296  }
297 }
298 template< int Degree,class Real>
299 void BSplineData<Degree,Real>::clearDotTables( int flags )
300 {
301  if( (flags & VV_DOT_FLAG) && vvDotTable ) DeletePointer( vvDotTable );
302  if( (flags & DV_DOT_FLAG) && dvDotTable ) DeletePointer( dvDotTable );
303  if( (flags & DD_DOT_FLAG) && ddDotTable ) DeletePointer( ddDotTable );
304 }
305 template< int Degree , class Real >
306 void BSplineData< Degree , Real >::setSampleSpan( int idx , int& start , int& end , double smooth ) const
307 {
308  int d , off , res;
309  BinaryNode< double >::DepthAndOffset( idx , d , off );
310  res = 1<<d;
311  double _start = ( off + 0.5 - 0.5*(Degree+1) ) / res - smooth;
312  double _end = ( off + 0.5 + 0.5*(Degree+1) ) / res + smooth;
313  // (start)/(sampleCount-1) >_start && (start-1)/(sampleCount-1)<=_start
314  // => start > _start * (sampleCount-1 ) && start <= _start*(sampleCount-1) + 1
315  // => _start * (sampleCount-1) + 1 >= start > _start * (sampleCount-1)
316  start = int( floor( _start * (sampleCount-1) + 1 ) );
317  if( start<0 ) start = 0;
318  // (end)/(sampleCount-1)<_end && (end+1)/(sampleCount-1)>=_end
319  // => end < _end * (sampleCount-1 ) && end >= _end*(sampleCount-1) - 1
320  // => _end * (sampleCount-1) > end >= _end * (sampleCount-1) - 1
321  end = int( ceil( _end * (sampleCount-1) - 1 ) );
322  if( end>=sampleCount ) end = sampleCount-1;
323 }
324 template< int Degree,class Real>
325 void BSplineData<Degree,Real>::setValueTables( int flags , double smooth )
326 {
327  clearValueTables();
328  if( flags & VALUE_FLAG ) valueTables = NewPointer< Real >( functionCount*sampleCount );
329  if( flags & D_VALUE_FLAG ) dValueTables = NewPointer< Real >( functionCount*sampleCount );
330  PPolynomial<Degree+1> function;
331  PPolynomial<Degree> dFunction;
332  for( int i=0 ; i<functionCount ; i++ )
333  {
334  if(smooth>0)
335  {
336  function = baseFunctions[i].MovingAverage(smooth);
337  dFunction = baseFunctions[i].derivative().MovingAverage(smooth);
338  }
339  else
340  {
341  function = baseFunctions[i];
342  dFunction = baseFunctions[i].derivative();
343  }
344  for( int j=0 ; j<sampleCount ; j++ )
345  {
346  double x=double(j)/(sampleCount-1);
347  if( flags & VALUE_FLAG ) valueTables[j*functionCount+i] = Real( function(x) );
348  if( flags & D_VALUE_FLAG ) dValueTables[j*functionCount+i] = Real( dFunction(x) );
349  }
350  }
351 }
352 template< int Degree,class Real>
353 void BSplineData<Degree,Real>::setValueTables( int flags , double valueSmooth , double derivativeSmooth )
354 {
355  clearValueTables();
356  if(flags & VALUE_FLAG) valueTables = NewPointer< Real >( functionCount*sampleCount );
357  if(flags & D_VALUE_FLAG) dValueTables = NewPointer< Real >( functionCount*sampleCount );
358  PPolynomial<Degree+1> function;
359  PPolynomial<Degree> dFunction;
360  for( int i=0 ; i<functionCount ; i++ )
361  {
362  if( valueSmooth>0 ) function=baseFunctions[i].MovingAverage( valueSmooth );
363  else function=baseFunctions[i];
364  if( derivativeSmooth>0 ) dFunction=baseFunctions[i].derivative().MovingAverage( derivativeSmooth );
365  else dFunction=baseFunctions[i].derivative();
366 
367  for( int j=0 ; j<sampleCount ; j++ )
368  {
369  double x=double(j)/(sampleCount-1);
370  if( flags & VALUE_FLAG ) valueTables[j*functionCount+i] = Real( function(x));
371  if( flags & D_VALUE_FLAG ) dValueTables[j*functionCount+i] = Real(dFunction(x));
372  }
373  }
374 }
375 
376 
377 template< int Degree,class Real>
378 void BSplineData<Degree,Real>::clearValueTables(void){
379  if( valueTables ) DeletePointer( valueTables );
380  if( dValueTables ) DeletePointer( dValueTables );
381 }
382 
383 template< int Degree,class Real>
384 inline int BSplineData<Degree,Real>::Index( int i1 , int i2 ) const { return i1*functionCount+i2; }
385 template< int Degree,class Real>
386 inline int BSplineData<Degree,Real>::SymmetricIndex( int i1 , int i2 )
387 {
388  if( i1>i2 ) return ((i1*i1+i1)>>1)+i2;
389  else return ((i2*i2+i2)>>1)+i1;
390 }
391 template< int Degree,class Real>
392 inline int BSplineData<Degree,Real>::SymmetricIndex( int i1 , int i2 , int& index )
393 {
394  if( i1<i2 )
395  {
396  index = ((i2*i2+i2)>>1)+i1;
397  return 1;
398  }
399  else
400  {
401  index = ((i1*i1+i1)>>1)+i2;
402  return 0;
403  }
404 }
405 
406 
408 // BSplineElements //
410 template< int Degree >
411 BSplineElements< Degree >::BSplineElements( int res , int offset , int boundary , int inset )
412 {
413  denominator = 1;
414  this->resize( res , BSplineElementCoefficients< Degree >() );
415 
416  for( int i=0 ; i<=Degree ; i++ )
417  {
418  int idx = -_off + offset + i;
419  if( idx>=0 && idx<res ) (*this)[idx][i] = 1;
420  }
421  if( boundary!=0 )
422  {
423  _addLeft( offset-2*res , boundary ) , _addRight( offset+2*res , boundary );
424  if( Degree&1 ) _addLeft( offset-res , boundary ) , _addRight( offset+res , boundary );
425  else _addLeft( -offset-1 , boundary ) , _addRight( -offset-1+2*res , boundary );
426  }
427  if( inset ) for( int i=0 ; i<inset && i<res ; i++ ) for( int j=0 ; j<=Degree ; j++ ) (*this)[i][j] = (*this)[res-1-i][j] = 0;
428 }
429 template< int Degree >
430 void BSplineElements< Degree >::_addLeft( int offset , int boundary )
431 {
432  int res = int( this->size() );
433  bool set = false;
434  for( int i=0 ; i<=Degree ; i++ )
435  {
436  int idx = -_off + offset + i;
437  if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set = true;
438  }
439  if( set ) _addLeft( offset-2*res , boundary );
440 }
441 template< int Degree >
442 void BSplineElements< Degree >::_addRight( int offset , int boundary )
443 {
444  int res = int( this->size() );
445  bool set = false;
446  for( int i=0 ; i<=Degree ; i++ )
447  {
448  int idx = -_off + offset + i;
449  if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set = true;
450  }
451  if( set ) _addRight( offset+2*res , boundary );
452 }
453 template< int Degree >
454 void BSplineElements< Degree >::upSample( BSplineElements< Degree >& high ) const
455 {
456  fprintf( stderr , "[ERROR] B-spline up-sampling not supported for degree %d\n" , Degree );
457  exit( 0 );
458 }
459 template<>
460 void BSplineElements< 1 >::upSample( BSplineElements< 1 >& high ) const
461 {
462  high.resize( size()*2 );
463  high.assign( high.size() , BSplineElementCoefficients<1>() );
464  for( int i=0 ; i<int(size()) ; i++ )
465  {
466  high[2*i+0][0] += 1 * (*this)[i][0];
467  high[2*i+0][1] += 0 * (*this)[i][0];
468  high[2*i+1][0] += 2 * (*this)[i][0];
469  high[2*i+1][1] += 1 * (*this)[i][0];
470 
471  high[2*i+0][0] += 1 * (*this)[i][1];
472  high[2*i+0][1] += 2 * (*this)[i][1];
473  high[2*i+1][0] += 0 * (*this)[i][1];
474  high[2*i+1][1] += 1 * (*this)[i][1];
475  }
476  high.denominator = denominator * 2;
477 }
478 template<>
479 void BSplineElements< 2 >::upSample( BSplineElements< 2 >& high ) const
480 {
481  /* /----\
482  / \
483  / \ = 1 /--\ +3 /--\ +3 /--\ +1 /--\
484  / \ / \ / \ / \ / \
485  |----------| |----------| |----------| |----------| |----------|
486  */
487 
488  high.resize( size()*2 );
489  high.assign( high.size() , BSplineElementCoefficients<2>() );
490  for( int i=0 ; i<int(size()) ; i++ )
491  {
492  high[2*i+0][0] += 1 * (*this)[i][0];
493  high[2*i+0][1] += 0 * (*this)[i][0];
494  high[2*i+0][2] += 0 * (*this)[i][0];
495  high[2*i+1][0] += 3 * (*this)[i][0];
496  high[2*i+1][1] += 1 * (*this)[i][0];
497  high[2*i+1][2] += 0 * (*this)[i][0];
498 
499  high[2*i+0][0] += 3 * (*this)[i][1];
500  high[2*i+0][1] += 3 * (*this)[i][1];
501  high[2*i+0][2] += 1 * (*this)[i][1];
502  high[2*i+1][0] += 1 * (*this)[i][1];
503  high[2*i+1][1] += 3 * (*this)[i][1];
504  high[2*i+1][2] += 3 * (*this)[i][1];
505 
506  high[2*i+0][0] += 0 * (*this)[i][2];
507  high[2*i+0][1] += 1 * (*this)[i][2];
508  high[2*i+0][2] += 3 * (*this)[i][2];
509  high[2*i+1][0] += 0 * (*this)[i][2];
510  high[2*i+1][1] += 0 * (*this)[i][2];
511  high[2*i+1][2] += 1 * (*this)[i][2];
512  }
513  high.denominator = denominator * 4;
514 }
515 
516 template< int Degree >
517 void BSplineElements< Degree >::differentiate( BSplineElements< Degree-1 >& d ) const
518 {
519  d.resize( this->size() );
520  d.assign( d.size() , BSplineElementCoefficients< Degree-1 >() );
521  for(int i=0 ; i<int(this->size()) ; i++ ) for(int j=0 ; j<=Degree ; j++ )
522  {
523  if( j>0 ) d[i][j-1] -= (*this)[i][j];
524  if( j<Degree ) d[i][j ] += (*this)[i][j];
525  }
526  d.denominator = denominator;
527 }
528 // If we were really good, we would implement this integral table to store
529 // rational values to improve precision...
530 template< int Degree1 , int Degree2 >
531 void SetBSplineElementIntegrals( double integrals[Degree1+1][Degree2+1] )
532 {
533  for(int i=0 ; i<=Degree1 ; i++ )
534  {
535  Polynomial< Degree1 > p1 = Polynomial< Degree1 >::BSplineComponent( i );
536  for(int j=0 ; j<=Degree2 ; j++ )
537  {
538  Polynomial< Degree2 > p2 = Polynomial< Degree2 >::BSplineComponent( j );
539  integrals[i][j] = ( p1 * p2 ).integral( 0 , 1 );
540  }
541  }
542 }