Developer Documentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
FunctionData.inl
1 /*
2 Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
3 All rights reserved.
4 
5 Redistribution and use in source and binary forms, with or without modification,
6 are permitted provided that the following conditions are met:
7 
8 Redistributions of source code must retain the above copyright notice, this list of
9 conditions and the following disclaimer. Redistributions in binary form must reproduce
10 the above copyright notice, this list of conditions and the following disclaimer
11 in the documentation and/or other materials provided with the distribution.
12 
13 Neither the name of the Johns Hopkins University nor the names of its contributors
14 may be used to endorse or promote products derived from this software without specific
15 prior written permission.
16 
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
18 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
19 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
20 SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
21 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
22 TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
23 BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
25 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
26 DAMAGE.
27 */
28 
30 // FunctionData //
32 template<int Degree,class Real>
33 FunctionData<Degree,Real>::FunctionData(void)
34 {
35  dotTable=dDotTable=d2DotTable=NULL;
36  valueTables=dValueTables=NULL;
37  res=0;
38 }
39 
40 template<int Degree,class Real>
41 FunctionData<Degree,Real>::~FunctionData(void)
42 {
43  if(res)
44  {
45  if( dotTable) delete[] dotTable;
46  if( dDotTable) delete[] dDotTable;
47  if(d2DotTable) delete[] d2DotTable;
48  if( valueTables) delete[] valueTables;
49  if(dValueTables) delete[] dValueTables;
50  }
51  dotTable=dDotTable=d2DotTable=NULL;
52  valueTables=dValueTables=NULL;
53  res=0;
54 }
55 
56 template<int Degree,class Real>
57 #if BOUNDARY_CONDITIONS
58 void FunctionData<Degree,Real>::set( const int& maxDepth , const PPolynomial<Degree>& F , const int& normalize , bool useDotRatios , bool reflectBoundary )
59 #else // !BOUNDARY_CONDITIONS
60 void FunctionData<Degree,Real>::set(const int& maxDepth,const PPolynomial<Degree>& F,const int& normalize , bool useDotRatios )
61 #endif // BOUNDARY_CONDITIONS
62 {
63  this->normalize = normalize;
64  this->useDotRatios = useDotRatios;
65 #if BOUNDARY_CONDITIONS
66  this->reflectBoundary = reflectBoundary;
67 #endif // BOUNDARY_CONDITIONS
68 
69  depth = maxDepth;
70  res = BinaryNode<double>::CumulativeCenterCount( depth );
71  res2 = (1<<(depth+1))+1;
72  baseFunctions = new PPolynomial<Degree+1>[res];
73  // Scale the function so that it has:
74  // 0] Value 1 at 0
75  // 1] Integral equal to 1
76  // 2] Square integral equal to 1
77  switch( normalize )
78  {
79  case 2:
80  baseFunction=F/sqrt((F*F).integral(F.polys[0].start,F.polys[F.polyCount-1].start));
81  break;
82  case 1:
83  baseFunction=F/F.integral(F.polys[0].start,F.polys[F.polyCount-1].start);
84  break;
85  default:
86  baseFunction=F/F(0);
87  }
88  dBaseFunction = baseFunction.derivative();
89 #if BOUNDARY_CONDITIONS
90  leftBaseFunction = baseFunction + baseFunction.shift( -1 );
91  rightBaseFunction = baseFunction + baseFunction.shift( 1 );
92  dLeftBaseFunction = leftBaseFunction.derivative();
93  dRightBaseFunction = rightBaseFunction.derivative();
94 #endif // BOUNDARY_CONDITIONS
95  double c1,w1;
96  for( int i=0 ; i<res ; i++ )
97  {
98  BinaryNode< double >::CenterAndWidth( i , c1 , w1 );
99 #if BOUNDARY_CONDITIONS
100  if( reflectBoundary )
101  {
102  int d , off;
103  BinaryNode< double >::DepthAndOffset( i , d , off );
104  if ( off==0 ) baseFunctions[i] = leftBaseFunction.scale( w1 ).shift( c1 );
105  else if( off==((1<<d)-1) ) baseFunctions[i] = rightBaseFunction.scale( w1 ).shift( c1 );
106  else baseFunctions[i] = baseFunction.scale( w1 ).shift( c1 );
107  }
108  else baseFunctions[i] = baseFunction.scale(w1).shift(c1);
109 #else // !BOUNDARY_CONDITIONS
110  baseFunctions[i] = baseFunction.scale(w1).shift(c1);
111 #endif // BOUNDARY_CONDITIONS
112  // Scale the function so that it has L2-norm equal to one
113  switch( normalize )
114  {
115  case 2:
116  baseFunctions[i]/=sqrt(w1);
117  break;
118  case 1:
119  baseFunctions[i]/=w1;
120  break;
121  }
122  }
123 }
124 template<int Degree,class Real>
125 void FunctionData<Degree,Real>::setDotTables( const int& flags )
126 {
127  clearDotTables( flags );
128  int size;
129  size = ( res*res + res )>>1;
130  if( flags & DOT_FLAG )
131  {
132  dotTable = new Real[size];
133  memset( dotTable , 0 , sizeof(Real)*size );
134  }
135  if( flags & D_DOT_FLAG )
136  {
137  dDotTable = new Real[size];
138  memset( dDotTable , 0 , sizeof(Real)*size );
139  }
140  if( flags & D2_DOT_FLAG )
141  {
142  d2DotTable = new Real[size];
143  memset( d2DotTable , 0 , sizeof(Real)*size );
144  }
145  double t1 , t2;
146  t1 = baseFunction.polys[0].start;
147  t2 = baseFunction.polys[baseFunction.polyCount-1].start;
148  for( int i=0 ; i<res ; i++ )
149  {
150  double c1 , c2 , w1 , w2;
151  BinaryNode<double>::CenterAndWidth( i , c1 , w1 );
152 #if BOUNDARY_CONDITIONS
153  int d1 , d2 , off1 , off2;
154  BinaryNode< double >::DepthAndOffset( i , d1 , off1 );
155  int boundary1 = 0;
156  if ( reflectBoundary && off1==0 ) boundary1 = -1;
157  else if( reflectBoundary && off1==( (1<<d1)-1 ) ) boundary1 = 1;
158 #endif // BOUNDARY_CONDITIONS
159  double start1 = t1 * w1 + c1;
160  double end1 = t2 * w1 + c1;
161  for( int j=0 ; j<=i ; j++ )
162  {
163  BinaryNode<double>::CenterAndWidth( j , c2 , w2 );
164 #if BOUNDARY_CONDITIONS
165  BinaryNode< double >::DepthAndOffset( j , d2 , off2 );
166  int boundary2 = 0;
167  if ( reflectBoundary && off2==0 ) boundary2 = -1;
168  else if( reflectBoundary && off2==( (1<<d2)-1 ) ) boundary2 = 1;
169 #endif // BOUNDARY_CONDITIONS
170  int idx = SymmetricIndex( i , j );
171 
172  double start = t1 * w2 + c2;
173  double end = t2 * w2 + c2;
174 #if BOUNDARY_CONDITIONS
175  if( reflectBoundary )
176  {
177  if( start<0 ) start = 0;
178  if( start>1 ) start = 1;
179  if( end <0 ) end = 0;
180  if( end >1 ) end = 1;
181  }
182 #endif // BOUNDARY_CONDITIONS
183 
184  if( start< start1 ) start = start1;
185  if( end > end1 ) end = end1;
186  if( start>= end ) continue;
187 
188 #if BOUNDARY_CONDITIONS
189  Real dot = dotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
190 #else // !BOUNDARY_CONDITIONS
191  Real dot = dotProduct( c1 , w1 , c2 , w2 );
192 #endif // BOUNDARY_CONDITIONS
193  if( fabs(dot)<1e-15 ) continue;
194  if( flags & DOT_FLAG ) dotTable[idx]=dot;
195  if( useDotRatios )
196  {
197 #if BOUNDARY_CONDITIONS
198  if( flags & D_DOT_FLAG ) dDotTable[idx] = -dDotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 ) / dot;
199  if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 ) / dot;
200 #else // !BOUNDARY_CONDITIONS
201  if( flags & D_DOT_FLAG ) dDotTable[idx] = -dDotProduct(c1,w1,c2,w2)/dot;
202  if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct(c1,w1,c2,w2)/dot;
203 #endif // BOUNDARY_CONDITIONS
204  }
205  else
206  {
207 #if BOUNDARY_CONDITIONS
208  if( flags & D_DOT_FLAG ) dDotTable[idx] = dDotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
209  if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
210 #else // !BOUNDARY_CONDTIONS
211  if( flags & D_DOT_FLAG ) dDotTable[idx] = dDotProduct(c1,w1,c2,w2);
212  if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct(c1,w1,c2,w2);
213 #endif // BOUNDARY_CONDITIONS
214  }
215  }
216  }
217 }
218 template<int Degree,class Real>
219 void FunctionData<Degree,Real>::clearDotTables( const int& flags )
220 {
221  if((flags & DOT_FLAG) && dotTable)
222  {
223  delete[] dotTable;
224  dotTable=NULL;
225  }
226  if((flags & D_DOT_FLAG) && dDotTable)
227  {
228  delete[] dDotTable;
229  dDotTable=NULL;
230  }
231  if((flags & D2_DOT_FLAG) && d2DotTable)
232  {
233  delete[] d2DotTable;
234  d2DotTable=NULL;
235  }
236 }
237 template<int Degree,class Real>
238 void FunctionData<Degree,Real>::setValueTables( const int& flags , const double& smooth )
239 {
240  clearValueTables();
241  if( flags & VALUE_FLAG ) valueTables = new Real[res*res2];
242  if( flags & D_VALUE_FLAG ) dValueTables = new Real[res*res2];
243  PPolynomial<Degree+1> function;
244  PPolynomial<Degree> dFunction;
245  for( int i=0 ; i<res ; i++ )
246  {
247  if(smooth>0)
248  {
249  function=baseFunctions[i].MovingAverage(smooth);
250  dFunction=baseFunctions[i].derivative().MovingAverage(smooth);
251  }
252  else
253  {
254  function=baseFunctions[i];
255  dFunction=baseFunctions[i].derivative();
256  }
257  for( int j=0 ; j<res2 ; j++ )
258  {
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));}
262  }
263  }
264 }
265 template<int Degree,class Real>
266 void FunctionData<Degree,Real>::setValueTables(const int& flags,const double& valueSmooth,const double& normalSmooth){
267  clearValueTables();
268  if(flags & VALUE_FLAG){ valueTables=new Real[res*res2];}
269  if(flags & D_VALUE_FLAG){dValueTables=new Real[res*res2];}
270  PPolynomial<Degree+1> function;
271  PPolynomial<Degree> dFunction;
272  for(int i=0;i<res;i++){
273  if(valueSmooth>0) { function=baseFunctions[i].MovingAverage(valueSmooth);}
274  else { function=baseFunctions[i];}
275  if(normalSmooth>0) {dFunction=baseFunctions[i].derivative().MovingAverage(normalSmooth);}
276  else {dFunction=baseFunctions[i].derivative();}
277 
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));}
282  }
283  }
284 }
285 
286 
287 template<int Degree,class Real>
288 void FunctionData<Degree,Real>::clearValueTables(void){
289  if( valueTables){delete[] valueTables;}
290  if(dValueTables){delete[] dValueTables;}
291  valueTables=dValueTables=NULL;
292 }
293 
294 #if BOUNDARY_CONDITIONS
295 template<int Degree,class Real>
296 Real FunctionData<Degree,Real>::dotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const
297 {
298  const PPolynomial< Degree > *b1 , *b2;
299  if ( boundary1==-1 ) b1 = & leftBaseFunction;
300  else if( boundary1== 0 ) b1 = & baseFunction;
301  else if( boundary1== 1 ) b1 = &rightBaseFunction;
302  if ( boundary2==-1 ) b2 = & leftBaseFunction;
303  else if( boundary2== 0 ) b2 = & baseFunction;
304  else if( boundary2== 1 ) b2 = &rightBaseFunction;
305  double r=fabs( baseFunction.polys[0].start );
306  switch( normalize )
307  {
308  case 2:
309  return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/sqrt(width1*width2));
310  case 1:
311  return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/(width1*width2));
312  default:
313  return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1);
314  }
315 }
316 template<int Degree,class Real>
317 Real FunctionData<Degree,Real>::dDotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const
318 {
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);
328  switch(normalize){
329  case 2:
330  return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/sqrt(width1*width2));
331  case 1:
332  return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/(width1*width2));
333  default:
334  return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r));
335  }
336 }
337 template<int Degree,class Real>
338 Real FunctionData<Degree,Real>::d2DotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const
339 {
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);
348  switch( normalize )
349  {
350  case 2:
351  return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/sqrt(width1*width2));
352  case 1:
353  return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/(width1*width2));
354  default:
355  return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2);
356  }
357 }
358 #else // !BOUNDARY_CONDITIONS
359 template<int Degree,class Real>
360 Real FunctionData<Degree,Real>::dotProduct(const double& center1,const double& width1,const double& center2,const double& width2) const{
361  double r=fabs(baseFunction.polys[0].start);
362  switch( normalize )
363  {
364  case 2:
365  return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/sqrt(width1*width2));
366  case 1:
367  return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/(width1*width2));
368  default:
369  return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1);
370  }
371 }
372 template<int Degree,class Real>
373 Real FunctionData<Degree,Real>::dDotProduct(const double& center1,const double& width1,const double& center2,const double& width2) const{
374  double r=fabs(baseFunction.polys[0].start);
375  switch(normalize){
376  case 2:
377  return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/sqrt(width1*width2));
378  case 1:
379  return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/(width1*width2));
380  default:
381  return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r));
382  }
383 }
384 template<int Degree,class Real>
385 Real FunctionData<Degree,Real>::d2DotProduct(const double& center1,const double& width1,const double& center2,const double& width2) const{
386  double r=fabs(baseFunction.polys[0].start);
387  switch(normalize){
388  case 2:
389  return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/sqrt(width1*width2));
390  case 1:
391  return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/(width1*width2));
392  default:
393  return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2);
394  }
395 }
396 #endif // BOUNDARY_CONDITIONS
397 template<int Degree,class Real>
398 inline int FunctionData<Degree,Real>::SymmetricIndex( const int& i1 , const int& i2 )
399 {
400  if( i1>i2 ) return ((i1*i1+i1)>>1)+i2;
401  else return ((i2*i2+i2)>>1)+i1;
402 }
403 template<int Degree,class Real>
404 inline int FunctionData<Degree,Real>::SymmetricIndex( const int& i1 , const int& i2 , int& index )
405 {
406  if( i1<i2 )
407  {
408  index = ((i2*i2+i2)>>1)+i1;
409  return 1;
410  }
411  else{
412  index = ((i1*i1+i1)>>1)+i2;
413  return 0;
414  }
415 }
OpenMesh::VectorT::dot
Scalar dot(const VectorT< Scalar, N > &_v1, const VectorT< Scalar, N > &_v2)
Definition: MeshNode2T.cc:263