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 "bspline_curve.h"
8
9namespace embree
10{
11 template<typename Vertex, typename Vertex_t = Vertex>
12 class __aligned(64) BSplinePatchT
13 {
14 typedef CatmullClark1RingT<Vertex,Vertex_t> CatmullClarkRing;
15 typedef CatmullClarkPatchT<Vertex,Vertex_t> CatmullClarkPatch;
16
17 public:
18
19 __forceinline BSplinePatchT () {}
20
21 __forceinline BSplinePatchT (const CatmullClarkPatch& patch) {
22 init(patch);
23 }
24
25 __forceinline BSplinePatchT(const CatmullClarkPatch& patch,
26 const BezierCurveT<Vertex>* border0,
27 const BezierCurveT<Vertex>* border1,
28 const BezierCurveT<Vertex>* border2,
29 const BezierCurveT<Vertex>* border3)
30 {
31 init(patch);
32 }
33
34 __forceinline BSplinePatchT (const HalfEdge* edge, const char* vertices, size_t stride) {
35 init(edge,vertices,stride);
36 }
37
38 __forceinline Vertex hard_corner(const Vertex& v01, const Vertex& v02,
39 const Vertex& v10, const Vertex& v11, const Vertex& v12,
40 const Vertex& v20, const Vertex& v21, const Vertex& v22)
41 {
42 return 4.0f*v11 - 2.0f*(v12+v21) + v22;
43 }
44
45 __forceinline Vertex soft_convex_corner( const Vertex& v01, const Vertex& v02,
46 const Vertex& v10, const Vertex& v11, const Vertex& v12,
47 const Vertex& v20, const Vertex& v21, const Vertex& v22)
48 {
49 return -8.0f*v11 + 4.0f*(v12+v21) + v22;
50 }
51
52 __forceinline Vertex convex_corner(const float vertex_crease_weight,
53 const Vertex& v01, const Vertex& v02,
54 const Vertex& v10, const Vertex& v11, const Vertex& v12,
55 const Vertex& v20, const Vertex& v21, const Vertex& v22)
56 {
57 if (std::isinf(vertex_crease_weight)) return hard_corner(v01,v02,v10,v11,v12,v20,v21,v22);
58 else return soft_convex_corner(v01,v02,v10,v11,v12,v20,v21,v22);
59 }
60
61 __forceinline Vertex load(const HalfEdge* edge, const char* vertices, size_t stride) {
62 return Vertex_t::loadu(vertices+edge->getStartVertexIndex()*stride);
63 }
64
65 __forceinline void init_border(const CatmullClarkRing& edge0,
66 Vertex& v01, Vertex& v02,
67 const Vertex& v11, const Vertex& v12,
68 const Vertex& v21, const Vertex& v22)
69 {
70 if (likely(edge0.has_opposite_back(0)))
71 {
72 v01 = edge0.back(2);
73 v02 = edge0.back(1);
74 } else {
75 v01 = 2.0f*v11-v21;
76 v02 = 2.0f*v12-v22;
77 }
78 }
79
80 __forceinline void init_corner(const CatmullClarkRing& edge0,
81 Vertex& v00, const Vertex& v01, const Vertex& v02,
82 const Vertex& v10, const Vertex& v11, const Vertex& v12,
83 const Vertex& v20, const Vertex& v21, const Vertex& v22)
84 {
85 const bool MAYBE_UNUSED has_back1 = edge0.has_opposite_back(1);
86 const bool has_back0 = edge0.has_opposite_back(0);
87 const bool has_front1 = edge0.has_opposite_front(1);
88 const bool MAYBE_UNUSED has_front2 = edge0.has_opposite_front(2);
89
90 if (likely(has_back0)) {
91 if (likely(has_front1)) { assert(has_back1 && has_front2); v00 = edge0.back(3); }
92 else { assert(!has_back1); v00 = 2.0f*v01-v02; }
93 }
94 else {
95 if (likely(has_front1)) { assert(!has_front2); v00 = 2.0f*v10-v20; }
96 else v00 = convex_corner(edge0.vertex_crease_weight,v01,v02,v10,v11,v12,v20,v21,v22);
97 }
98 }
99
100 void init(const CatmullClarkPatch& patch)
101 {
102 /* fill inner vertices */
103 const Vertex v11 = v[1][1] = patch.ring[0].vtx;
104 const Vertex v12 = v[1][2] = patch.ring[1].vtx;
105 const Vertex v22 = v[2][2] = patch.ring[2].vtx;
106 const Vertex v21 = v[2][1] = patch.ring[3].vtx;
107
108 /* fill border vertices */
109 init_border(patch.ring[0],v[0][1],v[0][2],v11,v12,v21,v22);
110 init_border(patch.ring[1],v[1][3],v[2][3],v12,v22,v11,v21);
111 init_border(patch.ring[2],v[3][2],v[3][1],v22,v21,v12,v11);
112 init_border(patch.ring[3],v[2][0],v[1][0],v21,v11,v22,v12);
113
114 /* fill corner vertices */
115 init_corner(patch.ring[0],v[0][0],v[0][1],v[0][2],v[1][0],v11,v12,v[2][0],v21,v22);
116 init_corner(patch.ring[1],v[0][3],v[1][3],v[2][3],v[0][2],v12,v22,v[0][1],v11,v21);
117 init_corner(patch.ring[2],v[3][3],v[3][2],v[3][1],v[2][3],v22,v21,v[1][3],v12,v11);
118 init_corner(patch.ring[3],v[3][0],v[2][0],v[1][0],v[3][1],v21,v11,v[3][2],v22,v12);
119 }
120
121 void init_border(const HalfEdge* edge0, const char* vertices, size_t stride,
122 Vertex& v01, Vertex& v02,
123 const Vertex& v11, const Vertex& v12,
124 const Vertex& v21, const Vertex& v22)
125 {
126 if (likely(edge0->hasOpposite()))
127 {
128 const HalfEdge* e = edge0->opposite()->next()->next();
129 v01 = load(e,vertices,stride);
130 v02 = load(e->next(),vertices,stride);
131 } else {
132 v01 = 2.0f*v11-v21;
133 v02 = 2.0f*v12-v22;
134 }
135 }
136
137 void init_corner(const HalfEdge* edge0, const char* vertices, size_t stride,
138 Vertex& v00, const Vertex& v01, const Vertex& v02,
139 const Vertex& v10, const Vertex& v11, const Vertex& v12,
140 const Vertex& v20, const Vertex& v21, const Vertex& v22)
141 {
142 const bool has_back0 = edge0->hasOpposite();
143 const bool has_front1 = edge0->prev()->hasOpposite();
144
145 if (likely(has_back0))
146 {
147 const HalfEdge* e = edge0->opposite()->next();
148 if (likely(has_front1))
149 {
150 assert(e->hasOpposite());
151 assert(edge0->prev()->opposite()->prev()->hasOpposite());
152 v00 = load(e->opposite()->prev(),vertices,stride);
153 }
154 else {
155 assert(!e->hasOpposite());
156 v00 = 2.0f*v01-v02;
157 }
158 }
159 else
160 {
161 if (likely(has_front1)) {
162 assert(!edge0->prev()->opposite()->prev()->hasOpposite());
163 v00 = 2.0f*v10-v20;
164 }
165 else {
166 assert(edge0->vertex_crease_weight == 0.0f || std::isinf(edge0->vertex_crease_weight));
167 v00 = convex_corner(edge0->vertex_crease_weight,v01,v02,v10,v11,v12,v20,v21,v22);
168 }
169 }
170 }
171
172 void init(const HalfEdge* edge0, const char* vertices, size_t stride)
173 {
174 assert( edge0->isRegularFace() );
175
176 /* fill inner vertices */
177 const Vertex v11 = v[1][1] = load(edge0,vertices,stride); const HalfEdge* edge1 = edge0->next();
178 const Vertex v12 = v[1][2] = load(edge1,vertices,stride); const HalfEdge* edge2 = edge1->next();
179 const Vertex v22 = v[2][2] = load(edge2,vertices,stride); const HalfEdge* edge3 = edge2->next();
180 const Vertex v21 = v[2][1] = load(edge3,vertices,stride); assert(edge0 == edge3->next());
181
182 /* fill border vertices */
183 init_border(edge0,vertices,stride,v[0][1],v[0][2],v11,v12,v21,v22);
184 init_border(edge1,vertices,stride,v[1][3],v[2][3],v12,v22,v11,v21);
185 init_border(edge2,vertices,stride,v[3][2],v[3][1],v22,v21,v12,v11);
186 init_border(edge3,vertices,stride,v[2][0],v[1][0],v21,v11,v22,v12);
187
188 /* fill corner vertices */
189 init_corner(edge0,vertices,stride,v[0][0],v[0][1],v[0][2],v[1][0],v11,v12,v[2][0],v21,v22);
190 init_corner(edge1,vertices,stride,v[0][3],v[1][3],v[2][3],v[0][2],v12,v22,v[0][1],v11,v21);
191 init_corner(edge2,vertices,stride,v[3][3],v[3][2],v[3][1],v[2][3],v22,v21,v[1][3],v12,v11);
192 init_corner(edge3,vertices,stride,v[3][0],v[2][0],v[1][0],v[3][1],v21,v11,v[3][2],v22,v12);
193 }
194
195 __forceinline BBox<Vertex> bounds() const
196 {
197 const Vertex* const cv = &v[0][0];
198 BBox<Vertex> bounds (cv[0]);
199 for (size_t i=1; i<16 ; i++)
200 bounds.extend( cv[i] );
201 return bounds;
202 }
203
204 __forceinline Vertex eval(const float uu, const float vv) const
205 {
206 const Vec4f v_n = BSplineBasis::eval(vv);
207 const Vertex_t curve0 = madd(v_n[0],v[0][0],madd(v_n[1],v[1][0],madd(v_n[2],v[2][0],v_n[3] * v[3][0])));
208 const Vertex_t curve1 = madd(v_n[0],v[0][1],madd(v_n[1],v[1][1],madd(v_n[2],v[2][1],v_n[3] * v[3][1])));
209 const Vertex_t curve2 = madd(v_n[0],v[0][2],madd(v_n[1],v[1][2],madd(v_n[2],v[2][2],v_n[3] * v[3][2])));
210 const Vertex_t curve3 = madd(v_n[0],v[0][3],madd(v_n[1],v[1][3],madd(v_n[2],v[2][3],v_n[3] * v[3][3])));
211
212 const Vec4f u_n = BSplineBasis::eval(uu);
213 return madd(u_n[0],curve0,madd(u_n[1],curve1,madd(u_n[2],curve2,u_n[3] * curve3)));
214 }
215
216 __forceinline Vertex eval_du(const float uu, const float vv) const
217 {
218 const Vec4f v_n = BSplineBasis::eval(vv);
219 const Vertex_t curve0 = madd(v_n[0],v[0][0],madd(v_n[1],v[1][0],madd(v_n[2],v[2][0],v_n[3] * v[3][0])));
220 const Vertex_t curve1 = madd(v_n[0],v[0][1],madd(v_n[1],v[1][1],madd(v_n[2],v[2][1],v_n[3] * v[3][1])));
221 const Vertex_t curve2 = madd(v_n[0],v[0][2],madd(v_n[1],v[1][2],madd(v_n[2],v[2][2],v_n[3] * v[3][2])));
222 const Vertex_t curve3 = madd(v_n[0],v[0][3],madd(v_n[1],v[1][3],madd(v_n[2],v[2][3],v_n[3] * v[3][3])));
223
224 const Vec4f u_n = BSplineBasis::derivative(uu);
225 return madd(u_n[0],curve0,madd(u_n[1],curve1,madd(u_n[2],curve2,u_n[3] * curve3)));
226 }
227
228 __forceinline Vertex eval_dv(const float uu, const float vv) const
229 {
230 const Vec4f v_n = BSplineBasis::derivative(vv);
231 const Vertex_t curve0 = madd(v_n[0],v[0][0],madd(v_n[1],v[1][0],madd(v_n[2],v[2][0],v_n[3] * v[3][0])));
232 const Vertex_t curve1 = madd(v_n[0],v[0][1],madd(v_n[1],v[1][1],madd(v_n[2],v[2][1],v_n[3] * v[3][1])));
233 const Vertex_t curve2 = madd(v_n[0],v[0][2],madd(v_n[1],v[1][2],madd(v_n[2],v[2][2],v_n[3] * v[3][2])));
234 const Vertex_t curve3 = madd(v_n[0],v[0][3],madd(v_n[1],v[1][3],madd(v_n[2],v[2][3],v_n[3] * v[3][3])));
235
236 const Vec4f u_n = BSplineBasis::eval(uu);
237 return madd(u_n[0],curve0,madd(u_n[1],curve1,madd(u_n[2],curve2,u_n[3] * curve3)));
238 }
239
240 __forceinline Vertex eval_dudu(const float uu, const float vv) const
241 {
242 const Vec4f v_n = BSplineBasis::eval(vv);
243 const Vertex_t curve0 = madd(v_n[0],v[0][0],madd(v_n[1],v[1][0],madd(v_n[2],v[2][0],v_n[3] * v[3][0])));
244 const Vertex_t curve1 = madd(v_n[0],v[0][1],madd(v_n[1],v[1][1],madd(v_n[2],v[2][1],v_n[3] * v[3][1])));
245 const Vertex_t curve2 = madd(v_n[0],v[0][2],madd(v_n[1],v[1][2],madd(v_n[2],v[2][2],v_n[3] * v[3][2])));
246 const Vertex_t curve3 = madd(v_n[0],v[0][3],madd(v_n[1],v[1][3],madd(v_n[2],v[2][3],v_n[3] * v[3][3])));
247
248 const Vec4f u_n = BSplineBasis::derivative2(uu);
249 return madd(u_n[0],curve0,madd(u_n[1],curve1,madd(u_n[2],curve2,u_n[3] * curve3)));
250 }
251
252 __forceinline Vertex eval_dvdv(const float uu, const float vv) const
253 {
254 const Vec4f v_n = BSplineBasis::derivative2(vv);
255 const Vertex_t curve0 = madd(v_n[0],v[0][0],madd(v_n[1],v[1][0],madd(v_n[2],v[2][0],v_n[3] * v[3][0])));
256 const Vertex_t curve1 = madd(v_n[0],v[0][1],madd(v_n[1],v[1][1],madd(v_n[2],v[2][1],v_n[3] * v[3][1])));
257 const Vertex_t curve2 = madd(v_n[0],v[0][2],madd(v_n[1],v[1][2],madd(v_n[2],v[2][2],v_n[3] * v[3][2])));
258 const Vertex_t curve3 = madd(v_n[0],v[0][3],madd(v_n[1],v[1][3],madd(v_n[2],v[2][3],v_n[3] * v[3][3])));
259
260 const Vec4f u_n = BSplineBasis::eval(uu);
261 return madd(u_n[0],curve0,madd(u_n[1],curve1,madd(u_n[2],curve2,u_n[3] * curve3)));
262 }
263
264 __forceinline Vertex eval_dudv(const float uu, const float vv) const
265 {
266 const Vec4f v_n = BSplineBasis::derivative(vv);
267 const Vertex_t curve0 = madd(v_n[0],v[0][0],madd(v_n[1],v[1][0],madd(v_n[2],v[2][0],v_n[3] * v[3][0])));
268 const Vertex_t curve1 = madd(v_n[0],v[0][1],madd(v_n[1],v[1][1],madd(v_n[2],v[2][1],v_n[3] * v[3][1])));
269 const Vertex_t curve2 = madd(v_n[0],v[0][2],madd(v_n[1],v[1][2],madd(v_n[2],v[2][2],v_n[3] * v[3][2])));
270 const Vertex_t curve3 = madd(v_n[0],v[0][3],madd(v_n[1],v[1][3],madd(v_n[2],v[2][3],v_n[3] * v[3][3])));
271
272 const Vec4f u_n = BSplineBasis::derivative(uu);
273 return madd(u_n[0],curve0,madd(u_n[1],curve1,madd(u_n[2],curve2,u_n[3] * curve3)));
274 }
275
276 __forceinline Vertex normal(const float uu, const float vv) const
277 {
278 const Vertex tu = eval_du(uu,vv);
279 const Vertex tv = eval_dv(uu,vv);
280 return cross(tu,tv);
281 }
282
283 template<typename T>
284 __forceinline Vec3<T> eval(const T& uu, const T& vv, const Vec4<T>& u_n, const Vec4<T>& v_n) const
285 {
286 const T curve0_x = madd(v_n[0],T(v[0][0].x),madd(v_n[1],T(v[1][0].x),madd(v_n[2],T(v[2][0].x),v_n[3] * T(v[3][0].x))));
287 const T curve1_x = madd(v_n[0],T(v[0][1].x),madd(v_n[1],T(v[1][1].x),madd(v_n[2],T(v[2][1].x),v_n[3] * T(v[3][1].x))));
288 const T curve2_x = madd(v_n[0],T(v[0][2].x),madd(v_n[1],T(v[1][2].x),madd(v_n[2],T(v[2][2].x),v_n[3] * T(v[3][2].x))));
289 const T curve3_x = madd(v_n[0],T(v[0][3].x),madd(v_n[1],T(v[1][3].x),madd(v_n[2],T(v[2][3].x),v_n[3] * T(v[3][3].x))));
290 const T x = madd(u_n[0],curve0_x,madd(u_n[1],curve1_x,madd(u_n[2],curve2_x,u_n[3] * curve3_x)));
291
292 const T curve0_y = madd(v_n[0],T(v[0][0].y),madd(v_n[1],T(v[1][0].y),madd(v_n[2],T(v[2][0].y),v_n[3] * T(v[3][0].y))));
293 const T curve1_y = madd(v_n[0],T(v[0][1].y),madd(v_n[1],T(v[1][1].y),madd(v_n[2],T(v[2][1].y),v_n[3] * T(v[3][1].y))));
294 const T curve2_y = madd(v_n[0],T(v[0][2].y),madd(v_n[1],T(v[1][2].y),madd(v_n[2],T(v[2][2].y),v_n[3] * T(v[3][2].y))));
295 const T curve3_y = madd(v_n[0],T(v[0][3].y),madd(v_n[1],T(v[1][3].y),madd(v_n[2],T(v[2][3].y),v_n[3] * T(v[3][3].y))));
296 const T y = madd(u_n[0],curve0_y,madd(u_n[1],curve1_y,madd(u_n[2],curve2_y,u_n[3] * curve3_y)));
297
298 const T curve0_z = madd(v_n[0],T(v[0][0].z),madd(v_n[1],T(v[1][0].z),madd(v_n[2],T(v[2][0].z),v_n[3] * T(v[3][0].z))));
299 const T curve1_z = madd(v_n[0],T(v[0][1].z),madd(v_n[1],T(v[1][1].z),madd(v_n[2],T(v[2][1].z),v_n[3] * T(v[3][1].z))));
300 const T curve2_z = madd(v_n[0],T(v[0][2].z),madd(v_n[1],T(v[1][2].z),madd(v_n[2],T(v[2][2].z),v_n[3] * T(v[3][2].z))));
301 const T curve3_z = madd(v_n[0],T(v[0][3].z),madd(v_n[1],T(v[1][3].z),madd(v_n[2],T(v[2][3].z),v_n[3] * T(v[3][3].z))));
302 const T z = madd(u_n[0],curve0_z,madd(u_n[1],curve1_z,madd(u_n[2],curve2_z,u_n[3] * curve3_z)));
303
304 return Vec3<T>(x,y,z);
305 }
306
307 template<typename T>
308 __forceinline Vec3<T> eval(const T& uu, const T& vv) const
309 {
310 const Vec4<T> u_n = BSplineBasis::eval(uu);
311 const Vec4<T> v_n = BSplineBasis::eval(vv);
312 return eval(uu,vv,u_n,v_n);
313 }
314
315 template<typename T>
316 __forceinline Vec3<T> eval_du(const T& uu, const T& vv) const
317 {
318 const Vec4<T> u_n = BSplineBasis::derivative(uu);
319 const Vec4<T> v_n = BSplineBasis::eval(vv);
320 return eval(uu,vv,u_n,v_n);
321 }
322
323 template<typename T>
324 __forceinline Vec3<T> eval_dv(const T& uu, const T& vv) const
325 {
326 const Vec4<T> u_n = BSplineBasis::eval(uu);
327 const Vec4<T> v_n = BSplineBasis::derivative(vv);
328 return eval(uu,vv,u_n,v_n);
329 }
330
331 template<typename T>
332 __forceinline Vec3<T> eval_dudu(const T& uu, const T& vv) const
333 {
334 const Vec4<T> u_n = BSplineBasis::derivative2(uu);
335 const Vec4<T> v_n = BSplineBasis::eval(vv);
336 return eval(uu,vv,u_n,v_n);
337 }
338
339 template<typename T>
340 __forceinline Vec3<T> eval_dvdv(const T& uu, const T& vv) const
341 {
342 const Vec4<T> u_n = BSplineBasis::eval(uu);
343 const Vec4<T> v_n = BSplineBasis::derivative2(vv);
344 return eval(uu,vv,u_n,v_n);
345 }
346
347 template<typename T>
348 __forceinline Vec3<T> eval_dudv(const T& uu, const T& vv) const
349 {
350 const Vec4<T> u_n = BSplineBasis::derivative(uu);
351 const Vec4<T> v_n = BSplineBasis::derivative(vv);
352 return eval(uu,vv,u_n,v_n);
353 }
354
355 template<typename T>
356 __forceinline Vec3<T> normal(const T& uu, const T& vv) const {
357 return cross(eval_du(uu,vv),eval_dv(uu,vv));
358 }
359
360 void eval(const float u, const float v,
361 Vertex* P, Vertex* dPdu, Vertex* dPdv, Vertex* ddPdudu, Vertex* ddPdvdv, Vertex* ddPdudv,
362 const float dscale = 1.0f) const
363 {
364 if (P) {
365 *P = eval(u,v);
366 }
367 if (dPdu) {
368 assert(dPdu); *dPdu = eval_du(u,v)*dscale;
369 assert(dPdv); *dPdv = eval_dv(u,v)*dscale;
370 }
371 if (ddPdudu) {
372 assert(ddPdudu); *ddPdudu = eval_dudu(u,v)*sqr(dscale);
373 assert(ddPdvdv); *ddPdvdv = eval_dvdv(u,v)*sqr(dscale);
374 assert(ddPdudv); *ddPdudv = eval_dudv(u,v)*sqr(dscale);
375 }
376 }
377
378 template<class vfloat>
379 __forceinline vfloat eval(const size_t i, const vfloat& uu, const vfloat& vv, const Vec4<vfloat>& u_n, const Vec4<vfloat>& v_n) const
380 {
381 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]))));
382 const vfloat curve1_x = madd(v_n[0],vfloat(v[0][1][i]),madd(v_n[1],vfloat(v[1][1][i]),madd(v_n[2],vfloat(v[2][1][i]),v_n[3] * vfloat(v[3][1][i]))));
383 const vfloat curve2_x = madd(v_n[0],vfloat(v[0][2][i]),madd(v_n[1],vfloat(v[1][2][i]),madd(v_n[2],vfloat(v[2][2][i]),v_n[3] * vfloat(v[3][2][i]))));
384 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]))));
385 return madd(u_n[0],curve0_x,madd(u_n[1],curve1_x,madd(u_n[2],curve2_x,u_n[3] * curve3_x)));
386 }
387
388 template<typename vbool, typename vfloat>
389 void eval(const vbool& valid, const vfloat& uu, const vfloat& vv,
390 float* P, float* dPdu, float* dPdv, float* ddPdudu, float* ddPdvdv, float* ddPdudv,
391 const float dscale, const size_t dstride, const size_t N) const
392 {
393 if (P) {
394 const Vec4<vfloat> u_n = BSplineBasis::eval(uu);
395 const Vec4<vfloat> v_n = BSplineBasis::eval(vv);
396 for (size_t i=0; i<N; i++) vfloat::store(valid,P+i*dstride,eval(i,uu,vv,u_n,v_n));
397 }
398 if (dPdu)
399 {
400 {
401 assert(dPdu);
402 const Vec4<vfloat> u_n = BSplineBasis::derivative(uu);
403 const Vec4<vfloat> v_n = BSplineBasis::eval(vv);
404 for (size_t i=0; i<N; i++) vfloat::store(valid,dPdu+i*dstride,eval(i,uu,vv,u_n,v_n)*dscale);
405 }
406 {
407 assert(dPdv);
408 const Vec4<vfloat> u_n = BSplineBasis::eval(uu);
409 const Vec4<vfloat> v_n = BSplineBasis::derivative(vv);
410 for (size_t i=0; i<N; i++) vfloat::store(valid,dPdv+i*dstride,eval(i,uu,vv,u_n,v_n)*dscale);
411 }
412 }
413 if (ddPdudu)
414 {
415 {
416 assert(ddPdudu);
417 const Vec4<vfloat> u_n = BSplineBasis::derivative2(uu);
418 const Vec4<vfloat> v_n = BSplineBasis::eval(vv);
419 for (size_t i=0; i<N; i++) vfloat::store(valid,ddPdudu+i*dstride,eval(i,uu,vv,u_n,v_n)*sqr(dscale));
420 }
421 {
422 assert(ddPdvdv);
423 const Vec4<vfloat> u_n = BSplineBasis::eval(uu);
424 const Vec4<vfloat> v_n = BSplineBasis::derivative2(vv);
425 for (size_t i=0; i<N; i++) vfloat::store(valid,ddPdvdv+i*dstride,eval(i,uu,vv,u_n,v_n)*sqr(dscale));
426 }
427 {
428 assert(ddPdudv);
429 const Vec4<vfloat> u_n = BSplineBasis::derivative(uu);
430 const Vec4<vfloat> v_n = BSplineBasis::derivative(vv);
431 for (size_t i=0; i<N; i++) vfloat::store(valid,ddPdudv+i*dstride,eval(i,uu,vv,u_n,v_n)*sqr(dscale));
432 }
433 }
434 }
435
436 friend __forceinline embree_ostream operator<<(embree_ostream o, const BSplinePatchT& p)
437 {
438 for (size_t y=0; y<4; y++)
439 for (size_t x=0; x<4; x++)
440 o << "[" << y << "][" << x << "] " << p.v[y][x] << embree_endl;
441 return o;
442 }
443
444 public:
445 Vertex v[4][4];
446 };
447
448 typedef BSplinePatchT<Vec3fa,Vec3fa_t> BSplinePatch3fa;
449}
450