cnq3/code/renderer/shaders/crp/common.hlsli

/*
===========================================================================
Copyright (C) 2023-2024 Gian 'myT' Schellenbaum

This file is part of Challenge Quake 3 (CNQ3).

Challenge Quake 3 is free software; you can redistribute it
and/or modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of the License,
or (at your option) any later version.

Challenge Quake 3 is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with Challenge Quake 3. If not, see <https://www.gnu.org/licenses/>.
===========================================================================
*/
// shared utilities


#pragma once


#include "../common/state_bits.h.hlsli"
#include "../common/blend.hlsli"


#define PI    3.1415926535897932384626433832795
#define PI_D2 (PI / 2.0)
#define PI_D4 (PI / 4.0)
#define PI_M2 (PI * 2.0)

#define INT8_MIN  0x80
#define INT16_MIN 0x8000
#define INT32_MIN 0x80000000
#define INT64_MIN 0x8000000000000000

#define INT8_MAX  0x7F
#define INT16_MAX 0x7FFF
#define INT32_MAX 0x7FFFFFFF
#define INT64_MAX 0x7FFFFFFFFFFFFFFF

#define UINT8_MAX  0xFF
#define UINT16_MAX 0xFFFF
#define UINT32_MAX 0xFFFFFFFF
#define UINT64_MAX 0xFFFFFFFFFFFFFFFF

#define FLT_INF asfloat(0x7F800000)
#define FLT_NAN asfloat(0xFFC00000)
#define FLT_MAX asfloat(0x7F7FFFFF)

#define HALF_MAX 65504.0


typedef RaytracingAccelerationStructure RTAS;


float DegToRad(float deg)
{
	return PI * (deg / 180.0);
}

float RadToDeg(float rad)
{
	return 180.0 * (rad / PI);
}

float Brightness(float3 color)
{
	float brightness = dot(color, float3(0.299, 0.587, 0.114));

	return brightness;
}

float3 ColorAtBrightness(float3 color, float targetBrightness)
{
	float brightness = Brightness(color);
	if(brightness <= 0.000001)
	{
		color = float3(0.5, 0.5, 0.5);
		brightness = Brightness(color);
	}
	float brightnessScale = targetBrightness / brightness;
	float3 result = color * brightnessScale;

	return result;
}

float4 MakeGreyscale(float4 input, float amount)
{
	float grey = dot(input.rgb, float3(0.299, 0.587, 0.114));
	float4 result = lerp(input, float4(grey, grey, grey, input.a), amount);

	return result;
}

float LinearDepth(float zwDepth, float zNear, float zFar)
{
	float n = zNear;
	float f = zFar;
	float zw = zwDepth;
	float zl = (f * n) / (n + zw * (f - n));

	return zl;
}

float LinearDepth(float zwDepth, float3 constants)
{
	return constants.x / (constants.y + zwDepth * constants.z);
}

float PostProjectionDepth(float viewDepth, float zNear, float zFar)
{
	float n = zNear;
	float f = zFar;
	float zv = viewDepth;
	float zw = (n * (f - zv)) / ((f - n) * zv);

	return zw;
}

float4 FSTrianglePosFromVertexId(uint id)
{
	return float4(
		(float)(id / 2) * 4.0 - 1.0,
		(float)(id % 2) * 4.0 - 1.0,
		0.0,
		1.0);
}

float2 FSTriangleTCFromVertexId(uint id)
{
	return float2(
		(float)(id / 2) * 2.0,
		1.0 - (float)(id % 2) * 2.0);
}

uint PackColor(float4 c)
{
	uint4 u = uint4(saturate(c) * 255.0);
	uint r = u.r | (u.g << 8) | (u.b << 16) | (u.a << 24);

	return r;
}

float4 UnpackColor(uint c)
{
	uint4 u = uint4(c & 0xFFu, (c >> 8) & 0xFFu, (c >> 16) & 0xFFu, (c >> 24) & 0xFFu);
	float4 r = float4(u) / 255.0;

	return r;
}

float EaseInCubic(float x)
{
	return x * x * x;
}

float EaseOutCubic(float x)
{
	float y = 1.0 - x;

	return 1.0 - y * y * y;
}

