CoCalc -- gregory

GitHub Repository: godotengine/godot
Path: blob/master/thirdparty/embree/kernels/subdiv/gregory_patch.h
⁹⁹¹³ views
1
// Copyright 2009-2021 Intel Corporation
2
// SPDX-License-Identifier: Apache-2.0
3

4
#pragma once
5

6
#include "catmullclark_patch.h"
7
#include "bezier_patch.h"
8
#include "bezier_curve.h"
9
#include "catmullclark_coefficients.h"
10

11
namespace embree
12
{  
13
  template<typename Vertex, typename Vertex_t = Vertex>
14
  class __aligned(64) GregoryPatchT
15
  {
16
    typedef CatmullClarkPatchT<Vertex,Vertex_t> CatmullClarkPatch;
17
    typedef GeneralCatmullClarkPatchT<Vertex,Vertex_t> GeneralCatmullClarkPatch;
18
    typedef CatmullClark1RingT<Vertex,Vertex_t> CatmullClark1Ring;
19
    typedef BezierCurveT<Vertex> BezierCurve;
20

21
  public:
22
    Vertex v[4][4];
23
    Vertex f[2][2];
24

25
    __forceinline GregoryPatchT() {}
26

27
    __forceinline GregoryPatchT(const CatmullClarkPatch& patch) {
28
      init(patch);
29
    }
30

31
    __forceinline GregoryPatchT(const CatmullClarkPatch& patch, 
32
                                const BezierCurve* border0, const BezierCurve* border1, const BezierCurve* border2, const BezierCurve* border3) 
33
    {
34
      init_crackfix(patch,border0,border1,border2,border3);
35
    }
36

37
    __forceinline GregoryPatchT (const HalfEdge* edge, const char* vertices, size_t stride) { 
38
      init(CatmullClarkPatch(edge,vertices,stride));
39
    }
40
      
41
    __forceinline Vertex& p0() { return v[0][0]; }
42
    __forceinline Vertex& p1() { return v[0][3]; }
43
    __forceinline Vertex& p2() { return v[3][3]; }
44
    __forceinline Vertex& p3() { return v[3][0]; }
45
    
46
    __forceinline Vertex& e0_p() { return v[0][1]; }
47
    __forceinline Vertex& e0_m() { return v[1][0]; }
48
    __forceinline Vertex& e1_p() { return v[1][3]; }
49
    __forceinline Vertex& e1_m() { return v[0][2]; }
50
    __forceinline Vertex& e2_p() { return v[3][2]; }
51
    __forceinline Vertex& e2_m() { return v[2][3]; }
52
    __forceinline Vertex& e3_p() { return v[2][0]; }
53
    __forceinline Vertex& e3_m() { return v[3][1]; }
54
    
55
    __forceinline Vertex& f0_p() { return v[1][1]; }
56
    __forceinline Vertex& f1_p() { return v[1][2]; }
57
    __forceinline Vertex& f2_p() { return v[2][2]; }
58
    __forceinline Vertex& f3_p() { return v[2][1]; }
59
    __forceinline Vertex& f0_m() { return f[0][0]; }
60
    __forceinline Vertex& f1_m() { return f[0][1]; }
61
    __forceinline Vertex& f2_m() { return f[1][1]; }
62
    __forceinline Vertex& f3_m() { return f[1][0]; }
63
    
64
    __forceinline const Vertex& p0() const { return v[0][0]; }
65
    __forceinline const Vertex& p1() const { return v[0][3]; }
66
    __forceinline const Vertex& p2() const { return v[3][3]; }
67
    __forceinline const Vertex& p3() const { return v[3][0]; }
68
    
69
    __forceinline const Vertex& e0_p() const { return v[0][1]; }
70
    __forceinline const Vertex& e0_m() const { return v[1][0]; }
71
    __forceinline const Vertex& e1_p() const { return v[1][3]; }
72
    __forceinline const Vertex& e1_m() const { return v[0][2]; }
73
    __forceinline const Vertex& e2_p() const { return v[3][2]; }
74
    __forceinline const Vertex& e2_m() const { return v[2][3]; }
75
    __forceinline const Vertex& e3_p() const { return v[2][0]; }
76
    __forceinline const Vertex& e3_m() const { return v[3][1]; }
77
    
78
    __forceinline const Vertex& f0_p() const { return v[1][1]; }
79
    __forceinline const Vertex& f1_p() const { return v[1][2]; }
80
    __forceinline const Vertex& f2_p() const { return v[2][2]; }
81
    __forceinline const Vertex& f3_p() const { return v[2][1]; }
82
    __forceinline const Vertex& f0_m() const { return f[0][0]; }
83
    __forceinline const Vertex& f1_m() const { return f[0][1]; }
84
    __forceinline const Vertex& f2_m() const { return f[1][1]; }
85
    __forceinline const Vertex& f3_m() const { return f[1][0]; }
86
    
87
    __forceinline Vertex initCornerVertex(const CatmullClarkPatch& irreg_patch, const size_t index) {
88
      return irreg_patch.ring[index].getLimitVertex();
89
    }
90
    
91
    __forceinline Vertex initPositiveEdgeVertex(const CatmullClarkPatch& irreg_patch, const size_t index, const Vertex& p_vtx) {
92
      return madd(1.0f/3.0f,irreg_patch.ring[index].getLimitTangent(),p_vtx);
93
    }
94
    
95
    __forceinline Vertex initNegativeEdgeVertex(const CatmullClarkPatch& irreg_patch, const size_t index, const Vertex& p_vtx) {
96
      return madd(1.0f/3.0f,irreg_patch.ring[index].getSecondLimitTangent(),p_vtx);
97
    }
98

99
    __forceinline Vertex initPositiveEdgeVertex2(const CatmullClarkPatch& irreg_patch, const size_t index, const Vertex& p_vtx) 
100
    {
101
      CatmullClark1Ring3fa r0,r1,r2;
102
      irreg_patch.ring[index].subdivide(r0);
103
      r0.subdivide(r1);
104
      r1.subdivide(r2);
105
      return madd(8.0f/3.0f,r2.getLimitTangent(),p_vtx);
106
    }
107
    
108
    __forceinline Vertex initNegativeEdgeVertex2(const CatmullClarkPatch& irreg_patch, const size_t index, const Vertex& p_vtx) 
109
    {
110
      CatmullClark1Ring3fa r0,r1,r2;
111
      irreg_patch.ring[index].subdivide(r0);
112
      r0.subdivide(r1);
113
      r1.subdivide(r2);
114
      return madd(8.0f/3.0f,r2.getSecondLimitTangent(),p_vtx);
115
    }
116
    
117
    void initFaceVertex(const CatmullClarkPatch& irreg_patch, 
118
			const size_t index, 
119
			const Vertex& p_vtx, 
120
                        const Vertex& e0_p_vtx, 
121
			const Vertex& e1_m_vtx, 
122
			const unsigned int face_valence_p1,
123
 			const Vertex& e0_m_vtx,	
124
			const Vertex& e3_p_vtx,	
125
			const unsigned int face_valence_p3,
126
			Vertex& f_p_vtx, 
127
			Vertex& f_m_vtx)
128
    {
129
      const unsigned int face_valence = irreg_patch.ring[index].face_valence;
130
      const unsigned int edge_valence = irreg_patch.ring[index].edge_valence;
131
      const unsigned int border_index = irreg_patch.ring[index].border_index;
132
      
133
      const Vertex& vtx     = irreg_patch.ring[index].vtx;
134
      const Vertex e_i      = irreg_patch.ring[index].getEdgeCenter(0);
135
      const Vertex c_i_m_1  = irreg_patch.ring[index].getQuadCenter(0);
136
      const Vertex e_i_m_1  = irreg_patch.ring[index].getEdgeCenter(1);
137
      
138
      Vertex c_i, e_i_p_1;
139
      const bool hasHardEdge0 =
140
        std::isinf(irreg_patch.ring[index].vertex_crease_weight) &&
141
        std::isinf(irreg_patch.ring[index].crease_weight[0]);
142
                
143
      if (unlikely((border_index == edge_valence-2) || hasHardEdge0))
144
      {
145
        /* mirror quad center and edge mid-point */
146
        c_i     = madd(2.0f, e_i - c_i_m_1, c_i_m_1);
147
        e_i_p_1 = madd(2.0f, vtx - e_i_m_1, e_i_m_1);
148
      }
149
      else
150
      {
151
        c_i     = irreg_patch.ring[index].getQuadCenter( face_valence-1 );
152
        e_i_p_1 = irreg_patch.ring[index].getEdgeCenter( face_valence-1 );
153
      }
154
      
155
      Vertex c_i_m_2, e_i_m_2;
156
      const bool hasHardEdge1 =
157
        std::isinf(irreg_patch.ring[index].vertex_crease_weight) &&
158
        std::isinf(irreg_patch.ring[index].crease_weight[1]);
159
      
160
      if (unlikely(border_index == 2 || hasHardEdge1))
161
      {
162
        /* mirror quad center and edge mid-point */
163
        c_i_m_2  = madd(2.0f, e_i_m_1 - c_i_m_1, c_i_m_1);
164
        e_i_m_2  = madd(2.0f, vtx - e_i, + e_i);
165
      }
166
      else
167
      {
168
        c_i_m_2  = irreg_patch.ring[index].getQuadCenter( 1 );
169
        e_i_m_2  = irreg_patch.ring[index].getEdgeCenter( 2 );
170
      }      
171
      
172
      const float d = 3.0f;
173
      //const float c     = cosf(2.0f*M_PI/(float)face_valence);
174
      //const float c_e_p = cosf(2.0f*M_PI/(float)face_valence_p1);
175
      //const float c_e_m = cosf(2.0f*M_PI/(float)face_valence_p3);
176
      
177
      const float c     = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence);
178
      const float c_e_p = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence_p1);
179
      const float c_e_m = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence_p3);
