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
29 #ifndef BSPLINE_DATA_INCLUDED
30 #define BSPLINE_DATA_INCLUDED
31
32 #define NEW_BSPLINE_CODE
33
34 #include "BinaryNode.h"
35 #include "PPolynomial.h"
36 #include "Array.h"
37
38 enum BoundaryType
39 {
40 BOUNDARY_FREE ,
41 BOUNDARY_DIRICHLET ,
42 BOUNDARY_NEUMANN ,
43 BOUNDARY_COUNT
44 };
45 const char* BoundaryNames[] = { "free" , "Dirichlet" , "Neumann" };
HasPartitionOfUnity(void)46 template< BoundaryType BType > inline bool HasPartitionOfUnity( void ){ return BType!=BOUNDARY_DIRICHLET; }
47
48 // This class represents a function that is a linear combination of B-spline elements.
49 // The coeff member indicating how much of each element is present.
50 // [WARNING] The ordering of B-spline elements is in the opposite order from that returned by Polynomial::BSplineComponent
51 template< int Degree >
52 struct BSplineElementCoefficients
53 {
54 int coeffs[Degree+1];
BSplineElementCoefficientsBSplineElementCoefficients55 BSplineElementCoefficients( void ){ memset( coeffs , 0 , sizeof( coeffs ) ); }
56 int& operator[]( int idx ){ return coeffs[idx]; }
57 const int& operator[]( int idx ) const { return coeffs[idx]; }
58 };
59
60 // This class represents a function on the the interval, partitioned into "res" blocks.
61 // On each block, the function is a degree-Degree polynomial, represented by the coefficients
62 // in the associated BSplineElementCoefficients.
63 // [NOTE] This representation of a function is agnostic to the type of boundary conditions (though the constructor is not).
64 template< int Degree >
65 struct BSplineElements : public std::vector< BSplineElementCoefficients< Degree > >
66 {
67 static const bool _Primal = (Degree&1);
68 static const int _Off = (Degree+1)/2;
69 static int _ReflectLeft ( int offset , int res );
70 static int _ReflectRight( int offset , int res );
71 static int _RotateLeft ( int offset , int res );
72 static int _RotateRight ( int offset , int res );
73 template< bool Left > void _addPeriodic( int offset , bool negate );
74 public:
75 // Coefficients are ordered as "/" "-" "\"
76 // [WARNING] This is the opposite of the order in Polynomial::BSplineComponent
77 int denominator;
78
BSplineElementsBSplineElements79 BSplineElements( void ) { denominator = 1; }
80 BSplineElements( int res , int offset , BoundaryType bType );
81
82 void upSample( BSplineElements& high ) const;
83 template< unsigned int D >
84 void differentiate( BSplineElements< Degree-D >& d ) const;
85
86 void print( FILE* fp=stdout ) const
87 {
88 for( int i=0 ; i<std::vector< BSplineElementCoefficients< Degree > >::size() ; i++ )
89 {
90 printf( "%d]" , i );
91 for( int j=0 ; j<=Degree ; j++ ) printf( " %d" , (*this)[i][j] );
92 printf( " (%d)\n" , denominator );
93 }
94 }
polynomialBSplineElements95 Polynomial< Degree > polynomial( int idx ) const
96 {
97 int res = (int)std::vector< BSplineElementCoefficients< Degree > >::size();
98 Polynomial< Degree > P;
99 if( idx>=0 && idx<res ) for( int d=0 ; d<=Degree ; d++ ) P += Polynomial< Degree >::BSplineComponent( Degree-d ).scale( 1./res ).shift( (idx+0.)/res ) * ( (*this)[idx][d] );
100 return P / denominator;
101 }
pPolynomialBSplineElements102 PPolynomial< Degree > pPolynomial( void ) const
103 {
104 int res = (int)std::vector< BSplineElementCoefficients< Degree > >::size();
105 PPolynomial< Degree > P;
106 P.polyCount = res + 1;
107 P.polys = AllocPointer< StartingPolynomial< Degree > >( P.polyCount );
108 for( int i=0 ; i<P.polyCount ; i++ ) P.polys[i].start = (i+0.) / res , P.polys[i].p = polynomial(i);
109 for( int i=res ; i>=1 ; i-- ) P.polys[i].p -= P.polys[i-1].p;
110 return P.compress(0);
111 }
112 };
113 template< int Degree , int DDegree > struct Differentiator { static void Differentiate( const BSplineElements< Degree >& bse , BSplineElements< DDegree >& dbse ); };
114 template< int Degree > struct Differentiator< Degree , Degree >{ static void Differentiate( const BSplineElements< Degree >& bse , BSplineElements< Degree >& dbse ); };
115 #define BSPLINE_SET_BOUNDS( name , s , e ) \
116 static const int name ## Start = (s); \
117 static const int name ## End = (e); \
118 static const int name ## Size = (e)-(s)+1
119
120 // Assumes that x is non-negative
121 #define _FLOOR_OF_HALF( x ) ( (x) >>1 )
122 #define _CEIL_OF_HALF( x ) ( ( (x)+1 )>>1 )
123 // Done with the assumption
124 #define FLOOR_OF_HALF( x ) ( (x)<0 ? - _CEIL_OF_HALF( -(x) ) : _FLOOR_OF_HALF( x ) )
125 #define CEIL_OF_HALF( x ) ( (x)<0 ? - _FLOOR_OF_HALF( -(x) ) : _CEIL_OF_HALF( x ) )
126 #define SMALLEST_INTEGER_LARGER_THAN_HALF( x ) ( CEIL_OF_HALF( (x)+1 ) )
127 #define LARGEST_INTEGER_SMALLER_THAN_HALF( x ) ( FLOOR_OF_HALF( (x)-1 ) )
128 #define SMALLEST_INTEGER_LARGER_THAN_OR_EQUAL_TO_HALF( x ) ( CEIL_OF_HALF( x ) )
129 #define LARGEST_INTEGER_SMALLER_THAN_OR_EQUAL_TO_HALF( x ) ( FLOOR_OF_HALF( x ) )
130
131 template< int Degree >
132 struct BSplineSupportSizes
133 {
134 inline static int Nodes( int depth ){ return ( 1<<depth ) + ( Degree&1 ); }
135 inline static bool OutOfBounds( int depth , int offset ){ return offset>=0 || offset<Nodes(depth); }
136 // An index is interiorly supported if its support is in the range [0,1<<depth)
137 inline static void InteriorSupportedSpan( int depth , int& begin , int& end ){ begin = -SupportStart , end = (1<<depth)-SupportEnd; }
138
139 // If the degree is even, we use a dual basis and functions are centered at the center of the interval
140 // It the degree is odd, we use a primal basis and functions are centered at the left end of the interval
141 // The function at index I is supported in:
142 // Support( I ) = [ I - (Degree+1-Inset)/2 , I + (Degree+1+Inset)/2 ]
143 // [NOTE] The value of ( Degree + 1 +/- Inset ) is always even
144 static const int Inset = (Degree&1) ? 0 : 1;
145 BSPLINE_SET_BOUNDS( Support , -( (Degree+1)/2 ) , Degree/2 );
146 BSPLINE_SET_BOUNDS( ChildSupport , 2*SupportStart , 2*(SupportEnd+1)-1 );
147 BSPLINE_SET_BOUNDS( Corner , SupportStart+1 , SupportEnd );
148 BSPLINE_SET_BOUNDS( ChildCorner , 2*SupportStart+1 , 2*SupportEnd + 1 );
149
150 // Setting I=0, we are looking for the smallest/largest integers J such that:
151 // Support( 0 ) CONTAINS Support( J )
152 // <=> [-(Degree+1-Inset) , (Degree+1+Inset) ] CONTAINS [ J-(Degree+1-Inset)/2 , J+(Degree+1+Inset)/2 ]
153 // Which is the same as the smallest/largest integers J such that:
154 // J - (Degree+1-Inset)/2 >= -(Degree+1-Inset) | J + (Degree+1+Inset)/2 <= (Degree+1+Inset)
155 // <=> J >= -(Degree+1-Inset)/2 | J <= (Degree+1+Inset)/2
156 BSPLINE_SET_BOUNDS( UpSample , - ( Degree + 1 - Inset ) / 2 , ( Degree + 1 + Inset ) /2 );
157
158 // Setting I=0/1, we are looking for the smallest/largest integers J such that:
159 // Support( J ) CONTAINS Support( 0/1 )
160 // <=> [ 2*J - (Degree+1-Inset) , 2*J + (Degree+1+Inset) ] CONTAINS [ 0/1 - (Degree+1-Inset)/2 , 0/1 + (Degree+1+Inset)/2 ]
161 // Which is the same as the smallest/largest integers J such that:
162 // 2*J + (Degree+1+Inset) >= 0/1 + (Degree+1+Inset)/2 | 2*J - (Degree+1-Inset) <= 0/1 - (Degree+1-Inset)/2
163 // <=> 2*J >= 0/1 - (Degree+1+Inset)/2 | 2*J <= 0/1 + (Degree+1-Inset)/2
164 BSPLINE_SET_BOUNDS( DownSample0 , SMALLEST_INTEGER_LARGER_THAN_OR_EQUAL_TO_HALF( 0 - ( Degree + 1 + Inset ) / 2 ) , LARGEST_INTEGER_SMALLER_THAN_OR_EQUAL_TO_HALF( 0 + ( Degree + 1 - Inset ) / 2 ) );
165 BSPLINE_SET_BOUNDS( DownSample1 , SMALLEST_INTEGER_LARGER_THAN_OR_EQUAL_TO_HALF( 1 - ( Degree + 1 + Inset ) / 2 ) , LARGEST_INTEGER_SMALLER_THAN_OR_EQUAL_TO_HALF( 1 + ( Degree + 1 - Inset ) / 2 ) );
166 static const int DownSampleStart[] , DownSampleEnd[] , DownSampleSize[];
167 };
168 template< int Degree > const int BSplineSupportSizes< Degree >::DownSampleStart[] = { DownSample0Start , DownSample1Start };
169 template< int Degree > const int BSplineSupportSizes< Degree >::DownSampleEnd [] = { DownSample0End , DownSample1End };
170 template< int Degree > const int BSplineSupportSizes< Degree >::DownSampleSize [] = { DownSample0Size , DownSample1Size };
171
172
173 // Given a B-Spline of degree Degree1 at position i, this gives the offsets of the B-splines of degree Degree2 that just overlap with it.
174 template< int Degree1 , int Degree2 >
175 struct BSplineOverlapSizes
176 {
177 typedef BSplineSupportSizes< Degree1 > EData1;
178 typedef BSplineSupportSizes< Degree2 > EData2;
179 BSPLINE_SET_BOUNDS( Overlap , EData1:: SupportStart - EData2::SupportEnd , EData1:: SupportEnd - EData2::SupportStart );
180 BSPLINE_SET_BOUNDS( ChildOverlap , EData1::ChildSupportStart - EData2::SupportEnd , EData1::ChildSupportEnd - EData2::SupportStart );
181 BSPLINE_SET_BOUNDS( OverlapSupport , OverlapStart + EData2::SupportStart , OverlapEnd + EData2::SupportEnd );
182 BSPLINE_SET_BOUNDS( ChildOverlapSupport , ChildOverlapStart + EData2::SupportStart , ChildOverlapEnd + EData2::SupportEnd );
183
184 // Setting I=0/1, we are looking for the smallest/largest integers J such that:
185 // Support( 2*J ) * 2 INTERSECTION Support( 0/1 ) NON-EMPTY
186 // <=> [ 2*J - (Degree2+1-Inset2) , 2*J + (Degree2+1+Inset2) ] INTERSECTION [ 0/1 - (Degree1+1-Inset1)/2 , 0/1 + (Degree1+1+Inset1)/2 ] NON-EMPTY
187 // Which is the same as the smallest/largest integers J such that:
188 // 0/1 - (Degree1+1-Inset1)/2 < 2*J + (Degree2+1+Inset2) | 0/1 + (Degree1+1+Inset1)/2 > 2*J - (Degree2+1-Inset2)
189 // <=> 2*J > 0/1 - ( 2*Degree2 + Degree1 + 3 + 2*Inset2 - Inset1 ) / 2 | 2*J < 0/1 + ( 2*Degree2 + Degree1 + 3 - 2*Inset2 + Inset1 ) / 2
190 BSPLINE_SET_BOUNDS( ParentOverlap0 , SMALLEST_INTEGER_LARGER_THAN_HALF( 0 - ( 2*Degree2 + Degree1 + 3 + 2*EData2::Inset - EData1::Inset ) / 2 ) , LARGEST_INTEGER_SMALLER_THAN_HALF( 0 + ( 2*Degree2 + Degree1 + 3 - 2*EData2::Inset + EData1::Inset ) / 2 ) );
191 BSPLINE_SET_BOUNDS( ParentOverlap1 , SMALLEST_INTEGER_LARGER_THAN_HALF( 1 - ( 2*Degree2 + Degree1 + 3 + 2*EData2::Inset - EData1::Inset ) / 2 ) , LARGEST_INTEGER_SMALLER_THAN_HALF( 1 + ( 2*Degree2 + Degree1 + 3 - 2*EData2::Inset + EData1::Inset ) / 2 ) );
192 static const int ParentOverlapStart[] , ParentOverlapEnd[] , ParentOverlapSize[];
193 };
194 template< int Degree1 , int Degree2 > const int BSplineOverlapSizes< Degree1 , Degree2 >::ParentOverlapStart[] = { ParentOverlap0Start , ParentOverlap1Start };
195 template< int Degree1 , int Degree2 > const int BSplineOverlapSizes< Degree1 , Degree2 >::ParentOverlapEnd [] = { ParentOverlap0End , ParentOverlap1End };
196 template< int Degree1 , int Degree2 > const int BSplineOverlapSizes< Degree1 , Degree2 >::ParentOverlapSize [] = { ParentOverlap0Size , ParentOverlap1Size };
197
198 template< int Degree , BoundaryType BType >
199 class BSplineEvaluationData
200 {
201 public:
202 static const int Pad = (BType==BOUNDARY_FREE ) ? BSplineSupportSizes< Degree >::SupportEnd : ( (Degree&1) && BType==BOUNDARY_DIRICHLET ) ? -1 : 0;
203 inline static int Begin( int depth ){ return -Pad; }
204 inline static int End ( int depth ){ return (1<<depth) + (Degree&1) + Pad; }
205 inline static bool OutOfBounds( int depth , int offset ){ return offset<Begin(depth) || offset>=End(depth); }
206
207 static const int OffsetStart = -BSplineSupportSizes< Degree >::SupportStart , OffsetStop = BSplineSupportSizes< Degree >::SupportEnd + ( Degree&1 ) , IndexSize = OffsetStart + OffsetStop + 1 + 2 * Pad;
208 static int OffsetToIndex( int depth , int offset )
209 {
210 int dim = BSplineSupportSizes< Degree >::Nodes( depth );
211 if ( offset<OffsetStart ) return Pad + offset;
212 else if( offset>=dim-OffsetStop ) return Pad + OffsetStart + 1 + offset - ( dim-OffsetStop );
213 else return Pad + OffsetStart;
214 }
215 static inline int IndexToOffset( int depth , int idx ){ return ( idx-Pad<=OffsetStart ? idx - Pad : ( BSplineSupportSizes< Degree >::Nodes(depth) + Pad - IndexSize + idx ) ); }
216
217 BSplineEvaluationData( void );
218
219 // [NOTE] The offset represents the node position, not the index of the function
220 static double Value( int depth , int off , double s , bool derivative );
221
222 // Note that this struct stores the components in left-to-right order
223 struct BSplineComponents
224 {
225 protected:
226 Polynomial< Degree > _polys[Degree+1];
227 public:
228 BSplineComponents( void ){ ; }
229 BSplineComponents( int depth , int offset );
230 const Polynomial< Degree >& operator[] ( int idx ) const { return _polys[idx]; }
231 BSplineComponents derivative( void ) const;
232 void printnl( void ) const { for( int d=0 ; d<=Degree ; d++ ) printf( "[%d] " , d ) , _polys[d].printnl(); }
233 };
234 struct BSplineUpSamplingCoefficients
235 {
236 protected:
237 int _coefficients[ BSplineSupportSizes< Degree >::UpSampleSize ];
238 public:
239 BSplineUpSamplingCoefficients( void ){ ; }
240 BSplineUpSamplingCoefficients( int depth , int offset );
241 double operator[] ( int idx ){ return (double)_coefficients[idx] / (1<<Degree); }
242 };
243
244 struct CenterEvaluator
245 {
246 struct Evaluator
247 {
248 protected:
249 friend BSplineEvaluationData;
250 int _depth;
251 double _ccValues[2][IndexSize][BSplineSupportSizes< Degree >::SupportSize];
252 public:
253 #ifdef BRUNO_LEVY_FIX
254 Evaluator( void ){ _depth = 0 ; memset( _ccValues , 0 , sizeof(_ccValues) ); }
255 #endif // BRUNO_LEVY_FIX
256 double value( int fIdx , int cIdx , bool d ) const;
257 int depth( void ) const { return _depth; }
258 };
259 struct ChildEvaluator
260 {
261 protected:
262 friend BSplineEvaluationData;
263 int _parentDepth;
264 double _pcValues[2][IndexSize][BSplineSupportSizes< Degree >::ChildSupportSize];
265 public:
266 #ifdef BRUNO_LEVY_FIX
267 ChildEvaluator( void ){ _parentDepth = 0 ; memset( _pcValues , 0 , sizeof(_pcValues) ); }
268 #endif // BRUNO_LEVY_FIX
269 double value( int fIdx , int cIdx , bool d ) const;
270 int parentDepth( void ) const { return _parentDepth; }
271 int childDepth( void ) const { return _parentDepth+1; }
272 };
273 };
274 static void SetCenterEvaluator( typename CenterEvaluator::Evaluator& evaluator , int depth );
275 static void SetChildCenterEvaluator( typename CenterEvaluator::ChildEvaluator& evaluator , int parentDepth );
276
277 struct CornerEvaluator
278 {
279 struct Evaluator
280 {
281 protected:
282 friend BSplineEvaluationData;
283 int _depth;
284 double _ccValues[2][IndexSize][BSplineSupportSizes< Degree >::CornerSize];
285 public:
286 #ifdef BRUNO_LEVY_FIX
287 Evaluator( void ){ _depth = 0 ; memset( _ccValues , 0 , sizeof( _ccValues ) ); }
288 #endif // BRUNO_LEVY_FIX
289 double value( int fIdx , int cIdx , bool d ) const;
290 int depth( void ) const { return _depth; }
291 };
292 struct ChildEvaluator
293 {
294 protected:
295 friend BSplineEvaluationData;
296 int _parentDepth;
297 double _pcValues[2][IndexSize][BSplineSupportSizes< Degree >::ChildCornerSize];
298 public:
299 #ifdef BRUNO_LEVY_FIX
300 ChildEvaluator( void ){ _parentDepth = 0 ; memset( _pcValues , 0 , sizeof( _pcValues ) ); }
301 #endif // BRUNO_LEVY_FIX
302 double value( int fIdx , int cIdx , bool d ) const;
303 int parentDepth( void ) const { return _parentDepth; }
304 int childDepth( void ) const { return _parentDepth+1; }
305 };
306 };
307 static void SetCornerEvaluator( typename CornerEvaluator::Evaluator& evaluator , int depth );
308 static void SetChildCornerEvaluator( typename CornerEvaluator::ChildEvaluator& evaluator , int parentDepth );
309
310 struct Evaluator
311 {
312 typename CenterEvaluator::Evaluator centerEvaluator;
313 typename CornerEvaluator::Evaluator cornerEvaluator;
314 double centerValue( int fIdx , int cIdx , bool d ) const { return centerEvaluator.value( fIdx , cIdx , d ); }
315 double cornerValue( int fIdx , int cIdx , bool d ) const { return cornerEvaluator.value( fIdx , cIdx , d ); }
316 };
317 static void SetEvaluator( Evaluator& evaluator , int depth ){ SetCenterEvaluator( evaluator.centerEvaluator , depth ) , SetCornerEvaluator( evaluator.cornerEvaluator , depth ); }
318
319 struct ChildEvaluator
320 {
321 typename CenterEvaluator::ChildEvaluator centerEvaluator;
322 typename CornerEvaluator::ChildEvaluator cornerEvaluator;
323 double centerValue( int fIdx , int cIdx , bool d ) const { return centerEvaluator.value( fIdx , cIdx , d ); }
324 double cornerValue( int fIdx , int cIdx , bool d ) const { return cornerEvaluator.value( fIdx , cIdx , d ); }
325 };
326 static void SetChildEvaluator( ChildEvaluator& evaluator , int depth ){ SetChildCenterEvaluator( evaluator.centerEvaluator , depth ) , SetChildCornerEvaluator( evaluator.cornerEvaluator , depth ); }
327
328 struct UpSampleEvaluator
329 {
330 protected:
331 friend BSplineEvaluationData;
332 int _lowDepth;
333 double _pcValues[IndexSize][BSplineSupportSizes< Degree >::UpSampleSize];
334 public:
335 #ifdef BRUNO_LEVY_FIX
336 UpSampleEvaluator( void ){ _lowDepth = 0 ; memset( _pcValues , 0 , sizeof( _pcValues ) ); }
337 #endif // BRUNO_LEVY_FIX
338 double value( int pIdx , int cIdx ) const;
339 int lowDepth( void ) const { return _lowDepth; }
340 };
341 static void SetUpSampleEvaluator( UpSampleEvaluator& evaluator , int lowDepth );
342 };
343
344 template< int Degree1 , BoundaryType BType1 , int Degree2 , BoundaryType BType2 >
345 class BSplineIntegrationData
346 {
347 public:
348 static const int OffsetStart = - BSplineOverlapSizes< Degree1 , Degree2 >::OverlapSupportStart , OffsetStop = BSplineOverlapSizes< Degree1 , Degree2 >::OverlapSupportEnd + ( Degree1&1 ) , IndexSize = OffsetStart + OffsetStop + 1 + 2 * BSplineEvaluationData< Degree1 , BType1 >::Pad;
349 static int OffsetToIndex( int depth , int offset )
350 {
351 int dim = BSplineSupportSizes< Degree1 >::Nodes( depth );
352 if ( offset<OffsetStart ) return BSplineEvaluationData< Degree1 , BType1 >::Pad + offset;
353 else if( offset>=dim-OffsetStop ) return BSplineEvaluationData< Degree1 , BType1 >::Pad + OffsetStart + 1 + offset - ( dim-OffsetStop );
354 else return BSplineEvaluationData< Degree1 , BType1 >::Pad + OffsetStart;
355 }
356 static inline int IndexToOffset( int depth , int idx ){ return ( idx-BSplineEvaluationData< Degree1 , BType1 >::Pad<=OffsetStart ? idx-BSplineEvaluationData< Degree1 , BType1 >::Pad : ( BSplineSupportSizes< Degree1 >::Nodes(depth) + BSplineEvaluationData< Degree1 , BType1 >::Pad - IndexSize + idx ) ); }
357
358 template< unsigned int D1 , unsigned int D2 > static double Dot( int depth1 , int off1 , int depth2 , int off2 );
359 // An index is interiorly overlapped if the support of its overlapping neighbors is in the range [0,1<<depth)
360 inline static void InteriorOverlappedSpan( int depth , int& begin , int& end ){ begin = -BSplineOverlapSizes< Degree1 , Degree2 >::OverlapStart-BSplineSupportSizes< Degree2 >::SupportStart , end = (1<<depth)-BSplineOverlapSizes< Degree1 , Degree2 >::OverlapEnd-BSplineSupportSizes< Degree2 >::SupportEnd; }
361
362 struct FunctionIntegrator
363 {
364 template< unsigned int D1 , unsigned int D2 >
365 struct Integrator
366 {
367 protected:
368 friend BSplineIntegrationData;
369 int _depth;
370 double _ccIntegrals[D1+1][D2+1][IndexSize][BSplineOverlapSizes< Degree1 , Degree2 >::OverlapSize];
371 public:
372 #ifdef BRUNO_LEVY_FIX
373 Integrator( void ){ _depth = 0 ; memset(_ccIntegrals, 0, sizeof(_ccIntegrals)); }
374 #endif // BRUNO_LEVY_FIX
375 double dot( int fIdx1 , int fidx2 , int d1 , int d2 ) const;
376 int depth( void ) const { return _depth; }
377 };
378 template< unsigned int D1 , unsigned int D2 >
379 struct ChildIntegrator
380 {
381 protected:
382 friend BSplineIntegrationData;
383 int _parentDepth;
384 double _pcIntegrals[D1+1][D2+1][IndexSize][BSplineOverlapSizes< Degree1 , Degree2 >::ChildOverlapSize];
385 public:
386 #ifdef BRUNO_LEVY_FIX
387 ChildIntegrator( void ){ _parentDepth = 0 ; memset( _pcIntegrals , 0 , sizeof( _pcIntegrals ) ); }
388 #endif // BRUNO_LEVY_FIX
389 double dot( int fIdx1 , int fidx2 , int d1 , int d2 ) const;
390 int parentDepth( void ) const { return _parentDepth; }
391 int childDepth( void ) const { return _parentDepth+1; }
392 };
393 };
394 // D1 and D2 indicate the number of derivatives that should be taken
395 template< unsigned int D1 , unsigned int D2 >
396 static void SetIntegrator( typename FunctionIntegrator::template Integrator< D1 , D2 >& integrator , int depth );
397 template< unsigned int D1 , unsigned int D2 >
398 static void SetChildIntegrator( typename FunctionIntegrator::template ChildIntegrator< D1 , D2 >& integrator , int parentDepth );
399
400 protected:
401 // _D1 and _D2 indicate the total number of derivatives the integrator will be storing
402 template< unsigned int D1 , unsigned int D2 , unsigned int _D1 , unsigned int _D2 >
403 struct _IntegratorSetter
404 {
405 static void Set( typename FunctionIntegrator::template Integrator< _D1 , _D2 >& integrator , int depth );
406 static void Set( typename FunctionIntegrator::template ChildIntegrator< _D1 , _D2 >& integrator , int depth );
407 };
408
409 template< unsigned int D1 , unsigned int D2 , unsigned int _D1 , unsigned int _D2 , class Integrator >
410 struct IntegratorSetter
411 {
412 static void Set2D( Integrator& integrator , int depth );
413 static void Set1D( Integrator& integrator , int depth );
414 };
415 template< unsigned int D1 , unsigned int _D1 , unsigned int _D2 , class Integrator >
416 struct IntegratorSetter< D1 , 0 , _D1 , _D2 , Integrator >
417 {
418 static void Set2D( Integrator& integrator , int depth );
419 static void Set1D( Integrator& integrator , int depth );
420 };
421 template< unsigned int D2 , unsigned int _D1 , unsigned int _D2 , class Integrator >
422 struct IntegratorSetter< 0 , D2 , _D1 , _D2 , Integrator >
423 {
424 static void Set2D( Integrator& integrator , int depth );
425 static void Set1D( Integrator& integrator , int depth );
426 };
427 template< unsigned int _D1 , unsigned int _D2 , class Integrator >
428 struct IntegratorSetter< 0 , 0 , _D1 , _D2 , Integrator >
429 {
430 static void Set2D( Integrator& integrator , int depth );
431 static void Set1D( Integrator& integrator , int depth );
432 };
433 };
434 #undef BSPLINE_SET_BOUNDS
435 #undef _FLOOR_OF_HALF
436 #undef _CEIL_OF_HALF
437 #undef FLOOR_OF_HALF
438 #undef CEIL_OF_HALF
439 #undef SMALLEST_INTEGER_LARGER_THAN_HALF
440 #undef LARGEST_INTEGER_SMALLER_THAN_HALF
441 #undef SMALLEST_INTEGER_LARGER_THAN_OR_EQUAL_TO_HALF
442 #undef LARGEST_INTEGER_SMALLER_THAN_OR_EQUAL_TO_HALF
443
444 template< int Degree , BoundaryType BType >
445 struct BSplineData
446 {
447 inline static int TotalFunctionCount( int depth ){ return depth<0 ? 0 : (1<<(depth+1)) - 1 + (depth+1) * ( (Degree&1) + 2 * BSplineEvaluationData< Degree , BType >::Pad ); }
448 inline static int FunctionIndex( int depth , int offset ){ return TotalFunctionCount( depth-1 ) + offset + BSplineEvaluationData< Degree , BType >::Pad; }
449 inline static void FactorFunctionIndex( int idx , int& depth , int& offset )
450 {
451 int dim;
452 depth = 0;
453 while( idx>=( dim = BSplineEvaluationData< Degree , BType >::End( depth ) - BSplineEvaluationData< Degree , BType >::Begin( depth ) ) ) idx -= dim , depth++;
454 offset = idx - BSplineEvaluationData< Degree , BType >::Pad;
455 }
456 inline static void FunctionSpan( int depth , int& fStart , int& fEnd ){ fStart = TotalFunctionCount( depth-1 ) , fEnd = TotalFunctionCount( depth ); }
457 inline static int RemapOffset( int depth , int idx , bool& reflect );
458
459 size_t functionCount;
460 Pointer( typename BSplineEvaluationData< Degree , BType >::BSplineComponents ) baseBSplines;
461 Pointer( typename BSplineEvaluationData< Degree , BType >::BSplineComponents ) dBaseBSplines;
462
463 BSplineData( int maxDepth );
464 ~BSplineData( void );
465 };
466
467 template< int Degree1 , int Degree2 > void SetBSplineElementIntegrals( double integrals[Degree1+1][Degree2+1] );
468
469
470 #include "BSplineData.inl"
471 #endif // BSPLINE_DATA_INCLUDED