float EaseInOutCubic(float x)
{
	if(x < 0.5)
	{
		return 4 * x * x * x;
	}

	float y = -2 * x + 2;

	return 1 - 0.5 * y * y * y;
}

float EaseInQuad(float x)
{
	return x * x;
}

float EaseOutQuad(float x)
{
	float y = 1.0 - x;

	return 1.0 - y * y;
}

float EaseInExp(float x)
{
	return x == 0.0 ? 0.0 : pow(2.0, 10.0 * x - 10.0);
}

float EaseOutExp(float x)
{
	return x == 1.0 ? 1.0 : 1.0 - pow(2.0, -10.0 * x);
}

float smoothstep01(float x)
{
	return smoothstep(0.0, 1.0, x);
}

// Oct*: octahedron normal vector encoding
// original code from "A Survey of Efficient Representations for Independent Unit Vectors"
// further improved by Krzysztof Narkowicz and Rune Stubbe

float2 OctWrap(float2 v)
{
	return (1.0 - abs(v.yx)) * select(v.xy >= 0.0, 1.0, -1.0);
}

float2 OctEncode(float3 n)
{
	n /= (abs(n.x) + abs(n.y) + abs(n.z));
	n.xy = n.z >= 0.0 ? n.xy : OctWrap(n.xy);
	n.xy = n.xy * 0.5 + 0.5;

	return n.xy;
}

float3 OctDecode(float2 f)
{
	f = f * 2.0 - 1.0;
	float3 n = float3(f.x, f.y, 1.0 - abs(f.x) - abs(f.y));
	float t = saturate(-n.z);
	n.xy += select(n.xy >= 0.0, -t, t);

	return normalize(n);
}

float3 GetPositionFromDepth(float2 tc01, float depthZW, matrix invMatrix)
{
	float x = tc01.x * 2.0 - 1.0;
	float y = (1.0 - tc01.y) * 2.0 - 1.0;
	float4 position = mul(float4(x, y, depthZW, 1.0), invMatrix);
	float3 result = position.xyz / position.w;

	return result;
}

float3 TransformNormal(float3 normal, matrix transform)
{
	return mul(transform, float4(normal, 0)).xyz;
}

float3 TransformPoint(float3 position, matrix transform)
{
	float4 result = mul(transform, float4(position, 1));
	return result.xyz / result.w;
}

float3 RandomColorFromUInt(uint id)
{
	float r = frac(0.420 + 1.337 * id);
	float g = frac(0.69 + 1.666 * id);
	float b = frac(0.13 + 1.777 * id);

	return float3(r, g, b);
}

float3 BiasPosition(float3 position, float3 normal)
{
	float3 result = position + sign(normal) * abs(position * 0.0000002);

	return result;
}

// from Mauricio Vives, https://gist.github.com/pixnblox/5e64b0724c186313bc7b6ce096b08820
// Projects the specified position (point) onto the plane with the specified origin and normal.
float3 ProjectPointOnPlane(float3 position, float3 planeOrigin, float3 planeNormal)
{
	return position - dot(position - planeOrigin, planeNormal) * planeNormal;
}

// from Mauricio Vives, https://gist.github.com/pixnblox/5e64b0724c186313bc7b6ce096b08820
// Computes the shading position of the specified geometric position and vertex positions and
// normals. For a triangle with normals describing a convex surface, this point will be slightly
// above the surface. For a concave surface, the geometry position is used directly.
// NOTE: The difference between the shading position and geometry position is significant when
// casting shadow rays. If the geometric position is used, a triangle may fully shadow itself when
// it should be partly lit based on the shading normals; this is the "shadow terminator" problem.
float3 GetShadingPosition(
	float3 geomPosition, float3 shadingNormal,
	float3 positions[3], float3 normals[3], float3 barycentrics)
{
	// Project the geometric position (inside the triangle) to the planes defined by the vertex
	// positions and normals.
	float3 p0 = ProjectPointOnPlane(geomPosition, positions[0], normals[0]);
	float3 p1 = ProjectPointOnPlane(geomPosition, positions[1], normals[1]);
	float3 p2 = ProjectPointOnPlane(geomPosition, positions[2], normals[2]);

	// Interpolate the projected positions using the barycentric coordinates, which gives the
	// shading position.
	float3 shadingPosition = p0 * barycentrics.x + p1 * barycentrics.y + p2 * barycentrics.z;

	// Return the shading position for a convex triangle, where the shading point is above the
	// triangle based on the shading normal. Otherwise use the geometric position.
	bool convex = dot(shadingPosition - geomPosition, shadingNormal) > 0.0;
	float3 result = convex ? shadingPosition : BiasPosition(geomPosition, shadingNormal);

	return result;
}

