GPURayTracer/shaders/ray_tracing/RayTracing.compute

777 lines
24 KiB
Plaintext
Raw Normal View History

2021-08-07 23:14:51 +02:00
#version 430 core
#define PI 3.1415926538f
layout(local_size_x = 8, local_size_y = 8) in;
layout(rgba32f) uniform image2D img_output;
uniform sampler2D u_skybox;
uniform sampler2D u_tex1;
uniform sampler2D u_tex2;
uniform sampler2D u_tex3;
uniform sampler2D u_tex4;
uniform sampler2D u_tex5;
precision highp float;
struct SceneObject{
vec4 position;
vec4 size;
vec4 orientation;
vec4 albedo;
vec4 specular;
vec4 emission;
ivec4 texture_id;
float refractive_index;
float transparency;
float smoothness;
int type;
int offset;
int vert_num;
int obj_id;
int padding;
};
struct Ray{
vec3 origin;
vec3 direction;
vec3 color;
};
layout(std140, binding=2) uniform shader_data{
vec4 _camera_position; // 4 4
vec4 _camera_rotation; // 4 8
vec2 _pixel_offset; // 2 10
ivec2 _screen_size; // 2 12
int _iterations; // 1 13
float _seed; // 1 14
int _object_count; // 1 15
int _sample; // 1 16
int _samples_per_frame; // 1 17
int _fov; // 1 18
vec2 padding; // 2 20
};
layout (std140, binding = 3) uniform object_data{
SceneObject objects[1];
};
layout (std430, binding = 4) buffer index_buffer{
int faces[];
};
layout (std430, binding = 5) buffer vert_buffer{
float vertices[];
};
layout (std430, binding = 6) buffer texture_vert_buffer{
float t_vertices[];
};
uniform float u_gamma;
float distlimit = 10000.0f;
float seed = _seed;
vec2 pixel_id;
vec3 DegToRad(vec3 vect){
return vect * float(PI)/180.0f;
}
float saturate(float value){
return clamp(value,0.0,1.0);
}
float sdot(vec3 x, vec3 y){
float f = 1.0f;
return saturate(dot(x, y) * f);
}
float sdot(vec3 x, vec3 y, float f){
return saturate(dot(x, y) * f);
}
// A single iteration of Bob Jenkins' One-At-A-Time hashing algorithm.
uint hash( uint x ) {
x += ( x << 10u );
x ^= ( x >> 6u );
x += ( x << 3u );
x ^= ( x >> 11u );
x += ( x << 15u );
return x;
}
uint hash( uvec3 v ) { return hash( v.x ^ hash(v.y) ^ hash(v.z) ); }
// Construct a float with half-open range [0:1] using low 23 bits.
// All zeroes yields 0.0, all ones yields the next smallest representable value below 1.0.
float floatConstruct( uint m ) {
const uint ieeeMantissa = 0x007FFFFFu; // binary32 mantissa bitmask
const uint ieeeOne = 0x3F800000u; // 1.0 in IEEE binary32
m &= ieeeMantissa; // Keep only mantissa bits (fractional part)
m |= ieeeOne; // Add fractional part to 1.0
float f = uintBitsToFloat( m ); // Range [1:2]
return f - 1.0; // Range [0:1]
}
// Pseudo-random value in half-open range [0:1].
float Rand() { seed += 1.0f; return floatConstruct(hash(floatBitsToUint(vec3(pixel_id, seed)))); }
float SmoothnessToPhongAlpha(float s){
return pow(1000.0f, s * s);
}
mat3x3 GetTangentSpace(vec3 n){
// Choose a helper vector for the cross product
vec3 helper = vec3(1, 0, 0);
if (abs(n.x) > 0.99f)
helper = vec3(0, 0, 1);
// Generate vectors
vec3 t = normalize(cross(n, helper));
vec3 b = normalize(cross(n, t));
return mat3x3(
t.x, b.x, n.x,
t.y, b.y, n.y,
t.z, b.z, n.z
);
}
vec3 SampleHemisphere(vec3 normal, float alpha){
// Sample the hemisphere, where alpha determines the kind of the sampling
float cosTheta = pow(Rand(), 1.0f / (alpha + 1.0f));
float sinTheta = sqrt(1.0f - cosTheta * cosTheta);
float phi = 2 * PI * Rand();
vec3 tangentSpaceDir = vec3(cos(phi) * sinTheta, sin(phi) * sinTheta, cosTheta);
// Transform direction to world space
return normalize(tangentSpaceDir * GetTangentSpace(normal));
}
vec3 RefractRay(vec3 direction, vec3 n, float n1, float n2){
vec3 new_direction;
direction = normalize(direction);
n = normalize(n);
float dotp = dot(n, -direction);
float ratio = n1/n2;
vec3 e1 = (ratio * direction); // n1/n2 * i
float e2 = (ratio * dotp); // n1/n2 * cosB
float e3 = (ratio * ratio) * (1.0f - (dotp * dotp)); // (n1/n2)^2 * (1 - cosB^2)
new_direction = e1 + ((e2 - sqrt(1.0f - e3)) * n); // n1/n2 * i + (n1/n2 * cosB - /(1 - sinT^2) * n
return normalize(new_direction);
}
float GetReflectance(vec3 direction, vec3 normal, float n1, float n2){
float raycast;
float cosI = dot(-direction, normal);
float ratio = n1/n2;
float sin2T = (ratio * ratio) * (1.0f - (cosI * cosI));
if(sin2T <= 1.0f){
// Reflection or Transmission
float cosT = sqrt(1.0f - sin2T);
float R1 = pow((n1 * cosI - n2 * cosT)/(n1 * cosI + n2 * cosT), 2);
float R2 = pow((n2 * cosI - n1 * cosT)/(n2 * cosI + n1 * cosT), 2);
float Ravg = (R1 + R2) / 2.0f;
raycast = Ravg;
}else{
// Total Internal Reflection
raycast = 1.0f;
}
return raycast;
}
vec3 RotateVector(vec3 vect, vec3 rotation){
rotation = DegToRad(rotation);
float xr = rotation.x;
float yr = rotation.y;
float zr = rotation.z;
mat3x3 Rx = mat3x3(
1, 0, 0,
0, cos(xr), -sin(xr),
0, sin(xr), cos(xr)
);
mat3x3 Ry = mat3x3(
cos(yr), 0, sin(yr),
0, 1, 0,
-sin(yr), 0, cos(yr)
);
mat3x3 Rz = mat3x3(
cos(zr), -sin(zr), 0,
sin(zr), cos(zr), 0,
0, 0, 1
);
return (((Ry*Rx)*Rz)*vect);
}
Ray CreateRay(vec3 origin, vec3 direction){
Ray ray;
ray.origin = origin;
ray.direction = normalize(direction);
return ray;
}
Ray CreateCameraRay(float x, float y){
float n = 0.1f;
float fov_v = (float(_fov)/2.0f) * (PI/180.0f);
float ratio = (float(_screen_size.x)/float(_screen_size.y));
float fov_h = atan( tan(fov_v) * ratio );
x = (x - 0.5f) * 2;
y = (y - 0.5f) * 2;
vec3 plane = vec3(x * tan(fov_h) * n,y * tan(fov_v) * n,n);
vec3 direction = RotateVector(normalize(plane),_camera_rotation.xyz);
return CreateRay(_camera_position.xyz,direction);
}
vec4 CheckIntersectionWithSphere(vec3 r1, vec3 r2, int i){
vec3 oc = r1 - objects[i].position.xyz;
float a = dot(r2,r2);
float b = 2.0f * dot(oc, r2);
float c = dot(oc, oc) - objects[i].size.x * objects[i].size.x;
float p = c/a;
float q = -b/a;
float d = b * b - (4 * a * c);
if(p < 0){
//inside
float val = (-b + sqrt(d)) / (2.0f*a);
vec3 hit = r1 + r2 * val;
return vec4(r1 + r2 * val, val);
}else{
if(q > 0){
//hit
float val = (-b - sqrt(d)) / (2.0f*a);
vec3 hit = r1 + r2 * val;
return vec4(r1 + r2 * val, val);
}else{
//not hit
return vec4(vec3(0,0,0),-1.0f);
}
}
}
vec4 CheckIntersectionWithPlane(vec3 r1, vec3 r2, int i){
vec3 normal = normalize(RotateVector(vec3(0,1,0), objects[i].orientation.xyz));
float dotP = dot(normal, r2);
vec3 pos = objects[i].position.xyz;
int ans = int(abs(dotP) >= 0.001f);
float d = dot((pos - r1),normal)/dotP;
return vec4(r1 + r2 * d, ans * (d + 1) - 1);
if(abs(dotP) >= 0.001f){
float d = dot((pos - r1),normal)/dotP;
if(d < 0)
return vec4(vec3(0,0,0),-1.0f);
return vec4(r1 + r2 * d, d);
}
return vec4(vec3(0,0,0),-1.0f);
}
vec3 GetTriangleVertex(int id, int object_id){
vec3 pos = objects[object_id].position.xyz;
vec3 size = objects[object_id].size.xyz;
vec3 rotation = objects[object_id].orientation.xyz;
id = id * 3;
vec3 vertex_from_array = vec3(vertices[id], vertices[id+1], vertices[id+2]);
vec3 vert = RotateVector(vertex_from_array * size, rotation);
return vert + pos;
}
vec2 GetTriangleUV(int id){
id = id * 2;
vec2 vertex_from_array = vec2(t_vertices[id], t_vertices[id+1]);
return vertex_from_array;
}
vec4 CheckIntersectionWithTriangleUV(vec3 r1, vec3 r2, int face_id, int object_id, inout vec2 uv, inout vec3 normal){
const float EPSILON = 0.0000001;
vec3 pos = objects[object_id].position.xyz;
vec3 size = objects[object_id].size.xyz;
vec3 vertex0 = GetTriangleVertex(faces[(face_id * 6) + 0], object_id);
vec3 vertex1 = GetTriangleVertex(faces[(face_id * 6) + 1], object_id);
vec3 vertex2 = GetTriangleVertex(faces[(face_id * 6) + 2], object_id);
vec3 edge1, edge2, h, s, q;
float a,f;
edge1 = vertex1 - vertex0;
edge2 = vertex2 - vertex0;
h = cross(r2, edge2);
vec3 N = cross(edge1,edge2); // N
bool reverse = dot(N, r2) > 0;
normal = (reverse ? -1 : 1) * normalize(N);
a = dot(edge1, h);
if (a > -EPSILON && a < EPSILON)
return vec4(0.0f, 0.0f, 0.0f, -1.0f); // This ray is parallel to this triangle.
f = 1.0/a;
s = r1 - vertex0;
uv.x = f * dot(s, h);
if (uv.x < 0.0 || uv.x > 1.0)
return vec4(0.0f, 0.0f, 0.0f, -1.0f); // This ray is parallel to this triangle.
q = cross(s, edge1);
uv.y = f * dot(r2, q);
if (uv.y < 0.0 || uv.x + uv.y > 1.0)
return vec4(0.0f, 0.0f, 0.0f, -1.0f); // This ray is parallel to this triangle.
// At this stage we can compute t to find out where the intersection point is on the line.
float t = f * dot(edge2, q);
if (t > EPSILON) // ray intersection
{
// map uv
float u = uv.x;
float v = uv.y;
float w = 1.0 - uv.x - uv.y;
// A(v0) - w
// B(v1) - u
// C(v2) - v
vec2 uv0 = GetTriangleUV(faces[(face_id * 6) + 3]);
vec2 uv1 = GetTriangleUV(faces[(face_id * 6) + 4]);
vec2 uv2 = GetTriangleUV(faces[(face_id * 6) + 5]);
uv = uv0 * w + uv1 * u + uv2 * v;
return vec4(r1 + r2 * t, t); // This ray is parallel to this triangle.
}
else // This means that there is a line intersection but not a ray intersection.
return vec4(0.0f, 0.0f, 0.0f, -1.0f); // This ray is parallel to this triangle.
}
vec4 CheckIntersectionWithMesh(vec3 r1, vec3 r2, int object_id, inout vec3 normal, inout vec2 uv){
vec4 raycast;
float val = 100000.0f;
vec3 hit_point;
vec3 temp_normal;
vec2 temp_uv;
int vert_num = objects[object_id].vert_num;
int offset = objects[object_id].offset;
for(int i = 0; i < vert_num; i++){
raycast = CheckIntersectionWithTriangleUV(r1, r2, offset + i, object_id, temp_uv, temp_normal);
if (raycast.w >= 0 && raycast.w < val){
hit_point = raycast.xyz;
val = raycast.w;
normal = temp_normal;
uv = temp_uv;
}
}
return vec4(hit_point, val);
}
float GetEnergy(vec3 color){
return dot(color, vec3(1.0f / 3.0f));
}
vec3 GetTextureColor(vec2 uv, int object_id, int map){
int id = objects[object_id].texture_id[map];
switch(id){
case 0:
return vec3(1.0f,1.0f,1.0f);
case 1:
return texture(u_tex1, uv).xyz;
case 2:
return texture(u_tex2, uv).xyz;
case 3:
return texture(u_tex3, uv).xyz;
case 4:
return texture(u_tex4, uv).xyz;
case 5:
return texture(u_tex5, uv).xyz;
default:
return vec3(1.0f,1.0f,1.0f);
}
}
vec3 GetNormal_old(vec3 pos, int object_id){
vec3 normal = vec3(0,1,0);
if(objects[object_id].type == 0) //sphere
normal = normalize(pos - objects[object_id].position.xyz);
else if(objects[object_id].type == 1) //plane
normal = normalize(RotateVector(vec3(0,1,0), objects[object_id].orientation.xyz));
else if(objects[object_id].type == 2){ //mesh
int v_offset = objects[object_id].offset;
vec3 v0 = GetTriangleVertex(faces[(v_offset * 6) + 0], object_id);
vec3 v1 = GetTriangleVertex(faces[(v_offset * 6) + 1], object_id);
vec3 v2 = GetTriangleVertex(faces[(v_offset * 6) + 2], object_id);
vec3 A = v1 - v0; // edge 0
vec3 B = v2 - v0; // edge 1
normal = normalize(cross(A, B));
}
return normal;
}
vec3 GetNormal(vec3 pos, vec3 dir, int object_id){
vec3 normal = vec3(0,1,0);
if(objects[object_id].type == 0) //sphere
normal = normalize(pos - objects[object_id].position.xyz);
else if(objects[object_id].type == 1) //plane
normal = normalize(RotateVector(vec3(0,1,0), objects[object_id].orientation.xyz));
else if(objects[object_id].type == 2){ //mesh
int v_offset = objects[object_id].offset;
vec3 v0 = GetTriangleVertex(faces[(v_offset * 6) + 0], object_id);
vec3 v1 = GetTriangleVertex(faces[(v_offset * 6) + 1], object_id);
vec3 v2 = GetTriangleVertex(faces[(v_offset * 6) + 2], object_id);
vec3 A = v1 - v0; // edge 0
vec3 B = v2 - v0; // edge 1
normal = normalize(cross(A, B));
}
bool reverse = dot(normal, dir) > 0;
normal = (reverse ? -1 : 1) * normalize(normal);
return normal;
}
vec3 GetNormal(vec3 pos, vec3 dir, int object_id, inout bool reversed){
vec3 normal = vec3(0,1,0);
if(objects[object_id].type == 0) //sphere
normal = normalize(pos - objects[object_id].position.xyz);
else if(objects[object_id].type == 1) //plane
normal = normalize(RotateVector(vec3(0,1,0), objects[object_id].orientation.xyz));
else if(objects[object_id].type == 2){ //mesh
int v_offset = objects[object_id].offset;
vec3 v0 = GetTriangleVertex(faces[(v_offset * 6) + 0], object_id);
vec3 v1 = GetTriangleVertex(faces[(v_offset * 6) + 1], object_id);
vec3 v2 = GetTriangleVertex(faces[(v_offset * 6) + 2], object_id);
vec3 A = v1 - v0; // edge 0
vec3 B = v2 - v0; // edge 1
normal = normalize(cross(A, B));
}
reversed = dot(normal, dir) > 0;
normal = (reversed ? -1 : 1) * normalize(normal);
return normal;
}
vec4 CastRay(Ray ray, inout vec3 normal, inout vec2 uv){
int object_id = -1;
vec4 raycast;
float val = distlimit;
vec3 hit_point;
vec3 temp_normal;
vec2 temp_uv;
for(int i = 0; i < _object_count; i++){
vec3 r1 = ray.origin;
vec3 r2 = ray.direction;
vec3 op = objects[i].position.xyz;
float r = objects[i].size.x;
if(objects[i].type == 0){
raycast = CheckIntersectionWithSphere(r1,r2,i);
if(raycast.w > 0){
bool reversed;
temp_normal = GetNormal(raycast.xyz,r2,i,reversed);
float theta = 1.0f - (acos((reversed ? -1 : 1) * temp_normal.y) / PI);
float phi = atan((reversed ? -1 : 1) * temp_normal.z, (reversed ? -1 : 1) * temp_normal.x) / (2 * PI) + 0.5f;
temp_uv.x = mod((phi + (objects[i].orientation[1] / (2.0 * PI))), 1);
temp_uv.y = theta;
}
}
else if(objects[i].type == 1){
raycast = CheckIntersectionWithPlane(r1,r2,i);
if(raycast.w > 0)
temp_normal = GetNormal(raycast.xyz,r2, i);
}
else if(objects[i].type == 2){
raycast = CheckIntersectionWithMesh(r1,r2,i, temp_normal, temp_uv);
}
else
raycast = vec4(0,0,0,distlimit);
if (raycast.w >= 0 && raycast.w < val){
hit_point = raycast.xyz;
val = raycast.w;
object_id = i;
normal = temp_normal;
uv = temp_uv;
}
}
return vec4(hit_point,object_id);
}
vec3 GetSkyColor(vec3 direction){
float theta = 1.0f - (acos(direction.y) / PI);
float phi = atan(direction.x, direction.z) / (2 * PI) + 0.5f;
return texture(u_skybox, vec2(phi, theta)).xyz;
}
void main(){
ivec2 id = ivec2(gl_GlobalInvocationID.xy);
vec4 pixel = vec4(0.0f, 0.0f, 0.0f, 1.0f);
vec4 raycast;
vec4 raycast_n;
vec3 color;
vec3 reflection;
vec3 normal;
vec3 albedo;
vec3 specular;
vec3 energy;
vec2 offset = _pixel_offset;
vec2 uv;
Ray reflection_ray;
Ray camera_ray;
float x, y;
float diffuse_chance;
float specular_chance;
float refraction_chance;
float sum;
float roulette;
float f;
float alpha;
float dist;
float air_transparency = 1.0f;
int face_id;
int object_id;
bool inside_object[100];
float refractive_indices_queue[100];
refractive_indices_queue[0] = air_transparency;
int refractive_indices_count = 1;
for(int i = 0; i < 100; i++){
inside_object[i] = false;
}
pixel_id = id;
for(int s = 0; s < _samples_per_frame; s++){
// map x,y to [0,1]
x = (((float(id.x) - 0.5f + offset.x)/_screen_size.x));
y = (((float(id.y) - 0.5f + offset.y)/_screen_size.y));
// new random offset
offset = vec2(1.0f * (abs(Rand()) - 0.5f) ,1.0f * (abs(Rand()) - 0.5f));
color = vec3(0.0f, 0.0f, 0.0f);
// get pixel for random generator
pixel_id = vec2(id) + offset;
// create a ray and cast it
camera_ray = CreateCameraRay(x,y);
raycast = CastRay(camera_ray, normal, uv);
object_id = int(raycast.w);
// if hit any object
if(object_id >= 0){
// calculate distance between rays
dist = length(_camera_position.xyz - raycast.xyz);
// initial energy values
energy = vec3(1.0f, 1.0f, 1.0f);
// add emmission from map
color += objects[object_id].emission.xyz * GetTextureColor(uv, object_id, 3);
//color += normal;
reflection_ray = camera_ray;
// recast ray multiple times
for(int i = 0; i < _iterations; i++){
roulette = Rand();
// get specular and albedo
specular = objects[object_id].specular.xyz * GetTextureColor(uv, object_id, 1);
albedo = min(1.0f - specular, objects[object_id].albedo.xyz * GetTextureColor(uv, object_id, 0));
// check if object is translucent
bool translucent = (objects[object_id].transparency > 0.01f);
if(translucent){
// refractive indicies
float n1;
float n2;
// 0 -> 1
if(!inside_object[object_id]){
// n1 is set to last index
n1 = refractive_indices_queue[refractive_indices_count-1];
// n2 is set to new object
n2 = objects[object_id].refractive_index;
}
else{ // 1 -> 0
// n1 is set to last index
n1 = refractive_indices_queue[refractive_indices_count-1];
// n2 is set to the one before that
n2 = refractive_indices_queue[refractive_indices_count-2];
}
f = 1.0f;
if(objects[object_id].smoothness <= 0.9999f){
alpha = SmoothnessToPhongAlpha(objects[object_id].smoothness);
normal = SampleHemisphere(normal, alpha);
f = (alpha + 2) / (alpha + 1);
}
// calculate reflection and transmission
float val_R = GetReflectance(reflection_ray.direction.xyz, normal, n1, n2);
float val_T = 1.0f - val_R;
bool refract = (roulette < val_T);
bool entering = !inside_object[object_id];
// update queue
if(refract){
inside_object[object_id] = !inside_object[object_id];
reflection = RefractRay(reflection_ray.direction, normal, n1, n2);
if(entering){
// add n2 to queue
refractive_indices_count++;
refractive_indices_queue[refractive_indices_count - 1] = n2;
}else{
// leaving (1 -> 0)
refractive_indices_count--;
}
}
else{
reflection = reflect(reflection_ray.direction, normal);
// no change to queue
}
vec3 ray_normal;
if(refract)
ray_normal = normal * -1;
else
ray_normal = normal;
// create and cast reflection ray
reflection_ray = CreateRay(raycast.xyz + ray_normal * 0.001f, reflection);
// calculate loss due to partial transparency
float transmittance = 1;
if(!entering)
transmittance = pow(10, -1 * (1.0f - objects[object_id].transparency) * dist);
// update energy
energy *= specular * transmittance;
}
else{
// calculate chances
diffuse_chance = GetEnergy(albedo.xyz) + 0.0001f;
specular_chance = GetEnergy(specular.xyz) + 0.0001f;
sum = diffuse_chance + specular_chance;
diffuse_chance /= sum;
specular_chance /= sum;
if(roulette < diffuse_chance){
//diffuse
reflection = SampleHemisphere(normal, 1.0f);
// update energy
energy *= albedo * 2 * sdot(normal, reflection);
}else{
// specular
reflection = normalize(reflection_ray.direction - 2 * dot(reflection_ray.direction, normal) * normal);
f = 1.0f;
if(objects[object_id].smoothness <= 0.9999f){
alpha = SmoothnessToPhongAlpha(objects[object_id].smoothness);
reflection = SampleHemisphere(reflection, alpha);
f = (alpha + 2) / (alpha + 1);
}
// update energy
energy *= specular * 2 * sdot(normal, reflection);
}
//create and cast reflected ray
reflection_ray = CreateRay(raycast.xyz + normal * 0.001f, reflection);
}
// cast new ray
raycast_n = CastRay(reflection_ray, normal, uv);
object_id = int(raycast_n.w);
//if hit sky, break and get color * energy
if(object_id == -1){
color += energy * GetSkyColor(reflection_ray.direction);
break;
}else{
color += energy * objects[object_id].emission.xyz * GetTextureColor(uv, object_id, 3);
dist = length(raycast_n.xyz - raycast.xyz);
}
// if energy is low, no need to keep going
if((energy.x + energy.y + energy.z <= 0.01f && i > 1)){
break;
}
raycast = raycast_n;
}
}
else{
// sample skybox
color = GetSkyColor(camera_ray.direction);
}
// add sample
pixel += vec4(color.xyz, 0.0f);
}
// get average over samples
pixel = vec4(pixel.xyz/float(_samples_per_frame), 1.0f);
// update the new average
vec4 old = imageLoad(img_output, id);
pixel = old + (vec4(pixel) - old)/_sample;
// update pixel in texture
imageStore(img_output, id, pixel);
}