180

181
      const Vertex r_e_p = 1.0f/3.0f * (e_i_m_1 - e_i_p_1) + 2.0f/3.0f * (c_i_m_1 - c_i);
182
      const Vertex r_e_m = 1.0f/3.0f * (e_i     - e_i_m_2) + 2.0f/3.0f * (c_i_m_1 - c_i_m_2);
183

184
      f_p_vtx = 1.0f / d * (c_e_p * p_vtx + (d - 2.0f*c - c_e_p) * e0_p_vtx + 2.0f*c* e1_m_vtx + r_e_p);      
185
      f_m_vtx = 1.0f / d * (c_e_m * p_vtx + (d - 2.0f*c - c_e_m) * e0_m_vtx + 2.0f*c* e3_p_vtx + r_e_m);     
186
    }
187

188
    __noinline void init(const CatmullClarkPatch& patch)
189
    {
190
      assert( patch.ring[0].hasValidPositions() );
191
      assert( patch.ring[1].hasValidPositions() );
192
      assert( patch.ring[2].hasValidPositions() );
193
      assert( patch.ring[3].hasValidPositions() );
194
      
195
      p0() = initCornerVertex(patch,0);
196
      p1() = initCornerVertex(patch,1);
197
      p2() = initCornerVertex(patch,2);
198
      p3() = initCornerVertex(patch,3);
199

200
      e0_p() = initPositiveEdgeVertex(patch,0, p0());
201
      e1_p() = initPositiveEdgeVertex(patch,1, p1());
202
      e2_p() = initPositiveEdgeVertex(patch,2, p2());
203
      e3_p() = initPositiveEdgeVertex(patch,3, p3());
204

205
      e0_m() = initNegativeEdgeVertex(patch,0, p0());
206
      e1_m() = initNegativeEdgeVertex(patch,1, p1());
207
      e2_m() = initNegativeEdgeVertex(patch,2, p2());
208
      e3_m() = initNegativeEdgeVertex(patch,3, p3());
209

210
      const unsigned int face_valence_p0 = patch.ring[0].face_valence;
211
      const unsigned int face_valence_p1 = patch.ring[1].face_valence;
212
      const unsigned int face_valence_p2 = patch.ring[2].face_valence;
213
      const unsigned int face_valence_p3 = patch.ring[3].face_valence;
214
      
215
      initFaceVertex(patch,0,p0(),e0_p(),e1_m(),face_valence_p1,e0_m(),e3_p(),face_valence_p3,f0_p(),f0_m() );
216
      initFaceVertex(patch,1,p1(),e1_p(),e2_m(),face_valence_p2,e1_m(),e0_p(),face_valence_p0,f1_p(),f1_m() );
217
      initFaceVertex(patch,2,p2(),e2_p(),e3_m(),face_valence_p3,e2_m(),e1_p(),face_valence_p1,f2_p(),f2_m() );
218
      initFaceVertex(patch,3,p3(),e3_p(),e0_m(),face_valence_p0,e3_m(),e2_p(),face_valence_p3,f3_p(),f3_m() );
219

220
    }
221

222
    __noinline void init_crackfix(const CatmullClarkPatch& patch, 
223
                                  const BezierCurve* border0, 
224
                                  const BezierCurve* border1,
225
                                  const BezierCurve* border2, 
226
                                  const BezierCurve* border3)
