CoCalc -- bezier

GitHub Repository: godotengine/godot
Path: blob/master/thirdparty/embree/kernels/subdiv/bezier_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_curve.h"
8

9
namespace embree
10
{  
11
  template<class T, class S>
12
    static __forceinline T deCasteljau(const S& uu, const T& v0, const T& v1, const T& v2, const T& v3)
13
  {
14
    const T v0_1 = lerp(v0,v1,uu);
15
    const T v1_1 = lerp(v1,v2,uu);
16
    const T v2_1 = lerp(v2,v3,uu);
17
    const T v0_2 = lerp(v0_1,v1_1,uu);
18
    const T v1_2 = lerp(v1_1,v2_1,uu);
19
    const T v0_3 = lerp(v0_2,v1_2,uu);
20
    return v0_3;
21
  }
22
  
23
  template<class T, class S>
24
    static __forceinline T deCasteljau_tangent(const S& uu, const T& v0, const T& v1, const T& v2, const T& v3)
25
  {
26
    const T v0_1 = lerp(v0,v1,uu);
27
    const T v1_1 = lerp(v1,v2,uu);
28
    const T v2_1 = lerp(v2,v3,uu);
29
    const T v0_2 = lerp(v0_1,v1_1,uu);
30
    const T v1_2 = lerp(v1_1,v2_1,uu);
31
    return S(3.0f)*(v1_2-v0_2);
32
  }
33

34
  template<typename Vertex>
35
    __forceinline Vertex computeInnerBezierControlPoint(const Vertex v[4][4], const size_t y, const size_t x) {
36
    return 1.0f / 36.0f * (16.0f * v[y][x] + 4.0f * (v[y-1][x] +  v[y+1][x] + v[y][x-1] + v[y][x+1]) + (v[y-1][x-1] + v[y+1][x+1] + v[y-1][x+1] + v[y+1][x-1]));
37
  }
38
  
39
  template<typename Vertex>
40
    __forceinline Vertex computeTopEdgeBezierControlPoint(const Vertex v[4][4], const size_t y, const size_t x) {
41
    return 1.0f / 18.0f * (8.0f * v[y][x] + 4.0f * v[y-1][x] + 2.0f * (v[y][x-1] + v[y][x+1]) + (v[y-1][x-1] + v[y-1][x+1]));
42
  }
43

44
  template<typename Vertex>
45
    __forceinline Vertex computeBottomEdgeBezierControlPoint(const Vertex v[4][4], const size_t y, const size_t x) {
46
    return 1.0f / 18.0f * (8.0f * v[y][x] + 4.0f * v[y+1][x] + 2.0f * (v[y][x-1] + v[y][x+1]) + v[y+1][x-1] + v[y+1][x+1]);
47
  }
48
  
49
  template<typename Vertex>
50
    __forceinline Vertex computeLeftEdgeBezierControlPoint(const Vertex v[4][4], const size_t y, const size_t x) {
51
    return 1.0f / 18.0f * (8.0f * v[y][x] + 4.0f * v[y][x-1] + 2.0f * (v[y-1][x] + v[y+1][x]) + v[y-1][x-1] + v[y+1][x-1]);
52
  }
53
  
54
  template<typename Vertex>
55
    __forceinline Vertex computeRightEdgeBezierControlPoint(const Vertex v[4][4], const size_t y, const size_t x) {
56
    return 1.0f / 18.0f * (8.0f * v[y][x] + 4.0f * v[y][x+1] + 2.0f * (v[y-1][x] + v[y+1][x]) + v[y-1][x+1] + v[y+1][x+1]);
57
  }
58
  
59
  template<typename Vertex>
60
    __forceinline Vertex computeCornerBezierControlPoint(const Vertex v[4][4], const size_t y, const size_t x, const ssize_t delta_y, const ssize_t delta_x)
61
  {
62
    return 1.0f / 9.0f * (4.0f * v[y][x] + 2.0f * (v[y+delta_y][x] + v[y][x+delta_x]) + v[y+delta_y][x+delta_x]);
63
  }
64

65
  template<typename Vertex, typename Vertex_t>
66
    class __aligned(64) BezierPatchT
67
  {
68
   public:
69
      Vertex matrix[4][4];
70
    
71
  public:
72

73
    __forceinline BezierPatchT() {}
74

75
    __forceinline BezierPatchT (const HalfEdge* edge, const char* vertices, size_t stride);
76

77
    __forceinline BezierPatchT(const CatmullClarkPatchT<Vertex,Vertex_t>& patch);
78

79
    __forceinline BezierPatchT(const CatmullClarkPatchT<Vertex,Vertex_t>& patch,
80
                               const BezierCurveT<Vertex>* border0,
81
                               const BezierCurveT<Vertex>* border1,
82
                               const BezierCurveT<Vertex>* border2,
83
                               const BezierCurveT<Vertex>* border3);
84
                               
85
    __forceinline BezierPatchT(const BSplinePatchT<Vertex,Vertex_t>& source)
86
    {
87
      /* compute inner bezier control points */
88
      matrix[0][0] = computeInnerBezierControlPoint(source.v,1,1);
89
      matrix[0][3] = computeInnerBezierControlPoint(source.v,1,2);
90
      matrix[3][3] = computeInnerBezierControlPoint(source.v,2,2);
91
      matrix[3][0] = computeInnerBezierControlPoint(source.v,2,1);
92
      
93
      /* compute top edge control points */
94
      matrix[0][1] = computeRightEdgeBezierControlPoint(source.v,1,1);
95
      matrix[0][2] = computeLeftEdgeBezierControlPoint(source.v,1,2); 
96
      
97
      /* compute bottom edge control points */
98
      matrix[3][1] = computeRightEdgeBezierControlPoint(source.v,2,1);
99
      matrix[3][2] = computeLeftEdgeBezierControlPoint(source.v,2,2);
100
      
101
      /* compute left edge control points */
102
      matrix[1][0] = computeBottomEdgeBezierControlPoint(source.v,1,1);
103
      matrix[2][0] = computeTopEdgeBezierControlPoint(source.v,2,1);
104
      
105
      /* compute right edge control points */
106
      matrix[1][3] = computeBottomEdgeBezierControlPoint(source.v,1,2);
107
      matrix[2][3] = computeTopEdgeBezierControlPoint(source.v,2,2);
108
      
109
      /* compute corner control points */
110
      matrix[1][1] = computeCornerBezierControlPoint(source.v,1,1, 1, 1);
111
      matrix[1][2] = computeCornerBezierControlPoint(source.v,1,2, 1,-1);
112
      matrix[2][2] = computeCornerBezierControlPoint(source.v,2,2,-1,-1);
113
      matrix[2][1] = computeCornerBezierControlPoint(source.v,2,1,-1, 1);      
114
    }
115

116
    static __forceinline Vertex_t bilinear(const Vec4f Bu, const Vertex matrix[4][4], const Vec4f Bv)
117
    {
118
      const Vertex_t M0 = madd(Bu.x,matrix[0][0],madd(Bu.y,matrix[0][1],madd(Bu.z,matrix[0][2],Bu.w * matrix[0][3]))); 
119
      const Vertex_t M1 = madd(Bu.x,matrix[1][0],madd(Bu.y,matrix[1][1],madd(Bu.z,matrix[1][2],Bu.w * matrix[1][3])));
120
      const Vertex_t M2 = madd(Bu.x,matrix[2][0],madd(Bu.y,matrix[2][1],madd(Bu.z,matrix[2][2],Bu.w * matrix[2][3])));
121
      const Vertex_t M3 = madd(Bu.x,matrix[3][0],madd(Bu.y,matrix[3][1],madd(Bu.z,matrix[3][2],Bu.w * matrix[3][3])));
122
      return madd(Bv.x,M0,madd(Bv.y,M1,madd(Bv.z,M2,Bv.w*M3)));
123
    }
124

125
    static __forceinline Vertex_t eval(const Vertex matrix[4][4], const float uu, const float vv) 
126
    {      
127
      const Vec4f Bu = BezierBasis::eval(uu);
128
      const Vec4f Bv = BezierBasis::eval(vv);
129
      return bilinear(Bu,matrix,Bv);
130
    }
131

132
    static __forceinline Vertex_t eval_du(const Vertex matrix[4][4], const float uu, const float vv) 
133
    {
134
      const Vec4f Bu = BezierBasis::derivative(uu);
135
      const Vec4f Bv = BezierBasis::eval(vv);
136
      return bilinear(Bu,matrix,Bv);
137
    }
138

139
    static __forceinline Vertex_t eval_dv(const Vertex matrix[4][4], const float uu, const float vv) 
140
    {
141
      const Vec4f Bu = BezierBasis::eval(uu);
142
      const Vec4f Bv = BezierBasis::derivative(vv);
143
      return bilinear(Bu,matrix,Bv);
144
    }
145

146
    static __forceinline Vertex_t eval_dudu(const Vertex matrix[4][4], const float uu, const float vv) 
147
    {
148
      const Vec4f Bu = BezierBasis::derivative2(uu);
149
      const Vec4f Bv = BezierBasis::eval(vv);
150
      return bilinear(Bu,matrix,Bv);
151
    }
152

153
    static __forceinline Vertex_t eval_dvdv(const Vertex matrix[4][4], const float uu, const float vv) 
154
    {
155
      const Vec4f Bu = BezierBasis::eval(uu);
156
      const Vec4f Bv = BezierBasis::derivative2(vv);
157
      return bilinear(Bu,matrix,Bv);
158
    }
159

160
    static __forceinline Vertex_t eval_dudv(const Vertex matrix[4][4], const float uu, const float vv) 
161
    {
162
      const Vec4f Bu = BezierBasis::derivative(uu);
163
      const Vec4f Bv = BezierBasis::derivative(vv);
164
      return bilinear(Bu,matrix,Bv);
165
    }
166

167
    static __forceinline Vertex_t normal(const Vertex matrix[4][4], const float uu, const float vv) 
168
    {
169
      const Vertex_t dPdu = eval_du(matrix,uu,vv);
170
      const Vertex_t dPdv = eval_dv(matrix,uu,vv);
171
      return cross(dPdu,dPdv);
172
    }
173

174
    __forceinline Vertex_t normal(const float uu, const float vv) 
175
    {
176
      const Vertex_t dPdu = eval_du(matrix,uu,vv);
177
      const Vertex_t dPdv = eval_dv(matrix,uu,vv);
178
      return cross(dPdu,dPdv);
179
    }
180

181
    __forceinline Vertex_t eval(const float uu, const float vv) const {
182
      return eval(matrix,uu,vv);
183
    }
184

185
    __forceinline Vertex_t eval_du(const float uu, const float vv) const { 
186
      return eval_du(matrix,uu,vv);
187
    }
188

189
    __forceinline Vertex_t eval_dv(const float uu, const float vv) const {
190
      return eval_dv(matrix,uu,vv);
191
    }
192

193
    __forceinline Vertex_t eval_dudu(const float uu, const float vv) const { 
194
      return eval_dudu(matrix,uu,vv);
195
    }
196
    
197
    __forceinline Vertex_t eval_dvdv(const float uu, const float vv) const { 
198
      return eval_dvdv(matrix,uu,vv);
199
    }
200

201
    __forceinline Vertex_t eval_dudv(const float uu, const float vv) const { 
202
      return eval_dudv(matrix,uu,vv);
203
    }
204

205
    __forceinline void eval(const float u, const float v, Vertex* P, Vertex* dPdu, Vertex* dPdv, Vertex* ddPdudu, Vertex* ddPdvdv, Vertex* ddPdudv, const float dscale = 1.0f) const
206
    {
207
      if (P) {
208
        *P = eval(u,v); 
209
      }
210
      if (dPdu) {
211
        assert(dPdu); *dPdu = eval_du(u,v)*dscale; 
212
        assert(dPdv); *dPdv = eval_dv(u,v)*dscale; 
213
      }
214
      if (ddPdudu) {
215
        assert(ddPdudu); *ddPdudu = eval_dudu(u,v)*sqr(dscale); 
216
        assert(ddPdvdv); *ddPdvdv = eval_dvdv(u,v)*sqr(dscale); 
217
        assert(ddPdudv); *ddPdudv = eval_dudv(u,v)*sqr(dscale); 
218
      }
219
    }
220

221
    template<class vfloat>
222
      __forceinline vfloat eval(const size_t i, const vfloat& uu, const vfloat& vv, const Vec4<vfloat>& u_n, const Vec4<vfloat>& v_n) const
223
      {
224
        const vfloat curve0_x = v_n[0] * vfloat(matrix[0][0][i]) + v_n[1] * vfloat(matrix[1][0][i]) + v_n[2] * vfloat(matrix[2][0][i]) + v_n[3] * vfloat(matrix[3][0][i]);
225
        const vfloat curve1_x = v_n[0] * vfloat(matrix[0][1][i]) + v_n[1] * vfloat(matrix[1][1][i]) + v_n[2] * vfloat(matrix[2][1][i]) + v_n[3] * vfloat(matrix[3][1][i]);
226
        const vfloat curve2_x = v_n[0] * vfloat(matrix[0][2][i]) + v_n[1] * vfloat(matrix[1][2][i]) + v_n[2] * vfloat(matrix[2][2][i]) + v_n[3] * vfloat(matrix[3][2][i]);
227
        const vfloat curve3_x = v_n[0] * vfloat(matrix[0][3][i]) + v_n[1] * vfloat(matrix[1][3][i]) + v_n[2] * vfloat(matrix[2][3][i]) + v_n[3] * vfloat(matrix[3][3][i]);
228
        return u_n[0] * curve0_x + u_n[1] * curve1_x + u_n[2] * curve2_x + u_n[3] * curve3_x;
229
      }
230

231
    template<typename vbool, typename vfloat>
232
      __forceinline void eval(const vbool& valid, const vfloat& uu, const vfloat& vv, 
233
                              float* P, float* dPdu, float* dPdv, float* ddPdudu, float* ddPdvdv, float* ddPdudv,
234
                              const float dscale, const size_t dstride, const size_t N) const
235
      {
236
        if (P) {
237
          const Vec4<vfloat> u_n = BezierBasis::eval(uu); 
238
          const Vec4<vfloat> v_n = BezierBasis::eval(vv); 
239
          for (size_t i=0; i<N; i++) vfloat::store(valid,P+i*dstride,eval(i,uu,vv,u_n,v_n));
240
        }
241
        if (dPdu) 
242
        {
243
          {
244
            assert(dPdu);
245
            const Vec4<vfloat> u_n = BezierBasis::derivative(uu);
246
            const Vec4<vfloat> v_n = BezierBasis::eval(vv); 
247
            for (size_t i=0; i<N; i++) vfloat::store(valid,dPdu+i*dstride,eval(i,uu,vv,u_n,v_n)*dscale);
248
          }
249
          {
250
            assert(dPdv);
251
            const Vec4<vfloat> u_n = BezierBasis::eval(uu);
252
            const Vec4<vfloat> v_n = BezierBasis::derivative(vv); 
253
            for (size_t i=0; i<N; i++) vfloat::store(valid,dPdv+i*dstride,eval(i,uu,vv,u_n,v_n)*dscale);
254
          }
255
        }
256
        if (ddPdudu) 
257
        {
258
          {
259
            assert(ddPdudu);
260
            const Vec4<vfloat> u_n = BezierBasis::derivative2(uu);
261
            const Vec4<vfloat> v_n = BezierBasis::eval(vv); 
262
            for (size_t i=0; i<N; i++) vfloat::store(valid,ddPdudu+i*dstride,eval(i,uu,vv,u_n,v_n)*sqr(dscale));
263
          }
264
          {
265
            assert(ddPdvdv);
266
            const Vec4<vfloat> u_n = BezierBasis::eval(uu);
267
            const Vec4<vfloat> v_n = BezierBasis::derivative2(vv); 
268
            for (size_t i=0; i<N; i++) vfloat::store(valid,ddPdvdv+i*dstride,eval(i,uu,vv,u_n,v_n)*sqr(dscale));
269
          }
270
          {
271
            assert(ddPdudv);
272
            const Vec4<vfloat> u_n = BezierBasis::derivative(uu);
273
            const Vec4<vfloat> v_n = BezierBasis::derivative(vv); 
274
            for (size_t i=0; i<N; i++) vfloat::store(valid,ddPdudv+i*dstride,eval(i,uu,vv,u_n,v_n)*sqr(dscale));
275
          }
276
        }
277
      }
278

279
    template<typename T>
280
      static __forceinline Vec3<T> eval(const Vertex matrix[4][4], const T& uu, const T& vv) 
281
    {      
282
      const T one_minus_uu = 1.0f - uu;
283
      const T one_minus_vv = 1.0f - vv;      
284

285
      const T B0_u = one_minus_uu * one_minus_uu * one_minus_uu;
286
      const T B0_v = one_minus_vv * one_minus_vv * one_minus_vv;
287
      const T B1_u = 3.0f * (one_minus_uu * uu * one_minus_uu);
288
      const T B1_v = 3.0f * (one_minus_vv * vv * one_minus_vv);
289
      const T B2_u = 3.0f * (uu * one_minus_uu * uu);
290
      const T B2_v = 3.0f * (vv * one_minus_vv * vv);
291
      const T B3_u = uu * uu * uu;
292
      const T B3_v = vv * vv * vv;
293
      
294
      const T x = 
295
        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))), 
