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
9namespace 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