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> eval(const T & u)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> derivative(const T & u)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> derivative2(const T & u)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: PrecomputedBSplineBasisPrecomputedBSplineBasis56 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 BSplineCurveTBSplineCurveT81 __forceinline BSplineCurveT() {} 82 BSplineCurveTBSplineCurveT83 __forceinline BSplineCurveT(const Vertex& v0, const Vertex& v1, const Vertex& v2, const Vertex& v3) 84 : v0(v0), v1(v1), v2(v2), v3(v3) {} 85 beginBSplineCurveT86 __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 endBSplineCurveT90 __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 centerBSplineCurveT94 __forceinline Vertex center() const { 95 return 0.25f*(v0+v1+v2+v3); 96 } 97 boundsBSplineCurveT98 __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 xfm_prBSplineCurveT106 __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 evalBSplineCurveT115 __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 eval_duBSplineCurveT121 __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 eval_duduBSplineCurveT127 __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 evalBSplineCurveT133 __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> vevalBSplineCurveT141 __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> veval_duBSplineCurveT148 __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> veval_duduBSplineCurveT155 __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> vevalBSplineCurveT162 __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> eval0BSplineCurveT169 __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> eval1BSplineCurveT180 __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> derivative0BSplineCurveT191 __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> derivative1BSplineCurveT202 __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 */ accurateRoundBoundsBSplineCurveT213 __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 */ accurateFlatBoundsBSplineCurveT238 __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> convert(const BezierCurveT<Vertex> & icurve,BezierCurveT<Vertex> & ocurve)280 __forceinline void convert(const BezierCurveT<Vertex>& icurve, BezierCurveT<Vertex>& ocurve) { 281 ocurve = icurve; 282 } 283 284 template<typename Vertex> convert(const BSplineCurveT<Vertex> & icurve,BSplineCurveT<Vertex> & ocurve)285 __forceinline void convert(const BSplineCurveT<Vertex>& icurve, BSplineCurveT<Vertex>& ocurve) { 286 ocurve = icurve; 287 } 288 289 template<typename Vertex> convert(const BezierCurveT<Vertex> & icurve,BSplineCurveT<Vertex> & ocurve)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> convert(const BSplineCurveT<Vertex> & icurve,BezierCurveT<Vertex> & ocurve)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> enlargeRadiusToMinWidth(const IntersectContext * context,const CurveGeometry * geom,const Vec3fa & ray_org,const BSplineCurveT<Vec3ff> & curve)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