296
        madd(B1_v,madd(B0_u,matrix[1][0].x,madd(B1_u,matrix[1][1].x,madd(B2_u,matrix[1][2].x,B3_u*matrix[1][3].x))),
297
        madd(B2_v,madd(B0_u,matrix[2][0].x,madd(B1_u,matrix[2][1].x,madd(B2_u,matrix[2][2].x,B3_u*matrix[2][3].x))),
298
             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)))))); 
299

300
      const T y = 
301
        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))), 
302
        madd(B1_v,madd(B0_u,matrix[1][0].y,madd(B1_u,matrix[1][1].y,madd(B2_u,matrix[1][2].y,B3_u*matrix[1][3].y))),
303
        madd(B2_v,madd(B0_u,matrix[2][0].y,madd(B1_u,matrix[2][1].y,madd(B2_u,matrix[2][2].y,B3_u*matrix[2][3].y))),
304
             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)))))); 
305
      
306
      const T z = 
307
        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))), 
308
        madd(B1_v,madd(B0_u,matrix[1][0].z,madd(B1_u,matrix[1][1].z,madd(B2_u,matrix[1][2].z,B3_u*matrix[1][3].z))),
309
        madd(B2_v,madd(B0_u,matrix[2][0].z,madd(B1_u,matrix[2][1].z,madd(B2_u,matrix[2][2].z,B3_u*matrix[2][3].z))),
310
             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)))))); 