227
    {
228
      assert( patch.ring[0].hasValidPositions() );
229
      assert( patch.ring[1].hasValidPositions() );
230
      assert( patch.ring[2].hasValidPositions() );
231
      assert( patch.ring[3].hasValidPositions() );
232
      
233
      p0() = initCornerVertex(patch,0);
234
      p1() = initCornerVertex(patch,1);
235
      p2() = initCornerVertex(patch,2);
236
      p3() = initCornerVertex(patch,3);
237

238
      e0_p() = initPositiveEdgeVertex(patch,0, p0());
239
      e1_p() = initPositiveEdgeVertex(patch,1, p1());
240
      e2_p() = initPositiveEdgeVertex(patch,2, p2());
241
      e3_p() = initPositiveEdgeVertex(patch,3, p3());
242

243
      e0_m() = initNegativeEdgeVertex(patch,0, p0());
244
      e1_m() = initNegativeEdgeVertex(patch,1, p1());
245
      e2_m() = initNegativeEdgeVertex(patch,2, p2());
246
      e3_m() = initNegativeEdgeVertex(patch,3, p3());
247

248
      if (unlikely(border0 != nullptr)) 
249
      {         
250
        p0()   = border0->v0;
251
        e0_p() = border0->v1; 
252
        e1_m() = border0->v2; 
253
        p1()   = border0->v3;
254
      }
255
      
256
      if (unlikely(border1 != nullptr))
257
      {          
258
        p1()   = border1->v0; 
259
        e1_p() = border1->v1; 
260
        e2_m() = border1->v2; 
261
        p2()   = border1->v3; 
262
      }
263

264
      if (unlikely(border2 != nullptr))
265
      {          
266
        p2()   = border2->v0; 
267
        e2_p() = border2->v1; 
268
        e3_m() = border2->v2; 
269
        p3()   = border2->v3; 
270
      }
271

272
      if (unlikely(border3 != nullptr))
273
      {          
274
        p3()   = border3->v0; 
275
        e3_p() = border3->v1; 
276
        e0_m() = border3->v2; 
277
        p0()   = border3->v3; 
278
      }
279

280
      const unsigned int face_valence_p0 = patch.ring[0].face_valence;
281
      const unsigned int face_valence_p1 = patch.ring[1].face_valence;
282
      const unsigned int face_valence_p2 = patch.ring[2].face_valence;
283
      const unsigned int face_valence_p3 = patch.ring[3].face_valence;
284
      
285
      initFaceVertex(patch,0,p0(),e0_p(),e1_m(),face_valence_p1,e0_m(),e3_p(),face_valence_p3,f0_p(),f0_m() );
286
      initFaceVertex(patch,1,p1(),e1_p(),e2_m(),face_valence_p2,e1_m(),e0_p(),face_valence_p0,f1_p(),f1_m() );
287
      initFaceVertex(patch,2,p2(),e2_p(),e3_m(),face_valence_p3,e2_m(),e1_p(),face_valence_p1,f2_p(),f2_m() );
288
      initFaceVertex(patch,3,p3(),e3_p(),e0_m(),face_valence_p0,e3_m(),e2_p(),face_valence_p3,f3_p(),f3_m() );
289
    }
290

291
    
292
    void computeGregoryPatchFacePoints(const unsigned int face_valence,
293
				       const Vertex& r_e_p, 
294
				       const Vertex& r_e_m, 					 
295
				       const Vertex& p_vtx, 
296
				       const Vertex& e0_p_vtx, 
297
				       const Vertex& e1_m_vtx, 
298
				       const unsigned int face_valence_p1,
299
				       const Vertex& e0_m_vtx,	
300
				       const Vertex& e3_p_vtx,	
301
				       const unsigned int face_valence_p3,
302
				       Vertex& f_p_vtx, 
303
				       Vertex& f_m_vtx,
304
                                       const float d = 3.0f)
305
    {
306
      //const float c     = cosf(2.0*M_PI/(float)face_valence);
307
      //const float c_e_p = cosf(2.0*M_PI/(float)face_valence_p1);
308
      //const float c_e_m = cosf(2.0*M_PI/(float)face_valence_p3);
309

310
      const float c     = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence);
311
      const float c_e_p = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence_p1);
312
      const float c_e_m = CatmullClarkPrecomputedCoefficients::table.cos_2PI_div_n(face_valence_p3);
313

314

315
      f_p_vtx = 1.0f / d * (c_e_p * p_vtx + (d - 2.0f*c - c_e_p) * e0_p_vtx + 2.0f*c* e1_m_vtx + r_e_p);      
316
      f_m_vtx = 1.0f / d * (c_e_m * p_vtx + (d - 2.0f*c - c_e_m) * e0_m_vtx + 2.0f*c* e3_p_vtx + r_e_m);      
317
      f_p_vtx = 1.0f / d * (c_e_p * p_vtx + (d - 2.0f*c - c_e_p) * e0_p_vtx + 2.0f*c* e1_m_vtx + r_e_p);      
318
      f_m_vtx = 1.0f / d * (c_e_m * p_vtx + (d - 2.0f*c - c_e_m) * e0_m_vtx + 2.0f*c* e3_p_vtx + r_e_m);
319
    }
320

321
    __noinline void init(const GeneralCatmullClarkPatch& patch)
322
    {
323
      assert(patch.size() == 4);
324
#if 0
325
      CatmullClarkPatch qpatch; patch.init(qpatch);
326
      init(qpatch);
327
#else
328
      const float face_valence_p0 = patch.ring[0].face_valence;
329
      const float face_valence_p1 = patch.ring[1].face_valence;
330
      const float face_valence_p2 = patch.ring[2].face_valence;
331
      const float face_valence_p3 = patch.ring[3].face_valence;
332

333
      Vertex p0_r_p, p0_r_m;
334
      patch.ring[0].computeGregoryPatchEdgePoints( p0(), e0_p(), e0_m(), p0_r_p, p0_r_m );
335

336
      Vertex p1_r_p, p1_r_m;
337
      patch.ring[1].computeGregoryPatchEdgePoints( p1(), e1_p(), e1_m(), p1_r_p, p1_r_m );
338
      
339
      Vertex p2_r_p, p2_r_m;
340
      patch.ring[2].computeGregoryPatchEdgePoints( p2(), e2_p(), e2_m(), p2_r_p, p2_r_m );
341

342
      Vertex p3_r_p, p3_r_m;
343
      patch.ring[3].computeGregoryPatchEdgePoints( p3(), e3_p(), e3_m(), p3_r_p, p3_r_m );
344

345
      computeGregoryPatchFacePoints(face_valence_p0, p0_r_p, p0_r_m, p0(), e0_p(), e1_m(), face_valence_p1, e0_m(), e3_p(), face_valence_p3, f0_p(), f0_m() );
346
      computeGregoryPatchFacePoints(face_valence_p1, p1_r_p, p1_r_m, p1(), e1_p(), e2_m(), face_valence_p2, e1_m(), e0_p(), face_valence_p0, f1_p(), f1_m() );
347
      computeGregoryPatchFacePoints(face_valence_p2, p2_r_p, p2_r_m, p2(), e2_p(), e3_m(), face_valence_p3, e2_m(), e1_p(), face_valence_p1, f2_p(), f2_m() );
348
      computeGregoryPatchFacePoints(face_valence_p3, p3_r_p, p3_r_m, p3(), e3_p(), e0_m(), face_valence_p0, e3_m(), e2_p(), face_valence_p3, f3_p(), f3_m() );
349

350
#endif
351
    }
352
   
353
    
354
    __forceinline void convert_to_bezier()
355
    {
356
      f0_p() = (f0_p() + f0_m()) * 0.5f;
357
      f1_p() = (f1_p() + f1_m()) * 0.5f;
358
      f2_p() = (f2_p() + f2_m()) * 0.5f;
359
      f3_p() = (f3_p() + f3_m()) * 0.5f;
360
      f0_m() = Vertex( zero );
361
      f1_m() = Vertex( zero );
362
      f2_m() = Vertex( zero );
363
      f3_m() = Vertex( zero );      
364
    }
365
    
