29#include "poisson_exceptions.h"
30#include "binary_node.h"
59 template<
int Degree >
inline bool LeftOverlap(
unsigned int,
int offset )
65 template<
int Degree >
inline bool RightOverlap(
unsigned int,
int offset )
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;
130 dBaseFunction = baseFunction.derivative();
133 for(
int i=0 ; i<
Degree+3 ; i++ )
137 if( i<=
Degree )
sPolys[i].p += baseBSpline[i ].shift( -1 );
138 if( i>=1 && i<=
Degree+1 )
sPolys[i].p += baseBSpline[i-1];
142 for(
int i=0 ; i<
Degree+3 ; i++ )
147 if( i>=1 && i<=
Degree+1 )
sPolys[i].p += baseBSpline[i-1].shift( 1 );
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];
185 if( flags & DV_DOT_FLAG )
190 if( flags & DD_DOT_FLAG )
192 ddDotTable =
new Real[size];
208 for(
int d1=0 ;
d1<=depth ;
d1++ )
214 b1.differentiate(
db1 );
219 for(
int i=0 ; i<
int(
b1.
size()) ; i++ )
for(
int j=0 ; j<=
Degree ; j++ )
222 if(
b1[i][j] )
end1 = i+1;
237 b2.differentiate(
db2 );
239 int idx = SymmetricIndex(
ii ,
jj );
260 vvDot /= (
b1.denominator *
b2.denominator );
261 dvDot /= (
b1.denominator *
b2.denominator );
262 vdDot /= (
b1.denominator *
b2.denominator );
263 ddDot /= (
b1.denominator *
b2.denominator );
265 if( flags & VV_DOT_FLAG ) vvDotTable [idx] =
Real(
vvDot );
270 if( flags & DD_DOT_FLAG ) ddDotTable[idx ] =
Real(
ddDot /
vvDot );
276 if( flags & DD_DOT_FLAG ) ddDotTable[idx ] =
Real(
ddDot );
282 b1.differentiate(
db1 );
284 for(
int i=0 ; i<
int(
b1.
size()) ; i++ )
for(
int j=0 ; j<=
Degree ; j++ )
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;
313 if( start<0 ) start = 0;
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++ )
332 function = baseFunctions[i].MovingAverage(smooth);
333 dFunction = baseFunctions[i].derivative().MovingAverage(smooth);
338 dFunction = baseFunctions[i].derivative();
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++){
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;
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];
467 d.denominator = denominator;
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 );
Iterator class for point clouds with or without given indices.
std::size_t size() const
Size of the range the iterator is going through.
void set(int maxDepth, bool useDotRatios=true, bool reflectBoundary=false)
static int SymmetricIndex(int i1, int i2)
virtual void setValueTables(int flags, double smooth=0)
virtual void clearDotTables(int flags)
virtual void setDotTables(int flags)
virtual void clearValueTables(void)
void setSampleSpan(int idx, int &start, int &end, double smooth=0) const
int Index(int i1, int i2) const
static void CenterAndWidth(int depth, int offset, Real ¢er, Real &width)
static int CumulativeCenterCount(int maxDepth)
static int CenterIndex(int depth, int offSet)
static int CornerCount(int depth)
static int CenterCount(int depth)
static void DepthAndOffset(int idx, int &depth, int &offset)
static PPolynomial BSpline(double radius=0.5)
An exception that is thrown when the arguments number or type is wrong/unhandled.
static Polynomial BSplineComponent(int i)
bool RightOverlap(unsigned int, int offset)
bool LeftOverlap(unsigned int, int offset)
int ReflectLeft(unsigned int, int offset)
void SetBSplineElementIntegrals(double integrals[Degree1+1][Degree2+1])
int ReflectRight(unsigned int depth, int offset)
void differentiate(BSplineElements< Degree-1 > &d) const
void _addRight(int offset, int boundary)
void _addLeft(int offset, int boundary)
void upSample(BSplineElements &high) const