From be918c5898dbb10b0a60bb2a7e31b23c12b654e3 Mon Sep 17 00:00:00 2001
From: Martin Opat <martinopat.opat@gmail.com>
Date: Mon, 30 Dec 2024 10:47:51 +0100
Subject: [PATCH] Exported ray-caster into its own file

---
 src/illumination/raycaster.h | 127 +++++++++++++++++++++++++++++++++++
 1 file changed, 127 insertions(+)
 create mode 100644 src/illumination/raycaster.h

diff --git a/src/illumination/raycaster.h b/src/illumination/raycaster.h
new file mode 100644
index 0000000..8677944
--- /dev/null
+++ b/src/illumination/raycaster.h
@@ -0,0 +1,127 @@
+#ifndef RAYCASTER_H
+#define RAYCASTER_H
+
+#include <cuda_runtime.h>
+#include "linalg/linalg.h"
+#include "consts.h"
+#include "shading.h"
+
+
+
+// Raycast + phong, TODO: Consider wrapping in a class
+__global__ void raycastKernel(float* volumeData, unsigned char* framebuffer, int d_volumeWidth, int d_volumeHeight, int d_volumeDepth) {
+    int px = blockIdx.x * blockDim.x + threadIdx.x;
+    int py = blockIdx.y * blockDim.y + threadIdx.y;
+    if (px >= IMAGE_WIDTH || py >= IMAGE_HEIGHT) return;
+
+    float accumR = 0.0f;
+    float accumG = 0.0f;
+    float accumB = 0.0f;
+
+    // Multiple samples per pixel
+    for (int s = 0; s < SAMPLES_PER_PIXEL; s++) {
+        // Map to [-1, 1]
+        float u = ((px + 0.5f) / IMAGE_WIDTH ) * 2.0f - 1.0f;
+        float v = ((py + 0.5f) / IMAGE_HEIGHT) * 2.0f - 1.0f;
+
+        // TODO: Move this (and all similar transformation code) to its own separate file
+        float tanHalfFov = tanf(fov * 0.5f);
+        u *= tanHalfFov;
+        v *= tanHalfFov;
+
+        // Find ray direction
+        Vec3 cameraRight = (cameraDir.cross(cameraUp)).normalize();
+        cameraUp = (cameraRight.cross(cameraDir)).normalize();
+        Vec3 rayDir = (cameraDir + cameraRight*u + cameraUp*v).normalize();
+
+        // Intersect (for simplicity just a 3D box from 0 to 1 in all dimensions) - TODO: Think about whether this is the best way to do this
+        float tNear = 0.0f;
+        float tFar  = 1e6f;
+        auto intersectAxis = [&](float start, float dirVal) {
+            if (fabsf(dirVal) < epsilon) {
+                if (start < 0.f || start > 1.f) {
+                    tNear = 1e9f;
+                    tFar  = -1e9f;
+                }
+            } else {
+                float t0 = (0.0f - start) / dirVal;
+                float t1 = (1.0f - start) / dirVal;
+                if (t0>t1) { 
+                    float tmp=t0; 
+                    t0=t1; 
+                    t1=tmp; 
+                }
+                if (t0>tNear) tNear = t0;
+                if (t1<tFar ) tFar  = t1;
+            }
+        };
+
+        intersectAxis(cameraPos.x, rayDir.x);
+        intersectAxis(cameraPos.y, rayDir.y);
+        intersectAxis(cameraPos.z, rayDir.z);
+
+        if (tNear > tFar) continue;  // No intersectionn
+        if (tNear < 0.0f) tNear = 0.0f;
+
+        float colorR = 0.0f, colorG = 0.0f, colorB = 0.0f;
+        float alphaAccum = 0.0f;
+
+        float tCurrent = tNear;
+        while (tCurrent < tFar && alphaAccum < alphaAcumLimit) {
+            Point3 pos = cameraPos + rayDir * tCurrent;
+
+            // Convert to volume indices
+            float fx = pos.x * (d_volumeWidth  - 1);
+            float fy = pos.y * (d_volumeHeight - 1);
+            float fz = pos.z * (d_volumeDepth  - 1);
+            int ix = (int)roundf(fx);
+            int iy = (int)roundf(fy);
+            int iz = (int)roundf(fz);
+
+            // Sample
+            float density = sampleVolumeNearest(volumeData, d_volumeWidth, d_volumeHeight, d_volumeDepth, ix, iy, iz);
+
+            // Basic transfer function. TODO: Move to a separate file, and then improve
+            float alphaSample = density * 0.1f;
+            // float alphaSample = 1.0f - expf(-density * 0.1f);
+            Color3 baseColor = Color3(density, 0.1f*density, 1.f - density);  // TODO: Implement a proper transfer function
+
+            // If density ~ 0, skip shading
+            if (density > minAllowedDensity) {
+                Vec3 grad = computeGradient(volumeData, d_volumeWidth, d_volumeHeight, d_volumeDepth, ix, iy, iz);
+                Vec3 normal = -grad.normalize();
+
+                Vec3 lightDir = (lightPos - pos).normalize();
+                Vec3 viewDir  = -rayDir.normalize();
+
+                // Apply Phong
+                Vec3 shadedColor = phongShading(normal, lightDir, viewDir, baseColor);
+
+                // Compose
+                colorR     += (1.0f - alphaAccum) * shadedColor.x * alphaSample;
+                colorG     += (1.0f - alphaAccum) * shadedColor.y * alphaSample;
+                colorB     += (1.0f - alphaAccum) * shadedColor.z * alphaSample;
+                alphaAccum += (1.0f - alphaAccum) * alphaSample;
+            }
+
+            tCurrent += stepSize;
+        }
+
+        accumR += colorR;
+        accumG += colorG;
+        accumB += colorB;
+    }
+
+    // Average samples
+    accumR /= (float)SAMPLES_PER_PIXEL;
+    accumG /= (float)SAMPLES_PER_PIXEL;
+    accumB /= (float)SAMPLES_PER_PIXEL;
+
+    // Final colour
+    int fbIndex = (py * IMAGE_WIDTH + px) * 3;
+    framebuffer[fbIndex + 0] = (unsigned char)(fminf(accumR, 1.f) * 255);
+    framebuffer[fbIndex + 1] = (unsigned char)(fminf(accumG, 1.f) * 255);
+    framebuffer[fbIndex + 2] = (unsigned char)(fminf(accumB, 1.f) * 255);
+}
+
+#endif // RAYCASTER_H
\ No newline at end of file