1
2/// π
3const PI: f32 = 3.141592653589793;
4/// 2π
5const PI_2: f32 = 6.283185307179586;
6/// π/2
7const HALF_PI: f32 = 1.57079632679;
8/// π/3
9const FRAC_PI_3: f32 = 1.0471975512;
10/// exp(1)
11const E: f32 = 2.718281828459045;
12
13fn affine2_to_square(affine: mat3x2<f32>) -> mat3x3<f32> {
14    return mat3x3<f32>(
15        vec3<f32>(affine[0].xy, 0.0),
16        vec3<f32>(affine[1].xy, 0.0),
17        vec3<f32>(affine[2].xy, 1.0),
18    );
19}
20
21fn affine3_to_square(affine: mat3x4<f32>) -> mat4x4<f32> {
22    return transpose(mat4x4<f32>(
23        affine[0],
24        affine[1],
25        affine[2],
26        vec4<f32>(0.0, 0.0, 0.0, 1.0),
27    ));
28}
29
30fn mat2x4_f32_to_mat3x3_unpack(
31    a: mat2x4<f32>,
32    b: f32,
33) -> mat3x3<f32> {
34    return mat3x3<f32>(
35        a[0].xyz,
36        vec3<f32>(a[0].w, a[1].xy),
37        vec3<f32>(a[1].zw, b),
38    );
39}
40
41/// Extracts the square portion of an affine matrix: i.e. discards the
42/// translation.
43fn affine3_to_mat3x3(affine: mat4x3<f32>) -> mat3x3<f32> {
44    return mat3x3<f32>(affine[0].xyz, affine[1].xyz, affine[2].xyz);
45}
46
47/// Returns the inverse of a 3x3 matrix.
48fn inverse_mat3x3(matrix: mat3x3<f32>) -> mat3x3<f32> {
49    let tmp0 = cross(matrix[1], matrix[2]);
50    let tmp1 = cross(matrix[2], matrix[0]);
51    let tmp2 = cross(matrix[0], matrix[1]);
52    let inv_det = 1.0 / dot(matrix[2], tmp2);
53    return transpose(mat3x3<f32>(tmp0 * inv_det, tmp1 * inv_det, tmp2 * inv_det));
54}
55
56/// Returns the inverse of an affine matrix.
57///
58/// https://en.wikipedia.org/wiki/Affine_transformation#Groups
59fn inverse_affine3(affine: mat4x3<f32>) -> mat4x3<f32> {
60    let matrix3 = affine3_to_mat3x3(affine);
61    let inv_matrix3 = inverse_mat3x3(matrix3);
62    return mat4x3<f32>(inv_matrix3[0], inv_matrix3[1], inv_matrix3[2], -(inv_matrix3 * affine[3]));
63}
64
65/// Extracts the upper 3x3 portion of a 4x4 matrix.
66fn mat4x4_to_mat3x3(m: mat4x4<f32>) -> mat3x3<f32> {
67    return mat3x3<f32>(m[0].xyz, m[1].xyz, m[2].xyz);
68}
69
70/// Copy the sign bit from B onto A.
71/// copysign allows proper handling of negative zero to match the rust implementation of orthonormalize
72fn copysign(a: f32, b: f32) -> f32 {
73    return bitcast<f32>((bitcast<u32>(a) & 0x7FFFFFFF) | (bitcast<u32>(b) & 0x80000000));
74}
75
76/// Constructs a right-handed orthonormal basis from a given unit Z vector.
77///
78/// NOTE: requires unit-length (normalized) input to function properly.
79///
80/// https://jcgt.org/published/0006/01/01/paper.pdf
81/// this method of constructing a basis from a vec3 is also used by `glam::Vec3::any_orthonormal_pair`
82/// the construction of the orthonormal basis up and right vectors here needs to precisely match the rust
83/// implementation in bevy_light/spot_light.rs:spot_light_world_from_view
84fn orthonormalize(z_basis: vec3<f32>) -> mat3x3<f32> {
85    let sign = copysign(1.0, z_basis.z);
86    let a = -1.0 / (sign + z_basis.z);
87    let b = z_basis.x * z_basis.y * a;
88    let x_basis = vec3(1.0 + sign * z_basis.x * z_basis.x * a, sign * b, -sign * z_basis.x);
89    let y_basis = vec3(b, sign + z_basis.y * z_basis.y * a, -z_basis.y);
90    return mat3x3(x_basis, y_basis, z_basis);
91}
92
93/// Rotates a direction by a rotation quaternion q.
94fn quat_rotate(q: vec4<f32>, dir: vec3<f32>) -> vec3<f32> {
95    let t = 2.0 * cross(q.xyz, dir);
96    return dir + q.w * t + cross(q.xyz, t);
97}
98
99/// Returns true if any part of a sphere is on the positive side of a plane.
100///
101/// `sphere_center.w` should be 1.0.
102///
103/// This is used for frustum culling.
104fn sphere_intersects_plane_half_space(
105    plane: vec4<f32>,
106    sphere_center: vec4<f32>,
107    sphere_radius: f32
108) -> bool {
109    return dot(plane, sphere_center) + sphere_radius > 0.0;
110}
111
112/// Returns the distances along the ray to its intersections with a sphere
113/// centered at the origin.
114///
115/// r: distance from the sphere center to the ray origin
116/// mu: cosine of the zenith angle
117/// sphere_radius: radius of the sphere
118///
119/// Returns vec2(t0, t1). If there is no intersection, returns vec2(-1.0).
120fn ray_sphere_intersect(r: f32, mu: f32, sphere_radius: f32) -> vec2<f32> {
121    let discriminant = r * r * (mu * mu - 1.0) + sphere_radius * sphere_radius;
122
123    // No intersection
124    if discriminant < 0.0 {
125        return vec2(-1.0);
126    }
127
128    let q = -r * mu;
129    let sqrt_discriminant = sqrt(discriminant);
130
131    // Return both intersection distances
132    return vec2(
133        q - sqrt_discriminant,
134        q + sqrt_discriminant
135    );
136}
137
138/// pow() but safe for NaNs/negatives
139fn powsafe(color: vec3<f32>, power: f32) -> vec3<f32> {
140    return pow(abs(color), vec3(power)) * sign(color);
141}
142
143/// https://en.wikipedia.org/wiki/Vector_projection#Vector_projection_2
144fn project_onto(lhs: vec3<f32>, rhs: vec3<f32>) -> vec3<f32> {
145    let other_len_sq_rcp = 1.0 / dot(rhs, rhs);
146    return rhs * dot(lhs, rhs) * other_len_sq_rcp;
147}
148
149// Below are fast approximations of common irrational and trig functions. These
150// are likely most useful when raymarching, for example, where complete numeric
151// accuracy can be sacrificed for greater sample count.
152
153/// Slightly less accurate than fast_acos_4, but much simpler.
154fn fast_acos(in_x: f32) -> f32 {
155    let x = abs(in_x);
156    var res = -0.156583 * x + HALF_PI;
157    res *= sqrt(1.0 - x);
158    return select(PI - res, res, in_x >= 0.0);
159}
160
161/// 4th order polynomial approximation
162/// 4 VGRP, 16 ALU Full Rate
163/// 7 * 10^-5 radians precision
164/// Reference : Handbook of Mathematical Functions (chapter : Elementary Transcendental Functions), M. Abramowitz and I.A. Stegun, Ed.
165fn fast_acos_4(x: f32) -> f32 {
166    let x1 = abs(x);
167    let x2 = x1 * x1;
168    let x3 = x2 * x1;
169    var s: f32;
170
171    s = -0.2121144 * x1 + 1.5707288;
172    s = 0.0742610 * x2 + s;
173    s = -0.0187293 * x3 + s;
174    s = sqrt(1.0 - x1) * s;
175
176	// acos function mirroring
177    return select(PI - s, s, x >= 0.0);
178}
179
180fn fast_atan2(y: f32, x: f32) -> f32 {
181    var t0 = max(abs(x), abs(y));
182    var t1 = min(abs(x), abs(y));
183    var t3 = t1 / t0;
184    var t4 = t3 * t3;
185
186    t0 = 0.0872929;
187    t0 = t0 * t4 - 0.301895;
188    t0 = t0 * t4 + 1.0;
189    t3 = t0 * t3;
190
191    t3 = select(t3, (0.5 * PI) - t3, abs(y) > abs(x));
192    t3 = select(t3, PI - t3, x < 0);
193    t3 = select(-t3, t3, y > 0);
194
195    return t3;
196}