366
    static __forceinline void computeInnerVertices(const Vertex matrix[4][4], const Vertex f_m[2][2], const float uu, const float vv,
367
						   Vertex_t& matrix_11, Vertex_t& matrix_12, Vertex_t& matrix_22, Vertex_t& matrix_21)
368
    {
369
      if (unlikely(uu == 0.0f || uu == 1.0f || vv == 0.0f || vv == 1.0f)) 
370
      {
371
	matrix_11 = matrix[1][1];
372
	matrix_12 = matrix[1][2];
373
	matrix_22 = matrix[2][2];
374
	matrix_21 = matrix[2][1];	 
375
      }
376
      else
377
      {
378
	const Vertex_t f0_p = matrix[1][1];
379
	const Vertex_t f1_p = matrix[1][2];
380
	const Vertex_t f2_p = matrix[2][2];
381
	const Vertex_t f3_p = matrix[2][1];
382
        
383
	const Vertex_t f0_m = f_m[0][0];
384
	const Vertex_t f1_m = f_m[0][1];
385
	const Vertex_t f2_m = f_m[1][1];
386
	const Vertex_t f3_m = f_m[1][0];
387
        
388
	matrix_11 = (      uu  * f0_p +       vv  * f0_m)*rcp(uu+vv);
389
	matrix_12 = ((1.0f-uu) * f1_m +       vv  * f1_p)*rcp(1.0f-uu+vv);
390
	matrix_22 = ((1.0f-uu) * f2_p + (1.0f-vv) * f2_m)*rcp(2.0f-uu-vv);
391
	matrix_21 = (      uu  * f3_m + (1.0f-vv) * f3_p)*rcp(1.0f+uu-vv);
392
      }
393
    } 
394

395
    template<typename vfloat>
396
    static __forceinline void computeInnerVertices(const Vertex v[4][4], const Vertex f[2][2], 
397
                                                   size_t i, const vfloat& uu, const vfloat& vv, vfloat& matrix_11, vfloat& matrix_12, vfloat& matrix_22, vfloat& matrix_21) 
398
    {
399
      const auto m_border = (uu == 0.0f) | (uu == 1.0f) | (vv == 0.0f) | (vv == 1.0f);
400

401
      const vfloat f0_p = v[1][1][i];
402
      const vfloat f1_p = v[1][2][i];
403
      const vfloat f2_p = v[2][2][i];
404
      const vfloat f3_p = v[2][1][i];
405
      
406
      const vfloat f0_m = f[0][0][i];
407
      const vfloat f1_m = f[0][1][i];
408
      const vfloat f2_m = f[1][1][i];
409
      const vfloat f3_m = f[1][0][i];
410
      
411
      const vfloat one_minus_uu = vfloat(1.0f) - uu;
412
      const vfloat one_minus_vv = vfloat(1.0f) - vv;      
413
      
414
      const vfloat f0_i = (          uu * f0_p +           vv * f0_m) * rcp(uu+vv);
415
      const vfloat f1_i = (one_minus_uu * f1_m +           vv * f1_p) * rcp(one_minus_uu+vv);
416
      const vfloat f2_i = (one_minus_uu * f2_p + one_minus_vv * f2_m) * rcp(one_minus_uu+one_minus_vv);
417
      const vfloat f3_i = (          uu * f3_m + one_minus_vv * f3_p) * rcp(uu+one_minus_vv);
418
      
419
      matrix_11 = select(m_border,f0_p,f0_i);
420
      matrix_12 = select(m_border,f1_p,f1_i);
421
      matrix_22 = select(m_border,f2_p,f2_i);
422
      matrix_21 = select(m_border,f3_p,f3_i);
423
    }
424

425
    static __forceinline Vertex eval(const Vertex matrix[4][4], const Vertex f[2][2], const float& uu, const float& vv) 
426
    {
427
      Vertex_t v_11, v_12, v_22, v_21;
428
      computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
429
      
430
      const Vec4<float> Bu = BezierBasis::eval(uu);
431
      const Vec4<float> Bv = BezierBasis::eval(vv);
432
      
433
      return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))), 
434
                  madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11        ,madd(Bu.z,v_12        ,Bu.w * matrix[1][3]))), 
435
                       madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21        ,madd(Bu.z,v_22        ,Bu.w * matrix[2][3]))), 
436
                            Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3])))))); 
437
    }
438

439
    static __forceinline Vertex eval_du(const Vertex matrix[4][4], const Vertex f[2][2], const float uu, const float vv) // approximative derivative
440
    {
441
      Vertex_t v_11, v_12, v_22, v_21;
442
      computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
443
      
444
      const Vec4<float> Bu = BezierBasis::derivative(uu);
445
      const Vec4<float> Bv = BezierBasis::eval(vv);
446

447
      return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))), 
448
                  madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11        ,madd(Bu.z,v_12        ,Bu.w * matrix[1][3]))), 
449
                       madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21        ,madd(Bu.z,v_22        ,Bu.w * matrix[2][3]))), 
450
                            Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3])))))); 
451
    }
452

453
    static __forceinline Vertex eval_dv(const Vertex matrix[4][4], const Vertex f[2][2], const float uu, const float vv) // approximative derivative
454
    {
455
      Vertex_t v_11, v_12, v_22, v_21;
456
      computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
457
      
458
      const Vec4<float> Bu = BezierBasis::eval(uu);
459
      const Vec4<float> Bv = BezierBasis::derivative(vv);
460
 
461
      return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))), 
462
                  madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11        ,madd(Bu.z,v_12        ,Bu.w * matrix[1][3]))), 
463
                       madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21        ,madd(Bu.z,v_22        ,Bu.w * matrix[2][3]))), 
464
                            Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3])))))); 
465
    }
466

467
    static __forceinline Vertex eval_dudu(const Vertex matrix[4][4], const Vertex f[2][2], const float uu, const float vv) // approximative derivative
468
    {
469
      Vertex_t v_11, v_12, v_22, v_21;
470
      computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
471
      
472
      const Vec4<float> Bu = BezierBasis::derivative2(uu);
473
      const Vec4<float> Bv = BezierBasis::eval(vv);
474
 
475
      return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))), 
476
                  madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11        ,madd(Bu.z,v_12        ,Bu.w * matrix[1][3]))), 
477
                       madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21        ,madd(Bu.z,v_22        ,Bu.w * matrix[2][3]))), 
478
                            Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3])))))); 
479
     }
480

481
    static __forceinline Vertex eval_dvdv(const Vertex matrix[4][4], const Vertex f[2][2], const float uu, const float vv) // approximative derivative
482
    {
483
      Vertex_t v_11, v_12, v_22, v_21;
484
      computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
485
      
486
      const Vec4<float> Bu = BezierBasis::eval(uu);
487
      const Vec4<float> Bv = BezierBasis::derivative2(vv);
488

489
      return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))), 
490
                  madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11        ,madd(Bu.z,v_12        ,Bu.w * matrix[1][3]))), 
491
                       madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21        ,madd(Bu.z,v_22        ,Bu.w * matrix[2][3]))), 
492
                            Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3])))))); 
493
    }
494

495
    static __forceinline Vertex eval_dudv(const Vertex matrix[4][4], const Vertex f[2][2], const float uu, const float vv) // approximative derivative
