29 #include "poisson_exceptions.h" 30 #include "binary_node.h" 59 template<
int Degree >
inline bool LeftOverlap(
unsigned int,
int offset )
62 if( Degree & 1 )
return (offset < 1+Degree) && (offset > -1-Degree );
63 else return (offset < Degree) && (offset > -2-Degree );
65 template<
int Degree >
inline bool RightOverlap(
unsigned int,
int offset )
68 if( Degree & 1 )
return (offset > 2-1-Degree) && (offset < 2+1+Degree );
69 else return (offset > 2-2-Degree) && (offset < 2+ Degree );
71 template<
int Degree >
inline int ReflectLeft(
unsigned int,
int offset )
73 if( Degree&1 )
return -offset;
74 else return -1-offset;
76 template<
int Degree >
inline int ReflectRight(
unsigned int depth ,
int offset )
79 if( Degree&1 )
return r -offset;
80 else return r-1-offset;
83 template<
int Degree ,
class Real >
86 vvDotTable = dvDotTable = ddDotTable = NULL;
87 valueTables = dValueTables = NULL;
90 functionCount = sampleCount = 0;
93 template<
int Degree ,
class Real >
102 delete[] valueTables;
103 delete[] dValueTables;
105 delete[] baseFunctions;
106 delete[] baseBSplines;
108 vvDotTable = dvDotTable = ddDotTable = NULL;
109 valueTables = dValueTables=NULL;
110 baseFunctions = NULL;
115 template<
int Degree,
class Real>
118 this->useDotRatios = useDotRatios;
119 this->reflectBoundary = reflectBoundary;
126 baseBSplines =
new BSplineComponents[functionCount];
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();
133 for(
int i=0 ; i<Degree+3 ; i++ )
135 sPolys[i].
start = double(-(Degree+1)/2) + i - 1.5;
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;
141 leftBaseFunction.set( sPolys , Degree+3 );
142 for(
int i=0 ; i<Degree+3 ; i++ )
144 sPolys[i].
start = double(-(Degree+1)/2) + i - 0.5;
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;
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;
157 for(
int i=0 ; i<functionCount ; i++ )
160 baseFunctions[i] = baseFunction.scale(w).shift(c);
161 baseBSplines[i] = baseBSpline.scale(w).shift(c);
162 if( reflectBoundary )
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);
174 template<
int Degree,
class Real>
177 clearDotTables( flags );
178 int size = ( functionCount*functionCount + functionCount )>>1;
179 int fullSize = functionCount*functionCount;
180 if( flags & VV_DOT_FLAG )
182 vvDotTable =
new Real[size];
183 memset( vvDotTable , 0 ,
sizeof(
Real)*size );
185 if( flags & DV_DOT_FLAG )
187 dvDotTable =
new Real[fullSize];
188 memset( dvDotTable , 0 ,
sizeof(
Real)*fullSize );
190 if( flags & DD_DOT_FLAG )
192 ddDotTable =
new Real[size];
193 memset( ddDotTable , 0 ,
sizeof(
Real)*size );
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 );
208 for(
int d1=0 ; d1<=depth ; d1++ )
209 for(
int off1=0 ; off1<(1<<d1) ; off1++ )
219 for(
int i=0 ; i<int(b1.size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
221 if( b1[i][j] && start1==-1 ) start1 = i;
222 if( b1[i][j] ) end1 = i+1;
224 for(
int d2=d1 ; d2<=depth ; d2++ )
226 for(
int off2=0 ; off2<(1<<d2) ; off2++ )
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;
239 int idx = SymmetricIndex( ii , jj );
240 int idx1 = Index( ii , jj ) , idx2 = Index( jj , ii );
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++ )
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];
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];
264 if( fabs(vvDot)<1e-15 )
continue;
265 if( flags & VV_DOT_FLAG ) vvDotTable [idx] =
Real( vvDot );
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 );
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 );
284 for(
int i=0 ; i<int(b1.size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
286 if( b1[i][j] && start1==-1 ) start1 = i;
287 if( b1[i][j] ) end1 = i+1;
292 template<
int Degree,
class Real>
295 if (flags & VV_DOT_FLAG) {
296 delete[] vvDotTable ; vvDotTable = NULL;
297 delete[] dvDotTable ; dvDotTable = NULL;
298 delete[] ddDotTable ; ddDotTable = NULL;
301 template<
int Degree ,
class Real >
307 double _start = ( off + 0.5 - 0.5*(Degree+1) ) / res - smooth;
308 double _end = ( off + 0.5 + 0.5*(Degree+1) ) / res + smooth;
312 start = int( floor( _start * (sampleCount-1) + 1 ) );
313 if( start<0 ) start = 0;
317 end = int( ceil( _end * (sampleCount-1) - 1 ) );
318 if( end>=sampleCount ) end = sampleCount-1;
320 template<
int Degree,
class Real>
324 if( flags & VALUE_FLAG ) valueTables =
new Real[functionCount*sampleCount];
325 if( flags & D_VALUE_FLAG ) dValueTables =
new Real[functionCount*sampleCount];
328 for(
int i=0 ; i<functionCount ; i++ )
337 function = baseFunctions[i];
340 for(
int j=0 ; j<sampleCount ; j++ )
342 double x=double(j)/(sampleCount-1);
343 if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=
Real(
function(x));}
344 if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=
Real(dFunction(x));}
348 template<
int Degree,
class Real>
351 if(flags & VALUE_FLAG){ valueTables=
new Real[functionCount*sampleCount];}
352 if(flags & D_VALUE_FLAG){dValueTables=
new Real[functionCount*sampleCount];}
355 for(
int i=0;i<functionCount;i++){
356 if(valueSmooth>0) {
function=baseFunctions[i].
MovingAverage(valueSmooth);}
357 else {
function=baseFunctions[i];}
359 else {dFunction=baseFunctions[i].
derivative();}
361 for(
int j=0;j<sampleCount;j++){
362 double x=double(j)/(sampleCount-1);
363 if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=
Real(
function(x));}
364 if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=
Real(dFunction(x));}
370 template<
int Degree,
class Real>
372 delete[] valueTables;
373 delete[] dValueTables;
374 valueTables=dValueTables=NULL;
377 template<
int Degree,
class Real>
379 template<
int Degree,
class Real>
382 if( i1>i2 )
return ((i1*i1+i1)>>1)+i2;
383 else return ((i2*i2+i2)>>1)+i1;
385 template<
int Degree,
class Real>
390 index = ((i2*i2+i2)>>1)+i1;
395 index = ((i1*i1+i1)>>1)+i2;
404 template<
int Degree >
410 for(
int i=0 ; i<=Degree ; i++ )
412 int idx = -_off + offset + i;
413 if( idx>=0 && idx<res ) (*this)[idx][i] = 1;
417 _addLeft( offset-2*res , boundary ) , _addRight( offset+2*res , boundary );
418 if( Degree&1 ) _addLeft( offset-res , boundary ) , _addRight( offset+res , boundary );
419 else _addLeft( -offset-1 , boundary ) , _addRight( -offset-1+2*res , boundary );
422 template<
int Degree >
425 int res = int( this->size() );
427 for(
int i=0 ; i<=Degree ; i++ )
429 int idx = -_off + offset + i;
430 if( idx>=0 && idx<res ) (*this)[idx][i] += boundary ,
set =
true;
432 if(
set ) _addLeft( offset-2*res , boundary );
434 template<
int Degree >
437 int res = int( this->size() );
439 for(
int i=0 ; i<=Degree ; i++ )
441 int idx = -_off + offset + i;
442 if( idx>=0 && idx<res ) (*this)[idx][i] += boundary ,
set =
true;
444 if(
set ) _addRight( offset+2*res , boundary );
446 template<
int Degree >
457 template<
int Degree >
460 d.resize( this->size() );
462 for(
int i=0 ; i<int(this->size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
464 if( j-1>=0 ) d[i][j-1] -= (*this)[i][j];
465 if( j<Degree ) d[i][j ] += (*this)[i][j];
471 template<
int Degree1 ,
int Degree2 >
474 for(
int i=0 ; i<=Degree1 ; i++ )
477 for(
int j=0 ; j<=Degree2 ; j++ )
480 integrals[i][j] = ( p1 * p2 ).integral( 0 , 1 );
PPolynomial< Degree-1 > derivative(void) const
virtual void clearDotTables(int flags)
static PPolynomial BSpline(double radius=0.5)
StartingPolynomial shift(double t) const
static int SymmetricIndex(int i1, int i2)
PPolynomial< Degree+1 > MovingAverage(double radius)
static int CumulativeCenterCount(int maxDepth)
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)
void _addLeft(int offset, int boundary)
bool RightOverlap(unsigned int, int offset)
void upSample(BSplineElements &high) const
virtual void setValueTables(int flags, double smooth=0)
void _addRight(int offset, int boundary)
virtual void clearValueTables(void)
static void CenterAndWidth(int depth, int offset, Real ¢er, Real &width)
bool LeftOverlap(unsigned int, int offset)
int Index(int i1, int i2) const
virtual void setDotTables(int flags)
int ReflectLeft(unsigned int, 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)
void SetBSplineElementIntegrals(double integrals[Degree1+1][Degree2+1])
static int CornerCount(int depth)
static Polynomial BSplineComponent(int i)
void setSampleSpan(int idx, int &start, int &end, double smooth=0) const
static void DepthAndOffset(int idx, int &depth, int &offset)