1 | // Copyright 2009-2021 Intel Corporation |
2 | // SPDX-License-Identifier: Apache-2.0 |
3 | |
4 | #pragma once |
5 | |
6 | #include "../common/default.h" |
7 | #include "bezier_curve.h" |
8 | |
9 | namespace embree |
10 | { |
11 | class BSplineBasis |
12 | { |
13 | public: |
14 | |
15 | template<typename T> |
16 | static __forceinline Vec4<T> eval(const T& u) |
17 | { |
18 | const T t = u; |
19 | const T s = T(1.0f) - u; |
20 | const T n0 = s*s*s; |
21 | const T n1 = (4.0f*(s*s*s)+(t*t*t)) + (12.0f*((s*t)*s) + 6.0f*((t*s)*t)); |
22 | const T n2 = (4.0f*(t*t*t)+(s*s*s)) + (12.0f*((t*s)*t) + 6.0f*((s*t)*s)); |
23 | const T n3 = t*t*t; |
24 | return T(1.0f/6.0f)*Vec4<T>(n0,n1,n2,n3); |
25 | } |
26 | |
27 | template<typename T> |
28 | static __forceinline Vec4<T> derivative(const T& u) |
29 | { |
30 | const T t = u; |
31 | const T s = 1.0f - u; |
32 | const T n0 = -s*s; |
33 | const T n1 = -t*t - 4.0f*(t*s); |
34 | const T n2 = s*s + 4.0f*(s*t); |
35 | const T n3 = t*t; |
36 | return T(0.5f)*Vec4<T>(n0,n1,n2,n3); |
37 | } |
38 | |
39 | template<typename T> |
40 | static __forceinline Vec4<T> derivative2(const T& u) |
41 | { |
42 | const T t = u; |
43 | const T s = 1.0f - u; |
44 | const T n0 = s; |
45 | const T n1 = t - 2.0f*s; |
46 | const T n2 = s - 2.0f*t; |
47 | const T n3 = t; |
48 | return Vec4<T>(n0,n1,n2,n3); |
49 | } |
50 | }; |
51 | |
52 | struct PrecomputedBSplineBasis |
53 | { |
54 | enum { N = 16 }; |
55 | public: |
56 | PrecomputedBSplineBasis() {} |
57 | PrecomputedBSplineBasis(int shift); |
58 | |
59 | /* basis for bspline evaluation */ |
60 | public: |
61 | float c0[N+1][N+1]; |
62 | float c1[N+1][N+1]; |
63 | float c2[N+1][N+1]; |
64 | float c3[N+1][N+1]; |
65 | |
66 | /* basis for bspline derivative evaluation */ |
67 | public: |
68 | float d0[N+1][N+1]; |
69 | float d1[N+1][N+1]; |
70 | float d2[N+1][N+1]; |
71 | float d3[N+1][N+1]; |
72 | }; |
73 | extern PrecomputedBSplineBasis bspline_basis0; |
74 | extern PrecomputedBSplineBasis bspline_basis1; |
75 | |
76 | template<typename Vertex> |
77 | struct BSplineCurveT |
78 | { |
79 | Vertex v0,v1,v2,v3; |
80 | |
81 | __forceinline BSplineCurveT() {} |
82 | |
83 | __forceinline BSplineCurveT(const Vertex& v0, const Vertex& v1, const Vertex& v2, const Vertex& v3) |
84 | : v0(v0), v1(v1), v2(v2), v3(v3) {} |
85 | |
86 | __forceinline Vertex begin() const { |
87 | return madd(1.0f/6.0f,v0,madd(2.0f/3.0f,v1,1.0f/6.0f*v2)); |
88 | } |
89 | |
90 | __forceinline Vertex end() const { |
91 | return madd(1.0f/6.0f,v1,madd(2.0f/3.0f,v2,1.0f/6.0f*v3)); |
92 | } |
93 | |
94 | __forceinline Vertex center() const { |
95 | return 0.25f*(v0+v1+v2+v3); |
96 | } |
97 | |
98 | __forceinline BBox<Vertex> bounds() const { |
99 | return merge(BBox<Vertex>(v0),BBox<Vertex>(v1),BBox<Vertex>(v2),BBox<Vertex>(v3)); |
100 | } |
101 | |
102 | __forceinline friend BSplineCurveT operator -( const BSplineCurveT& a, const Vertex& b ) { |
103 | return BSplineCurveT(a.v0-b,a.v1-b,a.v2-b,a.v3-b); |
104 | } |
105 | |
106 | __forceinline BSplineCurveT<Vec3ff> xfm_pr(const LinearSpace3fa& space, const Vec3fa& p) const |
107 | { |
108 | const Vec3ff q0(xfmVector(space,(Vec3fa)v0-p), v0.w); |
109 | const Vec3ff q1(xfmVector(space,(Vec3fa)v1-p), v1.w); |
110 | const Vec3ff q2(xfmVector(space,(Vec3fa)v2-p), v2.w); |
111 | const Vec3ff q3(xfmVector(space,(Vec3fa)v3-p), v3.w); |
112 | return BSplineCurveT<Vec3ff>(q0,q1,q2,q3); |
113 | } |
114 | |
115 | __forceinline Vertex eval(const float t) const |
116 | { |
117 | const Vec4<float> b = BSplineBasis::eval(t); |
118 | return madd(b.x,v0,madd(b.y,v1,madd(b.z,v2,b.w*v3))); |
119 | } |
120 | |
121 | __forceinline Vertex eval_du(const float t) const |
122 | { |
123 | const Vec4<float> b = BSplineBasis::derivative(t); |
124 | return madd(b.x,v0,madd(b.y,v1,madd(b.z,v2,b.w*v3))); |
125 | } |
126 | |
127 | __forceinline Vertex eval_dudu(const float t) const |
128 | { |
129 | const Vec4<float> b = BSplineBasis::derivative2(t); |
130 | return madd(b.x,v0,madd(b.y,v1,madd(b.z,v2,b.w*v3))); |
131 | } |
132 | |
133 | __forceinline void eval(const float t, Vertex& p, Vertex& dp, Vertex& ddp) const |
134 | { |
135 | p = eval(t); |
136 | dp = eval_du(t); |
137 | ddp = eval_dudu(t); |
138 | } |
139 | |
140 | template<int M> |
141 | __forceinline Vec4vf<M> veval(const vfloat<M>& t) const |
142 | { |
143 | const Vec4vf<M> b = BSplineBasis::eval(t); |
144 | return madd(b.x, Vec4vf<M>(v0), madd(b.y, Vec4vf<M>(v1), madd(b.z, Vec4vf<M>(v2), b.w * Vec4vf<M>(v3)))); |
145 | } |
146 | |
147 | template<int M> |
148 | __forceinline Vec4vf<M> veval_du(const vfloat<M>& t) const |
149 | { |
150 | const Vec4vf<M> b = BSplineBasis::derivative(t); |
151 | return madd(b.x, Vec4vf<M>(v0), madd(b.y, Vec4vf<M>(v1), madd(b.z, Vec4vf<M>(v2), b.w * Vec4vf<M>(v3)))); |
152 | } |
153 | |
154 | template<int M> |
155 | __forceinline Vec4vf<M> veval_dudu(const vfloat<M>& t) const |
156 | { |
157 | const Vec4vf<M> b = BSplineBasis::derivative2(t); |
158 | return madd(b.x, Vec4vf<M>(v0), madd(b.y, Vec4vf<M>(v1), madd(b.z, Vec4vf<M>(v2), b.w * Vec4vf<M>(v3)))); |
159 | } |
160 | |
161 | template<int M> |
162 | __forceinline void veval(const vfloat<M>& t, Vec4vf<M>& p, Vec4vf<M>& dp) const |
163 | { |
164 | p = veval<M>(t); |
165 | dp = veval_du<M>(t); |
166 | } |
167 | |
168 | template<int M> |
169 | __forceinline Vec4vf<M> eval0(const int ofs, const int size) const |
170 | { |
171 | assert(size <= PrecomputedBSplineBasis::N); |
172 | assert(ofs <= size); |
173 | return madd(vfloat<M>::loadu(&bspline_basis0.c0[size][ofs]), Vec4vf<M>(v0), |
174 | madd(vfloat<M>::loadu(&bspline_basis0.c1[size][ofs]), Vec4vf<M>(v1), |
175 | madd(vfloat<M>::loadu(&bspline_basis0.c2[size][ofs]), Vec4vf<M>(v2), |
176 | vfloat<M>::loadu(&bspline_basis0.c3[size][ofs]) * Vec4vf<M>(v3)))); |
177 | } |
178 | |
179 | template<int M> |
180 | __forceinline Vec4vf<M> eval1(const int ofs, const int size) const |
181 | { |
182 | assert(size <= PrecomputedBSplineBasis::N); |
183 | assert(ofs <= size); |
184 | return madd(vfloat<M>::loadu(&bspline_basis1.c0[size][ofs]), Vec4vf<M>(v0), |
185 | madd(vfloat<M>::loadu(&bspline_basis1.c1[size][ofs]), Vec4vf<M>(v1), |
186 | madd(vfloat<M>::loadu(&bspline_basis1.c2[size][ofs]), Vec4vf<M>(v2), |
187 | vfloat<M>::loadu(&bspline_basis1.c3[size][ofs]) * Vec4vf<M>(v3)))); |
188 | } |
189 | |
190 | template<int M> |
191 | __forceinline Vec4vf<M> derivative0(const int ofs, const int size) const |
192 | { |
193 | assert(size <= PrecomputedBSplineBasis::N); |
194 | assert(ofs <= size); |
195 | return madd(vfloat<M>::loadu(&bspline_basis0.d0[size][ofs]), Vec4vf<M>(v0), |
196 | madd(vfloat<M>::loadu(&bspline_basis0.d1[size][ofs]), Vec4vf<M>(v1), |
197 | madd(vfloat<M>::loadu(&bspline_basis0.d2[size][ofs]), Vec4vf<M>(v2), |
198 | vfloat<M>::loadu(&bspline_basis0.d3[size][ofs]) * Vec4vf<M>(v3)))); |
199 | } |
200 | |
201 | template<int M> |
202 | __forceinline Vec4vf<M> derivative1(const int ofs, const int size) const |
203 | { |
204 | assert(size <= PrecomputedBSplineBasis::N); |
205 | assert(ofs <= size); |
206 | return madd(vfloat<M>::loadu(&bspline_basis1.d0[size][ofs]), Vec4vf<M>(v0), |
207 | madd(vfloat<M>::loadu(&bspline_basis1.d1[size][ofs]), Vec4vf<M>(v1), |
208 | madd(vfloat<M>::loadu(&bspline_basis1.d2[size][ofs]), Vec4vf<M>(v2), |
209 | vfloat<M>::loadu(&bspline_basis1.d3[size][ofs]) * Vec4vf<M>(v3)))); |
210 | } |
211 | |
212 | /* calculates bounds of bspline curve geometry */ |
213 | __forceinline BBox3fa accurateRoundBounds() const |
214 | { |
215 | const int N = 7; |
216 | const float scale = 1.0f/(3.0f*(N-1)); |
217 | Vec4vfx pl(pos_inf), pu(neg_inf); |
218 | for (int i=0; i<=N; i+=VSIZEX) |
219 | { |
220 | vintx vi = vintx(i)+vintx(step); |
221 | vboolx valid = vi <= vintx(N); |
222 | const Vec4vfx p = eval0<VSIZEX>(i,N); |
223 | const Vec4vfx dp = derivative0<VSIZEX>(i,N); |
224 | const Vec4vfx pm = p-Vec4vfx(scale)*select(vi!=vintx(0),dp,Vec4vfx(zero)); |
225 | const Vec4vfx pp = p+Vec4vfx(scale)*select(vi!=vintx(N),dp,Vec4vfx(zero)); |
226 | pl = select(valid,min(pl,p,pm,pp),pl); // FIXME: use masked min |
227 | pu = select(valid,max(pu,p,pm,pp),pu); // FIXME: use masked min |
228 | } |
229 | const Vec3fa lower(reduce_min(pl.x),reduce_min(pl.y),reduce_min(pl.z)); |
230 | const Vec3fa upper(reduce_max(pu.x),reduce_max(pu.y),reduce_max(pu.z)); |
231 | const float r_min = reduce_min(pl.w); |
232 | const float r_max = reduce_max(pu.w); |
233 | const Vec3fa upper_r = Vec3fa(max(abs(r_min),abs(r_max))); |
234 | return enlarge(BBox3fa(lower,upper),upper_r); |
235 | } |
236 | |
237 | /* calculates bounds when tessellated into N line segments */ |
238 | __forceinline BBox3fa accurateFlatBounds(int N) const |
239 | { |
240 | if (likely(N == 4)) |
241 | { |
242 | const Vec4vf4 pi = eval0<4>(0,4); |
243 | const Vec3fa lower(reduce_min(pi.x),reduce_min(pi.y),reduce_min(pi.z)); |
244 | const Vec3fa upper(reduce_max(pi.x),reduce_max(pi.y),reduce_max(pi.z)); |
245 | const Vec3fa upper_r = Vec3fa(reduce_max(abs(pi.w))); |
246 | const Vec3ff pe = end(); |
247 | return enlarge(BBox3fa(min(lower,pe),max(upper,pe)),max(upper_r,Vec3fa(abs(pe.w)))); |
248 | } |
249 | else |
250 | { |
251 | Vec3vfx pl(pos_inf), pu(neg_inf); vfloatx ru(0.0f); |
252 | for (int i=0; i<=N; i+=VSIZEX) |
253 | { |
254 | vboolx valid = vintx(i)+vintx(step) <= vintx(N); |
255 | const Vec4vfx pi = eval0<VSIZEX>(i,N); |
256 | |
257 | pl.x = select(valid,min(pl.x,pi.x),pl.x); // FIXME: use masked min |
258 | pl.y = select(valid,min(pl.y,pi.y),pl.y); |
259 | pl.z = select(valid,min(pl.z,pi.z),pl.z); |
260 | |
261 | pu.x = select(valid,max(pu.x,pi.x),pu.x); // FIXME: use masked min |
262 | pu.y = select(valid,max(pu.y,pi.y),pu.y); |
263 | pu.z = select(valid,max(pu.z,pi.z),pu.z); |
264 | |
265 | ru = select(valid,max(ru,abs(pi.w)),ru); |
266 | } |
267 | const Vec3fa lower(reduce_min(pl.x),reduce_min(pl.y),reduce_min(pl.z)); |
268 | const Vec3fa upper(reduce_max(pu.x),reduce_max(pu.y),reduce_max(pu.z)); |
269 | const Vec3fa upper_r(reduce_max(ru)); |
270 | return enlarge(BBox3fa(lower,upper),upper_r); |
271 | } |
272 | } |
273 | |
274 | friend __forceinline embree_ostream operator<<(embree_ostream cout, const BSplineCurveT& curve) { |
275 | return cout << "BSplineCurve { v0 = " << curve.v0 << ", v1 = " << curve.v1 << ", v2 = " << curve.v2 << ", v3 = " << curve.v3 << " }" ; |
276 | } |
277 | }; |
278 | |
279 | template<typename Vertex> |
280 | __forceinline void convert(const BezierCurveT<Vertex>& icurve, BezierCurveT<Vertex>& ocurve) { |
281 | ocurve = icurve; |
282 | } |
283 | |
284 | template<typename Vertex> |
285 | __forceinline void convert(const BSplineCurveT<Vertex>& icurve, BSplineCurveT<Vertex>& ocurve) { |
286 | ocurve = icurve; |
287 | } |
288 | |
289 | template<typename Vertex> |
290 | __forceinline void convert(const BezierCurveT<Vertex>& icurve, BSplineCurveT<Vertex>& ocurve) |
291 | { |
292 | const Vertex v0 = madd(6.0f,icurve.v0,madd(-7.0f,icurve.v1,2.0f*icurve.v2)); |
293 | const Vertex v1 = msub(2.0f,icurve.v1,icurve.v2); |
294 | const Vertex v2 = msub(2.0f,icurve.v2,icurve.v1); |
295 | const Vertex v3 = madd(2.0f,icurve.v1,madd(-7.0f,icurve.v2,6.0f*icurve.v3)); |
296 | ocurve = BSplineCurveT<Vertex>(v0,v1,v2,v3); |
297 | } |
298 | |
299 | template<typename Vertex> |
300 | __forceinline void convert(const BSplineCurveT<Vertex>& icurve, BezierCurveT<Vertex>& ocurve) |
301 | { |
302 | const Vertex v0 = madd(1.0f/6.0f,icurve.v0,madd(2.0f/3.0f,icurve.v1,1.0f/6.0f*icurve.v2)); |
303 | const Vertex v1 = madd(2.0f/3.0f,icurve.v1,1.0f/3.0f*icurve.v2); |
304 | const Vertex v2 = madd(1.0f/3.0f,icurve.v1,2.0f/3.0f*icurve.v2); |
305 | const Vertex v3 = madd(1.0f/6.0f,icurve.v1,madd(2.0f/3.0f,icurve.v2,1.0f/6.0f*icurve.v3)); |
306 | ocurve = BezierCurveT<Vertex>(v0,v1,v2,v3); |
307 | } |
308 | |
309 | template<typename CurveGeometry> |
310 | __forceinline BSplineCurveT<Vec3ff> enlargeRadiusToMinWidth(const IntersectContext* context, const CurveGeometry* geom, const Vec3fa& ray_org, const BSplineCurveT<Vec3ff>& curve) |
311 | { |
312 | return BSplineCurveT<Vec3ff>(enlargeRadiusToMinWidth(context,geom,ray_org,curve.v0), |
313 | enlargeRadiusToMinWidth(context,geom,ray_org,curve.v1), |
314 | enlargeRadiusToMinWidth(context,geom,ray_org,curve.v2), |
315 | enlargeRadiusToMinWidth(context,geom,ray_org,curve.v3)); |
316 | } |
317 | |
318 | typedef BSplineCurveT<Vec3fa> BSplineCurve3fa; |
319 | } |
320 | |
321 | |