Add Reflection, Clamp, and Merge Header Files

I added 3D Vector Reflection which reflects a vector along a normal, a function that clamps a variable, and merged header files so that it will be possible to do math with random vectors.

vectormath.h → custommath.h

@@ -1,12 +1,41 @@
 #include <math.h>
 
+// Classic "Unsafe" Minimum and Maximum Functions
+// People on Stack Exchange said it was bad but it works Fine for me...
+#define min(a, b) (a < b ? a : b)
+#define max(a, b) (a > b ? a : b)
+
+// Clamp Function
+#define clamp(n, a, b) max(a, min(b, n))
+
+// High-Quality Random Numbers from Michael0884, ported from: https://www.shadertoy.com/view/wltcRS
+unsigned int ns;
+#define INIT_RNG ns = 185730U*sample+(x+y*resolutionx);
+
+void pcg(){
+    unsigned int state = ns*747796405U+2891336453U;
+    unsigned int word = ((state >> ((state >> 28U) + 4U)) ^ state)*277803737U;
+    ns = (word >> 22U) ^ word;
+}
+
+float random(){
+    pcg();
+    return (float)ns/(float)0xffffffffU;
+}
+
 // 3D Vector
 typedef struct{
     float x;
     float y;
     float z;
 } vec3;
 
+// 2D Vector
+typedef struct{
+    float x;
+    float y;
+} vec2;
+
 // Convert 3 Floating-Point Numbers to a 3D Vector
 vec3 float3(float x, float y, float z){
     vec3 vec;
@@ -16,6 +45,14 @@ vec3 float3(float x, float y, float z){
     return vec;
 };
 
+// Convert 3 Floating-Point Numbers to a 3D Vector
+vec2 float2(float x, float y){
+    vec2 vec;
+    vec.x = x;
+    vec.y = y;
+    return vec;
+};
+
 // Convert a Floating-Point Number to a 3D Vector
 vec3 floatf3(float x){
     vec3 vec;
@@ -76,4 +113,12 @@ float vector3Length(vec3 vec){
 vec3 normalize(vec3 vec){
     float length = vector3Length(vec);
     return float3(vec.x/length, vec.y/length, vec.z/length);
-}
+}
+
+vec3 reflect(vec3 vector, vec3 normal){
+    return vec3Sub(vector, vec3Mult(floatf3(2.0*vec3dotp(vector, normal)), normal));
+}
+
+// TODO: Sphereically Uniform Random (https://www.shadertoy.com/view/WttyWX)
+//vec2 nrand2(float sigma, vec2 mean){vec2 Z = rand2(); return mean+sigma*sqrt(-2.0*log(Z.x))*vec2(cos(twopi*Z.y),sin(twopi*Z.y));}
+//vec3 nrand3(float sigma, vec3 mean){vec4 Z = rand4(); return mean+sigma*sqrt(-2.0*log(Z.xxy))*vec3(cos(twopi*Z.z), sin(twopi*Z.z), cos(twopi*Z.w));}

pathtrace.c

@@ -5,21 +5,17 @@
 
 // Include Custom Utilities
 #include "targa.h"
-#include "rng.h"
-#include "vectormath.h"
+#include "custommath.h"
 
-#define samples     (unsigned int) 1024U
-#define resolutionx (unsigned int) 320U
-#define resolutiony (unsigned int) 240U
-#define maxmarches  (unsigned int) 128U
+// Parameters
+#define samples     256U
+#define resolutionx 320U
+#define resolutiony 240U
+#define maxmarches  128U
 #define collisiondist 0.01
 #define scenesize 32.0
 #define camfov 1.0
 
-// Classic "Unsafe" Minimum and Maximum Functions
-#define MIN(a, b) (a < b ? a : b)
-#define MAX(a, b) (a > b ? a : b)
-
 float distestimator(vec3 pos){
     return vector3Length(pos)-0.5;
 }
@@ -29,7 +25,7 @@ float light(vec3 pos){
 }
 
 vec3 pathtrace(vec3 raydir, vec3 rayori){
-    vec3 lastpos = rayori, raypos = rayori;
+    vec3 raypos = rayori, lastpos = rayori, normal;
     float distest;
     for(unsigned int i = 0U; i < maxmarches; i++){
         if(vector3Length(raypos) > scenesize){break;}
@@ -38,21 +34,22 @@ vec3 pathtrace(vec3 raydir, vec3 rayori){
         }
         distest = distestimator(raypos);
         if(distest < collisiondist){
+            normal = normalize(raypos);
             // TODO: Spherically Uniform Random- Right now this is Basically a Cube which is Incorrect
-            raydir = normalize(float3(rand()-0.5, rand()-0.5, rand()-0.5));
+            raydir = reflect(raydir, normalize(vec3Add(normal, vec3Mult(floatf3(2.0), float3(random()-0.5, random()-0.5, random()-0.5)))));
             raypos = lastpos;
         }
         lastpos = raypos;
-        raypos = vec3Add(raypos, vec3Mult(raydir, floatf3(0.5*MIN(distest, light(raypos)))));
+        raypos = vec3Add(raypos, vec3Mult(raydir, floatf3(min(distest, light(raypos)))));
     }
     return floatf3(0.0);
 }
 
 int main(){
     // Set-Up Variables
-    unsigned int resolutionmax = MAX(resolutionx, resolutiony), byte = 0U;
-    vec3 raydir, outCol;
-    float uvx, uvy;
+    unsigned int resolutionmax = max(resolutionx, resolutiony), byte = 0U;
+    vec3 raydir, normal, outCol;
+    vec2 uv;
     uint8_t imageBuffer[resolutionx * resolutiony * 3U] = {0U};
 
     // Execute Rendering
@@ -61,16 +58,15 @@ int main(){
         // Monte-Carlo Rendering
         for(unsigned int sample = 0U; sample < samples; sample++){
             INIT_RNG;
-            uvx = 2.0*((x+0.5)-0.5*resolutionx)/resolutionmax;
-            uvy = 2.0*((y+0.5)-0.5*resolutiony)/resolutionmax;
-            raydir = normalize(float3(uvx*camfov, uvy*camfov, 1.0));
+            uv = float2(2.0*((x+0.5)-0.5*resolutionx)/resolutionmax, 2.0*((y+0.5)-0.5*resolutiony)/resolutionmax);
+            raydir = normalize(float3(uv.x*camfov, uv.y*camfov, 1.0));
             outCol = vec3Add(outCol, vec3Mult(pathtrace(raydir, float3(0.0, 0.0, -2.0)), floatf3(32.0)));
         }
         outCol = vec3Div(outCol, floatf3(samples));
-        imageBuffer[byte   ] = MAX(0U, MIN(255U, round(outCol.x*255U)));
-        imageBuffer[byte+1U] = MAX(0U, MIN(255U, round(outCol.y*255U)));
-        imageBuffer[byte+2U] = MAX(0U, MIN(255U, round(outCol.z*255U)));
-        byte += 3;
+        imageBuffer[byte   ] = clamp(round(outCol.z*255U), 0U, 255U);
+        imageBuffer[byte+1U] = clamp(round(outCol.y*255U), 0U, 255U);
+        imageBuffer[byte+2U] = clamp(round(outCol.x*255U), 0U, 255U);
+        byte += 3U;
     }
     }