496
    {
497
      Vertex_t v_11, v_12, v_22, v_21;
498
      computeInnerVertices(matrix,f,uu,vv,v_11, v_12, v_22, v_21);
499
      
500
      const Vec4<float> Bu = BezierBasis::derivative(uu);
501
      const Vec4<float> Bv = BezierBasis::derivative(vv);
502

503
      return madd(Bv.x,madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))), 
504
                  madd(Bv.y,madd(Bu.x,matrix[1][0],madd(Bu.y,v_11        ,madd(Bu.z,v_12        ,Bu.w * matrix[1][3]))), 
505
                       madd(Bv.z,madd(Bu.x,matrix[2][0],madd(Bu.y,v_21        ,madd(Bu.z,v_22        ,Bu.w * matrix[2][3]))), 
506
                            Bv.w*madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3])))))); 
507
    }
508

509
    __forceinline Vertex eval(const float uu, const float vv) const {
510
      return eval(v,f,uu,vv);
511
    }
512

513
    __forceinline Vertex eval_du( const float uu, const float vv) const {
514
      return eval_du(v,f,uu,vv);
515
    }
516

517
    __forceinline Vertex eval_dv( const float uu, const float vv) const {
518
      return eval_dv(v,f,uu,vv);
519
    }
520

521
    __forceinline Vertex eval_dudu( const float uu, const float vv) const {
522
      return eval_dudu(v,f,uu,vv);
523
    }
524

525
    __forceinline Vertex eval_dvdv( const float uu, const float vv) const {
526
      return eval_dvdv(v,f,uu,vv);
527
    }
528

529
    __forceinline Vertex eval_dudv( const float uu, const float vv) const {
530
      return eval_dudv(v,f,uu,vv);
531
    }
532

533
    static __forceinline Vertex normal(const Vertex matrix[4][4], const Vertex f_m[2][2], const float uu, const float vv)  // FIXME: why not using basis functions
534
    {
535
      /* interpolate inner vertices */
536
      Vertex_t matrix_11, matrix_12, matrix_22, matrix_21;
537
      computeInnerVertices(matrix,f_m,uu,vv,matrix_11, matrix_12, matrix_22, matrix_21);
538
      
539
      /* tangentU */
540
      const Vertex_t col0 = deCasteljau(vv, (Vertex_t)matrix[0][0], (Vertex_t)matrix[1][0], (Vertex_t)matrix[2][0], (Vertex_t)matrix[3][0]);
541
      const Vertex_t col1 = deCasteljau(vv, (Vertex_t)matrix[0][1], (Vertex_t)matrix_11   , (Vertex_t)matrix_21   , (Vertex_t)matrix[3][1]);
542
      const Vertex_t col2 = deCasteljau(vv, (Vertex_t)matrix[0][2], (Vertex_t)matrix_12   , (Vertex_t)matrix_22   , (Vertex_t)matrix[3][2]);
543
      const Vertex_t col3 = deCasteljau(vv, (Vertex_t)matrix[0][3], (Vertex_t)matrix[1][3], (Vertex_t)matrix[2][3], (Vertex_t)matrix[3][3]);
544
      
545
      const Vertex_t tangentU = deCasteljau_tangent(uu, col0, col1, col2, col3);
546
      
547
      /* tangentV */
548
      const Vertex_t row0 = deCasteljau(uu, (Vertex_t)matrix[0][0], (Vertex_t)matrix[0][1], (Vertex_t)matrix[0][2], (Vertex_t)matrix[0][3]);
549
      const Vertex_t row1 = deCasteljau(uu, (Vertex_t)matrix[1][0], (Vertex_t)matrix_11   , (Vertex_t)matrix_12   , (Vertex_t)matrix[1][3]);
550
      const Vertex_t row2 = deCasteljau(uu, (Vertex_t)matrix[2][0], (Vertex_t)matrix_21   , (Vertex_t)matrix_22   , (Vertex_t)matrix[2][3]);
551
      const Vertex_t row3 = deCasteljau(uu, (Vertex_t)matrix[3][0], (Vertex_t)matrix[3][1], (Vertex_t)matrix[3][2], (Vertex_t)matrix[3][3]);
552
      
553
      const Vertex_t tangentV = deCasteljau_tangent(vv, row0, row1, row2, row3);
554
      
555
      /* normal = tangentU x tangentV */
556
      const Vertex_t n = cross(tangentU,tangentV);
557
      
558
      return n;     
559
    }
560
   
561
    __forceinline Vertex normal( const float uu, const float vv) const {
562
      return normal(v,f,uu,vv);
563
    }    
564
    
565
    __forceinline void eval(const float u, const float v, 
566
                            Vertex* P, Vertex* dPdu, Vertex* dPdv, 
567
                            Vertex* ddPdudu, Vertex* ddPdvdv, Vertex* ddPdudv,
568
                            const float dscale = 1.0f) const
569
    {
570
      if (P) {
571
        *P = eval(u,v); 
572
      }
573
      if (dPdu) {
574
        assert(dPdu); *dPdu = eval_du(u,v)*dscale; 
575
        assert(dPdv); *dPdv = eval_dv(u,v)*dscale; 
576
      }
577
      if (ddPdudu) {
578
        assert(ddPdudu); *ddPdudu = eval_dudu(u,v)*sqr(dscale); 
579
        assert(ddPdvdv); *ddPdvdv = eval_dvdv(u,v)*sqr(dscale); 
580
        assert(ddPdudv); *ddPdudv = eval_dudv(u,v)*sqr(dscale); 
581
      }
582
    }
583

584
    template<class vfloat>
585
    static __forceinline vfloat eval(const Vertex v[4][4], const Vertex f[2][2], 
586
                                     const size_t i, const vfloat& uu, const vfloat& vv, const Vec4<vfloat>& u_n, const Vec4<vfloat>& v_n,
587
                                     vfloat& matrix_11, vfloat& matrix_12, vfloat& matrix_22, vfloat& matrix_21)
588
    {
589
      const vfloat curve0_x = madd(v_n[0],vfloat(v[0][0][i]),madd(v_n[1],vfloat(v[1][0][i]),madd(v_n[2],vfloat(v[2][0][i]),v_n[3] * vfloat(v[3][0][i]))));
590
      const vfloat curve1_x = madd(v_n[0],vfloat(v[0][1][i]),madd(v_n[1],vfloat(matrix_11 ),madd(v_n[2],vfloat(matrix_21 ),v_n[3] * vfloat(v[3][1][i]))));
591
      const vfloat curve2_x = madd(v_n[0],vfloat(v[0][2][i]),madd(v_n[1],vfloat(matrix_12 ),madd(v_n[2],vfloat(matrix_22 ),v_n[3] * vfloat(v[3][2][i]))));
592
      const vfloat curve3_x = madd(v_n[0],vfloat(v[0][3][i]),madd(v_n[1],vfloat(v[1][3][i]),madd(v_n[2],vfloat(v[2][3][i]),v_n[3] * vfloat(v[3][3][i]))));
593
      return madd(u_n[0],curve0_x,madd(u_n[1],curve1_x,madd(u_n[2],curve2_x,u_n[3] * curve3_x)));
594
    }
