38 template<
int Degree,
class Real>
41 dotTable=dDotTable=d2DotTable=NULL;
42 valueTables=dValueTables=NULL;
46 template<
int Degree,
class Real>
51 if( dotTable)
delete[] dotTable;
52 if( dDotTable)
delete[] dDotTable;
53 if(d2DotTable)
delete[] d2DotTable;
54 if( valueTables)
delete[] valueTables;
55 if(dValueTables)
delete[] dValueTables;
57 dotTable=dDotTable=d2DotTable=NULL;
58 valueTables=dValueTables=NULL;
62 template<
int Degree,
class Real>
63 #if BOUNDARY_CONDITIONS
65 #else // !BOUNDARY_CONDITIONS
67 #endif // BOUNDARY_CONDITIONS
69 this->normalize = normalize;
70 this->useDotRatios = useDotRatios;
71 #if BOUNDARY_CONDITIONS
72 this->reflectBoundary = reflectBoundary;
73 #endif // BOUNDARY_CONDITIONS
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();
100 #endif // BOUNDARY_CONDITIONS
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);
115 #else // !BOUNDARY_CONDITIONS
116 baseFunctions[i] = baseFunction.scale(w1).shift(c1);
117 #endif // BOUNDARY_CONDITIONS
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];
139 memset( dotTable , 0 ,
sizeof(
Real)*size );
141 if( flags & D_DOT_FLAG )
143 dDotTable =
new Real[size];
144 memset( dDotTable , 0 ,
sizeof(
Real)*size );
146 if( flags & D2_DOT_FLAG )
148 d2DotTable =
new Real[size];
149 memset( d2DotTable , 0 ,
sizeof(
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
159 int d1 , d2 , off1 , off2;
162 if ( reflectBoundary && off1==0 ) boundary1 = -1;
163 else if( reflectBoundary && off1==( (1<<d1)-1 ) ) boundary1 = 1;
164 #endif // BOUNDARY_CONDITIONS
165 double start1 = t1 * w1 + c1;
166 double end1 = t2 * w1 + c1;
167 for(
int j=0 ; j<=i ; j++ )
170 #if BOUNDARY_CONDITIONS
173 if ( reflectBoundary && off2==0 ) boundary2 = -1;
174 else if( reflectBoundary && off2==( (1<<d2)-1 ) ) boundary2 = 1;
175 #endif // 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;
188 #endif // BOUNDARY_CONDITIONS
190 if( start< start1 ) start = start1;
191 if( end > end1 ) end = end1;
192 if( start>= end )
continue;
194 #if BOUNDARY_CONDITIONS
195 Real dot = dotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
196 #else // !BOUNDARY_CONDITIONS
197 Real dot = dotProduct( c1 , w1 , c2 , w2 );
198 #endif // BOUNDARY_CONDITIONS
199 if( fabs(dot)<1e-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;
206 #else // !BOUNDARY_CONDITIONS
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;
209 #endif // BOUNDARY_CONDITIONS
213 #if BOUNDARY_CONDITIONS
214 if( flags & D_DOT_FLAG ) dDotTable[idx] = dDotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
215 if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
216 #else // !BOUNDARY_CONDTIONS
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);
219 #endif // BOUNDARY_CONDITIONS
224 template<
int Degree,
class Real>
227 if((flags & DOT_FLAG) && dotTable)
232 if((flags & D_DOT_FLAG) && dDotTable)
237 if((flags & D2_DOT_FLAG) && d2DotTable)
243 template<
int Degree,
class Real>
247 if( flags & VALUE_FLAG ) valueTables =
new Real[res*res2];
248 if( flags & D_VALUE_FLAG ) dValueTables =
new Real[res*res2];
251 for(
int i=0 ; i<res ; i++ )
260 function=baseFunctions[i];
263 for(
int j=0 ; j<res2 ; j++ )
265 double x=double(j)/(res2-1);
266 if(flags & VALUE_FLAG){ valueTables[j*res+i]=
Real(
function(x));}
267 if(flags & D_VALUE_FLAG){dValueTables[j*res+i]=
Real(dFunction(x));}
271 template<
int Degree,
class Real>
274 if(flags & VALUE_FLAG){ valueTables=
new Real[res*res2];}
275 if(flags & D_VALUE_FLAG){dValueTables=
new Real[res*res2];}
278 for(
int i=0;i<res;i++){
279 if(valueSmooth>0) {
function=baseFunctions[i].
MovingAverage(valueSmooth);}
280 else {
function=baseFunctions[i];}
282 else {dFunction=baseFunctions[i].
derivative();}
284 for(
int j=0;j<res2;j++){
285 double x=double(j)/(res2-1);
286 if(flags & VALUE_FLAG){ valueTables[j*res+i]=
Real(
function(x));}
287 if(flags & D_VALUE_FLAG){dValueTables[j*res+i]=
Real(dFunction(x));}
293 template<
int Degree,
class Real>
295 if( valueTables){
delete[] valueTables;}
296 if(dValueTables){
delete[] dValueTables;}
297 valueTables=dValueTables=NULL;
300 #if BOUNDARY_CONDITIONS
301 template<
int Degree,
class Real>
305 if ( boundary1==-1 ) b1 = & leftBaseFunction;
306 else if( boundary1== 0 ) b1 = & baseFunction;
307 else if( boundary1== 1 ) b1 = &rightBaseFunction;
308 if ( boundary2==-1 ) b2 = & leftBaseFunction;
309 else if( boundary2== 0 ) b2 = & baseFunction;
310 else if( boundary2== 1 ) b2 = &rightBaseFunction;
311 double r=fabs( baseFunction.polys[0].start );
315 return Real(((*b1)*b2->
scale(width2/width1).
shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/sqrt(width1*width2));
317 return Real(((*b1)*b2->
scale(width2/width1).
shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/(width1*width2));
319 return Real(((*b1)*b2->
scale(width2/width1).
shift((center2-center1)/width1)).integral(-2*r,2*r)*width1);
322 template<
int Degree,
class Real>
325 const PPolynomial< Degree-1 > *b1;
326 const PPolynomial< Degree > *b2;
327 if ( boundary1==-1 ) b1 = & dLeftBaseFunction;
328 else if( boundary1== 0 ) b1 = & dBaseFunction;
329 else if( boundary1== 1 ) b1 = &dRightBaseFunction;
330 if ( boundary2==-1 ) b2 = & leftBaseFunction;
331 else if( boundary2== 0 ) b2 = & baseFunction;
332 else if( boundary2== 1 ) b2 = & rightBaseFunction;
333 double r=fabs(baseFunction.polys[0].start);
336 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/sqrt(width1*width2));
338 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/(width1*width2));
340 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r));
343 template<
int Degree,
class Real>
346 const PPolynomial< Degree-1 > *b1 , *b2;
347 if ( boundary1==-1 ) b1 = & dLeftBaseFunction;
348 else if( boundary1== 0 ) b1 = & dBaseFunction;
349 else if( boundary1== 1 ) b1 = &dRightBaseFunction;
350 if ( boundary2==-1 ) b2 = & dLeftBaseFunction;
351 else if( boundary2== 0 ) b2 = & dBaseFunction;
352 else if( boundary2== 1 ) b2 = &dRightBaseFunction;
353 double r=fabs(baseFunction.polys[0].start);
357 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/sqrt(width1*width2));
359 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/(width1*width2));
361 return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2);
364 #else // !BOUNDARY_CONDITIONS
365 template<
int Degree,
class Real>
367 double r=fabs(baseFunction.polys[0].start);
371 return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/sqrt(width1*width2));
373 return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/(width1*width2));
375 return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1);
378 template<
int Degree,
class Real>
380 double r=fabs(baseFunction.polys[0].start);
383 return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/sqrt(width1*width2));
385 return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/(width1*width2));
387 return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r));
390 template<
int Degree,
class Real>
392 double r=fabs(baseFunction.polys[0].start);
395 return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/sqrt(width1*width2));
397 return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/(width1*width2));
399 return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2);
402 #endif // BOUNDARY_CONDITIONS
403 template<
int Degree,
class Real>
406 if( i1>i2 )
return ((i1*i1+i1)>>1)+i2;
407 else return ((i2*i2+i2)>>1)+i1;
409 template<
int Degree,
class Real>
414 index = ((i2*i2+i2)>>1)+i1;
418 index = ((i1*i1+i1)>>1)+i2;
static int SymmetricIndex(const int &i1, const int &i2)
virtual void clearValueTables(void)
Real dotProduct(const double ¢er1, const double &width1, const double ¢er2, const double &width2, int boundary1, int boundary2) const
PPolynomial< Degree+1 > MovingAverage(double radius)
PPolynomial scale(double s) const
static int CumulativeCenterCount(int maxDepth)
virtual void setDotTables(const int &flags)
virtual void setValueTables(const int &flags, const double &smooth=0)
PPolynomial< Degree-1 > derivative(void) const
virtual void clearDotTables(const int &flags)
static void CenterAndWidth(int depth, int offset, Real ¢er, Real &width)
Real dDotProduct(const double ¢er1, const double &width1, const double ¢er2, const double &width2, int boundary1, int boundary2) const
PPolynomial shift(double t) const
Real d2DotProduct(const double ¢er1, const double &width1, const double ¢er2, const double &width2, int boundary1, int boundary2) const
void set(const int &maxDepth, const PPolynomial< Degree > &F, const int &normalize, bool useDotRatios=true, bool reflectBoundary=false)
static void DepthAndOffset(int idx, int &depth, int &offset)