// based on "Hacking the Shadow Terminator" by Johannes Hanika in "Ray Tracing Gems II"
float3 GetShadingPositionV2(float3 geomPosition, float3 positions[3], float3 normals[3], float3 barycentrics)
{
	float3 tmpu = geomPosition - positions[0];
	float3 tmpv = geomPosition - positions[1];
	float3 tmpw = geomPosition - positions[2];
	float dotu = min(0.0, dot(tmpu, normals[0]));
	float dotv = min(0.0, dot(tmpv, normals[1]));
	float dotw = min(0.0, dot(tmpw, normals[2]));
	tmpu -= dotu * normals[0];
	tmpv -= dotv * normals[1];
	tmpw -= dotw * normals[2];
	float3 shadingPosition = geomPosition + 1.0 * (barycentrics.x * tmpu + barycentrics.y * tmpv + barycentrics.z * tmpw);

	return shadingPosition;
}

template<typename T>
T trilerp(T v0, T v1, T v2, float3 barycentrics)
{
	return
		barycentrics.x * v0 +
		barycentrics.y * v1 +
		barycentrics.z * v2;
}

template<>
float trilerp(float v0, float v1, float v2, float3 barycentrics)
{
	return dot(float3(v0, v1, v2), barycentrics);
}

// Interleaved Gradient Noise by Jorge Jimenez
// from "Next Generation Post Processing in Call of Duty: Advanced Warfare"
float InterleavedGradientNoise(float2 uv)
{
	float3 magic = float3(0.06711056, 0.00583715, 52.9829189);
	return frac(magic.z * frac(dot(uv, magic.xy)));
}

template<typename T>
bool IsInRange(T p, T min, T max)
{
	return all(p >= min) && all(p <= max);
}

bool Is01(float2 p)
{
	return IsInRange(p, float2(0, 0), float2(1, 1));
}

bool Is01(float3 p)
{
	return IsInRange(p, float3(0, 0, 0), float3(1, 1, 1));
}

bool Is01(float4 p)
{
	return IsInRange(p, float4(0, 0, 0, 0), float4(1, 1, 1, 1));
}

bool IsInTexture(int2 tc, int2 textureSize)
{
	return all(tc >= int2(0, 0)) && all(tc < textureSize);
}

bool IsInTexture(int3 tc, int3 textureSize)
{
	return all(tc >= int3(0, 0, 0)) && all(tc < textureSize);
}

template<typename T>
uint2 GetTextureSize(Texture2D<T> texture0)
{
	uint2 size;
	texture0.GetDimensions(size.x, size.y);

	return size;
}

template<typename T>
uint2 GetTextureSize(RWTexture2D<T> texture0)
{
	uint2 size;
	texture0.GetDimensions(size.x, size.y);

	return size;
}

template<typename T>
uint3 GetTextureSize(Texture3D<T> texture0)
{
	uint3 size;
	texture0.GetDimensions(size.x, size.y, size.z);

	return size;
}

template<typename T>
uint3 GetTextureSize(RWTexture3D<T> texture0)
{
	uint3 size;
	texture0.GetDimensions(size.x, size.y, size.z);

	return size;
}

// by Sakib Saikia, https://sakibsaikia.github.io/graphics/2022/01/04/Nan-Checks-In-HLSL.html
bool IsNan(float x)
{
	return (asuint(x) & 0x7FFFFFFFu) > 0x7F800000u;
}

bool isnan(float x)
{
	return IsNan(x);
}