311
      
312
      return Vec3<T>(x,y,z);
313
    }
314

315
    template<typename vfloat>
316
      __forceinline Vec3<vfloat> eval(const vfloat& uu, const vfloat& vv) const {     
317
      return eval(matrix,uu,vv);
318
    }
319

320
    template<class T>
321
      static __forceinline Vec3<T> normal(const Vertex matrix[4][4], const T& uu, const T& vv) 
322
    {
323
      
324
      const Vec3<T> matrix_00 = Vec3<T>(matrix[0][0].x,matrix[0][0].y,matrix[0][0].z);
325
      const Vec3<T> matrix_01 = Vec3<T>(matrix[0][1].x,matrix[0][1].y,matrix[0][1].z);
326
      const Vec3<T> matrix_02 = Vec3<T>(matrix[0][2].x,matrix[0][2].y,matrix[0][2].z);
327
      const Vec3<T> matrix_03 = Vec3<T>(matrix[0][3].x,matrix[0][3].y,matrix[0][3].z);
328

329
      const Vec3<T> matrix_10 = Vec3<T>(matrix[1][0].x,matrix[1][0].y,matrix[1][0].z);
330
      const Vec3<T> matrix_11 = Vec3<T>(matrix[1][1].x,matrix[1][1].y,matrix[1][1].z);
331
      const Vec3<T> matrix_12 = Vec3<T>(matrix[1][2].x,matrix[1][2].y,matrix[1][2].z);
332
      const Vec3<T> matrix_13 = Vec3<T>(matrix[1][3].x,matrix[1][3].y,matrix[1][3].z);
333

334
      const Vec3<T> matrix_20 = Vec3<T>(matrix[2][0].x,matrix[2][0].y,matrix[2][0].z);
335
      const Vec3<T> matrix_21 = Vec3<T>(matrix[2][1].x,matrix[2][1].y,matrix[2][1].z);
336
      const Vec3<T> matrix_22 = Vec3<T>(matrix[2][2].x,matrix[2][2].y,matrix[2][2].z);
337
      const Vec3<T> matrix_23 = Vec3<T>(matrix[2][3].x,matrix[2][3].y,matrix[2][3].z);
338

339
      const Vec3<T> matrix_30 = Vec3<T>(matrix[3][0].x,matrix[3][0].y,matrix[3][0].z);
340
      const Vec3<T> matrix_31 = Vec3<T>(matrix[3][1].x,matrix[3][1].y,matrix[3][1].z);
341
      const Vec3<T> matrix_32 = Vec3<T>(matrix[3][2].x,matrix[3][2].y,matrix[3][2].z);
342
      const Vec3<T> matrix_33 = Vec3<T>(matrix[3][3].x,matrix[3][3].y,matrix[3][3].z);
343
            
344
      /* tangentU */
345
      const Vec3<T> col0 = deCasteljau(vv, matrix_00, matrix_10, matrix_20, matrix_30);
346
      const Vec3<T> col1 = deCasteljau(vv, matrix_01, matrix_11, matrix_21, matrix_31);
347
      const Vec3<T> col2 = deCasteljau(vv, matrix_02, matrix_12, matrix_22, matrix_32);
348
      const Vec3<T> col3 = deCasteljau(vv, matrix_03, matrix_13, matrix_23, matrix_33);
349
      
350
      const Vec3<T> tangentU = deCasteljau_tangent(uu, col0, col1, col2, col3);
351
      
352
      /* tangentV */
353
      const Vec3<T> row0 = deCasteljau(uu, matrix_00, matrix_01, matrix_02, matrix_03);
354
      const Vec3<T> row1 = deCasteljau(uu, matrix_10, matrix_11, matrix_12, matrix_13);
355
      const Vec3<T> row2 = deCasteljau(uu, matrix_20, matrix_21, matrix_22, matrix_23);
356
      const Vec3<T> row3 = deCasteljau(uu, matrix_30, matrix_31, matrix_32, matrix_33);
357
      
358
      const Vec3<T> tangentV = deCasteljau_tangent(vv, row0, row1, row2, row3);
359
      
360
      /* normal = tangentU x tangentV */
361
      const Vec3<T> n = cross(tangentU,tangentV);
362
      return n;
363
    }
364

365
    template<typename vfloat>
366
      __forceinline Vec3<vfloat> normal(const vfloat& uu, const vfloat& vv) const {     
367
      return normal(matrix,uu,vv);
368
    }
369
  };
370

371
  typedef BezierPatchT<Vec3fa,Vec3fa_t> BezierPatch3fa;
372
}
373

374
Product

Resources

Company
