57 template<
int Degree >
inline bool LeftOverlap(
unsigned int depth ,
int offset )
60 if( Degree & 1 )
return (offset < 1+Degree) && (offset > -1-Degree );
61 else return (offset < Degree) && (offset > -2-Degree );
63 template<
int Degree >
inline bool RightOverlap(
unsigned int depth ,
int offset )
67 if( Degree & 1 )
return (offset > 2-1-Degree) && (offset < 2+1+Degree );
68 else return (offset > 2-2-Degree) && (offset < 2+ Degree );
70 template<
int Degree >
inline int ReflectLeft(
unsigned int depth ,
int offset )
72 if( Degree&1 )
return -offset;
73 else return -1-offset;
75 template<
int Degree >
inline int ReflectRight(
unsigned int depth ,
int offset )
78 if( Degree&1 )
return r -offset;
79 else return r-1-offset;
82 template<
int Degree ,
class Real >
85 vvDotTable = dvDotTable = ddDotTable = NULL;
86 valueTables = dValueTables = NULL;
87 functionCount = sampleCount = 0;
90 template<
int Degree ,
class Real >
95 if( vvDotTable )
delete[] vvDotTable;
96 if( dvDotTable )
delete[] dvDotTable;
97 if( ddDotTable )
delete[] ddDotTable;
99 if( valueTables )
delete[] valueTables;
100 if( dValueTables )
delete[] dValueTables;
102 vvDotTable = dvDotTable = ddDotTable = NULL;
103 valueTables = dValueTables=NULL;
107 template<
int Degree,
class Real>
110 this->useDotRatios = useDotRatios;
111 this->reflectBoundary = reflectBoundary;
121 for(
int i=0 ; i<=Degree ; i++ ) baseBSpline[i] = Polynomial< Degree >::BSplineComponent( i ).shift(
double(-(Degree+1)/2) + i - 0.5 );
122 dBaseFunction = baseFunction.derivative();
125 for(
int i=0 ; i<Degree+3 ; i++ )
127 sPolys[i].
start = double(-(Degree+1)/2) + i - 1.5;
129 if( i<=Degree ) sPolys[i].
p += baseBSpline[i ].shift( -1 );
130 if( i>=1 && i<=Degree+1 ) sPolys[i].
p += baseBSpline[i-1];
131 for(
int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
133 leftBaseFunction.set( sPolys , Degree+3 );
134 for(
int i=0 ; i<Degree+3 ; i++ )
136 sPolys[i].
start = double(-(Degree+1)/2) + i - 0.5;
138 if( i<=Degree ) sPolys[i].
p += baseBSpline[i ];
139 if( i>=1 && i<=Degree+1 ) sPolys[i].
p += baseBSpline[i-1].shift( 1 );
140 for(
int j=0 ; j<i ; j++ ) sPolys[i].p -= sPolys[j].p;
142 rightBaseFunction.set( sPolys , Degree+3 );
143 dLeftBaseFunction = leftBaseFunction.derivative();
144 dRightBaseFunction = rightBaseFunction.derivative();
145 leftBSpline = rightBSpline = baseBSpline;
146 leftBSpline [1] += leftBSpline[2].
shift( -1 ) , leftBSpline[0] *= 0;
147 rightBSpline[1] += rightBSpline[0].
shift( 1 ) , rightBSpline[2] *= 0;
149 for(
int i=0 ; i<functionCount ; i++ )
152 baseFunctions[i] = baseFunction.scale(w).shift(c);
153 baseBSplines[i] = baseBSpline.scale(w).shift(c);
154 if( reflectBoundary )
159 if ( off==0 ) baseFunctions[i] = leftBaseFunction.scale(w).shift(c);
160 else if( off==r-1 ) baseFunctions[i] = rightBaseFunction.scale(w).shift(c);
161 if ( off==0 ) baseBSplines[i] = leftBSpline.scale(w).shift(c);
162 else if( off==r-1 ) baseBSplines[i] = rightBSpline.scale(w).shift(c);
166 template<
int Degree,
class Real>
169 clearDotTables( flags );
170 int size = ( functionCount*functionCount + functionCount )>>1;
171 int fullSize = functionCount*functionCount;
172 if( flags & VV_DOT_FLAG )
174 vvDotTable =
new Real[size];
175 memset( vvDotTable , 0 ,
sizeof(
Real)*size );
177 if( flags & DV_DOT_FLAG )
179 dvDotTable =
new Real[fullSize];
180 memset( dvDotTable , 0 ,
sizeof(
Real)*fullSize );
182 if( flags & DD_DOT_FLAG )
184 ddDotTable =
new Real[size];
185 memset( ddDotTable , 0 ,
sizeof(
Real)*size );
187 double vvIntegrals[Degree+1][Degree+1];
188 double vdIntegrals[Degree+1][Degree ];
189 double dvIntegrals[Degree ][Degree+1];
190 double ddIntegrals[Degree ][Degree ];
191 int vvSums[Degree+1][Degree+1];
192 int vdSums[Degree+1][Degree ];
193 int dvSums[Degree ][Degree+1];
194 int ddSums[Degree ][Degree ];
195 SetBSplineElementIntegrals< Degree , Degree >( vvIntegrals );
200 for(
int d1=0 ; d1<=depth ; d1++ )
201 for(
int off1=0 ; off1<(1<<d1) ; off1++ )
211 for(
int i=0 ; i<int(b1.size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
213 if( b1[i][j] && start1==-1 ) start1 = i;
214 if( b1[i][j] ) end1 = i+1;
216 for(
int d2=d1 ; d2<=depth ; d2++ )
218 for(
int off2=0 ; off2<(1<<d2) ; off2++ )
220 int start2 = off2-Degree;
221 int end2 = off2+Degree+1;
222 if( start2>=end1 || start1>=end2 )
continue;
223 start2 = std::max< int >( start1 , start2 );
224 end2 = std::min< int >( end1 , end2 );
225 if( d1==d2 && off2<off1 )
continue;
231 int idx = SymmetricIndex( ii , jj );
232 int idx1 = Index( ii , jj ) , idx2 = Index( jj , ii );
234 memset( vvSums , 0 ,
sizeof(
int ) * ( Degree+1 ) * ( Degree+1 ) );
235 memset( vdSums , 0 ,
sizeof(
int ) * ( Degree+1 ) * ( Degree ) );
236 memset( dvSums , 0 ,
sizeof(
int ) * ( Degree ) * ( Degree+1 ) );
237 memset( ddSums , 0 ,
sizeof(
int ) * ( Degree ) * ( Degree ) );
238 for(
int i=start2 ; i<end2 ; i++ )
240 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) vvSums[j][k] += b1[i][j] * b2[i][k];
241 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) vdSums[j][k] += b1[i][j] * db2[i][k];
242 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) dvSums[j][k] += db1[i][j] * b2[i][k];
243 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) ddSums[j][k] += db1[i][j] * db2[i][k];
245 double vvDot = 0 , dvDot = 0 , vdDot = 0 , ddDot = 0;
246 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) vvDot += vvIntegrals[j][k] * vvSums[j][k];
247 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) vdDot += vdIntegrals[j][k] * vdSums[j][k];
248 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) dvDot += dvIntegrals[j][k] * dvSums[j][k];
249 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) ddDot += ddIntegrals[j][k] * ddSums[j][k];
256 if( fabs(vvDot)<1e-15 )
continue;
257 if( flags & VV_DOT_FLAG ) vvDotTable [idx] =
Real( vvDot );
260 if( flags & DV_DOT_FLAG ) dvDotTable[idx1] =
Real( dvDot / vvDot );
261 if( flags & DV_DOT_FLAG ) dvDotTable[idx2] =
Real( vdDot / vvDot );
262 if( flags & DD_DOT_FLAG ) ddDotTable[idx ] =
Real( ddDot / vvDot );
266 if( flags & DV_DOT_FLAG ) dvDotTable[idx1] =
Real( dvDot );
267 if( flags & DV_DOT_FLAG ) dvDotTable[idx2] =
Real( dvDot );
268 if( flags & DD_DOT_FLAG ) ddDotTable[idx ] =
Real( ddDot );
276 for(
int i=0 ; i<int(b1.size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
278 if( b1[i][j] && start1==-1 ) start1 = i;
279 if( b1[i][j] ) end1 = i+1;
284 template<
int Degree,
class Real>
287 if( (flags & VV_DOT_FLAG) && vvDotTable )
delete[] vvDotTable , vvDotTable = NULL;
288 if( (flags & DV_DOT_FLAG) && dvDotTable )
delete[] dvDotTable , dvDotTable = NULL;
289 if( (flags & DD_DOT_FLAG) && ddDotTable )
delete[] ddDotTable , ddDotTable = NULL;
291 template<
int Degree ,
class Real >
297 double _start = ( off + 0.5 - 0.5*(Degree+1) ) / res - smooth;
298 double _end = ( off + 0.5 + 0.5*(Degree+1) ) / res + smooth;
302 start = int( floor( _start * (sampleCount-1) + 1 ) );
303 if( start<0 ) start = 0;
307 end = int( ceil( _end * (sampleCount-1) - 1 ) );
308 if( end>=sampleCount ) end = sampleCount-1;
310 template<
int Degree,
class Real>
314 if( flags & VALUE_FLAG ) valueTables =
new Real[functionCount*sampleCount];
315 if( flags & D_VALUE_FLAG ) dValueTables =
new Real[functionCount*sampleCount];
318 for(
int i=0 ; i<functionCount ; i++ )
327 function = baseFunctions[i];
330 for(
int j=0 ; j<sampleCount ; j++ )
332 double x=double(j)/(sampleCount-1);
333 if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=
Real(
function(x));}
334 if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=
Real(dFunction(x));}
338 template<
int Degree,
class Real>
341 if(flags & VALUE_FLAG){ valueTables=
new Real[functionCount*sampleCount];}
342 if(flags & D_VALUE_FLAG){dValueTables=
new Real[functionCount*sampleCount];}
345 for(
int i=0;i<functionCount;i++){
346 if(valueSmooth>0) {
function=baseFunctions[i].
MovingAverage(valueSmooth);}
347 else {
function=baseFunctions[i];}
349 else {dFunction=baseFunctions[i].
derivative();}
351 for(
int j=0;j<sampleCount;j++){
352 double x=double(j)/(sampleCount-1);
353 if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=
Real(
function(x));}
354 if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=
Real(dFunction(x));}
360 template<
int Degree,
class Real>
362 if( valueTables){
delete[] valueTables;}
363 if(dValueTables){
delete[] dValueTables;}
364 valueTables=dValueTables=NULL;
367 template<
int Degree,
class Real>
369 template<
int Degree,
class Real>
372 if( i1>i2 )
return ((i1*i1+i1)>>1)+i2;
373 else return ((i2*i2+i2)>>1)+i1;
375 template<
int Degree,
class Real>
380 index = ((i2*i2+i2)>>1)+i1;
385 index = ((i1*i1+i1)>>1)+i2;
394 template<
int Degree >
400 for(
int i=0 ; i<=Degree ; i++ )
402 int idx = -_off + offset + i;
403 if( idx>=0 && idx<res ) (*this)[idx][i] = 1;
407 _addLeft( offset-2*res , boundary ) , _addRight( offset+2*res , boundary );
408 if( Degree&1 ) _addLeft( offset-res , boundary ) , _addRight( offset+res , boundary );
409 else _addLeft( -offset-1 , boundary ) , _addRight( -offset-1+2*res , boundary );
412 template<
int Degree >
415 int res = int( this->size() );
417 for(
int i=0 ; i<=Degree ; i++ )
419 int idx = -_off + offset + i;
420 if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set =
true;
422 if( set ) _addLeft( offset-2*res , boundary );
424 template<
int Degree >
427 int res = int( this->size() );
429 for(
int i=0 ; i<=Degree ; i++ )
431 int idx = -_off + offset + i;
432 if( idx>=0 && idx<res ) (*this)[idx][i] += boundary , set =
true;
434 if( set ) _addRight( offset+2*res , boundary );
436 template<
int Degree >
439 fprintf( stderr ,
"[ERROR] B-spline up-sampling not supported for degree %d\n" , Degree );
445 high.resize( size()*2 );
447 for(
int i=0 ; i<int(size()) ; i++ )
449 high[2*i+0][0] += 1 * (*this)[i][0];
450 high[2*i+0][1] += 0 * (*this)[i][0];
451 high[2*i+1][0] += 2 * (*this)[i][0];
452 high[2*i+1][1] += 1 * (*this)[i][0];
454 high[2*i+0][0] += 1 * (*this)[i][1];
455 high[2*i+0][1] += 2 * (*this)[i][1];
456 high[2*i+1][0] += 0 * (*this)[i][1];
457 high[2*i+1][1] += 1 * (*this)[i][1];
470 high.resize( size()*2 );
472 for(
int i=0 ; i<int(size()) ; i++ )
474 high[2*i+0][0] += 1 * (*this)[i][0];
475 high[2*i+0][1] += 0 * (*this)[i][0];
476 high[2*i+0][2] += 0 * (*this)[i][0];
477 high[2*i+1][0] += 3 * (*this)[i][0];
478 high[2*i+1][1] += 1 * (*this)[i][0];
479 high[2*i+1][2] += 0 * (*this)[i][0];
481 high[2*i+0][0] += 3 * (*this)[i][1];
482 high[2*i+0][1] += 3 * (*this)[i][1];
483 high[2*i+0][2] += 1 * (*this)[i][1];
484 high[2*i+1][0] += 1 * (*this)[i][1];
485 high[2*i+1][1] += 3 * (*this)[i][1];
486 high[2*i+1][2] += 3 * (*this)[i][1];
488 high[2*i+0][0] += 0 * (*this)[i][2];
489 high[2*i+0][1] += 1 * (*this)[i][2];
490 high[2*i+0][2] += 3 * (*this)[i][2];
491 high[2*i+1][0] += 0 * (*this)[i][2];
492 high[2*i+1][1] += 0 * (*this)[i][2];
493 high[2*i+1][2] += 1 * (*this)[i][2];
498 template<
int Degree >
501 d.resize( this->size() );
503 for(
int i=0 ; i<int(this->size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
505 if( j-1>=0 ) d[i][j-1] -= (*this)[i][j];
506 if( j<Degree ) d[i][j ] += (*this)[i][j];
512 template<
int Degree1 ,
int Degree2 >
515 for(
int i=0 ; i<=Degree1 ; i++ )
518 for(
int j=0 ; j<=Degree2 ; j++ )
521 integrals[i][j] = ( p1 * p2 ).integral( 0 , 1 );
virtual void clearDotTables(int flags)
StartingPolynomial shift(double t) const
bool RightOverlap(unsigned int depth, int offset)
static int CenterIndex(int depth, int offSet)
static int SymmetricIndex(int i1, int i2)
virtual void clearValueTables(void)
void _addLeft(int offset, int boundary)
static int CenterCount(int depth)
static int CumulativeCenterCount(int maxDepth)
PPolynomial< Degree-1 > derivative(void) const
bool LeftOverlap(unsigned int depth, int offset)
virtual void setDotTables(int flags)
void upSample(BSplineElements &high) const
void differentiate(BSplineElements< Degree-1 > &d) const
void _addRight(int offset, int boundary)
int Index(int i1, int i2) const
static void DepthAndOffset(int idx, int &depth, int &offset)
void set(int maxDepth, bool useDotRatios=true, bool reflectBoundary=false)
static Polynomial BSplineComponent(int i)
void setSampleSpan(int idx, int &start, int &end, double smooth=0) const
void SetBSplineElementIntegrals(double integrals[Degree1+1][Degree2+1])
virtual void setValueTables(int flags, double smooth=0)
int ReflectRight(unsigned int depth, int offset)
int ReflectLeft(unsigned int depth, int offset)
static PPolynomial BSpline(double radius=0.5)
static int CornerCount(int depth)
PPolynomial< Degree+1 > MovingAverage(double radius)
static void CenterAndWidth(int depth, int offset, Real ¢er, Real &width)