yum-archive/2ner

A toon shader for Unity's BIRP.

git clone https://git.yummers.dev/yum-archive/2ner

yumRemove worldPos interpolator4885fb4

master
branchesmasterall refs →
10.5 KiB369 linesraw
1#ifndef __MATH_INC
2#define __MATH_INC
3
4#include "pema99.cginc"
5
6#define PI              3.14159265358979323846264f
7#define TAU             (2.0f * PI)
8#define HALF_PI         (PI * 0.5f)
9#define RCP_PI          (1.0f / PI)
10#define RCP_TAU         (1.0f / TAU)
11#define PHI             1.618033989f
12#define RCP_PHI         0.618033989f
13#define SQRT_2          1.414213562f
14#define SQRT_2_RCP      0.707106781f
15#define RCP_SQRT_2      0.707106781f
16#define RCP_SQRT_3      0.577350269f
17#define TWO_OVER_THREE  0.6666666666666666f
18#define SQRT_3_OVER_2   0.8660254037844386f
19#define EULERS_CONSTANT 2.718281828f
20
21
22float pow5(float x)
23{
24  float tmp = x * x;
25  return (tmp * tmp) * x;
26}
27
28// Wrap NoL. Assume it's already clamped.
29// At k=0, you get standard lambertian shading.
30// At k=0.5, you get half-lambertian shading.
31// At k=1.0, you get flat shading.
32// k must be on [0, 1].
33// Energy preserving, within some small bound.
34float wrapNoL(float NoL, float k) {
35  float lambertian = NoL;
36  float half_lambertian = pow(max(1e-4, (NoL + 0.5f) / (1.0f + 0.5f)), 2);
37  float flat = RCP_PI;
38
39  if (k < 0.5) {
40    return lerp(lambertian, half_lambertian, k * 2.0f);
41  } else {
42    return lerp(half_lambertian, flat, k * 2.0f - 1.0f);
43  }
44}
45
46float halfLambertianNoL(float NoL) {
47		// https://www.iro.umontreal.ca/~derek/files/jgt_wrap_final.pdf
48    float tmp = (NoL + 1)  * 0.5;
49    return tmp * tmp;
50}
51
52float rand1(float p)
53{
54  return frac(sin(p) * 43758.5453123);
55}
56
57float rand2(float2 p)
58{
59  return frac(sin(dot(p, float2(12.9898, 78.233))) * 43758.5453123);
60}
61
62inline float rand3_dot(float3 p)
63{
64  return dot(p, float3(151.0, 157.0, 163.0));
65}
66
67float3 rand3_hash(float3 p)
68{
69    // Improved Murmurhash3 by Squirrel Eiserloh (GDC 2017)
70    p = float3(dot(p, float3(127.1, 311.7, 74.7)),
71               dot(p, float3(269.5, 183.3, 246.1)),
72               dot(p, float3(113.5, 271.9, 124.6)));
73    return -1.0 + 2.0 * frac(sin(p) * 43758.5453123);
74}
75
76float rand3(float3 p)
77{
78    return frac(rand3_hash(p).x);
79}
80
81float2 domainWarp1(float x, uint octaves, float strength, float scale, float speed)
82{
83  [loop]
84  for (uint i = 0; i < octaves; i++) {
85    x += strength * frac(sin(float2(
86      dot(x * scale, float2(12.9898, 78.233)),
87      dot(x * scale + 1, float2(12.9898, 78.233))) * 43758.5453123));
88  }
89  return x;
90}
91
92float2 domainWarp2(float2 uv, uint octaves, float strength, float scale, float speed)
93{
94  uv *= 0.001;
95  [loop]
96  for (uint i = 0; i < octaves; i++) {
97    uv += strength * frac(sin(float2(
98      dot(uv * scale, float2(12.9898, 78.233)),
99      dot(uv * scale, float2(36.7539, 50.3658)))) * 43758.5453123);
100  }
101  uv *= 1000;
102  return uv;
103}
104
105float determinant(float3x3 m)
106{
107  return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
108            - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]))
109            + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
110}
111
112float3x3 inverse(float3x3 m)
113{
114  float det = determinant(m);
115
116  float3x3 adj;
117  adj[0][0] =  (m[1][1] * m[2][2] - m[1][2] * m[2][1]);
118  adj[0][1] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]);
119  adj[0][2] =  (m[0][1] * m[1][2] - m[0][2] * m[1][1]);
120  
121  adj[1][0] = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]);
122  adj[1][1] =  (m[0][0] * m[2][2] - m[0][2] * m[2][0]);
123  adj[1][2] = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]);
124  
125  adj[2][0] =  (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
126  adj[2][1] = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]);
127  adj[2][2] =  (m[0][0] * m[1][1] - m[0][1] * m[1][0]);
128
129  return adj * (1.0 / det);
130}
131
132float3 domainWarp3(float3 pos, uint octaves, float strength, float scale, float offset)
133{
134  [loop]
135  for (uint i = 0; i < octaves; i++) {
136    pos += strength * frac(sin(float3(
137      rand3_dot(pos * scale + offset),
138      rand3_dot(pos * scale + offset + 1),
139      rand3_dot(pos * scale + offset + 2)) * 43758.5453123));
140  }
141  return pos;
142}
143
144void domainWarp3Normals(inout float3 normal, inout float3 tangent, float3 basePos, uint octaves, float strength, float scale, float offset)
145{
146    // Use the actual vertex position for correct derivative evaluation.
147    float3 p = basePos;
148    
149    // Start with the identity matrix for the total Jacobian.
150    float3x3 J = float3x3(
151        1.0, 0.0, 0.0,
152        0.0, 1.0, 0.0,
153        0.0, 0.0, 1.0
154    );
155    
156    const float k = 43758.5453123;
157    // Updated constant vector to match that of rand3_dot (used in domainWarp3)
158    const float3 c = float3(151.0, 157.0, 163.0);
159    
160    for (uint i = 0; i < octaves; i++)
161    {
162        // Compute the vector v using the same offsetting as in domainWarp3.
163        float3 v = float3(
164            dot(p * scale + float3(offset, offset, offset), c),
165            dot(p * scale + float3(offset + 1.0, offset + 1.0, offset + 1.0), c),
166            dot(p * scale + float3(offset + 2.0, offset + 2.0, offset + 2.0), c)
167        );
168        
169        // Compute the warp offset with frac.
170        float3 f_val = frac(sin(v) * k);
171        float3 warpOffset = strength * f_val;
172        
173        // Compute the derivative (Jacobian) of the offset.
174        float3 cos_v = cos(v);
175        float3x3 D = float3x3(
176            strength * k * scale * cos_v.x * c.x, strength * k * scale * cos_v.x * c.y, strength * k * scale * cos_v.x * c.z,
177            strength * k * scale * cos_v.y * c.x, strength * k * scale * cos_v.y * c.y, strength * k * scale * cos_v.y * c.z,
178            strength * k * scale * cos_v.z * c.x, strength * k * scale * cos_v.z * c.y, strength * k * scale * cos_v.z * c.z
179        );
180        
181        // The per–octave Jacobian is I + D.
182        float3x3 iterJacobian = float3x3(
183            1.0 + D[0][0],        D[0][1],        D[0][2],
184                   D[1][0], 1.0 + D[1][1],        D[1][2],
185                   D[2][0],        D[2][1], 1.0 + D[2][2]
186        );
187        
188        // Chain this iteration's Jacobian.
189        J = mul(iterJacobian, J);
190        
191        // Update p for the next iteration.
192        p += warpOffset;
193    }
194    
195    // Transform the normal via the inverse-transpose of the total Jacobian.
196    float3x3 invTransJ = transpose(inverse(J));
197    normal = normalize(mul(invTransJ, normal));
198    
199    // Transform the tangent via the forward total Jacobian.
200    tangent = normalize(mul(J, tangent));
201}
202
203// Alpha blend `dst` onto `src`.
204// Imagine two transparent planes. We're rendering a situation where you're
205// looking through `front` at `behind`.
206float4 alphaBlend(float4 behind, float4 front) {
207  return float4(front.rgb * front.a + behind.rgb * (1 - front.a), front.a + behind.a * (1 - front.a));
208}
209
210// Reoriented normal mapping
211// https://blog.selfshadow.com/publications/blending-in-detail/
212// Inputs are in tangent space.
213float3 blendNormalsHill12(float3 n0, float3 n1) {
214  n0.z += 1.0;
215  n1.xy = -n1.xy;
216  
217  return normalize(n0 * dot(n0, n1) - n1 * n0.z);
218}
219
220float luminance(float3 color) {
221  return dot(color, float3(0.2126, 0.7152, 0.0722));
222}
223
224float median(float3 x) {
225  // Get the min and max.
226  float x_min= min(min(x.r, x.g), x.b);
227  float x_max = max(max(x.r, x.g), x.b);
228  
229  // Compute (x.r + x.g + x.b) - (x_min + x_max). This gives us the median.
230  return (x.r + x.g + x.b) - (x_min + x_max);
231}
232
233// Quaternions
234float4 qmul(float4 q1, float4 q2)
235{
236  return float4(
237      q2.xyz * q1.w + q1.xyz * q2.w + cross(q1.xyz, q2.xyz),
238      q1.w * q2.w - dot(q1.xyz, q2.xyz));
239}
240
241// Vector rotation with a quaternion
242// https://blog.molecular-matters.com/2013/05/24/a-faster-quaternion-vector-multiplication/
243float3 rotate_vector(float3 v, float4 q)
244{
245  float3 t = 2.0 * cross(q.xyz, v);
246  return v + q.w * t + cross(q.xyz, t);
247}
248
249float4 get_quaternion(float3 axis_normal, float theta) {
250  return float4(axis_normal * sin(theta / 2), cos(theta / 2));
251}
252
253void calcNormalInScreenSpace(inout float3 normal, float3 objPos) {
254  normal = normalize(cross(ddy(objPos), ddx(objPos)));
255}
256
257// Formulae from here: https://www.rapidtables.com/convert/color/rgb-to-cmyk.html
258float4 rgbToCmyk(float3 rgb) {
259  float4 cmyk;
260  cmyk[3] = 1 - max(rgb.r, max(rgb.g, rgb.b));
261  cmyk[0] = (1 - rgb.r - cmyk[3]) / (1 - cmyk[3]);
262  cmyk[1] = (1 - rgb.g - cmyk[3]) / (1 - cmyk[3]);
263  cmyk[2] = (1 - rgb.b - cmyk[3]) / (1 - cmyk[3]);
264  return cmyk;
265}
266
267float rgbToCmyk_C(float3 rgb) {
268  float k = 1 - max(rgb.r, max(rgb.g, rgb.b));
269  float c = (1 - rgb.r - k) / (1 - k);
270  return c;
271}
272float rgbToCmyk_M(float3 rgb) {
273  float k = 1 - max(rgb.r, max(rgb.g, rgb.b));
274  float m = (1 - rgb.g - k) / (1 - k);
275  return m;
276}
277float rgbToCmyk_Y(float3 rgb) {
278  float k = 1 - max(rgb.r, max(rgb.g, rgb.b));
279  float y = (1 - rgb.b - k) / (1 - k);
280  return y;
281}
282float rgbToCmyk_K(float3 rgb) {
283  float k = 1 - max(rgb.r, max(rgb.g, rgb.b));
284  return k;
285}
286
287float3 cmykToRgb(float4 cmyk) {
288  return float3(
289    (1 - cmyk[0]) * (1 - cmyk[3]),
290    (1 - cmyk[1]) * (1 - cmyk[3]),
291    (1 - cmyk[2]) * (1 - cmyk[3]));
292}
293
294// Cartesian to cube hexagonal coordinates.
295// Based on this: https://backdrifting.net/post/064_hex_grids
296float3 cart_to_hex(float2 cart) {
297  float p = cart.x;
298  float q = dot(cart, float2(0.5f, SQRT_3_OVER_2));
299  float r = dot(cart, float2(0.5f, -SQRT_3_OVER_2));
300
301  return float3(p, q, r) * TWO_OVER_THREE;
302}
303
304float2 hex_to_cart(float3 cart) {
305  return float2(
306      cart[0] + (cart[1] + cart[2]) * 0.5f,
307      (cart[1] - cart[2]) * SQRT_3_OVER_2);
308}
309
310// Rotate 45 degrees.
311float2 rot45(float2 v) { return float2(v.x - v.y, v.x + v.y) * RCP_SQRT_2; }
312
313// p = point to get noise for
314float valueNoise2D(
315    float2 p) {
316  // quantized part
317  float2 q = floor(p);
318  // fractional part
319  float2 f = frac(p);
320
321  float l00 = rand2(q);
322  float l01 = rand2(q + float2(0, 1));
323  float l10 = rand2(q + float2(1, 0));
324  float l11 = rand2(q + float2(1, 1));
325
326  // Cubic interpolation.
327  f = f * f * (3.0f - 2.0f * f);
328
329  float l0 = lerp(l00, l01, f.y);
330  float l1 = lerp(l10, l11, f.y);
331  return lerp(l0, l1, f.x);
332}
333
334// p = point to get noise for
335float valueNoise3D(
336    float3 p) {
337  // quantized part
338  float3 q = floor(p);
339  // fractional part
340  float3 f = frac(p);
341
342  float l000 = rand3(q);
343  float l001 = rand3(q + float3(0, 0, 1));
344  float l010 = rand3(q + float3(0, 1, 0));
345  float l011 = rand3(q + float3(0, 1, 1));
346  float l100 = rand3(q + float3(1, 0, 0));
347  float l101 = rand3(q + float3(1, 0, 1));
348  float l110 = rand3(q + float3(1, 1, 0));
349  float l111 = rand3(q + float3(1, 1, 1));
350
351  // Cubic interpolation.
352  f = f * f * (3.0f - 2.0f * f);
353
354  float l00 = lerp(l000, l001, f.z);
355  float l01 = lerp(l010, l011, f.z);
356  float l10 = lerp(l100, l101, f.z);
357  float l11 = lerp(l110, l111, f.z);
358  float l0 = lerp(l00, l01, f.y);
359  float l1 = lerp(l10, l11, f.y);
360  return lerp(l0, l1, f.x);
361}
362
363// Fixed version of quilez's `tone` here:
364// https://iquilezles.org/articles/functions/
365float tone(float x, float k) {
366  return (x * (k + 1)) / (k * x + 1);
367}
368
369#endif  // __MATH_INC