Playing around with values to actually get something not fully opaque

src/consts.cu

@@ -5,10 +5,11 @@ __device__ Vec3 d_cameraDir;
 __device__ Vec3 d_cameraUp;
 __device__ Point3 d_lightPos;
 
-Point3 h_cameraPos = Point3::init(-300.0f, 200.0f, -300.0f);
+// Point3 h_cameraPos = Point3::init(300.0f, 200.0f, -700.0f);  // Camera for full data set
+Point3 h_cameraPos = Point3::init(50.0f, -50.0f, -75.0f);  // Camera for partially trimmed data set (TODO: Probably upside down atm)
 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();
+Vec3 h_cameraUp = Vec3::init(0.0, 0.0, 1.0).normalize();
 Point3 h_lightPos = Point3::init(1.5, 2.0, -1.0);
 
 void copyConstantsToDevice() {

src/consts.h

@@ -5,9 +5,12 @@
 #include <cmath>
 
 // --------------------------- Basic Constants ---------------------------
-const int VOLUME_WIDTH  = 576;
+// TODO: Right now this corresponds to the data set resolution (i.e., voxel per block in volume), however, we can separate this to allow for higher-resolution rendering
-const int VOLUME_HEIGHT = 361;
+// const int VOLUME_WIDTH  = 576;  // lon
-const int VOLUME_DEPTH  = 42;
+const int VOLUME_WIDTH  = 97;  // lon
+// const int VOLUME_HEIGHT = 361;  // lat
+const int VOLUME_HEIGHT = 71;  // lat
+const int VOLUME_DEPTH  = 42;  // lev
 
 const int IMAGE_WIDTH   = 1600;
 const int IMAGE_HEIGHT  = 1200;
@@ -15,14 +18,18 @@ const int IMAGE_HEIGHT  = 1200;
 const double epsilon = 1e-10f;
 const double infty   = 1e15f;  // This value is used to represent missing values in data
 
+// --------------------------- Dataset Constants ---------------------------
+const float MIN_TEMP = 184.0f;
+const float MAX_TEMP = 313.0f;
+
 
 // --------------------------- Raycasting Constants ---------------------------
-const int SAMPLES_PER_PIXEL = 8;  // TODO: Right now uses simple variance, consider using something more advanced (e.g., some commonly-used noise map)
+const int SAMPLES_PER_PIXEL = 1;  // TODO: Right now uses simple variance, consider using something more advanced (e.g., some commonly-used noise map)
 
-const float alphaAcumLimit = 0.65f;   // TODO: Idk what a good accumulation value is
+const float alphaAcumLimit = 0.4f;   // TODO: Idk what a good accumulation value is  <--- This finally does something when using alpha in both places at least
 const float minAllowedDensity = 0.001f;
 
-const float stepSize = 0.002f;
+const float stepSize = 0.02f;
 
 
 // --------------------------- Illumination Constants ---------------------------

src/gui/MainWindow.cpp

@@ -6,6 +6,17 @@
 
 #include "Shader.h"
 
+
+// TODO: Delete
+void saveImage2(const char* filename, unsigned char* framebuffer, int width, int height) {  // TODO: Figure out a better way to do this
+    std::ofstream imageFile(filename, std::ios::out | std::ios::binary);
+    imageFile << "P6\n" << width << " " << height << "\n255\n";
+    for (int i = 0; i < width * height * 3; i++) {
+        imageFile << framebuffer[i];
+    }
+    imageFile.close();
+}
+
 Window::Window(unsigned int w, unsigned int h) {
   this->w = w;
   this->h = h;
@@ -53,9 +64,17 @@ int Window::init(float* data) {
   if (init_quad(data)) return -1;
 	this->last_frame = std::chrono::steady_clock::now();
 
-  while (!glfwWindowShouldClose(window)) {
+  // TODO: These changes are temporary
+  // while (!glfwWindowShouldClose(window)) {
+  //   Window::tick();
+  // }
   Window::tick();
-  }
+  Window::tick();
+  // Save the image
+  unsigned char* pixels = new unsigned char[this->w * this->h * 3];
+  glReadPixels(0, 0, this->w, this->h, GL_RGB, GL_UNSIGNED_BYTE, pixels);
+  saveImage2("output.ppm", pixels, this->w, this->h);
+  delete[] pixels;
 
   Window::free(data);
   return 0;

src/illumination/Raycaster.cu

@@ -45,42 +45,92 @@ __device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const
     return c0 * (1.0f - dz) + c1 * dz;
 }
 
-// Function to map a temperature to an RGB color
+__device__ float opacityFromGradient(const Vec3 &grad) {
-__device__ Color3 temperatureToRGB(float temperature) {
+    float gradMag = grad.length();  // magnitude
-    // atm, the scalar field is normalized
+    float k = 1e-6f;               // tweak (the smaller the value, the less opacity)  // TODO: What should be the value of this?
-    const float minTemp = 0.0f; // coldest == deep blue
+    float alpha = 1.0f - expf(-k * gradMag);
-    const float maxTemp = 1.0f;  // hottest temperature == deep red
+    return alpha;
-
-    temperature = clamp(temperature, minTemp, maxTemp);
-    float t = normalize(temperature, minTemp, maxTemp);
-
-    float r, g, b;
-
-    if (t < 0.5f) { // From blue to green
-        t *= 2.0f; // Scale to [0, 1]
-        r = interpolate(0.0f, 0.0f, t);
-        g = interpolate(0.0f, 1.0f, t);
-        b = interpolate(1.0f, 0.0f, t);
-    } else { // From green to red
-        t = (t - 0.5f) * 2.0f; // Scale to [0, 1]
-        r = interpolate(0.0f, 1.0f, t);
-        g = interpolate(1.0f, 0.0f, t);
-        b = interpolate(0.0f, 0.0f, t);
 }
 
-    return Color3::init(r, g, b);
+struct ColorStop
+{
+    float pos;       // in [0,1]  
+    Color3 color;    // R,G,B in [0,1]
+};
+
+// TODO: Rename probably
+__device__ Color3 colorMapViridis(float normalizedT) {
+    // Here we redefine the color stops to go from deep blue (0.0) to purple (0.5) to deep red (1.0)
+    ColorStop tempStops[] = {
+        { 0.0f, Color3::init(0.0f, 0.0f, 1.0f) },   // deep blue
+        { 0.5f, Color3::init(0.5f, 0.0f, 0.5f) },   // purple
+        { 1.0f, Color3::init(1.0f, 0.0f, 0.0f) }    // deep red
+    };
+
+    // Clamp to [0,1]
+    if (normalizedT < 0.0f) normalizedT = 0.0f;
+    if (normalizedT > 1.0f) normalizedT = 1.0f;
+
+    // We have 3 stops => 2 intervals
+    const int N = 3;
+    for (int i = 0; i < N - 1; ++i)
+    {
+        float start = tempStops[i].pos;
+        float end   = tempStops[i + 1].pos;
+
+        if (normalizedT >= start && normalizedT <= end)
+        {
+            float localT = (normalizedT - start) / (end - start);
+            return interpolate(tempStops[i].color, tempStops[i + 1].color, localT);
         }
+    }
+    // Fallback if something goes out of [0,1] or numerical issues
+    return tempStops[N - 1].color;
+}
+
+// TODO: This is the old colour map, probably delete
+// // Function to map a temperature to an RGB color
+// __device__ Color3 temperatureToRGB(float temperature) {
+//     // atm, the scalar field is normalized
+//     const float minTemp = 184.f; // coldest == deep blue
+//     const float maxTemp = 312.f;  // hottest temperature == deep red
+
+//     if (temperature < minTemp) {
+//       return Color3::init(1.f, 1.f, 1.f);
+//     }
+//     temperature = clamp(temperature, minTemp, maxTemp);
+//     float t = normalize(temperature, minTemp, maxTemp);
+
+//     float r, g, b;
+
+//     if (t < 0.5f) { // From blue to green
+//         t *= 2.0f; // Scale to [0, 1]
+//         r = interpolate(0.0f, 0.0f, t);
+//         g = interpolate(0.0f, 1.0f, t);
+//         b = interpolate(1.0f, 0.0f, t);
+//     } else { // From green to red
+//         t = (t - 0.5f) * 2.0f; // Scale to [0, 1]
+//         r = interpolate(0.0f, 1.0f, t);
+//         g = interpolate(1.0f, 0.0f, t);
+//         b = interpolate(0.0f, 0.0f, t);
+//     }
+
+//     return Color3::init(r, g, b);
+// }
 
 
 // Transfer function
-__device__ float4 transferFunction(float density, Vec3 grad, Point3 pos, Vec3 rayDir) {
+__device__ float4 transferFunction(float density, const Vec3& grad, const Point3& pos, const Vec3& rayDir) {
-  float4 result;
   // Basic transfer function. TODO: Move to a separate file, and then improve
-  float alphaSample = density * 0.1f;
-  // result.w = 1.0f - expf(-density * 0.1f);
 
   // Color3 baseColor = Color3::init(density, 0.1f*density, 1.f - density);  // TODO: Implement a proper transfer function
-  Color3 baseColor = temperatureToRGB(density);
+  // Color3 baseColor = temperatureToRGB(density);
+  
+  float normDensity = (density - MIN_TEMP) / (MAX_TEMP - MIN_TEMP);
+  normDensity = clamp(normDensity, 0.0f, 1.0f);
+  Color3 baseColor = colorMapViridis(normDensity);
+  float alpha = opacityFromGradient(grad);
+  float alphaSample = density * alpha;  // TODO: Decide whether to keep alpha here or not
 
   Vec3 normal = -grad.normalize();
 
@@ -91,13 +141,14 @@ __device__ float4 transferFunction(float density, Vec3 grad, Point3 pos, Vec3 ra
   Vec3 shadedColor = phongShading(normal, lightDir, viewDir, baseColor);
 
   // Compose
-  result.x = (1.0f - alphaSample) * shadedColor.x * alphaSample;
+  float4 result;
-  result.y = (1.0f - alphaSample) * shadedColor.y * alphaSample;
+  result.x = shadedColor.x * alphaSample;
-  result.z = (1.0f - alphaSample) * shadedColor.z * alphaSample;
+  result.y = shadedColor.y * alphaSample;
-  result.w = (1.0f - alphaSample) * alphaSample;
+  result.z = shadedColor.z * alphaSample;
+  result.w = alpha;  // TODO: Again, decide if alpha here is correct or not
 
-  // TODO: Add silhouette - Take gradient of volume dot with view direction (if small then this is a silhouette)
+  // TODO: This is the black silhouette, technically if we are doing alpha based on gradient then it's kind of redundant (?) ... but could also be used for even sharper edges
-  if (grad.dot(viewDir) < epsilon / 100000.0f) {
+  if (grad.length() > epsilon && fabs(grad.normalize().dot(viewDir)) < 0.2f) {
     result.x = 0.0f;
     result.y = 0.0f;
     result.z = 0.0f;
@@ -163,10 +214,10 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) {
         intersectAxis(d_cameraPos.z, rayDir.z, 0.0f, (float)VOLUME_DEPTH);
 
         if (tNear > tFar){
-          // No intersectionn
+          // No intersection -> Set to brackground color (multiply by SAMPLES_PER_PIXEL because we divide by it later)
-          accumR = 0.9f;
+          accumR = 0.1f * (float)SAMPLES_PER_PIXEL;
-          accumG = 0.9f;
+          accumG = 0.1f * (float)SAMPLES_PER_PIXEL;
-          accumB = 0.9f;
+          accumB = 0.1f * (float)SAMPLES_PER_PIXEL;
         } else {
           if (tNear < 0.0f) tNear = 0.0f;
 
@@ -190,10 +241,10 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) {
               if (density > minAllowedDensity) {
                 Vec3 grad = computeGradient(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz);
                 float4 color = transferFunction(density, grad, pos, rayDir);
-                colorR += color.x;
+                colorR += color.x * (alphaAcumLimit - alphaAccum);
-                colorG += color.y;
+                colorG += color.y * (alphaAcumLimit - alphaAccum);
-                colorB += color.z;
+                colorB += color.z * (alphaAcumLimit - alphaAccum);
-                alphaAccum += color.w;
+                alphaAccum += color.w * (alphaAcumLimit - alphaAccum);
               }
 
 
@@ -204,13 +255,14 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) {
           accumG += colorG;
           accumB += colorB;
 
-          // float leftover = 1.0f - alphaAccum;
+          float leftover = 1.0 - alphaAccum;
-          // accumR = accumR + leftover * 0.9f;
+          accumR = accumR + leftover * 0.1f;
-          // accumG = accumG + leftover * 0.9f;
+          accumG = accumG + leftover * 0.1f;
-          // accumB = accumB + leftover * 0.9f;
+          accumB = accumB + leftover * 0.1f;
         }
     }
 
+
     // Average samples
     accumR /= (float)SAMPLES_PER_PIXEL;
     accumG /= (float)SAMPLES_PER_PIXEL;

src/main.cu

@@ -21,8 +21,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 path = "data";
     std::string variable = "T";
     DataReader dataReader(path, variable);
     size_t dataLength = dataReader.fileLength(idx);
@@ -31,8 +31,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 path = "data";
     std::string varU = "U";
     std::string varV = "V";
 
@@ -64,23 +64,29 @@ int main() {
   // Generate debug volume data
   float* hostVolume = new float[VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH];
   // generateVolume(hostVolume, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH);
+  int inftyCount=0;
   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 + 3*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH];
-    if (data[i + 3*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH] + epsilon >= infty) hostVolume[i] = 0.0f;
+    if (data[i + 3*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH] + epsilon >= infty) {hostVolume[i] = -infty; inftyCount++;}
   }
+  std::cout << "inftyCount: " << inftyCount << std::endl;
 
-  // Min-max normalization
+  float minVal = *std::min_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH, [](float a, float b) {
-  float minVal = *std::min_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH);
+    if (a <= epsilon) return false;
+    if (b <= epsilon) return true;
+    return a < b;
+  });
   float maxVal = *std::max_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH);
   std::cout << "minVal: " << minVal << " maxVal: " << maxVal << std::endl;
 
+  // Min-max normalization TODO: Decide whether to keep the normalization here but probably not
   // Normalize to [0, 1]
   // Temperature: min: 0 max: 1 avg: 0.776319 median: 0.790567
   // Speed: min: 0 max: 1 avg: 0.132117 median: 0.0837869
-  for (int i = 0; i < VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH; i++) {
+  // for (int i = 0; i < VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH; i++) {
-    hostVolume[i] = (hostVolume[i] - minVal) / (maxVal - minVal);
+  //   hostVolume[i] = (hostVolume[i] - minVal) / (maxVal - minVal);
-  }
+  // }
 
   // // print min, max, avg., and median values <--- the code actually does not work when this snippet is enabled so probably TODO: Delete this later
   // std::sort(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH);