Forgot to do atomic commits so ... booyah

src/consts.cu

@@ -5,8 +5,9 @@ __device__ Vec3 d_cameraDir;
 __device__ Vec3 d_cameraUp;
 __device__ Point3 d_lightPos;
 
-Point3 h_cameraPos = Point3::init(-0.7, -1.0, -2.0);
+Point3 h_cameraPos = Point3::init(-500.0f, 200.0f, 100.0f);
-Vec3 h_cameraDir = Vec3::init(0.4, 0.6, 1.0).normalize();
+Vec3 center = Vec3::init((float)VOLUME_WIDTH/2.0f, (float)VOLUME_HEIGHT/2.0f, (float)VOLUME_DEPTH/2.0f);
+Vec3 h_cameraDir = (center - h_cameraPos).normalize();
 Vec3 h_cameraUp = Vec3::init(0.0, 1.0, 0.0).normalize();
 Point3 h_lightPos = Point3::init(1.5, 2.0, -1.0);
 

src/consts.h

@@ -5,8 +5,8 @@
 #include <cmath>
 
 // --------------------------- Basic Constants ---------------------------
-const int VOLUME_WIDTH  = 49;
+const int VOLUME_WIDTH  = 576;
-const int VOLUME_HEIGHT = 51;
+const int VOLUME_HEIGHT = 361;
 const int VOLUME_DEPTH  = 42;
 
 const int IMAGE_WIDTH   = 800;

src/illumination/Raycaster.cu

@@ -40,15 +40,16 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) {
         // 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) {
+        auto intersectAxis = [&](float start, float dir, float minV, float maxV) {
-            if (fabsf(dirVal) < epsilon) {
+            if (fabsf(dir) < epsilon) {
-                if (start < 0.f || start > 1.f) {
+                // Ray parallel to axis. If outside min..max, no intersection.
+                if (start < minV || start > maxV) {
                     tNear = 1e9f;
                     tFar  = -1e9f;
                 }
             } else {
-                float t0 = (0.0f - start) / dirVal;
+                float t0 = (minV - start) / dir;
-                float t1 = (1.0f - start) / dirVal;
+                float t1 = (maxV - start) / dir;
                 if (t0 > t1) {
                     float tmp = t0;
                     t0 = t1;
@@ -59,9 +60,9 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) {
             }
         };
 
-        intersectAxis(d_cameraPos.x, rayDir.x);
+        intersectAxis(d_cameraPos.x, rayDir.x, 0.0f, (float)VOLUME_WIDTH);
-        intersectAxis(d_cameraPos.y, rayDir.y);
+        intersectAxis(d_cameraPos.y, rayDir.y, 0.0f, (float)VOLUME_HEIGHT);
-        intersectAxis(d_cameraPos.z, rayDir.z);
+        intersectAxis(d_cameraPos.z, rayDir.z, 0.0f, (float)VOLUME_DEPTH);
 
         if (tNear > tFar) continue;  // No intersectionn
         if (tNear < 0.0f) tNear = 0.0f;
@@ -74,12 +75,12 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) {
             Point3 pos = d_cameraPos + rayDir * tCurrent;
 
             // Convert to volume indices
-            float fx = pos.x * (VOLUME_WIDTH  - 1);
+            // float fx = pos.x * (VOLUME_WIDTH  - 1);
-            float fy = pos.y * (VOLUME_HEIGHT - 1);
+            // float fy = pos.y * (VOLUME_HEIGHT - 1);
-            float fz = pos.z * (VOLUME_DEPTH  - 1);
+            // float fz = pos.z * (VOLUME_DEPTH  - 1);
-            int ix = (int)roundf(fx);
+            int ix = (int)roundf(pos.x);
-            int iy = (int)roundf(fy);
+            int iy = (int)roundf(pos.y);
-            int iz = (int)roundf(fz);
+            int iz = (int)roundf(pos.z);
 
             // Sample
             float density = sampleVolumeNearest(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz);
@@ -107,6 +108,9 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) {
                 alphaAccum += (1.0f - alphaAccum) * alphaSample;
             }
 
+            // TODO: Add silhouette - Take gradient of volume dot with view direction (if small then this is a silhouette)
+            // TODO: Scale volume correctly
+
             tCurrent += stepSize;
         }
 

src/main.cu

@@ -20,7 +20,8 @@ static float* d_volume = nullptr;
 // * very similarly - actual code for loading new data as the simulation progresses - right now its effectively a static image loader
 
 void getTemperature(std::vector<float>& temperatureData, int idx = 0) {
-    std::string path = "data/trimmed";
+    // std::string path = "data/trimmed";
+    std::string path = "data";
     std::string variable = "T";
     DataReader dataReader(path, variable);
     size_t dataLength = dataReader.fileLength(idx);
@@ -29,7 +30,8 @@ void getTemperature(std::vector<float>& temperatureData, int idx = 0) {
 }
 
 void getSpeed(std::vector<float>& speedData, int idx = 0) {
-    std::string path = "data/trimmed";
+    // std::string path = "data/trimmed";
+    std::string path = "data";
     std::string varU = "U";
     std::string varV = "V";
 
@@ -52,9 +54,10 @@ void getSpeed(std::vector<float>& speedData, int idx = 0) {
 
 int main() {
   std::vector<float> data;
-  // getTemperature(data);
+  getTemperature(data);
-  getSpeed(data);
+  // getSpeed(data);
 
+  std::cout << "DATA size: " << data.size() << std::endl;
 
   // TODO: Eveontually remove debug below (i.e., eliminate for-loop etc.)
   // Generate debug volume data
@@ -62,13 +65,15 @@ int main() {
   // generateVolume(hostVolume, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH);
   for (int i = 0; i < VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH; i++) {  // TODO: This is technically an unnecessary artifact of the old code taking in a float* instead of a std::vector
     // Discard temperatures above a small star (supposedly, missing temperature values)
-    hostVolume[i] = data[i];
+    hostVolume[i] = data[i + 3*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH];
-    if (data[i] + epsilon >= infty) hostVolume[i] = 0.0f;
+    if (data[i + 3*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH] + epsilon >= infty) hostVolume[i] = 0.0f;
   }
 
   // Min-max normalization
   float minVal = *std::min_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH);
   float maxVal = *std::max_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH);
+  std::cout << "minVal: " << minVal << " maxVal: " << maxVal << std::endl;
+
   for (int i = 0; i < VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH; i++) {
     hostVolume[i] = (hostVolume[i] - minVal) / (maxVal - minVal);
   }