1 // Copyright Contributors to the Open Shading Language project.
2 // SPDX-License-Identifier: BSD-3-Clause
3 // https://github.com/AcademySoftwareFoundation/OpenShadingLanguage
4
5 #include <limits>
6
7
8 #include "oslexec_pvt.h"
9 #include <OSL/oslnoise.h>
10 #include <OSL/dual_vec.h>
11 #include <OSL/Imathx/Imathx.h>
12
13 #include <OpenImageIO/fmath.h>
14
15 OSL_NAMESPACE_ENTER
16
17 namespace pvt {
18
19 // TODO: It would be preferable to use the Imath versions of these functions in
20 // all cases, but these templates should suffice until a more complete
21 // device-friendly version of Imath is available.
22 namespace hostdevice {
23 template <typename T> inline OSL_HOSTDEVICE T clamp (T x, T lo, T hi);
24 #ifndef __CUDA_ARCH__
25 template <> inline OSL_HOSTDEVICE double clamp<double> (double x, double lo, double hi) { return Imath::clamp (x, lo, hi); }
26 template <> inline OSL_HOSTDEVICE float clamp<float> (float x, float lo, float hi) { return Imath::clamp (x, lo, hi); }
27 #else
28 template <> inline OSL_HOSTDEVICE double clamp<double> (double x, double lo, double hi) { return (x < lo) ? lo : ((x > hi) ? hi : x); }
29 template <> inline OSL_HOSTDEVICE float clamp<float> (float x, float lo, float hi) { return (x < lo) ? lo : ((x > hi) ? hi : x); }
30 #endif
31 }
32
33 static OSL_DEVICE constexpr float Gabor_Frequency = 2.0;
34 static OSL_DEVICE constexpr float Gabor_Impulse_Weight = 1.0f;
35
36 // The Gabor kernel in theory has infinite support (its envelope is
37 // a Gaussian). To restrict the distance at which we must sum the
38 // kernels, we only consider those whose Gaussian envelopes are
39 // above the truncation threshold, as a portion of the Gaussian's
40 // peak value.
41 static OSL_DEVICE const float Gabor_Truncate = 0.02f;
42
43
44
45 // Very fast random number generator based on [Borosh & Niederreiter 1983]
46 // linear congruential generator.
47 class fast_rng {
48 public:
49 // seed based on the cell containing P
50 OSL_DEVICE
51 fast_rng (const Vec3 &p, int seed=0) {
52 // Use guts of cellnoise
53 m_seed = inthash(unsigned(OIIO::ifloor(p.x)),
54 unsigned(OIIO::ifloor(p.y)),
55 unsigned(OIIO::ifloor(p.z)),
56 unsigned(seed));
57 if (! m_seed)
58 m_seed = 1;
59 }
60 // Return uniform on [0,1)
61 OSL_HOSTDEVICE
operator()62 float operator() () {
63 return (m_seed *= 3039177861u) / float(UINT_MAX);
64 }
65 // Return poisson distribution with the given mean
66 OSL_HOSTDEVICE
poisson(float mean)67 int poisson (float mean) {
68 float g = expf (-mean);
69 unsigned int em = 0;
70 float t = (*this)();
71 while (t > g) {
72 ++em;
73 t *= (*this)();
74 }
75 return em;
76 }
77 private:
78 unsigned int m_seed;
79 };
80
81 // The Gabor kernel is a harmonic (cosine) modulated by a Gaussian
82 // envelope. This version is augmented with a phase, per [Lagae2011].
83 // \param weight magnitude of the pulse
84 // \param omega orientation of the harmonic
85 // \param phi phase of the harmonic.
86 // \param bandwidth width of the gaussian envelope (called 'a'
87 // in [Lagae09].
88 // \param x the position being sampled
89 template <class VEC> // VEC should be Vec3 or Vec2
90 inline OSL_HOSTDEVICE Dual2<float>
gabor_kernel(const Dual2<float> & weight,const VEC & omega,const Dual2<float> & phi,float bandwidth,const Dual2<VEC> & x)91 gabor_kernel (const Dual2<float> &weight, const VEC &omega,
92 const Dual2<float> &phi, float bandwidth, const Dual2<VEC> &x)
93 {
94 // see Equation 1
95 Dual2<float> g = exp (float(-M_PI) * (bandwidth * bandwidth) * dot(x,x));
96 Dual2<float> h = cos (float(M_TWO_PI) * dot(omega,x) + phi);
97 return weight * g * h;
98 }
99
100
101
102 inline OSL_HOSTDEVICE void
slice_gabor_kernel_3d(const Dual2<float> & d,float w,float a,const Vec3 & omega,float phi,Dual2<float> & w_s,Vec2 & omega_s,Dual2<float> & phi_s)103 slice_gabor_kernel_3d (const Dual2<float> &d, float w, float a,
104 const Vec3 &omega, float phi,
105 Dual2<float> &w_s, Vec2 &omega_s, Dual2<float> &phi_s)
106 {
107 // Equation 6
108 w_s = w * exp(float(-M_PI) * (a*a)*(d*d));
109 //omega_s[0] = omega[0];
110 //omega_s[1] = omega[1];
111 //phi_s = phi - float(M_TWO_PI) * d * omega[2];
112 omega_s.x = omega.x;
113 omega_s.y = omega.y;
114 // A.W. think this was a bug, supposed to be omega.z not omega.x;
115 //phi_s = phi - float(M_TWO_PI) * d * omega.x;
116 phi_s = phi - float(M_TWO_PI) * d * omega.z;
117 }
118
119
120 namespace {
121 inline OSL_HOSTDEVICE Vec2
gabor_mul_m22_v2(const Matrix22 & m,const Vec2 & v)122 gabor_mul_m22_v2(const Matrix22& m, const Vec2& v)
123 {
124 float a = v.x * m.x[0][0] + v.y * m.x[1][0];
125 float b = v.x * m.x[0][1] + v.y * m.x[1][1];
126 return Vec2(a, b);
127 }
128
129 }
130
131 static OSL_HOSTDEVICE void
filter_gabor_kernel_2d(const Matrix22 & filter,const Dual2<float> & w,float a,const Vec2 & omega,const Dual2<float> & phi,Dual2<float> & w_f,float & a_f,Vec2 & omega_f,Dual2<float> & phi_f)132 filter_gabor_kernel_2d (const Matrix22 &filter, const Dual2<float> &w, float a,
133 const Vec2 &omega, const Dual2<float> &phi,
134 Dual2<float> &w_f, float &a_f,
135 Vec2 &omega_f, Dual2<float> &phi_f)
136 {
137 // Equation 10
138 Matrix22 Sigma_f = filter;
139 Dual2<float> c_G = w;
140 Vec2 mu_G = omega;
141 Matrix22 Sigma_G = (a * a / float(M_TWO_PI)) * Matrix22();
142 float c_F = 1.0f / (float(M_TWO_PI) * sqrtf(determinant(Sigma_f)));
143 Matrix22 Sigma_F = float(1.0 / (4.0 * M_PI * M_PI)) * Sigma_f.inverse();
144 Matrix22 Sigma_G_Sigma_F = Sigma_G + Sigma_F;
145 Dual2<float> c_GF = c_F * c_G
146 * (1.0f / (float(M_TWO_PI) * sqrtf(determinant(Sigma_G_Sigma_F))))
147 * expf(-0.5f * dot(gabor_mul_m22_v2(Sigma_G_Sigma_F.inverse(),mu_G), mu_G));
148 Matrix22 Sigma_G_i = Sigma_G.inverse();
149 Matrix22 Sigma_GF = (Sigma_F.inverse() + Sigma_G_i).inverse();
150 Matrix22 Sigma_GF_Gi = Sigma_GF * Sigma_G_i;
151 Vec2 mu_GF = gabor_mul_m22_v2(Sigma_GF_Gi, mu_G);
152 w_f = c_GF;
153 a_f = sqrtf(M_TWO_PI * sqrtf(determinant(Sigma_GF)));
154 omega_f = mu_GF;
155 phi_f = phi;
156 }
157
158
159 OSL_FORCEINLINE OSL_HOSTDEVICE float
wrap(float s,float period)160 wrap (float s, float period)
161 {
162 period = floorf (period);
163 if (period < 1.0f)
164 period = 1.0f;
165 return s - period * floorf (s / period);
166 }
167
168
169 // avoid aliasing issues
170 static OSL_FORCEINLINE OSL_HOSTDEVICE Vec3
wrap(const Vec3 & s,const Vec3 & period)171 wrap (const Vec3 &s, const Vec3 &period)
172 {
173 return Vec3 (wrap (s.x, period.x),
174 wrap (s.y, period.y),
175 wrap (s.z, period.z));
176 }
177
178
179 // Normalize v and set a and b to be unit vectors (any two unit vectors)
180 // that are orthogonal to v and each other. We get the first
181 // orthonormal by taking the cross product of v and (1,0,0), unless v
182 // points roughly toward (1,0,0), in which case we cross with (0,1,0).
183 // Either way, we get something orthogonal. Then cross(v,a) is mutually
184 // orthogonal to the other two.
185 inline OSL_HOSTDEVICE void
make_orthonormals(Vec3 & v,Vec3 & a,Vec3 & b)186 make_orthonormals (Vec3 &v, Vec3 &a, Vec3 &b)
187 {
188 // avoid aliasing issues by not using the [] operator
189 v.normalize();
190 if (fabsf(v.x) < 0.9f)
191 a.setValue (0.0f, v.z, -v.y); // v X (1,0,0)
192 else
193 a.setValue (-v.z, 0.0f, v.x); // v X (0,1,0)
194 a.normalize ();
195 b = v.cross (a);
196 // b.normalize (); // note: not necessary since v is unit length
197 }
198
199
200
201 // Helper function: per-component 'floor' of a Dual2<Vec3>.
202 inline OSL_HOSTDEVICE Vec3
floor(const Dual2<Vec3> & vd)203 floor (const Dual2<Vec3> &vd)
204 {
205 // avoid aliasing issues by not using the [] operator
206 const Vec3 &v (vd.val());
207 return Vec3 (floorf(v.x), floorf(v.y), floorf(v.z));
208 }
209
210 } // namespace pvt
211
212 OSL_NAMESPACE_EXIT
213