595
    
596
    template<typename vbool, typename vfloat>
597
    static __forceinline void eval(const Vertex v[4][4], const Vertex f[2][2], 
598
                                   const vbool& valid, const vfloat& uu, const vfloat& vv, 
599
                                   float* P, float* dPdu, float* dPdv, float* ddPdudu, float* ddPdvdv, float* ddPdudv,
600
                                   const float dscale, const size_t dstride, const size_t N) 
601
    {
602
      if (P) {
603
        const Vec4<vfloat> u_n = BezierBasis::eval(uu); 
604
        const Vec4<vfloat> v_n = BezierBasis::eval(vv); 
605
        for (size_t i=0; i<N; i++) {
606
          vfloat matrix_11, matrix_12, matrix_22, matrix_21;
607
          computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21); // FIXME: calculated multiple times
608
          vfloat::store(valid,P+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21));
609
        }
610
      }
611
      if (dPdu)
612
      {
613
        {
614
          assert(dPdu);
615
          const Vec4<vfloat> u_n = BezierBasis::derivative(uu); 
616
          const Vec4<vfloat> v_n = BezierBasis::eval(vv);
617
          for (size_t i=0; i<N; i++) {
618
            vfloat matrix_11, matrix_12, matrix_22, matrix_21;
619
            computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21);  // FIXME: calculated multiple times
620
            vfloat::store(valid,dPdu+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21)*dscale);
621
          }
622
        }
623
        {
624
          assert(dPdv);
625
          const Vec4<vfloat> u_n = BezierBasis::eval(uu); 
626
          const Vec4<vfloat> v_n = BezierBasis::derivative(vv);
627
          for (size_t i=0; i<N; i++) {
628
            vfloat matrix_11, matrix_12, matrix_22, matrix_21;
629
            computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21);  // FIXME: calculated multiple times
630
            vfloat::store(valid,dPdv+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21)*dscale);
631
          }
632
        }
633
      }
634
      if (ddPdudu)
635
      {
636
        {
637
          assert(ddPdudu);
638
          const Vec4<vfloat> u_n = BezierBasis::derivative2(uu); 
639
          const Vec4<vfloat> v_n = BezierBasis::eval(vv);
640
          for (size_t i=0; i<N; i++) {
641
            vfloat matrix_11, matrix_12, matrix_22, matrix_21;
642
            computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21);  // FIXME: calculated multiple times
643
            vfloat::store(valid,ddPdudu+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21)*sqr(dscale));
644
          }
645
        }
646
        {
647
          assert(ddPdvdv);
648
          const Vec4<vfloat> u_n = BezierBasis::eval(uu); 
649
          const Vec4<vfloat> v_n = BezierBasis::derivative2(vv);
650
          for (size_t i=0; i<N; i++) {
651
            vfloat matrix_11, matrix_12, matrix_22, matrix_21;
652
            computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21);  // FIXME: calculated multiple times
653
            vfloat::store(valid,ddPdvdv+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21)*sqr(dscale));
654
          }
655
        }
656
        {
657
          assert(ddPdudv);
658
          const Vec4<vfloat> u_n = BezierBasis::derivative(uu); 
659
          const Vec4<vfloat> v_n = BezierBasis::derivative(vv);
660
          for (size_t i=0; i<N; i++) {
661
            vfloat matrix_11, matrix_12, matrix_22, matrix_21;
662
            computeInnerVertices(v,f,i,uu,vv,matrix_11,matrix_12,matrix_22,matrix_21);  // FIXME: calculated multiple times
663
            vfloat::store(valid,ddPdudv+i*dstride,eval(v,f,i,uu,vv,u_n,v_n,matrix_11,matrix_12,matrix_22,matrix_21)*sqr(dscale));
664
          }
665
        }
666
      }
667
    }
668

669
    template<typename vbool, typename vfloat>
670
    __forceinline void eval(const vbool& valid, const vfloat& uu, const vfloat& vv, 
671
                            float* P, float* dPdu, float* dPdv, float* ddPdudu, float* ddPdvdv, float* ddPdudv,
672
                            const float dscale, const size_t dstride, const size_t N) const {
673
      eval(v,f,valid,uu,vv,P,dPdu,dPdv,ddPdudu,ddPdvdv,ddPdudv,dscale,dstride,N);
674
    }
675

676
    template<class T>
677
      static __forceinline Vec3<T> eval_t(const Vertex matrix[4][4], const Vec3<T> f[2][2], const T& uu, const T& vv) 