// from "Using Blue Noise For Raytraced Soft Shadows" by Alan Wolfe in "Ray Tracing Gems II"
// this turns the blue noise into a low discrepancy additive recurrence
float AnimateBlueNoise(float blueNoise, uint frameIndex)
{
	return frac(blueNoise + float(frameIndex % 32) * 0.61803399);
}

// credit: David Hoskins
float Hash1To1(float p)
{
	p = frac(p * .1031);
	p *= p + 33.33;
	p *= p + p;
	return frac(p);
}

// credit: David Hoskins
float Hash2To1(float2 p)
{
	float3 p3 = frac(float3(p.xyx) * .1031);
	p3 += dot(p3, p3.yzx + 33.33);
	return frac((p3.x + p3.y) * p3.z);
}

// credit: David Hoskins
float Hash3To1(float3 p3)
{
	p3 = frac(p3 * .1031);
	p3 += dot(p3, p3.yzx + 33.33);
	return frac((p3.x + p3.y) * p3.z);
}

// credit: David Hoskins
float2 Hash1To2(float p)
{
	float3 p3 = frac(float3(p, p, p) * float3(.1031, .1030, .0973));
	p3 += dot(p3, p3.yzx + 33.33);
	return frac((p3.xx + p3.yz) * p3.zy);
}

// credit: David Hoskins
float2 Hash2To2(float2 p)
{
	float3 p3 = frac(float3(p.xyx) * float3(.1031, .1030, .0973));
	p3 += dot(p3, p3.yzx + 33.33);
	return frac((p3.xx + p3.yz) * p3.zy);
}

// credit: David Hoskins
float2 Hash3To2(float3 p3)
{
	p3 = frac(p3 * float3(.1031, .1030, .0973));
	p3 += dot(p3, p3.yzx + 33.33);
	return frac((p3.xx + p3.yz) * p3.zy);
}

// credit: David Hoskins
float3 Hash1To3(float p)
{
	float3 p3 = frac(float3(p, p, p) * float3(.1031, .1030, .0973));
	p3 += dot(p3, p3.yzx + 33.33);
	return frac((p3.xxy + p3.yzz) * p3.zyx);
}

// credit: David Hoskins
float3 Hash2To3(float2 p)
{
	float3 p3 = frac(float3(p.xyx) * float3(.1031, .1030, .0973));
	p3 += dot(p3, p3.yxz + 33.33);
	return frac((p3.xxy + p3.yzz) * p3.zyx);
}

// credit: David Hoskins
float3 Hash3To3(float3 p3)
{
	p3 = frac(p3 * float3(.1031, .1030, .0973));
	p3 += dot(p3, p3.yxz + 33.33);
	return frac((p3.xxy + p3.yxx) * p3.zyx);
}

// credit: David Hoskins
float4 Hash1To4(float p)
{
	float4 p4 = frac(float4(p, p, p, p) * float4(.1031, .1030, .0973, .1099));
	p4 += dot(p4, p4.wzxy + 33.33);
	return frac((p4.xxyz + p4.yzzw) * p4.zywx);
}

// credit: David Hoskins
float4 Hash2To4(float2 p)
{
	float4 p4 = frac(float4(p.xyxy) * float4(.1031, .1030, .0973, .1099));
	p4 += dot(p4, p4.wzxy + 33.33);
	return frac((p4.xxyz + p4.yzzw) * p4.zywx);
}

// credit: David Hoskins
float4 Hash3To4(float3 p)
{
	float4 p4 = frac(float4(p.xyzx) * float4(.1031, .1030, .0973, .1099));
	p4 += dot(p4, p4.wzxy + 33.33);
	return frac((p4.xxyz + p4.yzzw) * p4.zywx);
}

// credit: David Hoskins
float4 Hash4To4(float4 p4)
{
	p4 = frac(p4 * float4(.1031, .1030, .0973, .1099));
	p4 += dot(p4, p4.wzxy + 33.33);
	return frac((p4.xxyz + p4.yzzw) * p4.zywx);
}

float2 NDCToTC(float2 ndc)
{
	float2 tc = ndc * float2(0.5, -0.5) + float2(0.5, 0.5);

	return tc;
}

float3 NDCToTC(float3 ndc)
{
	float3 tc = ndc * float3(0.5, -0.5, 0.5) + float3(0.5, 0.5, 0.5);

	return tc;
}

float2 TCToNDC(float2 tc)
{
	float2 ndc = (2.0 * tc - 1.0) * float2(1, -1);

	return ndc;
}

float3 TCToNDC(float3 tc)
{
	float3 ndc = (2.0 * tc - 1.0) * float3(1, -1, 1);

	return ndc;
}

// returns the longest vector
float2 vmax(float2 a, float2 b)
{
	float2 result = dot(a, a) > dot(b, b) ? a : b;

	return result;
}

uint PackHalf2(float2 input)
{
	uint2 d = f32tof16(input);
	uint result = d.x | (d.y << 16u);

	return result;
}

float2 UnpackHalf2(uint input)
{
	uint2 d = uint2(input & 0xFFFFu, input >> 16u);
	float2 result = f16tof32(d);

	return result;
}

float2 CartesianToPolar(float2 cartesian)
{
	float radius = length(cartesian);
	float angle = atan2(cartesian.y, cartesian.x);
	float2 polar = float2(radius, angle);

	return polar;
}

float2 PolarToCartesian(float2 polar)
{
	float sinAngle, cosAngle;
	sincos(polar.y, sinAngle, cosAngle);
	float2 cartesian = polar.x * float2(cosAngle, sinAngle);

	return cartesian;
}

// Beer-Lambert law
float Transmittance(float distance, float extinction)
{
	float transmittance = exp(-distance * extinction);

	return transmittance;
}

// phase function for Mie scattering
// g is in the range [-1;1]
// -1: backward scattering
//  0: isotropic
//  1: forward scattering
float HenyeyGreenstein(float cosTheta, float g)
{
	float g2 = g * g;
	float num = 1.0 - g2;
	float denom = 4.0 * PI * pow(1.0 + g2 - 2.0 * g * cosTheta, 1.5);
	float result = num / denom;

	return result;
}

uint FlattenIndex(uint3 tileIndex, uint3 tileResolution)
{
	return
		tileIndex.x +
		tileIndex.y * tileResolution.x +
		tileIndex.z * tileResolution.x * tileResolution.y;
}

int FlattenIndex(int3 tileIndex, int3 tileResolution)
{
	return
		tileIndex.x +
		tileIndex.y * tileResolution.x +
		tileIndex.z * tileResolution.x * tileResolution.y;
}

uint3 UnflattenIndex(uint flatIndex, uint3 tileResolution)
{
	uint w = tileResolution.x;
	uint h = tileResolution.y;
	uint wh = w * h;
	uint z = flatIndex / wh;
	flatIndex -= z * wh;
	uint y = flatIndex / w;
	uint x = flatIndex - y * w;
	uint3 result = uint3(x, y, z);

	return result;
}

int3 UnflattenIndex(int flatIndex, int3 tileResolution)
{
	int w = tileResolution.x;
	int h = tileResolution.y;
	int wh = w * h;
	int z = flatIndex / wh;
	flatIndex -= z * wh;
	int y = flatIndex / w;
	int x = flatIndex - y * w;
	int3 result = int3(x, y, z);

	return result;
}

void ClearBoundingBox(out int3 boxMin, out int3 boxMax)
{
	boxMin = int3(INT32_MAX, INT32_MAX, INT32_MAX);
	boxMax = int3(INT32_MIN, INT32_MIN, INT32_MIN);
}

void ClearBoundingBox(out float3 boxMin, out float3 boxMax)
{
	boxMin = float3(FLT_MAX, FLT_MAX, FLT_MAX);
	boxMax = float3(-FLT_MAX, -FLT_MAX, -FLT_MAX);
}

template<typename T>
void ExpandBoundingBox(inout T boxMin, inout T boxMax, T newPoint)
{
	boxMin = min(boxMin, newPoint);
	boxMax = max(boxMax, newPoint);
}

