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 | |
9 | namespace 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(x: 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(vertex_crease_weight: 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(edge: e,vertices,stride); |
130 | v02 = load(edge: 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(edge: 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(vertex_crease_weight: 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(edge: edge0,vertices,stride); const HalfEdge* edge1 = edge0->next(); |
178 | const Vertex v12 = v[1][2] = load(edge: edge1,vertices,stride); const HalfEdge* edge2 = edge1->next(); |
179 | const Vertex v22 = v[2][2] = load(edge: edge2,vertices,stride); const HalfEdge* edge3 = edge2->next(); |
180 | const Vertex v21 = v[2][1] = load(edge: 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(u: 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(u: 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(u: 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(u: 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(u: 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(u: 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(u: 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(u: 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(u: 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(u: 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(u: 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(u: 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(x: dscale); |
373 | assert(ddPdvdv); *ddPdvdv = eval_dvdv(u,v)*sqr(x: dscale); |
374 | assert(ddPdudv); *ddPdudv = eval_dudv(u,v)*sqr(x: 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(x: 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(x: 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(x: 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 | |