118 this->useDotRatios = useDotRatios;
119 this->reflectBoundary = reflectBoundary;
129 for(
int i=0 ; i<=Degree ; i++ ) baseBSpline[i] = Polynomial< Degree >::BSplineComponent( i ).shift(
static_cast<double>(-(Degree+1)/2) + i - 0.5 );
130 dBaseFunction = baseFunction.derivative();
133 for(
int i=0 ; i<Degree+3 ; i++ )
135 sPolys[i].
start =
static_cast<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 =
static_cast<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);
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]{};
184 if( flags & DV_DOT_FLAG )
186 dvDotTable =
new Real[fullSize]{};
188 if( flags & DD_DOT_FLAG )
190 ddDotTable =
new Real[size]{};
192 double vvIntegrals[Degree+1][Degree+1];
193 double vdIntegrals[Degree+1][Degree ];
194 double dvIntegrals[Degree ][Degree+1];
195 double ddIntegrals[Degree ][Degree ];
196 int vvSums[Degree+1][Degree+1];
197 int vdSums[Degree+1][Degree ];
198 int dvSums[Degree ][Degree+1];
199 int ddSums[Degree ][Degree ];
200 SetBSplineElementIntegrals< Degree , Degree >( vvIntegrals );
205 for(
int d1=0 ; d1<=depth ; d1++ )
206 for(
int off1=0 ; off1<(1<<d1) ; off1++ )
216 for(
int i=0 ; i<int(b1.size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
218 if( b1[i][j] && start1==-1 ) start1 = i;
219 if( b1[i][j] ) end1 = i+1;
221 for(
int d2=d1 ; d2<=depth ; d2++ )
223 for(
int off2=0 ; off2<(1<<d2) ; off2++ )
225 int start2 = off2-Degree;
226 int end2 = off2+Degree+1;
227 if( start2>=end1 || start1>=end2 )
continue;
228 start2 = std::max< int >( start1 , start2 );
229 end2 = std::min< int >( end1 , end2 );
230 if( d1==d2 && off2<off1 )
continue;
236 int idx = SymmetricIndex( ii , jj );
237 int idx1 = Index( ii , jj ) , idx2 = Index( jj , ii );
239 memset( vvSums , 0 ,
sizeof(
int ) * ( Degree+1 ) * ( Degree+1 ) );
240 memset( vdSums , 0 ,
sizeof(
int ) * ( Degree+1 ) * ( Degree ) );
241 memset( dvSums , 0 ,
sizeof(
int ) * ( Degree ) * ( Degree+1 ) );
242 memset( ddSums , 0 ,
sizeof(
int ) * ( Degree ) * ( Degree ) );
243 for(
int i=start2 ; i<end2 ; i++ )
245 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) vvSums[j][k] += b1[i][j] * b2[i][k];
246 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) vdSums[j][k] += b1[i][j] * db2[i][k];
247 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) dvSums[j][k] += db1[i][j] * b2[i][k];
248 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) ddSums[j][k] += db1[i][j] * db2[i][k];
250 double vvDot = 0 , dvDot = 0 , vdDot = 0 , ddDot = 0;
251 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) vvDot += vvIntegrals[j][k] * vvSums[j][k];
252 for(
int j=0 ; j<=Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) vdDot += vdIntegrals[j][k] * vdSums[j][k];
253 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k<=Degree ; k++ ) dvDot += dvIntegrals[j][k] * dvSums[j][k];
254 for(
int j=0 ; j< Degree ; j++ )
for(
int k=0 ; k< Degree ; k++ ) ddDot += ddIntegrals[j][k] * ddSums[j][k];
261 if( fabs(vvDot)<1e-15 )
continue;
262 if( flags & VV_DOT_FLAG ) vvDotTable [idx] =
Real( vvDot );
265 if( flags & DV_DOT_FLAG ) dvDotTable[idx1] =
Real( dvDot / vvDot );
266 if( flags & DV_DOT_FLAG ) dvDotTable[idx2] =
Real( vdDot / vvDot );
267 if( flags & DD_DOT_FLAG ) ddDotTable[idx ] =
Real( ddDot / vvDot );
271 if( flags & DV_DOT_FLAG ) dvDotTable[idx1] =
Real( dvDot );
272 if( flags & DV_DOT_FLAG ) dvDotTable[idx2] =
Real( dvDot );
273 if( flags & DD_DOT_FLAG ) ddDotTable[idx ] =
Real( ddDot );
281 for(
int i=0 ; i<int(b1.size()) ; i++ )
for(
int j=0 ; j<=Degree ; j++ )
283 if( b1[i][j] && start1==-1 ) start1 = i;
284 if( b1[i][j] ) end1 = i+1;
325 if( flags & VALUE_FLAG ) valueTables =
new Real[functionCount*sampleCount];
326 if( flags & D_VALUE_FLAG ) dValueTables =
new Real[functionCount*sampleCount];
329 for(
int i=0 ; i<functionCount ; i++ )
338 function = baseFunctions[i];
341 for(
int j=0 ; j<sampleCount ; j++ )
343 double x=
static_cast<double>(j)/(sampleCount-1);
344 if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=
Real( function(x));}
345 if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=
Real(dFunction(x));}
352 if(flags & VALUE_FLAG){ valueTables=
new Real[functionCount*sampleCount];}
353 if(flags & D_VALUE_FLAG){dValueTables=
new Real[functionCount*sampleCount];}
356 for(
int i=0;i<functionCount;i++){
357 if(valueSmooth>0) { function=baseFunctions[i].
MovingAverage(valueSmooth);}
358 else { function=baseFunctions[i];}
360 else {dFunction=baseFunctions[i].
derivative();}
362 for(
int j=0;j<sampleCount;j++){
363 double x=
static_cast<double>(j)/(sampleCount-1);
364 if(flags & VALUE_FLAG){ valueTables[j*functionCount+i]=
Real( function(x));}
365 if(flags & D_VALUE_FLAG){dValueTables[j*functionCount+i]=
Real(dFunction(x));}