678
    {
679
      typedef typename T::Bool M;
680
      const M m_border = (uu == 0.0f) | (uu == 1.0f) | (vv == 0.0f) | (vv == 1.0f);
681

682
      const Vec3<T> f0_p = Vec3<T>(matrix[1][1].x,matrix[1][1].y,matrix[1][1].z);
683
      const Vec3<T> f1_p = Vec3<T>(matrix[1][2].x,matrix[1][2].y,matrix[1][2].z);
684
      const Vec3<T> f2_p = Vec3<T>(matrix[2][2].x,matrix[2][2].y,matrix[2][2].z);
685
      const Vec3<T> f3_p = Vec3<T>(matrix[2][1].x,matrix[2][1].y,matrix[2][1].z);
686
      
687
      const Vec3<T> f0_m = f[0][0];
688
      const Vec3<T> f1_m = f[0][1];
689
      const Vec3<T> f2_m = f[1][1];
690
      const Vec3<T> f3_m = f[1][0];
691
      
692
      const T one_minus_uu = T(1.0f) - uu;
693
      const T one_minus_vv = T(1.0f) - vv;      
694
      
695
      const Vec3<T> f0_i = (          uu * f0_p +           vv * f0_m) * rcp(uu+vv);
696
      const Vec3<T> f1_i = (one_minus_uu * f1_m +           vv * f1_p) * rcp(one_minus_uu+vv);
697
      const Vec3<T> f2_i = (one_minus_uu * f2_p + one_minus_vv * f2_m) * rcp(one_minus_uu+one_minus_vv);
698
      const Vec3<T> f3_i = (          uu * f3_m + one_minus_vv * f3_p) * rcp(uu+one_minus_vv);
699
      
700
      const Vec3<T> F0( select(m_border,f0_p.x,f0_i.x), select(m_border,f0_p.y,f0_i.y), select(m_border,f0_p.z,f0_i.z) );
701
      const Vec3<T> F1( select(m_border,f1_p.x,f1_i.x), select(m_border,f1_p.y,f1_i.y), select(m_border,f1_p.z,f1_i.z) );
702
      const Vec3<T> F2( select(m_border,f2_p.x,f2_i.x), select(m_border,f2_p.y,f2_i.y), select(m_border,f2_p.z,f2_i.z) );
703
      const Vec3<T> F3( select(m_border,f3_p.x,f3_i.x), select(m_border,f3_p.y,f3_i.y), select(m_border,f3_p.z,f3_i.z) );
704

705
      const T B0_u = one_minus_uu * one_minus_uu * one_minus_uu;
706
      const T B0_v = one_minus_vv * one_minus_vv * one_minus_vv;
707
      const T B1_u = 3.0f * (one_minus_uu * uu * one_minus_uu);
708
      const T B1_v = 3.0f * (one_minus_vv * vv * one_minus_vv);
709
      const T B2_u = 3.0f * (uu * one_minus_uu * uu);
710
      const T B2_v = 3.0f * (vv * one_minus_vv * vv);
711
      const T B3_u = uu * uu * uu;
712
      const T B3_v = vv * vv * vv;
713

714
      const T x = madd(B0_v,madd(B0_u,matrix[0][0].x,madd(B1_u,matrix[0][1].x,madd(B2_u,matrix[0][2].x,B3_u * matrix[0][3].x))), 
715
                  madd(B1_v,madd(B0_u,matrix[1][0].x,madd(B1_u,F0.x          ,madd(B2_u,F1.x          ,B3_u * matrix[1][3].x))), 
716
                  madd(B2_v,madd(B0_u,matrix[2][0].x,madd(B1_u,F3.x          ,madd(B2_u,F2.x          ,B3_u * matrix[2][3].x))), 
717
                       B3_v*madd(B0_u,matrix[3][0].x,madd(B1_u,matrix[3][1].x,madd(B2_u,matrix[3][2].x,B3_u * matrix[3][3].x)))))); 
718

719
      const T y = madd(B0_v,madd(B0_u,matrix[0][0].y,madd(B1_u,matrix[0][1].y,madd(B2_u,matrix[0][2].y,B3_u * matrix[0][3].y))),
720
                  madd(B1_v,madd(B0_u,matrix[1][0].y,madd(B1_u,F0.y          ,madd(B2_u,F1.y          ,B3_u * matrix[1][3].y))),
721
                  madd(B2_v,madd(B0_u,matrix[2][0].y,madd(B1_u,F3.y          ,madd(B2_u,F2.y          ,B3_u * matrix[2][3].y))),
722
                       B3_v*madd(B0_u,matrix[3][0].y,madd(B1_u,matrix[3][1].y,madd(B2_u,matrix[3][2].y,B3_u * matrix[3][3].y))))));
723
      
724
      const T z = madd(B0_v,madd(B0_u,matrix[0][0].z,madd(B1_u,matrix[0][1].z,madd(B2_u,matrix[0][2].z,B3_u * matrix[0][3].z))),
725
                  madd(B1_v,madd(B0_u,matrix[1][0].z,madd(B1_u,F0.z          ,madd(B2_u,F1.z          ,B3_u * matrix[1][3].z))),
726
                  madd(B2_v,madd(B0_u,matrix[2][0].z,madd(B1_u,F3.z          ,madd(B2_u,F2.z          ,B3_u * matrix[2][3].z))),
727
                       B3_v*madd(B0_u,matrix[3][0].z,madd(B1_u,matrix[3][1].z,madd(B2_u,matrix[3][2].z,B3_u * matrix[3][3].z))))));
728
      
729
      return Vec3<T>(x,y,z);
730
    }
731

732
    template<class T>
733
    __forceinline Vec3<T> eval(const T& uu, const T& vv) const 
734
    {
735
      Vec3<T> ff[2][2];
736
      ff[0][0] = Vec3<T>(f[0][0]);
737
      ff[0][1] = Vec3<T>(f[0][1]);
738
      ff[1][1] = Vec3<T>(f[1][1]);
739
      ff[1][0] = Vec3<T>(f[1][0]);
740
      return eval_t(v,ff,uu,vv);
741
    }
742

743
    template<class T>
744
      static __forceinline Vec3<T> normal_t(const Vertex matrix[4][4], const Vec3<T> f[2][2], const T& uu, const T& vv) 
745
    {
746
      typedef typename T::Bool M;
747
      
748
      const Vec3<T> f0_p = Vec3<T>(matrix[1][1].x,matrix[1][1].y,matrix[1][1].z);
749
      const Vec3<T> f1_p = Vec3<T>(matrix[1][2].x,matrix[1][2].y,matrix[1][2].z);
750
      const Vec3<T> f2_p = Vec3<T>(matrix[2][2].x,matrix[2][2].y,matrix[2][2].z);
751
      const Vec3<T> f3_p = Vec3<T>(matrix[2][1].x,matrix[2][1].y,matrix[2][1].z);
752

753
      const Vec3<T> f0_m = f[0][0];
754
      const Vec3<T> f1_m = f[0][1];
755
      const Vec3<T> f2_m = f[1][1];
756
      const Vec3<T> f3_m = f[1][0];
757
      
758
      const T one_minus_uu = T(1.0f) - uu;
759
      const T one_minus_vv = T(1.0f) - vv;      
760
      
761
      const Vec3<T> f0_i = (          uu * f0_p +           vv * f0_m) * rcp(uu+vv);
762
      const Vec3<T> f1_i = (one_minus_uu * f1_m +           vv * f1_p) * rcp(one_minus_uu+vv);
763
      const Vec3<T> f2_i = (one_minus_uu * f2_p + one_minus_vv * f2_m) * rcp(one_minus_uu+one_minus_vv);
764
      const Vec3<T> f3_i = (          uu * f3_m + one_minus_vv * f3_p) * rcp(uu+one_minus_vv);
765

766
#if 1
767
      const M m_corner0 = (uu == 0.0f) & (vv == 0.0f);
768
      const M m_corner1 = (uu == 1.0f) & (vv == 0.0f);
769
      const M m_corner2 = (uu == 1.0f) & (vv == 1.0f);
770
      const M m_corner3 = (uu == 0.0f) & (vv == 1.0f);      
771
      const Vec3<T> matrix_11( select(m_corner0,f0_p.x,f0_i.x), select(m_corner0,f0_p.y,f0_i.y), select(m_corner0,f0_p.z,f0_i.z) );
772
      const Vec3<T> matrix_12( select(m_corner1,f1_p.x,f1_i.x), select(m_corner1,f1_p.y,f1_i.y), select(m_corner1,f1_p.z,f1_i.z) );
773
      const Vec3<T> matrix_22( select(m_corner2,f2_p.x,f2_i.x), select(m_corner2,f2_p.y,f2_i.y), select(m_corner2,f2_p.z,f2_i.z) );
774
      const Vec3<T> matrix_21( select(m_corner3,f3_p.x,f3_i.x), select(m_corner3,f3_p.y,f3_i.y), select(m_corner3,f3_p.z,f3_i.z) );
775
#else
776
      const M m_border = (uu == 0.0f) | (uu == 1.0f) | (vv == 0.0f) | (vv == 1.0f);
777
      const Vec3<T> matrix_11( select(m_border,f0_p.x,f0_i.x), select(m_border,f0_p.y,f0_i.y), select(m_border,f0_p.z,f0_i.z) );
778
      const Vec3<T> matrix_12( select(m_border,f1_p.x,f1_i.x), select(m_border,f1_p.y,f1_i.y), select(m_border,f1_p.z,f1_i.z) );
779
      const Vec3<T> matrix_22( select(m_border,f2_p.x,f2_i.x), select(m_border,f2_p.y,f2_i.y), select(m_border,f2_p.z,f2_i.z) );
780
      const Vec3<T> matrix_21( select(m_border,f3_p.x,f3_i.x), select(m_border,f3_p.y,f3_i.y), select(m_border,f3_p.z,f3_i.z) );
781
#endif
782
      
783
      const Vec3<T> matrix_00 = Vec3<T>(matrix[0][0].x,matrix[0][0].y,matrix[0][0].z);
784
      const Vec3<T> matrix_10 = Vec3<T>(matrix[1][0].x,matrix[1][0].y,matrix[1][0].z);
785
      const Vec3<T> matrix_20 = Vec3<T>(matrix[2][0].x,matrix[2][0].y,matrix[2][0].z);
786
      const Vec3<T> matrix_30 = Vec3<T>(matrix[3][0].x,matrix[3][0].y,matrix[3][0].z);
787
      
788
      const Vec3<T> matrix_01 = Vec3<T>(matrix[0][1].x,matrix[0][1].y,matrix[0][1].z);
789
      const Vec3<T> matrix_02 = Vec3<T>(matrix[0][2].x,matrix[0][2].y,matrix[0][2].z);
790
      const Vec3<T> matrix_03 = Vec3<T>(matrix[0][3].x,matrix[0][3].y,matrix[0][3].z);
791
      
792
      const Vec3<T> matrix_31 = Vec3<T>(matrix[3][1].x,matrix[3][1].y,matrix[3][1].z);
793
      const Vec3<T> matrix_32 = Vec3<T>(matrix[3][2].x,matrix[3][2].y,matrix[3][2].z);
794
      const Vec3<T> matrix_33 = Vec3<T>(matrix[3][3].x,matrix[3][3].y,matrix[3][3].z);
795
      
796
      const Vec3<T> matrix_13 = Vec3<T>(matrix[1][3].x,matrix[1][3].y,matrix[1][3].z);
797
      const Vec3<T> matrix_23 = Vec3<T>(matrix[2][3].x,matrix[2][3].y,matrix[2][3].z);
798
      
799
      /* tangentU */
800
      const Vec3<T> col0 = deCasteljau(vv, matrix_00, matrix_10, matrix_20, matrix_30);
801
      const Vec3<T> col1 = deCasteljau(vv, matrix_01, matrix_11, matrix_21, matrix_31);
802
      const Vec3<T> col2 = deCasteljau(vv, matrix_02, matrix_12, matrix_22, matrix_32);
803
      const Vec3<T> col3 = deCasteljau(vv, matrix_03, matrix_13, matrix_23, matrix_33);
804
      
805
      const Vec3<T> tangentU = deCasteljau_tangent(uu, col0, col1, col2, col3);
806
      
807
      /* tangentV */
808
      const Vec3<T> row0 = deCasteljau(uu, matrix_00, matrix_01, matrix_02, matrix_03);
809
      const Vec3<T> row1 = deCasteljau(uu, matrix_10, matrix_11, matrix_12, matrix_13);
810
      const Vec3<T> row2 = deCasteljau(uu, matrix_20, matrix_21, matrix_22, matrix_23);
811
      const Vec3<T> row3 = deCasteljau(uu, matrix_30, matrix_31, matrix_32, matrix_33);
812
      
813
      const Vec3<T> tangentV = deCasteljau_tangent(vv, row0, row1, row2, row3);
814
      
815
      /* normal = tangentU x tangentV */
816
      const Vec3<T> n = cross(tangentU,tangentV);
817
      return n;
818
    }
