From 7c046a099a564a43af12b2882e542b8c70ca9596 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 23:42:33 +0100 Subject: [PATCH] Playing around with values to actually get something not fully opaque --- src/consts.cu | 5 +- src/consts.h | 19 +++-- src/gui/MainWindow.cpp | 25 +++++- src/illumination/Raycaster.cu | 146 +++++++++++++++++++++++----------- src/main.cu | 26 +++--- 5 files changed, 153 insertions(+), 68 deletions(-) diff --git a/src/consts.cu b/src/consts.cu index 5fa57f3..c55069f 100644 --- a/src/consts.cu +++ b/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() { diff --git a/src/consts.h b/src/consts.h index adf961b..fb4a318 100644 --- a/src/consts.h +++ b/src/consts.h @@ -5,9 +5,12 @@ #include // --------------------------- Basic Constants --------------------------- -const int VOLUME_WIDTH = 576; -const int VOLUME_HEIGHT = 361; -const int VOLUME_DEPTH = 42; +// 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_WIDTH = 576; // lon +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 --------------------------- diff --git a/src/gui/MainWindow.cpp b/src/gui/MainWindow.cpp index a487328..3ab2dc6 100644 --- a/src/gui/MainWindow.cpp +++ b/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)) { - Window::tick(); - } + // 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; diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index df30958..e69d753 100644 --- a/src/illumination/Raycaster.cu +++ b/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__ Color3 temperatureToRGB(float temperature) { - // atm, the scalar field is normalized - const float minTemp = 0.0f; // coldest == deep blue - const float maxTemp = 1.0f; // hottest temperature == deep red - - 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); +__device__ float opacityFromGradient(const Vec3 &grad) { + float gradMag = grad.length(); // magnitude + float k = 1e-6f; // tweak (the smaller the value, the less opacity) // TODO: What should be the value of this? + float alpha = 1.0f - expf(-k * gradMag); + return alpha; } +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) { - float4 result; +__device__ float4 transferFunction(float density, const Vec3& grad, const Point3& pos, const Vec3& rayDir) { // 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; - result.y = (1.0f - alphaSample) * shadedColor.y * alphaSample; - result.z = (1.0f - alphaSample) * shadedColor.z * alphaSample; - result.w = (1.0f - alphaSample) * alphaSample; + float4 result; + result.x = shadedColor.x * alphaSample; + result.y = shadedColor.y * 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) - if (grad.dot(viewDir) < epsilon / 100000.0f) { + // 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.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 - accumR = 0.9f; - accumG = 0.9f; - accumB = 0.9f; + // No intersection -> Set to brackground color (multiply by SAMPLES_PER_PIXEL because we divide by it later) + accumR = 0.1f * (float)SAMPLES_PER_PIXEL; + accumG = 0.1f * (float)SAMPLES_PER_PIXEL; + 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; - colorG += color.y; - colorB += color.z; - alphaAccum += color.w; + colorR += color.x * (alphaAcumLimit - alphaAccum); + colorG += color.y * (alphaAcumLimit - alphaAccum); + colorB += color.z * (alphaAcumLimit - alphaAccum); + 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; - // accumR = accumR + leftover * 0.9f; - // accumG = accumG + leftover * 0.9f; - // accumB = accumB + leftover * 0.9f; + float leftover = 1.0 - alphaAccum; + accumR = accumR + leftover * 0.1f; + accumG = accumG + leftover * 0.1f; + accumB = accumB + leftover * 0.1f; } } + // Average samples accumR /= (float)SAMPLES_PER_PIXEL; accumG /= (float)SAMPLES_PER_PIXEL; diff --git a/src/main.cu b/src/main.cu index 2c14d13..43194f5 100644 --- a/src/main.cu +++ b/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& temperatureData, int idx = 0) { - // std::string path = "data/trimmed"; - std::string path = "data"; + std::string path = "data/trimmed"; + // 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& temperatureData, int idx = 0) { } void getSpeed(std::vector& speedData, int idx = 0) { - // std::string path = "data/trimmed"; - std::string path = "data"; + std::string path = "data/trimmed"; + // 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 minVal = *std::min_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH, [](float a, float b) { + 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++) { - hostVolume[i] = (hostVolume[i] - minVal) / (maxVal - minVal); - } + // for (int i = 0; i < VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH; i++) { + // 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);