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