// Credit: Riku Salminen
// dispatch the draw call with 36 indices
float3 CubeFromVertexID(uint vertexId)
{
	int tri = int(vertexId) / 3;
	int idx = int(vertexId) % 3;
	int face = tri / 2;
	int top = tri % 2;

	int dir = face % 3;
	int pos = face / 3;

	int nz = dir >> 1;
	int ny = dir & 1;
	int nx = 1 ^ (ny | nz);

	float3 d = float3(nx, ny, nz);
	float flip = 1 - 2 * pos;

	float3 n = flip * d;
	float3 u = -d.yzx;
	float3 v = flip * d.zxy;

	float mirror = -1 + 2 * top;
	float3 xyz = n + mirror * (1 - 2 * (idx & 1)) * u + mirror * (1 - 2 * (idx >> 1)) * v;

	return xyz;
}

// dispatch the draw call with 6 indices
float2 QuadFromVertexID(uint vertexId)
{
	float2 position;
	position.x = (vertexId >= 1 && vertexId <= 3) ? 1.0 : -1.0;
	position.y = (vertexId >= 2 && vertexId <= 4) ? -1.0 : 1.0;

	return position;
}

float3 AABoxIndexToWorldSpace(int3 index, float3 centerPosition, float3 textureSize, float3 worldScale)
{
	float3 position = centerPosition + worldScale * (float3(index)+float3(0.5, 0.5, 0.5) - 0.5 * textureSize);

	return position;
}

float3 AABoxTCToWorldSpace(float3 tc, float3 centerPosition, float3 textureSize, float3 worldScale)
{
	float3 position = centerPosition + worldScale * textureSize * (tc - float3(0.5, 0.5, 0.5));

	return position;
}

float3 AABoxWorldSpaceToTC(float3 position, float3 centerPosition, float3 textureSize, float3 worldScale)
{
	float3 boxSize = worldScale * textureSize;
	float3 boxMin = centerPosition - 0.5 * boxSize;
	float3 tc = (position - boxMin) / boxSize;

	return tc;
}

int3 AABoxWorldSpaceToIndex(float3 position, float3 centerPosition, float3 textureSize, float3 worldScale)
{
	float3 boxSize = worldScale * textureSize;
	float3 boxMin = centerPosition - 0.5 * boxSize;
	float3 indexF = (position - boxMin) / worldScale;
	int3 index = int3(indexF);

	return index;
}

float3 AABoxIndexToWorldSpace(int3 index, float3 centerPosition, float3 textureSize, float worldScale)
{
	return AABoxIndexToWorldSpace(index, centerPosition, textureSize, worldScale.xxx);
}

float3 AABoxTCToWorldSpace(float3 tc, float3 centerPosition, float3 textureSize, float worldScale)
{
	return AABoxTCToWorldSpace(tc, centerPosition, textureSize, worldScale.xxx);
}

float3 AABoxWorldSpaceToTC(float3 position, float3 centerPosition, float3 textureSize, float worldScale)
{
	return AABoxWorldSpaceToTC(position, centerPosition, textureSize, worldScale.xxx);
}

int3 AABoxWorldSpaceToIndex(float3 position, float3 centerPosition, float3 textureSize, float worldScale)
{
	return AABoxWorldSpaceToIndex(position, centerPosition, textureSize, worldScale.xxx);
}

template<typename T>
T min3(T v0, T v1, T v2)
{
	return min(v0, min(v1, v2));
}

template<typename T>
T max3(T v0, T v1, T v2)
{
	return max(v0, max(v1, v2));
}

template<typename T>
T min4(T v0, T v1, T v2, T v3)
{
	return min(min(v0, v1), min(v2, v3));
}

template<typename T>
T max4(T v0, T v1, T v2, T v3)
{
	return max(max(v0, v1), max(v2, v3));
}

// credit: Inigo Quilez
// returns t == -1 when nothing was hit
float RaytraceSphere(float3 rayOrigin, float3 rayDir, float3 spherePos, float sphereRadius)
{
	float3 oc = rayOrigin - spherePos;
	float b = dot(oc, rayDir);
	float c = dot(oc, oc) - sphereRadius * sphereRadius;
	float h = b * b - c;
	if(h < 0.0)
	{
		return -1.0;
	}
	h = sqrt(h);
	float t = -b - h;

	return t;
}

