1 // Copyright Contributors to the Open Shading Language project. 2 // SPDX-License-Identifier: BSD-3-Clause 3 // https://github.com/AcademySoftwareFoundation/OpenShadingLanguage 4 5 #pragma once 6 7 OSL_NAMESPACE_ENTER 8 9 namespace pvt { 10 11 // declaring a Spline namespace to avoid polluting the existing 12 // namespaces with all these templated helper functions. 13 namespace Spline { 14 15 struct SplineBasis { 16 int basis_step; 17 float basis[4][4]; 18 }; 19 20 // ======================================================== 21 // 22 // Interpolation bases for splines 23 // 24 // The order here is very important for the SplineInterp::create 25 // constructor below. Any additional modes should be added to 26 // the end, or SplineInterp::create updated as well. 27 // 28 // ======================================================== 29 30 enum { 31 kCatmullRom, 32 kBezier, 33 kBSpline, 34 kHermite, 35 kLinear, 36 kConstant, 37 kNumSplineTypes 38 }; 39 40 OSL_CONSTANT_DATA const static SplineBasis gBasisSet[kNumSplineTypes] = { 41 // 42 // catmullrom 43 // 44 { 1, { {(-1.0f/2.0f), ( 3.0f/2.0f), (-3.0f/2.0f), ( 1.0f/2.0f)}, 45 {( 2.0f/2.0f), (-5.0f/2.0f), ( 4.0f/2.0f), (-1.0f/2.0f)}, 46 {(-1.0f/2.0f), ( 0.0f/2.0f), ( 1.0f/2.0f), ( 0.0f/2.0f)}, 47 {( 0.0f/2.0f), ( 2.0f/2.0f), ( 0.0f/2.0f), ( 0.0f/2.0f)} } }, 48 // 49 // bezier 50 // 51 { 3, { {-1, 3, -3, 1}, 52 { 3, -6, 3, 0}, 53 {-3, 3, 0, 0}, 54 { 1, 0, 0, 0} } }, 55 // 56 // bspline 57 // 58 { 1, { {(-1.0f/6.0f), ( 3.0f/6.0f), (-3.0f/6.0f), (1.0f/6.0f)}, 59 {( 3.0f/6.0f), (-6.0f/6.0f), ( 3.0f/6.0f), (0.0f/6.0f)}, 60 {(-3.0f/6.0f), ( 0.0f/6.0f), ( 3.0f/6.0f), (0.0f/6.0f)}, 61 {( 1.0f/6.0f), ( 4.0f/6.0f), ( 1.0f/6.0f), (0.0f/6.0f)} } }, 62 // 63 // hermite 64 // 65 { 2, { { 2, 1, -2, 1}, 66 {-3, -2, 3, -1}, 67 { 0, 1, 0, 0}, 68 { 1, 0, 0, 0} } }, 69 // 70 // linear 71 // 72 { 1, { {0, 0, 0, 0}, 73 {0, 0, 0, 0}, 74 {0, -1, 1, 0}, 75 {0, 1, 0, 0} } }, 76 // 77 // constant 78 // 79 { 1, { {0, 0, 0, 0}, 80 {0, 0, 0, 0}, 81 {0, 0, 0, 0}, 82 {0, 0, 0, 0} } } 83 }; 84 85 struct SplineInterp { 86 const SplineBasis& spline; 87 const bool constant; 88 createSplineInterp89 OSL_HOSTDEVICE static SplineInterp create(StringParam basis_name) 90 { 91 if (basis_name == StringParams::catmullrom) 92 return { gBasisSet[kCatmullRom], false }; 93 if (basis_name == StringParams::bezier) 94 return { gBasisSet[kBezier], false }; 95 if (basis_name == StringParams::bspline) 96 return { gBasisSet[kBSpline], false }; 97 if (basis_name == StringParams::hermite) 98 return { gBasisSet[kHermite], false }; 99 if (basis_name == StringParams::constant) 100 return { gBasisSet[kConstant], true }; 101 102 // Default to linear 103 return { gBasisSet[kLinear], false }; 104 } 105 106 107 // We need to know explicitly whether the knots have 108 // derivatives associated with them because of the way 109 // Dual2<T> forms of arrays are stored.. Arrays with 110 // derivatives are stored: 111 // T T T... T.dx T.dx T.dx... T.dy T.dy T.dy... 112 // This means, we need to explicitly construct the Dual2<T> 113 // form of the knots on the fly. 114 // if 'is_dual' == true, then OUTTYPE == Dual2<INTYPE> 115 // if 'is_dual' == false, then OUTTYPE == INTYPE 116 117 // This functor will extract a T or a Dual2<T> type from a VaryingRef array 118 template <class OUTTYPE, class INTYPE, bool is_dual> 119 struct extractValueFromArray 120 { 121 OSL_HOSTDEVICE OUTTYPE operator()(const INTYPE *value, int array_length, int idx); 122 }; 123 124 template <class OUTTYPE, class INTYPE> 125 struct extractValueFromArray<OUTTYPE, INTYPE, true> 126 { 127 OSL_HOSTDEVICE OUTTYPE operator()(const INTYPE *value, int array_length, int idx) 128 { 129 return OUTTYPE( value[idx + 0*array_length], 130 value[idx + 1*array_length], 131 value[idx + 2*array_length] ); 132 } 133 }; 134 135 template <class OUTTYPE, class INTYPE> 136 struct extractValueFromArray<OUTTYPE, INTYPE, false> 137 { 138 OSL_HOSTDEVICE OUTTYPE operator()(const INTYPE *value, int /*array_length*/, int idx) 139 { 140 return OUTTYPE( value[idx] ); 141 } 142 }; 143 144 // Spline functor for use with the inverse function 145 template <class RTYPE, class XTYPE> 146 struct SplineFunctor { 147 OSL_HOSTDEVICE SplineFunctor (const SplineInterp& spline_, const float *knots_, 148 int knot_count_, int knot_arraylen_) 149 : spline(spline_), knots(knots_), knot_count(knot_count_), 150 knot_arraylen(knot_arraylen_) { } 151 152 OSL_HOSTDEVICE RTYPE operator() (XTYPE x) { 153 RTYPE v; 154 spline.evaluate<RTYPE,XTYPE,float,float,false> (v, x, knots, knot_count, knot_arraylen); 155 return v; 156 } 157 private: 158 const SplineInterp& spline; 159 const float *knots; 160 int knot_count, knot_arraylen; 161 }; 162 163 template <class RTYPE, class XTYPE, class CTYPE, class KTYPE, bool knot_derivs > 164 OSL_HOSTDEVICE void 165 evaluate(RTYPE &result, XTYPE &xval, const KTYPE *knots, 166 int knot_count, int knot_arraylen) const 167 { 168 using OIIO::clamp; 169 XTYPE x = clamp(xval, XTYPE(0.0), XTYPE(1.0)); 170 int nsegs = ((knot_count - 4) / spline.basis_step) + 1; 171 x = x*(float)nsegs; 172 float seg_x = removeDerivatives(x); 173 int segnum = (int)seg_x; 174 if (segnum < 0) 175 segnum = 0; 176 if (segnum > (nsegs-1)) 177 segnum = nsegs-1; 178 179 if (constant) { 180 // Special case for "constant" basis 181 RTYPE P = removeDerivatives (knots[segnum+1]); 182 assignment (result, P); 183 return; 184 } 185 186 // x is the position along segment 'segnum' 187 x = x - float(segnum); 188 int s = segnum * spline.basis_step; 189 190 // create a functor so we can cleanly(!) extract 191 // the knot elements 192 extractValueFromArray<CTYPE, KTYPE, knot_derivs> myExtract; 193 CTYPE P[4]; 194 for (int k = 0; k < 4; k++) { 195 P[k] = myExtract(knots, knot_arraylen, s + k); 196 } 197 198 CTYPE tk[4]; 199 for (int k = 0; k < 4; k++) { 200 tk[k] = spline.basis[k][0] * P[0] + 201 spline.basis[k][1] * P[1] + 202 spline.basis[k][2] * P[2] + 203 spline.basis[k][3] * P[3]; 204 } 205 206 RTYPE tresult; 207 // The following is what we want, but this gives me template errors 208 // which I'm too lazy to decipher: 209 // tresult = ((tk[0]*x + tk[1])*x + tk[2])*x + tk[3]; 210 tresult = (tk[0] * x + tk[1]); 211 tresult = (tresult * x + tk[2]); 212 tresult = (tresult * x + tk[3]); 213 assignment(result, tresult); 214 } 215 216 // Evaluate the inverse of a spline, i.e., solve for the x for which 217 // spline_evaluate(x) == y. 218 template <class YTYPE> 219 OSL_HOSTDEVICE void 220 inverse (YTYPE &x, YTYPE y, const float *knots, 221 int knot_count, int knot_arraylen) const 222 { 223 // account for out-of-range inputs, just clamp to the values we have 224 int lowindex = spline.basis_step == 1 ? 1 : 0; 225 int highindex = spline.basis_step == 1 ? knot_count-2 : knot_count-1; 226 bool increasing = knots[1] < knots[knot_count-2]; 227 if (increasing) { 228 if (y <= knots[lowindex]) { 229 x = YTYPE(0); 230 return; 231 } 232 if (y >= knots[highindex]) { 233 x = YTYPE(1); 234 return; 235 } 236 } else { 237 if (y >= knots[lowindex]) { 238 x = YTYPE(0); 239 return; 240 } 241 if (y <= knots[highindex]) { 242 x = YTYPE(1); 243 return; 244 } 245 } 246 247 248 SplineFunctor<YTYPE,YTYPE> S (*this, knots, knot_count, knot_arraylen); 249 // Because of the nature of spline interpolation, monotonic knots 250 // can still lead to a non-monotonic curve. To deal with this, 251 // search separately on each spline segment and hope for the best. 252 int nsegs = (knot_count - 4) / spline.basis_step + 1; 253 float nseginv = 1.0f / nsegs; 254 YTYPE r0 = 0.0; 255 x = 0; 256 for (int s = 0; s < nsegs; ++s) { // Search each interval 257 YTYPE r1 = nseginv * (s+1); 258 bool brack; 259 x = OIIO::invert (S, y, r0, r1, 32, YTYPE(1.0e-6), &brack); 260 if (brack) 261 return; 262 r0 = r1; // Start of next interval is end of this one 263 } 264 } 265 }; 266 267 268 }; // namespace Spline 269 }; // namespace pvt 270 OSL_NAMESPACE_EXIT 271 272