yum-archive/Tooner

A toon shader for Unity's BIRP.

git clone https://git.yummers.dev/yum-archive/Tooner

yumAdd bayer matrices & grabpass-based lens shader4bd9c2b

master
branchesmasterall refs →
7.8 KiB322 linesraw
1#include "pema99.cginc"
2
3#ifndef __MATH_INC
4#define __MATH_INC
5
6#define PI 3.14159265
7#define TAU PI * 2.0
8#define PHI 1.618033989
9
10static const float BayerM4x4[16] = {
11    0.0/15.0,   8.0/15.0,   2.0/15.0,   10.0/15.0,
12    12.0/15.0,  4.0/15.0,   14.0/15.0,  6.0/15.0,
13    3.0/15.0,   11.0/15.0,  1.0/15.0,   9.0/15.0,
14    15.0/15.0,  7.0/15.0,   13.0/15.0,  5.0/15.0
15};
16
17// xdd
18static const float BayerM8x8[64] = {
19    0.0/63.0, 32.0/63.0,  8.0/63.0, 40.0/63.0,  2.0/63.0, 34.0/63.0, 10.0/63.0, 42.0/63.0,
20   48.0/63.0, 16.0/63.0, 56.0/63.0, 24.0/63.0, 50.0/63.0, 18.0/63.0, 58.0/63.0, 26.0/63.0,
21   12.0/63.0, 44.0/63.0,  4.0/63.0, 36.0/63.0, 14.0/63.0, 46.0/63.0,  6.0/63.0, 38.0/63.0,
22   60.0/63.0, 28.0/63.0, 52.0/63.0, 20.0/63.0, 62.0/63.0, 30.0/63.0, 54.0/63.0, 22.0/63.0,
23    3.0/63.0, 35.0/63.0, 11.0/63.0, 43.0/63.0,  1.0/63.0, 33.0/63.0,  9.0/63.0, 41.0/63.0,
24   51.0/63.0, 19.0/63.0, 59.0/63.0, 27.0/63.0, 49.0/63.0, 17.0/63.0, 57.0/63.0, 25.0/63.0,
25   15.0/63.0, 47.0/63.0,  7.0/63.0, 39.0/63.0, 13.0/63.0, 45.0/63.0,  5.0/63.0, 37.0/63.0,
26   63.0/63.0, 31.0/63.0, 55.0/63.0, 23.0/63.0, 61.0/63.0, 29.0/63.0, 53.0/63.0, 21.0/63.0
27};
28
29// Hacky parameterizable whiteout blending. Probably some big mistakes but it
30// passes the eyeball test.
31// At w=0.5, this looks kinda like whiteout blending.
32// At w=0, this returns n0.
33// At w=1, this returns n1.
34#define MY_BLEND_NORMALS(n0, n1, w) normalize(float3((n0.xy * (1 - w) + n1.xy * w), lerp(1, n0.z, (1-w)) * lerp(1, n1.z, w)))
35
36float golden_lds(uint i)
37{
38  return glsl_mod(1.61803398875 * float(i), 1);
39}
40
41// Complex numbers
42typedef float2 complex;
43
44float creal(complex z0)
45{
46  return z0.x;
47}
48
49float cimag(complex z0)
50{
51  return z0.y;
52}
53
54float cnorm2(complex z0)
55{
56  return z0.x * z0.x + z0.y * z0.y;
57}
58
59float cnorm(complex z0)
60{
61  return sqrt(cnorm2(z0));
62}
63
64complex cconjugate(complex z)
65{
66  return float2(z.x, -z.y);
67}
68
69complex cmul(complex z0, complex z1)
70{
71  return float2(z0.x * z1.x - z0.y * z1.y, z0.x * z1.y + z0.y * z1.x);
72}
73
74complex cdiv(complex z0, complex z1)
75{
76  float re = creal(z0) * creal(z1) + cimag(z0) * cimag(z1);
77  float im = cimag(z0) * cimag(z1) - creal(z0) * creal(z1);
78  float z1_norm2 = cnorm2(z1);
79  return float2(re, im) / z1_norm2;
80}
81
82// Evaluates z0**n.
83// Uses Euler's identity to support fractional values of `n`.
84// Expensive.
85complex cpow_fractional(complex z0, float n)
86{
87  float r = sqrt(z0.x * z0.x + z0.y * z0.y);
88  float t = atan(z0.y / z0.x);
89  return pow(r, n) * float2(cos(t * n), sin(t * n));
90}
91
92// Evaluates z0**n.
93// Utilizes recursive squaring to support high values of `n`.
94// Cheap.
95complex cpow(complex z0, uint n)
96{
97  if (n == 0) {
98    return 1;
99  }
100
101  complex z = z0;
102  while (n > 1) {
103    if (n % 2 == 0) {
104      z = cmul(z, z);
105      n /= 2;
106    } else {
107      z = cmul(z, z0);
108      n -= 1;
109    }
110  }
111  return z;
112}
113
114// Quaternions
115float4 qmul(float4 q1, float4 q2)
116{
117  return float4(
118      q2.xyz * q1.w + q1.xyz * q2.w + cross(q1.xyz, q2.xyz),
119      q1.w * q2.w - dot(q1.xyz, q2.xyz)
120      );
121}
122
123// Vector rotation with a quaternion
124// http://mathworld.wolfram.com/Quaternion.html
125float3 rotate_vector(float3 v, float4 r)
126{
127  float4 r_c = r * float4(-1, -1, -1, 1);
128  return qmul(r, qmul(float4(v, 0), r_c)).xyz;
129}
130
131float4 get_quaternion(float3 axis_normal, float theta) {
132  return float4(
133      axis_normal * sin(theta / 2), cos(theta / 2));
134}
135
136// Differentiable approximation of the standard `max` function.
137float dmax(float a, float b, float k)
138{
139  return log2(exp2(k * a) + exp2(k * b)) / k;
140}
141
142// Differentiable approximation of the standard `min` function.
143float dmin(float a, float b, float k)
144{
145  return -1.0 * dmax(-1.0 * a, -1.0 * b, k);
146}
147
148float dabs(float a, float k)
149{
150  return log2(exp2(k * a) + exp2(-1.0 * k * a));
151}
152
153float dsaturate(float x, float k)
154{
155  return dmin(dmax(x, 0, k), 1, k);
156}
157
158float dclamp(float x, float lo, float hi, float k)
159{
160  return dmin(dmax(x, lo, k), hi, k);
161}
162
163float rand(uint seed) {
164  seed = seed * 747796405 + 2891336453;
165  uint result = ((seed >> ((seed >> 28) + 4)) ^ seed) * 277803737;
166  result = (result >> 22) ^ result;
167  return result / 4294967295.0;
168}
169
170// Generate a random number on [0, 1].
171float rand2(float2 p)
172{
173  return frac(sin(dot(p,
174          float2(12.9898, 78.233)))
175      * 43758.5453123);
176}
177
178// Generate a random number on [0, 1].
179float rand3(float3 p)
180{
181  return glsl_mod(sin(dot(p, float3(151.0, 157.0, 163.0))) * 997.0, 1.0);
182}
183
184float length2(float2 p)
185{
186  return p.x * p.x + p.y * p.y;
187}
188
189// 3 dimensional value noise. `p` is assumed to be a point inside a unit cube.
190// Theory: https://en.wikipedia.org/wiki/Value_noise
191float vnoise3d(float3 p)
192{
193  float3 pu = floor(p);
194  float3 pv = glsl_mod(frac(p), 1.0);
195
196  // Assign random numbers to the corner of a cube.
197  float n000 = rand3(pu + float3(0,0,0));
198  float n001 = rand3(pu + float3(0,0,1));
199  float n010 = rand3(pu + float3(0,1,0));
200  float n011 = rand3(pu + float3(0,1,1));
201  float n100 = rand3(pu + float3(1,0,0));
202  float n101 = rand3(pu + float3(1,0,1));
203  float n110 = rand3(pu + float3(1,1,0));
204  float n111 = rand3(pu + float3(1,1,1));
205
206  float n00 = lerp(n000, n001, pv.z);
207  float n01 = lerp(n010, n011, pv.z);
208  float n10 = lerp(n100, n101, pv.z);
209  float n11 = lerp(n110, n111, pv.z);
210
211  float n0 = lerp(n00, n01, pv.y);
212  float n1 = lerp(n10, n11, pv.y);
213  
214  float n = lerp(n0, n1, pv.x);
215
216  return n;
217}
218
219float fbm(float3 p, const int n_octaves, float w)
220{
221  float g = exp2(-w);
222  float a = 1.0;
223  float p_scale = 1.0;
224
225  float res = 0.0;
226  for (int i = 0; i < n_octaves; i++) {
227    res += a * vnoise3d(p * p_scale);
228
229    p_scale /= w;
230    a *= g;
231  }
232  return res;
233}
234
235float median(float x, float y, float z)
236{
237  return max(min(x, y), min(max(x, y), z));
238}
239
240float median(float3 x)
241{
242  return median(x.x, x.y, x.z);
243}
244
245// Yoinked from here
246//   https://www.scratchapixel.com/lessons/3d-basic-rendering/minimal-ray-tracer-rendering-simple-shapes/ray-sphere-intersection.html
247bool solveQuadratic(float a, float b, float c, out float x0, out float x1)
248{
249  float discriminant = b * b - 4 * a * c;
250  if (discriminant < 0) {
251    return false;
252  } else if (discriminant == 0) {
253    x0 = -0.5 * b / a;
254    x1 = x0;
255  } else {
256    float q = (b > 0) ?
257      -0.5 * (b + sqrt(discriminant)) :
258      -0.5 * (b - sqrt(discriminant));
259    x0 = q/a;
260    x1 = c/q;
261  }
262  float tmp_min = min(x0, x1);
263  float tmp_max = max(x0, x1);
264  x0 = tmp_min;
265  x1 = tmp_max;
266  return true;
267}
268
269float determinant(float3x3 m)
270{
271  return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
272            - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]))
273            + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
274}
275
276float3x3 invert(float3x3 m)
277{
278  float det = determinant(m);
279
280  float3x3 adj;
281  adj[0][0] =  (m[1][1] * m[2][2] - m[1][2] * m[2][1]);
282  adj[0][1] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]);
283  adj[0][2] =  (m[0][1] * m[1][2] - m[0][2] * m[1][1]);
284  
285  adj[1][0] = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]);
286  adj[1][1] =  (m[0][0] * m[2][2] - m[0][2] * m[2][0]);
287  adj[1][2] = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]);
288  
289  adj[2][0] =  (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
290  adj[2][1] = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]);
291  adj[2][2] =  (m[0][0] * m[1][1] - m[0][1] * m[1][0]);
292
293  return adj * (1.0 / det);
294}
295
296// Return largest number which divides into both 'a' and 'b'.
297// Uses the Euclidean algorithm: repeatedly divide larger by smaller number.
298uint gcd(uint a, uint b)
299{
300    #define GCD_MAX_ITER 24
301    for (uint i = 0; i < GCD_MAX_ITER; i++) {
302        if (b == 0) {
303            return a;
304        }
305        a = a % b;
306        // Swap a and b
307        uint tmp = a;
308        a = b;
309        b = tmp;
310    }
311    return 1;
312}
313
314float wrapNoL(float NoL, float factor) {
315    // Apply wrapped lighting correction
316		// https://www.iro.umontreal.ca/~derek/files/jgt_wrap_final.pdf
317    //float4 wrapped = (NoL + 1) * (NoL + 1) * .25;
318    return pow(max(1E-4, (NoL + factor) / (1 + factor)), 1 + factor);
319}
320
321#endif  // __MATH_INC
322