819

820
     template<class T>
821
    __forceinline Vec3<T> normal(const T& uu, const T& vv) const 
822
    {
823
      Vec3<T> ff[2][2];
824
      ff[0][0] = Vec3<T>(f[0][0]);
825
      ff[0][1] = Vec3<T>(f[0][1]);
826
      ff[1][1] = Vec3<T>(f[1][1]);
827
      ff[1][0] = Vec3<T>(f[1][0]);
828
      return normal_t(v,ff,uu,vv);
829
    }
830

831
    __forceinline BBox<Vertex> bounds() const
832
    {
833
      const Vertex *const cv = &v[0][0];
834
      BBox<Vertex> bounds (cv[0]);
835
      for (size_t i=1; i<16; i++) 
836
        bounds.extend( cv[i] );
837
      bounds.extend(f[0][0]);
838
      bounds.extend(f[1][0]);
839
      bounds.extend(f[1][1]);
840
      bounds.extend(f[1][1]);
841
      return bounds;
842
    }
843
    
844
    friend embree_ostream operator<<(embree_ostream o, const GregoryPatchT& p)
845
    {
846
      for (size_t y=0; y<4; y++)
847
	for (size_t x=0; x<4; x++)
848
	  o << "v[" << y << "][" << x << "] " << p.v[y][x] << embree_endl;
849
      
850
      for (size_t y=0; y<2; y++)
851
	for (size_t x=0; x<2; x++)
852
	  o << "f[" << y << "][" << x << "] " << p.f[y][x] << embree_endl;
853
      return o;
854
    } 
855
  };
856

857
  typedef GregoryPatchT<Vec3fa,Vec3fa_t> GregoryPatch3fa;
858

859
  template<typename Vertex, typename Vertex_t>
860
    __forceinline  BezierPatchT<Vertex,Vertex_t>::BezierPatchT (const HalfEdge* edge, const char* vertices, size_t stride) 
861
  {
862
    CatmullClarkPatchT<Vertex,Vertex_t> patch(edge,vertices,stride);
863
    GregoryPatchT<Vertex,Vertex_t> gpatch(patch); 
864
    gpatch.convert_to_bezier(); 
865
    for (size_t y=0; y<4; y++)
866
      for (size_t x=0; x<4; x++)
867
        matrix[y][x] = (Vertex_t)gpatch.v[y][x];
868
  }
869
  
870
   template<typename Vertex, typename Vertex_t>
871
    __forceinline BezierPatchT<Vertex,Vertex_t>::BezierPatchT(const CatmullClarkPatchT<Vertex,Vertex_t>& patch) 
872
    {
873
      GregoryPatchT<Vertex,Vertex_t> gpatch(patch); 
874
      gpatch.convert_to_bezier(); 
875
      for (size_t y=0; y<4; y++)
876
	for (size_t x=0; x<4; x++)
877
	  matrix[y][x] = (Vertex_t)gpatch.v[y][x];
878
    }
879

880
   template<typename Vertex, typename Vertex_t>
881
     __forceinline BezierPatchT<Vertex,Vertex_t>::BezierPatchT(const CatmullClarkPatchT<Vertex,Vertex_t>& patch, 
882
                                                               const BezierCurveT<Vertex>* border0,
883
                                                               const BezierCurveT<Vertex>* border1,
884
                                                               const BezierCurveT<Vertex>* border2,
885
                                                               const BezierCurveT<Vertex>* border3) 
886
    {
887
      GregoryPatchT<Vertex,Vertex_t> gpatch(patch,border0,border1,border2,border3); 
888
      gpatch.convert_to_bezier(); 
889
      for (size_t y=0; y<4; y++)
890
	for (size_t x=0; x<4; x++)
891
	  matrix[y][x] = (Vertex_t)gpatch.v[y][x];
892
    }
893
}
894

895
Product

Resources

Company