float3 DirectionFromLongLat(float longitude01, float latitude01)
{
	float lon = longitude01 * PI_M2;
	float lat = latitude01 * PI_M2;
	float sinLat, cosLat;
	sincos(lat, sinLat, cosLat);
	float sinLon, cosLon;
	sincos(lon, sinLon, cosLon);
	float3 direction = float3(cosLat * sinLon, sinLat * sinLon, cosLon);

	return direction;
}

float SphereVolume(float radius)
{
	float volume = ((4.0 / 3.0) * PI) * radius * radius * radius;

	return volume;
}

// "2D Polyhedral Bounds of a Clipped, Perspective-Projected 3D Sphere" by Mara and McGuire
float2 ProjectedSphereExtentsAxis(float xy, float z, float r)
{
	float t = sqrt(xy * xy + z * z - r * r);
	float min = (t * xy - r * z) / (t * z + r * xy);
	float max = (t * xy + r * z) / (t * z - r * xy);
	float2 result = float2(min, max);

	return result;
}

// "2D Polyhedral Bounds of a Clipped, Perspective-Projected 3D Sphere" by Mara and McGuire
float4 ProjectedSphereExtentsNDC(float3 spherePositionWS, float sphereRadius, matrix viewMatrix, matrix projMatrix)
{
	float4 spherePosVSw = mul(viewMatrix, float4(spherePositionWS, 1));
	float3 spherePosVS = spherePosVSw.xyz / spherePosVSw.w;
	float2 extentsX = ProjectedSphereExtentsAxis(spherePosVS.x, -spherePosVS.z, sphereRadius) * projMatrix[0][0] + projMatrix[2][0];
	float2 extentsY = ProjectedSphereExtentsAxis(spherePosVS.y,  spherePosVS.z, sphereRadius) * projMatrix[1][1] + projMatrix[2][1];
	float4 result = float4(extentsX.x, extentsY.x, extentsX.y, extentsY.y);

	return result;
}

float3 Project(float3 P, matrix m)
{
	float4 Qw = mul(m, float4(P, 1));
	float3 Q = Qw.xyz / Qw.w;

	return Q;
}

float VoxelStepSize(float3 dir, float3 voxelSize)
{
	float3 stepSize3 = voxelSize / max(abs(dir), float(0.000001).xxx);
	float stepSize = min3(stepSize3.x, stepSize3.y, stepSize3.z);

	return stepSize;
}

float2x2 RandomRotationMatrix2x2(float3 position)
{
	float angle = Hash3To1(position) * 2.0 * PI;
	float sin, cos;
	sincos(angle, sin, cos);
	float2x2 result = float2x2(cos, -sin, sin, cos);

	return result;
}

float3x3 RandomRotationMatrix3x3(float3 position)
{
	float3 angles = Hash3To3(position) * 2.0 * PI;
	float3 sin, cos;
	sincos(angles.x, sin.x, cos.x);
	sincos(angles.y, sin.y, cos.y);
	sincos(angles.z, sin.z, cos.z);
	float r0 = cos.x * cos.y;
	float r1 = cos.x * sin.y * sin.z - sin.x * cos.z;
	float r2 = cos.x * sin.y * cos.z + sin.x * sin.z;
	float r3 = sin.x * cos.y;
	float r4 = sin.x * sin.y * sin.z + cos.x * cos.z;
	float r5 = sin.x * sin.y * cos.z - cos.x * sin.z;
	float r6 = -sin.y;
	float r7 = cos.y * sin.z;
	float r8 = cos.y * cos.z;
	float3x3 result = float3x3(r0, r1, r2, r3, r4, r5, r6, r7, r8);

	return result;
}

float2x2 IdentityMatrix2x2()
{
	return float2x2(
		1, 0,
		0, 1);
}

float3x3 IdentityMatrix3x3()
{
	return float3x3(
		1, 0, 0,
		0, 1, 0,
		0, 0, 1);
}

matrix IdentityMatrix4x4()
{
	return matrix(
		1, 0, 0, 0,
		0, 1, 0, 0,
		0, 0, 1, 0,
		0, 0, 0, 1);
}