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