38 template<
int Degree,
class Real>
41 dotTable=dDotTable=d2DotTable=
NULL;
42 valueTables=dValueTables=
NULL;
46 template<
int Degree,
class Real>
55 delete[] dValueTables;
57 dotTable=dDotTable=d2DotTable=
NULL;
58 valueTables=dValueTables=
NULL;
62 template<
int Degree,
class Real>
63#if BOUNDARY_CONDITIONS
69 this->normalize = normalize;
70 this->useDotRatios = useDotRatios;
71#if BOUNDARY_CONDITIONS
72 this->reflectBoundary = reflectBoundary;
77 res2 = (1<<(depth+1))+1;
86 baseFunction=
F/
sqrt((
F*
F).integral(
F.polys[0].start,
F.polys[
F.polyCount-1].start));
89 baseFunction=
F/
F.integral(
F.polys[0].start,
F.polys[
F.polyCount-1].start);
94 dBaseFunction = baseFunction.derivative();
95#if BOUNDARY_CONDITIONS
96 leftBaseFunction = baseFunction + baseFunction.shift( -1 );
97 rightBaseFunction = baseFunction + baseFunction.shift( 1 );
98 dLeftBaseFunction = leftBaseFunction.derivative();
99 dRightBaseFunction = rightBaseFunction.derivative();
102 for(
int i=0 ; i<res ; i++ )
105#if BOUNDARY_CONDITIONS
106 if( reflectBoundary )
110 if ( off==0 ) baseFunctions[i] = leftBaseFunction.scale(
w1 ).shift( c1 );
111 else if( off==((1<<d)-1) ) baseFunctions[i] = rightBaseFunction.scale(
w1 ).shift( c1 );
112 else baseFunctions[i] = baseFunction.scale(
w1 ).shift( c1 );
114 else baseFunctions[i] = baseFunction.scale(
w1).shift(c1);
116 baseFunctions[i] = baseFunction.scale(
w1).shift(c1);
122 baseFunctions[i]/=
sqrt(
w1);
125 baseFunctions[i]/=
w1;
130 template<
int Degree,
class Real>
133 clearDotTables( flags );
135 size = ( res*res + res )>>1;
136 if( flags & DOT_FLAG )
138 dotTable =
new Real[size];
141 if( flags & D_DOT_FLAG )
143 dDotTable =
new Real[size];
146 if( flags & D2_DOT_FLAG )
148 d2DotTable =
new Real[size];
152 t1 = baseFunction.polys[0].start;
153 t2 = baseFunction.polys[baseFunction.polyCount-1].start;
154 for(
int i=0 ; i<res ; i++ )
156 double c1 , c2 ,
w1 ,
w2;
158#if BOUNDARY_CONDITIONS
167 for(
int j=0 ; j<=i ; j++ )
170#if BOUNDARY_CONDITIONS
176 int idx = SymmetricIndex( i , j );
178 double start =
t1 *
w2 + c2;
179 double end =
t2 *
w2 + c2;
180#if BOUNDARY_CONDITIONS
181 if( reflectBoundary )
183 if( start<0 ) start = 0;
184 if( start>1 ) start = 1;
185 if( end <0 ) end = 0;
186 if( end >1 ) end = 1;
192 if( start>= end )
continue;
194#if BOUNDARY_CONDITIONS
197 Real dot = dotProduct( c1 ,
w1 , c2 ,
w2 );
199 if(
fabs(dot)<1
e-15 )
continue;
200 if( flags & DOT_FLAG ) dotTable[idx]=dot;
203#if BOUNDARY_CONDITIONS
204 if( flags & D_DOT_FLAG ) dDotTable[idx] = -dDotProduct( c1 ,
w1 , c2 ,
w2 ,
boundary1 ,
boundary2 ) / dot;
205 if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct( c1 ,
w1 , c2 ,
w2 ,
boundary1 ,
boundary2 ) / dot;
207 if( flags & D_DOT_FLAG ) dDotTable[idx] = -dDotProduct(c1,
w1,c2,
w2)/dot;
208 if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct(c1,
w1,c2,
w2)/dot;
213#if BOUNDARY_CONDITIONS
217 if( flags & D_DOT_FLAG ) dDotTable[idx] = dDotProduct(c1,
w1,c2,
w2);
218 if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct(c1,
w1,c2,
w2);
224 template<
int Degree,
class Real>
227 if (flags & DOT_FLAG)
237 template<
int Degree,
class Real>
241 if( flags & VALUE_FLAG ) valueTables =
new Real[res*res2];
242 if( flags & D_VALUE_FLAG ) dValueTables =
new Real[res*res2];
245 for(
int i=0 ; i<res ; i++ )
249 function=baseFunctions[i].MovingAverage(smooth);
250 dFunction=baseFunctions[i].derivative().MovingAverage(smooth);
257 for(
int j=0 ; j<res2 ; j++ )
259 double x=
double(j)/(res2-1);
260 if(flags & VALUE_FLAG){ valueTables[j*res+i]=
Real(
function(x));}
261 if(flags & D_VALUE_FLAG){dValueTables[j*res+i]=
Real(
dFunction(x));}
265 template<
int Degree,
class Real>
268 if(flags & VALUE_FLAG){ valueTables=
new Real[res*res2];}
269 if(flags & D_VALUE_FLAG){dValueTables=
new Real[res*res2];}
272 for(
int i=0;i<res;i++){
276 else {
dFunction=baseFunctions[i].derivative();}
278 for(
int j=0;j<res2;j++){
279 double x=
double(j)/(res2-1);
280 if(flags & VALUE_FLAG){ valueTables[j*res+i]=
Real(
function(x));}
281 if(flags & D_VALUE_FLAG){dValueTables[j*res+i]=
Real(
dFunction(x));}
287 template<
int Degree,
class Real>
289 delete[] valueTables;
290 delete[] dValueTables;
291 valueTables=dValueTables=
NULL;
294#if BOUNDARY_CONDITIONS
295 template<
int Degree,
class Real>
305 double r=
fabs( baseFunction.polys[0].start );
316 template<
int Degree,
class Real>
319 const PPolynomial< Degree-1 > *b1;
320 const PPolynomial< Degree > *b2;
321 if ( boundary1==-1 ) b1 = & dLeftBaseFunction;
322 else if( boundary1== 0 ) b1 = & dBaseFunction;
323 else if( boundary1== 1 ) b1 = &dRightBaseFunction;
324 if ( boundary2==-1 ) b2 = & leftBaseFunction;
325 else if( boundary2== 0 ) b2 = & baseFunction;
326 else if( boundary2== 1 ) b2 = & rightBaseFunction;
327 double r=fabs(baseFunction.polys[0].start);
330 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/sqrt(width1*width2));
332 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/(width1*width2));
334 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r));
337 template<
int Degree,
class Real>
340 const PPolynomial< Degree-1 > *b1 , *b2;
341 if ( boundary1==-1 ) b1 = & dLeftBaseFunction;
342 else if( boundary1== 0 ) b1 = & dBaseFunction;
343 else if( boundary1== 1 ) b1 = &dRightBaseFunction;
344 if ( boundary2==-1 ) b2 = & dLeftBaseFunction;
345 else if( boundary2== 0 ) b2 = & dBaseFunction;
346 else if( boundary2== 1 ) b2 = &dRightBaseFunction;
347 double r=fabs(baseFunction.polys[0].start);
351 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/sqrt(width1*width2));
353 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/(width1*width2));
355 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2);
359 template<
int Degree,
class Real>
361 double r=
fabs(baseFunction.polys[0].start);
372 template<
int Degree,
class Real>
374 double r=
fabs(baseFunction.polys[0].start);
384 template<
int Degree,
class Real>
386 double r=
fabs(baseFunction.polys[0].start);
397 template<
int Degree,
class Real>
400 if( i1>i2 )
return ((i1*i1+i1)>>1)+i2;
401 else return ((i2*i2+i2)>>1)+i1;
403 template<
int Degree,
class Real>
408 index = ((i2*i2+i2)>>1)+i1;
412 index = ((i1*i1+i1)>>1)+i2;
Iterator class for point clouds with or without given indices.
static void CenterAndWidth(int depth, int offset, Real ¢er, Real &width)
static int CumulativeCenterCount(int maxDepth)
static void DepthAndOffset(int idx, int &depth, int &offset)
virtual void setDotTables(const int &flags)
Real dotProduct(const double ¢er1, const double &width1, const double ¢er2, const double &width2) const
Real d2DotProduct(const double ¢er1, const double &width1, const double ¢er2, const double &width2) const
virtual void clearValueTables(void)
static int SymmetricIndex(const int &i1, const int &i2)
virtual void setValueTables(const int &flags, const double &smooth=0)
void set(const int &maxDepth, const PPolynomial< Degree > &F, const int &normalize, bool useDotRatios=true)
virtual void clearDotTables(const int &flags)
Real dDotProduct(const double ¢er1, const double &width1, const double ¢er2, const double &width2) const