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