From d3b88aed302c7174740718b4adcdffd6bc45a33f Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Wed, 8 Jan 2025 22:43:30 +0100 Subject: [PATCH 01/25] Forgot to do atomic commits so ... booyah --- src/consts.cu | 5 ++-- src/consts.h | 4 +-- src/illumination/Raycaster.cu | 46 +++++++++++++++++++---------------- src/main.cu | 17 ++++++++----- 4 files changed, 41 insertions(+), 31 deletions(-) diff --git a/src/consts.cu b/src/consts.cu index b432f6d..27803b4 100644 --- a/src/consts.cu +++ b/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); -Vec3 h_cameraDir = Vec3::init(0.4, 0.6, 1.0).normalize(); +Point3 h_cameraPos = Point3::init(-500.0f, 200.0f, 100.0f); +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); diff --git a/src/consts.h b/src/consts.h index 0c1e8d4..1c6fbd2 100644 --- a/src/consts.h +++ b/src/consts.h @@ -5,8 +5,8 @@ #include // --------------------------- Basic Constants --------------------------- -const int VOLUME_WIDTH = 49; -const int VOLUME_HEIGHT = 51; +const int VOLUME_WIDTH = 576; +const int VOLUME_HEIGHT = 361; const int VOLUME_DEPTH = 42; const int IMAGE_WIDTH = 800; diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index a405e25..5b97e0d 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -40,28 +40,29 @@ __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) { - if (fabsf(dirVal) < epsilon) { - if (start < 0.f || start > 1.f) { - tNear = 1e9f; + auto intersectAxis = [&](float start, float dir, float minV, float maxV) { + if (fabsf(dir) < epsilon) { + // 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 t1 = (1.0f - start) / dirVal; - if (t0>t1) { - float tmp=t0; - t0=t1; - t1=tmp; + float t0 = (minV - start) / dir; + float t1 = (maxV - start) / dir; + if (t0 > t1) { + float tmp = t0; + t0 = t1; + t1 = tmp; } - if (t0>tNear) tNear = t0; - if (t1 tNear) tNear = t0; + if (t1 < tFar ) tFar = t1; } }; - intersectAxis(d_cameraPos.x, rayDir.x); - intersectAxis(d_cameraPos.y, rayDir.y); - intersectAxis(d_cameraPos.z, rayDir.z); + intersectAxis(d_cameraPos.x, rayDir.x, 0.0f, (float)VOLUME_WIDTH); + intersectAxis(d_cameraPos.y, rayDir.y, 0.0f, (float)VOLUME_HEIGHT); + 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 fy = pos.y * (VOLUME_HEIGHT - 1); - float fz = pos.z * (VOLUME_DEPTH - 1); - int ix = (int)roundf(fx); - int iy = (int)roundf(fy); - int iz = (int)roundf(fz); + // float fx = pos.x * (VOLUME_WIDTH - 1); + // float fy = pos.y * (VOLUME_HEIGHT - 1); + // float fz = pos.z * (VOLUME_DEPTH - 1); + int ix = (int)roundf(pos.x); + int iy = (int)roundf(pos.y); + 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; } diff --git a/src/main.cu b/src/main.cu index 10a3dfc..293cd9b 100644 --- a/src/main.cu +++ b/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& 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& temperatureData, int idx = 0) { } void getSpeed(std::vector& 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& speedData, int idx = 0) { int main() { std::vector data; - // getTemperature(data); - getSpeed(data); + getTemperature(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]; - if (data[i] + epsilon >= infty) hostVolume[i] = 0.0f; + 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; } // 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); } From c5b8b6cff0f384abfbff2237e3d11dd24eb774d2 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Wed, 8 Jan 2025 22:57:06 +0100 Subject: [PATCH 02/25] Position the camera better so that the scaled volume is actually visible --- src/consts.cu | 2 +- src/illumination/Raycaster.cu | 13 +++++++++++++ src/main.cu | 4 ++-- src/objs/sphere.h | 13 ------------- 4 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/consts.cu b/src/consts.cu index 27803b4..5fa57f3 100644 --- a/src/consts.cu +++ b/src/consts.cu @@ -5,7 +5,7 @@ __device__ Vec3 d_cameraDir; __device__ Vec3 d_cameraUp; __device__ Point3 d_lightPos; -Point3 h_cameraPos = Point3::init(-500.0f, 200.0f, 100.0f); +Point3 h_cameraPos = Point3::init(-300.0f, 200.0f, -300.0f); 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(); diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index 5b97e0d..705d882 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -10,6 +10,19 @@ #include #include "objs/sphere.h" +// Samples the voxel nearest to the given coordinates. TODO: Can be re-used in other places so move +__device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz) { + if (vx < 0) vx = 0; + if (vy < 0) vy = 0; + if (vz < 0) vz = 0; + if (vx >= volW) vx = volW - 1; + if (vy >= volH) vy = volH - 1; + if (vz >= volD) vz = volD - 1; + + int idx = vz * volW * volH + vy * volD + vx; + return volumeData[idx]; +} + // TODO: instead of IMAGEWIDTH and IMAGEHEIGHT this should reflect the windowSize; __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { diff --git a/src/main.cu b/src/main.cu index 293cd9b..968be31 100644 --- a/src/main.cu +++ b/src/main.cu @@ -54,8 +54,8 @@ void getSpeed(std::vector& speedData, int idx = 0) { int main() { std::vector data; - getTemperature(data); - // getSpeed(data); + // getTemperature(data); + getSpeed(data); std::cout << "DATA size: " << data.size() << std::endl; diff --git a/src/objs/sphere.h b/src/objs/sphere.h index dbac251..27e042d 100644 --- a/src/objs/sphere.h +++ b/src/objs/sphere.h @@ -78,17 +78,4 @@ __host__ void generateVolume(float* volumeData, int volW, int volH, int volD) { } } } -} - -// Samples the voxel nearest to the given coordinates. TODO: Can be re-used in other places so move -__device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz) { - if (vx < 0) vx = 0; - if (vy < 0) vy = 0; - if (vz < 0) vz = 0; - if (vx >= volW) vx = volW - 1; - if (vy >= volH) vy = volH - 1; - if (vz >= volD) vz = volD - 1; - - int idx = vz * volW * volH + vy * volD + vx; - return volumeData[idx]; } \ No newline at end of file From cc8ef125089a2928497aec1936a697424a996cda Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 11:43:58 +0100 Subject: [PATCH 03/25] Started doing basic silhoutte shading + changed background to gray --- src/illumination/Raycaster.cu | 137 +++++++++++++++++++++++----------- 1 file changed, 95 insertions(+), 42 deletions(-) diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index 705d882..d119dbc 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -23,6 +23,62 @@ __device__ float sampleVolumeNearest(float* volumeData, const int volW, const in return volumeData[idx]; } +// tri-linear interpolation - ready if necessary (but no visible improvement for full volume) +__device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz) { + int ix = (int)fx; + int iy = (int)fy; + int iz = (int)fz; + + float dx = fx - ix; + float dy = fy - iy; + float dz = fz - iz; + + float c00 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz) * (1.0f - dx) + sampleVolumeNearest(volumeData, volW, volH, volD, ix + 1, iy, iz) * dx; + float c10 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy + 1, iz) * (1.0f - dx) + sampleVolumeNearest(volumeData, volW, volH, volD, ix + 1, iy + 1, iz) * dx; + float c01 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz + 1) * (1.0f - dx) + sampleVolumeNearest(volumeData, volW, volH, volD, ix + 1, iy, iz + 1) * dx; + float c11 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy + 1, iz + 1) * (1.0f - dx) + sampleVolumeNearest(volumeData, volW, volH, volD, ix + 1, iy + 1, iz + 1) * dx; + + float c0 = c00 * (1.0f - dy) + c10 * dy; + float c1 = c01 * (1.0f - dy) + c11 * dy; + + return c0 * (1.0f - dz) + c1 * dz; +} + + +// Transfer function +__device__ float4 transferFunction(float density, Vec3 grad, Point3 pos, 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 + + Vec3 normal = -grad.normalize(); + + Vec3 lightDir = (d_lightPos - pos).normalize(); + Vec3 viewDir = -rayDir.normalize(); + + // Apply Phong + 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; + + // 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) { + result.x = 0.0f; + result.y = 0.0f; + result.z = 0.0f; + result.w = 1.0f; + } + + return result; +} + // TODO: instead of IMAGEWIDTH and IMAGEHEIGHT this should reflect the windowSize; __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { @@ -77,59 +133,56 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { intersectAxis(d_cameraPos.y, rayDir.y, 0.0f, (float)VOLUME_HEIGHT); intersectAxis(d_cameraPos.z, rayDir.z, 0.0f, (float)VOLUME_DEPTH); - if (tNear > tFar) continue; // No intersectionn - if (tNear < 0.0f) tNear = 0.0f; + if (tNear > tFar){ + // No intersectionn + accumR = 0.9f; + accumG = 0.9f; + accumB = 0.9f; + } else { + if (tNear < 0.0f) tNear = 0.0f; - float colorR = 0.0f, colorG = 0.0f, colorB = 0.0f; - float alphaAccum = 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 = d_cameraPos + rayDir * tCurrent; + float tCurrent = tNear; + while (tCurrent < tFar && alphaAccum < alphaAcumLimit) { + Point3 pos = d_cameraPos + rayDir * tCurrent; - // Convert to volume indices - // float fx = pos.x * (VOLUME_WIDTH - 1); - // float fy = pos.y * (VOLUME_HEIGHT - 1); - // float fz = pos.z * (VOLUME_DEPTH - 1); - int ix = (int)roundf(pos.x); - int iy = (int)roundf(pos.y); - int iz = (int)roundf(pos.z); + // Convert to volume indices + // float fx = pos.x * (VOLUME_WIDTH - 1); + // float fy = pos.y * (VOLUME_HEIGHT - 1); + // float fz = pos.z * (VOLUME_DEPTH - 1); + int ix = (int)roundf(pos.x); + int iy = (int)roundf(pos.y); + int iz = (int)roundf(pos.z); - // Sample - float density = sampleVolumeNearest(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz); + // Sample (pick appropriate method based on volume size) + float density = sampleVolumeNearest(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz); + // float density = sampleVolumeTrilinear(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, pos.x, pos.y, pos.z); - // 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::init(density, 0.1f*density, 1.f - density); // TODO: Implement a proper transfer function - - // If density ~ 0, skip shading - if (density > minAllowedDensity) { + // If density ~ 0, skip shading + if (density > minAllowedDensity) { Vec3 grad = computeGradient(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz); - Vec3 normal = -grad.normalize(); + float4 color = transferFunction(density, grad, pos, rayDir); + colorR += color.x; + colorG += color.y; + colorB += color.z; + alphaAccum += color.w; + } - Vec3 lightDir = (d_lightPos - pos).normalize(); - Vec3 viewDir = -rayDir.normalize(); - // Apply Phong - Vec3 shadedColor = phongShading(normal, lightDir, viewDir, baseColor); + tCurrent += stepSize; + } - // 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; - } + accumR += colorR; + accumG += colorG; + accumB += colorB; - // TODO: Add silhouette - Take gradient of volume dot with view direction (if small then this is a silhouette) - // TODO: Scale volume correctly - - tCurrent += stepSize; + // float leftover = 1.0f - alphaAccum; + // accumR = accumR + leftover * 0.9f; + // accumG = accumG + leftover * 0.9f; + // accumB = accumB + leftover * 0.9f; } - - accumR += colorR; - accumG += colorG; - accumB += colorB; } // Average samples From b107a2a6def926b477f683866aed0cac25a142dc Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 11:53:28 +0100 Subject: [PATCH 04/25] Added some more useful math functions --- src/linalg/mat.cu | 16 ++++++++++++++++ src/linalg/mat.h | 7 +++++++ 2 files changed, 23 insertions(+) diff --git a/src/linalg/mat.cu b/src/linalg/mat.cu index c5bdd4a..c42a8b3 100644 --- a/src/linalg/mat.cu +++ b/src/linalg/mat.cu @@ -44,3 +44,19 @@ __device__ unsigned int packUnorm4x8(float r, float g, float b, float a) { return u.out; } + +// Clamp a value between a min and max value +__device__ float clamp(float value, float min, float max) { + return fmaxf(min, fminf(value, max)); +} + +// Normalize a float to the range [0, 1] +__device__ float normalize(float value, float min, float max) { + return (value - min) / (max - min); +} + +// Interpolate between two values +template +__device__ T interpolate(T start, T end, float t) { + return start + t * (end - start); +} diff --git a/src/linalg/mat.h b/src/linalg/mat.h index cd070d7..559b8c5 100644 --- a/src/linalg/mat.h +++ b/src/linalg/mat.h @@ -7,4 +7,11 @@ __device__ Vec3 computeGradient(float* volumeData, const int volW, const int vol __device__ unsigned int packUnorm4x8(float r, float g, float b, float a); +__device__ float clamp(float value, float min, float max); +__device__ float normalize(float value, float min, float max); + +template +__device__ float interpolate(T start, T end, float t); + + #endif // MAT_H From 7cf829a20a8ad54a5ba17944397a67b785e0d36e Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 11:54:46 +0100 Subject: [PATCH 05/25] Used new normalize functions since we already have it --- src/main.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.cu b/src/main.cu index 968be31..41848b7 100644 --- a/src/main.cu +++ b/src/main.cu @@ -75,7 +75,7 @@ int main() { 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); + hostVolume[i] = normalize(hostVolume[i], minVal, maxVal); } // Allocate + copy data to GPU From ec31f3f897169378626c33411935b72dd4550630 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 12:00:14 +0100 Subject: [PATCH 06/25] Nevermind, main is not a host function so the new normalize function does not work there --- src/main.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.cu b/src/main.cu index 41848b7..968be31 100644 --- a/src/main.cu +++ b/src/main.cu @@ -75,7 +75,7 @@ int main() { std::cout << "minVal: " << minVal << " maxVal: " << maxVal << std::endl; for (int i = 0; i < VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH; i++) { - hostVolume[i] = normalize(hostVolume[i], minVal, maxVal); + hostVolume[i] = (hostVolume[i] - minVal) / (maxVal - minVal); } // Allocate + copy data to GPU From 52a47ed302d636bf875e5038438f53d19a268ae4 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 12:16:18 +0100 Subject: [PATCH 07/25] Fixed the newly defined math function to actually work with cuda --- src/linalg/mat.cu | 8 +------- src/linalg/mat.h | 6 +++++- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/src/linalg/mat.cu b/src/linalg/mat.cu index c42a8b3..c206cb1 100644 --- a/src/linalg/mat.cu +++ b/src/linalg/mat.cu @@ -53,10 +53,4 @@ __device__ float clamp(float value, float min, float max) { // Normalize a float to the range [0, 1] __device__ float normalize(float value, float min, float max) { return (value - min) / (max - min); -} - -// Interpolate between two values -template -__device__ T interpolate(T start, T end, float t) { - return start + t * (end - start); -} +} \ No newline at end of file diff --git a/src/linalg/mat.h b/src/linalg/mat.h index 559b8c5..4d2a1d5 100644 --- a/src/linalg/mat.h +++ b/src/linalg/mat.h @@ -10,8 +10,12 @@ __device__ unsigned int packUnorm4x8(float r, float g, float b, float a); __device__ float clamp(float value, float min, float max); __device__ float normalize(float value, float min, float max); +// Interpolate between two values template -__device__ float interpolate(T start, T end, float t); +__device__ T interpolate(T start, T end, float t) { + return start + t * (end - start); +} + #endif // MAT_H From 10324eb1ca1127411919258b369f51405f5942db Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 12:16:55 +0100 Subject: [PATCH 08/25] Made image size bigger cause I am blind --- src/consts.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/consts.h b/src/consts.h index 1c6fbd2..adf961b 100644 --- a/src/consts.h +++ b/src/consts.h @@ -9,8 +9,8 @@ const int VOLUME_WIDTH = 576; const int VOLUME_HEIGHT = 361; const int VOLUME_DEPTH = 42; -const int IMAGE_WIDTH = 800; -const int IMAGE_HEIGHT = 600; +const int IMAGE_WIDTH = 1600; +const int IMAGE_HEIGHT = 1200; const double epsilon = 1e-10f; const double infty = 1e15f; // This value is used to represent missing values in data From d978630f26fe9b2db11eb8768ec6c537916a69e1 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 12:17:16 +0100 Subject: [PATCH 09/25] Analyzed data for transfer-function purposes --- src/main.cu | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/main.cu b/src/main.cu index 968be31..2c14d13 100644 --- a/src/main.cu +++ b/src/main.cu @@ -10,6 +10,7 @@ #include #include "linalg/linalg.h" #include +#include static float* d_volume = nullptr; @@ -54,8 +55,8 @@ void getSpeed(std::vector& speedData, int idx = 0) { int main() { std::vector data; - // getTemperature(data); - getSpeed(data); + getTemperature(data); + // getSpeed(data); std::cout << "DATA size: " << data.size() << std::endl; @@ -74,10 +75,19 @@ int main() { float maxVal = *std::max_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH); std::cout << "minVal: " << minVal << " maxVal: " << maxVal << std::endl; + // 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); } + // // 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); + // float sum = std::accumulate(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH, 0.0f); + // float avg = sum / (VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH); + // std::cout << "min: " << hostVolume[0] << " max: " << hostVolume[VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH - 1] << " avg: " << avg << " median: " << hostVolume[VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH / 2] << std::endl; + // Allocate + copy data to GPU size_t volumeSize = sizeof(float) * VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH; cudaMalloc((void**)&d_volume, volumeSize); From fc2488c5780bcc3597253b7445c16eaff0cd5045 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 12:17:30 +0100 Subject: [PATCH 10/25] Improved transfer function business --- src/illumination/Raycaster.cu | 38 +++++++++++++++++++++++++++++------ 1 file changed, 32 insertions(+), 6 deletions(-) diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index d119dbc..df30958 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -10,6 +10,7 @@ #include #include "objs/sphere.h" +// TODO: Probbably move this transfer function business into a different file // Samples the voxel nearest to the given coordinates. TODO: Can be re-used in other places so move __device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz) { if (vx < 0) vx = 0; @@ -44,6 +45,32 @@ __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); +} + // Transfer function __device__ float4 transferFunction(float density, Vec3 grad, Point3 pos, Vec3 rayDir) { @@ -52,7 +79,8 @@ __device__ float4 transferFunction(float density, Vec3 grad, Point3 pos, Vec3 ra 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 = Color3::init(density, 0.1f*density, 1.f - density); // TODO: Implement a proper transfer function + Color3 baseColor = temperatureToRGB(density); Vec3 normal = -grad.normalize(); @@ -80,6 +108,7 @@ __device__ float4 transferFunction(float density, Vec3 grad, Point3 pos, Vec3 ra } + // TODO: instead of IMAGEWIDTH and IMAGEHEIGHT this should reflect the windowSize; __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { int px = blockIdx.x * blockDim.x + threadIdx.x; @@ -149,16 +178,13 @@ __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 fy = pos.y * (VOLUME_HEIGHT - 1); - // float fz = pos.z * (VOLUME_DEPTH - 1); int ix = (int)roundf(pos.x); int iy = (int)roundf(pos.y); int iz = (int)roundf(pos.z); // Sample (pick appropriate method based on volume size) - float density = sampleVolumeNearest(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz); - // float density = sampleVolumeTrilinear(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, pos.x, pos.y, pos.z); + // float density = sampleVolumeNearest(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz); + float density = sampleVolumeTrilinear(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, pos.x, pos.y, pos.z); // If density ~ 0, skip shading if (density > minAllowedDensity) { From 1f335dffc9a08edebacd42a28fb4b8fe994e8fcb Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 20:36:14 +0100 Subject: [PATCH 11/25] Added support for float multiplication of vec3 --- src/linalg/vec.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/linalg/vec.h b/src/linalg/vec.h index f4bf332..06a1b38 100644 --- a/src/linalg/vec.h +++ b/src/linalg/vec.h @@ -20,6 +20,12 @@ struct Vec3 { // TODO: Maybe make this into a class ... maybe __host__ __device__ Vec3 operator*(double b) const { return Vec3::init(x * b, y * b, z * b); } __host__ __device__ Vec3& operator*=(double b) { x *= b; y *= b; z *= b; return *this; } + __host__ __device__ Vec3 operator*(float b) const { return Vec3::init(x * b, y * b, z * b); } + __host__ __device__ Vec3& operator*=(float b) { x *= b; y *= b; z *= b; return *this; } + friend __host__ __device__ Vec3 operator*(float a, const Vec3& b) { + return Vec3::init(a * b.x, a * b.y, a * b.z); + } + __host__ __device__ double dot(const Vec3& b) const { return x * b.x + y * b.y + z * b.z; } __host__ __device__ Vec3 cross(const Vec3& b) const { return Vec3::init(y * b.z - z * b.y, z * b.x - x * b.z, x * b.y - y * b.x); } From 09220ced2586acbd6bb38fd1bf68b25a55f440f0 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 20:36:43 +0100 Subject: [PATCH 12/25] Added func. for length of vec3 --- src/linalg/vec.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/linalg/vec.h b/src/linalg/vec.h index 06a1b38..281709a 100644 --- a/src/linalg/vec.h +++ b/src/linalg/vec.h @@ -30,6 +30,7 @@ struct Vec3 { // TODO: Maybe make this into a class ... maybe __host__ __device__ double dot(const Vec3& b) const { return x * b.x + y * b.y + z * b.z; } __host__ __device__ Vec3 cross(const Vec3& b) const { return Vec3::init(y * b.z - z * b.y, z * b.x - x * b.z, x * b.y - y * b.x); } __host__ __device__ Vec3 normalize() const { double len = sqrt(x * x + y * y + z * z); return Vec3::init(x / len, y / len, z / len); } + __host__ __device__ double length() const { return sqrt(x * x + y * y + z * z); } }; typedef Vec3 Point3; From 7c046a099a564a43af12b2882e542b8c70ca9596 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Thu, 9 Jan 2025 23:42:33 +0100 Subject: [PATCH 13/25] 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); From 6f2ed9e4df0040b99b6974841d949b2ee59ad634 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Fri, 10 Jan 2025 21:02:18 +0100 Subject: [PATCH 14/25] It finally looks about right but the code urgently needs to be tidied up --- src/consts.h | 9 +- src/illumination/FrameBuffer.cu | 4 +- src/illumination/FrameBuffer.h | 2 +- src/illumination/Raycaster.cu | 249 +++++++++++++++++++++++++------- src/main.cu | 33 ++++- 5 files changed, 231 insertions(+), 66 deletions(-) diff --git a/src/consts.h b/src/consts.h index fb4a318..dd2e4f6 100644 --- a/src/consts.h +++ b/src/consts.h @@ -19,14 +19,17 @@ 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; +const float MIN_TEMP = 210.0f; +const float MAX_TEMP = 240.0f; + +const float MIN_SPEED = 0.0F; +const float MAX_SPEED = 14.0f; // --------------------------- Raycasting Constants --------------------------- 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.4f; // TODO: Idk what a good accumulation value is <--- This finally does something when using alpha in both places at least +const float alphaAcumLimit = 1.0f; // 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.02f; diff --git a/src/illumination/FrameBuffer.cu b/src/illumination/FrameBuffer.cu index df4fe99..becfb86 100644 --- a/src/illumination/FrameBuffer.cu +++ b/src/illumination/FrameBuffer.cu @@ -4,9 +4,9 @@ __host__ FrameBuffer::FrameBuffer(unsigned int w, unsigned int h) : w(w), h(h) {} -__device__ void FrameBuffer::writePixel(int x, int y, float r, float g, float b) { +__device__ void FrameBuffer::writePixel(int x, int y, float r, float g, float b, float a) { int i = y * this->w + x; // the opengl buffer uses BGRA format; dunno why - this->buffer[i] = packUnorm4x8(b, g, r, 1.0f); + this->buffer[i] = packUnorm4x8(b, g, r, a); } diff --git a/src/illumination/FrameBuffer.h b/src/illumination/FrameBuffer.h index c7e83ad..fde7b1f 100644 --- a/src/illumination/FrameBuffer.h +++ b/src/illumination/FrameBuffer.h @@ -14,7 +14,7 @@ public: unsigned int h; __host__ FrameBuffer(unsigned int w, unsigned int h); - __device__ void writePixel(int x, int y, float r, float g, float b); + __device__ void writePixel(int x, int y, float r, float g, float b, float a); }; #endif // FRAMEBUFFER_H diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index e69d753..8376a5b 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -9,35 +9,110 @@ #include "shading.h" #include #include "objs/sphere.h" +#include // TODO: Probbably move this transfer function business into a different file // Samples the voxel nearest to the given coordinates. TODO: Can be re-used in other places so move __device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz) { + // TODO: I actually do not know which screen coordinates corresponds to which volume coordinates ... so let's try all options + + // x <-> width, y <-> height, z <-> depth + // if (vx < 0) vx = 0; + // if (vy < 0) vy = 0; + // if (vz < 0) vz = 0; + // if (vx >= volW) vx = volW - 1; + // if (vy >= volH) vy = volH - 1; + // if (vz >= volD) vz = volD - 1; + + // int idx = vz * volW * volH + vy * volD + vx; + // return volumeData[idx]; + + // x <-> width, y <-> depth, z <-> height + // if (vx < 0) vx = 0; + // if (vy < 0) vy = 0; + // if (vz < 0) vz = 0; + // if (vx >= volW) vx = volW - 1; + // if (vy >= volH) vy = volH - 1; + // if (vz >= volD) vz = volD - 1; + + // int idx = vy * volW * volD + vz * volW + vx; + // return volumeData[idx]; + + // x <-> height, y <-> width, z <-> depth <--- So far this is the best one if (vx < 0) vx = 0; if (vy < 0) vy = 0; if (vz < 0) vz = 0; - if (vx >= volW) vx = volW - 1; - if (vy >= volH) vy = volH - 1; + if (vx >= volH) vx = volH - 1; + if (vy >= volW) vy = volW - 1; if (vz >= volD) vz = volD - 1; - int idx = vz * volW * volH + vy * volD + vx; + int idx = vz * volW * volH + vx * volW + vy; return volumeData[idx]; + + // x <-> height, y <-> depth, z <-> width + // if (vx < 0) vx = 0; + // if (vy < 0) vy = 0; + // if (vz < 0) vz = 0; + // if (vx >= volH) vx = volH - 1; + // if (vy >= volD) vy = volD - 1; + // if (vz >= volW) vz = volW - 1; + + // int idx = vy * volW * volH + vz * volW + vx; + // return volumeData[idx]; + + + // x <-> depth, y <-> width, z <-> height + // if (vx < 0) vx = 0; + // if (vy < 0) vy = 0; + // if (vz < 0) vz = 0; + // if (vx >= volD) vx = volD - 1; + // if (vy >= volW) vy = volW - 1; + // if (vz >= volH) vz = volH - 1; + + // int idx = vz * volD * volW + vx * volW + vy; + // return volumeData[idx]; + + // x <-> depth, y <-> height, z <-> width + // if (vx < 0) vx = 0; + // if (vy < 0) vy = 0; + // if (vz < 0) vz = 0; + // if (vx >= volD) vx = volD - 1; + // if (vy >= volH) vy = volH - 1; + // if (vz >= volW) vz = volW - 1; + + // int idx = vy * volD * volW + vz * volD + vx; + // return volumeData[idx]; + } // tri-linear interpolation - ready if necessary (but no visible improvement for full volume) __device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz) { - int ix = (int)fx; - int iy = (int)fy; - int iz = (int)fz; + int ix = (int)floorf(fx); + int iy = (int)floorf(fy); + int iz = (int)floorf(fz); + // Clamp indices to valid range + int ix1 = min(ix + 1, volW - 1); + int iy1 = min(iy + 1, volH - 1); + int iz1 = min(iz + 1, volD - 1); + ix = max(ix, 0); + iy = max(iy, 0); + iz = max(iz, 0); + + // Compute weights float dx = fx - ix; float dy = fy - iy; float dz = fz - iz; - float c00 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz) * (1.0f - dx) + sampleVolumeNearest(volumeData, volW, volH, volD, ix + 1, iy, iz) * dx; - float c10 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy + 1, iz) * (1.0f - dx) + sampleVolumeNearest(volumeData, volW, volH, volD, ix + 1, iy + 1, iz) * dx; - float c01 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz + 1) * (1.0f - dx) + sampleVolumeNearest(volumeData, volW, volH, volD, ix + 1, iy, iz + 1) * dx; - float c11 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy + 1, iz + 1) * (1.0f - dx) + sampleVolumeNearest(volumeData, volW, volH, volD, ix + 1, iy + 1, iz + 1) * dx; + // Sample values + float c00 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz) * (1.0f - dx) + + sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy, iz) * dx; + float c10 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy1, iz) * (1.0f - dx) + + sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy1, iz) * dx; + float c01 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz1) * (1.0f - dx) + + sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy, iz1) * dx; + float c11 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy1, iz1) * (1.0f - dx) + + sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy1, iz1) * dx; float c0 = c00 * (1.0f - dy) + c10 * dy; float c1 = c01 * (1.0f - dy) + c11 * dy; @@ -45,9 +120,10 @@ __device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const return c0 * (1.0f - dz) + c1 * dz; } + __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 k = 1e-4f; // 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; } @@ -88,35 +164,65 @@ __device__ Color3 colorMapViridis(float normalizedT) { 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 +// Monochromatic colormap for speed +__device__ Color3 colorMapMonochrome(float normalizedSpeed) { + // Define the color stops: black (0.0) to white (1.0) + ColorStop speedStops[] = { + { 0.0f, Color3::init(0.0f, 0.0f, 0.0f) }, // No colour + { 1.0f, Color3::init(1.0f, 1.0f, 1.0f) } // White + }; -// if (temperature < minTemp) { -// return Color3::init(1.f, 1.f, 1.f); -// } -// temperature = clamp(temperature, minTemp, maxTemp); -// float t = normalize(temperature, minTemp, maxTemp); + // Clamp to [0,1] + if (normalizedSpeed < 0.0f) normalizedSpeed = 0.0f; + if (normalizedSpeed > 1.0f) normalizedSpeed = 1.0f; -// float r, g, b; + // Single interval (N = 2 for black to white) + const int N = 2; + for (int i = 0; i < N - 1; ++i) + { + float start = speedStops[i].pos; + float end = speedStops[i + 1].pos; -// 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); -// } + if (normalizedSpeed >= start && normalizedSpeed <= end) + { + float localT = (normalizedSpeed - start) / (end - start); + return interpolate(speedStops[i].color, speedStops[i + 1].color, localT); + } + } -// return Color3::init(r, g, b); -// } + // Fallback + return speedStops[N - 1].color; +} + +// Colormap function for gradient from blue to white to red +__device__ Color3 colorMapPythonLike(float normalizedValue) { + // Define control points for the colormap (approximation of coolwarm) + ColorStop coolwarmStops[] = { + { 0.0f, Color3::init(0.2298057f, 0.29871797f, 0.75368315f) }, // Dark Blue + { 0.25f, Color3::init(0.23437708f, 0.30554173f, 0.75967953f) }, // Mid Blue + { 0.5f, Color3::init(0.27582712f, 0.36671692f, 0.81255294f) }, // White + { 0.75f, Color3::init(0.79606387f, 0.84869321f, 0.93347147f) }, // Light Orange + { 1.0f, Color3::init(0.70567316f, 0.01555616f, 0.15023281f) } // Red + }; + + // Clamp the scalar value to [0, 1] + normalizedValue = fminf(fmaxf(normalizedValue, 0.0f), 1.0f); + + // Interpolate between the defined color stops + const int N = 5; // Number of control points + for (int i = 0; i < N - 1; ++i) { + float start = coolwarmStops[i].pos; + float end = coolwarmStops[i + 1].pos; + + if (normalizedValue >= start && normalizedValue <= end) { + float localT = (normalizedValue - start) / (end - start); + return interpolate(coolwarmStops[i].color, coolwarmStops[i + 1].color, localT); + } + } + + // Fallback (shouldn't reach here due to clamping) + return coolwarmStops[N - 1].color; +} // Transfer function @@ -127,9 +233,18 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3 // Color3 baseColor = temperatureToRGB(density); float normDensity = (density - MIN_TEMP) / (MAX_TEMP - MIN_TEMP); + // float normDensity = (density - MIN_SPEED) / (MAX_SPEED - MIN_SPEED); + normDensity = clamp(normDensity, 0.0f, 1.0f); - Color3 baseColor = colorMapViridis(normDensity); + + // Color3 baseColor = colorMapViridis(normDensity); + // Color3 baseColor = colorMapMonochrome(normDensity); + Color3 baseColor = colorMapPythonLike(normDensity); + + float alpha = opacityFromGradient(grad); + // alpha = 0.1f; + alpha = 1.0f / (1.0f + expf(-250.f * (normDensity - 0.5f))); // This is also quite nice, but the exponent should be parameterized float alphaSample = density * alpha; // TODO: Decide whether to keep alpha here or not Vec3 normal = -grad.normalize(); @@ -169,12 +284,19 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { float accumR = 0.0f; float accumG = 0.0f; float accumB = 0.0f; + float accumA = 1.0f; + // Initialize random state + curandState randState; + curand_init(1234, px + py * IMAGE_WIDTH, 0, &randState); + // Multiple samples per pixel // 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; + float jitterU = (curand_uniform(&randState) - 0.5f) / IMAGE_WIDTH; + float jitterV = (curand_uniform(&randState) - 0.5f) / IMAGE_HEIGHT; + float u = ((px + 0.5f + jitterU) / IMAGE_WIDTH ) * 2.0f - 1.0f; + float v = ((py + 0.5f + jitterV) / 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); @@ -186,7 +308,7 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { d_cameraUp = (cameraRight.cross(d_cameraDir)).normalize(); Vec3 rayDir = (d_cameraDir + cameraRight*u + d_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 + // Intersect float tNear = 0.0f; float tFar = 1e6f; auto intersectAxis = [&](float start, float dir, float minV, float maxV) { @@ -209,24 +331,25 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { } }; - intersectAxis(d_cameraPos.x, rayDir.x, 0.0f, (float)VOLUME_WIDTH); - intersectAxis(d_cameraPos.y, rayDir.y, 0.0f, (float)VOLUME_HEIGHT); + intersectAxis(d_cameraPos.x, rayDir.x, 0.0f, (float)VOLUME_HEIGHT); + intersectAxis(d_cameraPos.y, rayDir.y, 0.0f, (float)VOLUME_WIDTH); intersectAxis(d_cameraPos.z, rayDir.z, 0.0f, (float)VOLUME_DEPTH); - if (tNear > tFar){ + if (tNear > tFar) { // 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; + accumA = 1.0f * (float)SAMPLES_PER_PIXEL; } else { 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 = d_cameraPos + rayDir * tCurrent; + float t = tNear; + while (t < tFar && alphaAccum < alphaAcumLimit) { + Point3 pos = d_cameraPos + rayDir * t; // Convert to volume indices int ix = (int)roundf(pos.x); @@ -234,31 +357,46 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { int iz = (int)roundf(pos.z); // Sample (pick appropriate method based on volume size) - // float density = sampleVolumeNearest(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz); - float density = sampleVolumeTrilinear(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, pos.x, pos.y, pos.z); + float density = sampleVolumeNearest(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz); + // float density = sampleVolumeTrilinear(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, pos.x, pos.y, pos.z); // If density ~ 0, skip shading 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 * (alphaAcumLimit - alphaAccum); - colorG += color.y * (alphaAcumLimit - alphaAccum); - colorB += color.z * (alphaAcumLimit - alphaAccum); - alphaAccum += color.w * (alphaAcumLimit - alphaAccum); + float4 color = transferFunction(density, grad, pos, rayDir); // This already returns the alpha-weighted color + + // // Accumulate color, respecting current transparency + // colorR += (1.0f - alphaAccum) * color.x * color.w; + // colorG += (1.0f - alphaAccum) * color.y * color.w; + // colorB += (1.0f - alphaAccum) * color.z * color.w; + + // // Update the total accumulated alpha + // alphaAccum += (1.0f - alphaAccum) * color.w; + + colorR = (1.0f - alphaAccum) * color.x + colorR; + colorG = (1.0f - alphaAccum) * color.y + colorG; + colorB = (1.0f - alphaAccum) * color.z + colorB; + alphaAccum = (1 - alphaAccum) * color.w + alphaAccum; + } - tCurrent += stepSize; + t += stepSize; } + + // Calculate final colour accumR += colorR; accumG += colorG; accumB += colorB; + accumA += alphaAccum; + // Blend with background (for transparency) // TODO: Put background colour in a constant float leftover = 1.0 - alphaAccum; accumR = accumR + leftover * 0.1f; accumG = accumG + leftover * 0.1f; accumB = accumB + leftover * 0.1f; + // accumA = accumA + leftover * 0.0f; // This is set to 0.9 instead of 1.0 to be able to see where the volume is a little bit } } @@ -267,9 +405,10 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { accumR /= (float)SAMPLES_PER_PIXEL; accumG /= (float)SAMPLES_PER_PIXEL; accumB /= (float)SAMPLES_PER_PIXEL; + accumA /= (float)SAMPLES_PER_PIXEL; // Final colour - framebuffer.writePixel(px, py, accumR, accumG, accumB); + framebuffer.writePixel(px, py, accumR, accumG, accumB, accumA); // 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); diff --git a/src/main.cu b/src/main.cu index 43194f5..942e4ae 100644 --- a/src/main.cu +++ b/src/main.cu @@ -55,8 +55,8 @@ void getSpeed(std::vector& speedData, int idx = 0) { int main() { std::vector data; - getTemperature(data); - // getSpeed(data); + getTemperature(data, 0); + // getSpeed(data, 294); std::cout << "DATA size: " << data.size() << std::endl; @@ -65,13 +65,36 @@ int main() { 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 + for (int i = 0; i < VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH; i++) { // 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] = -infty; inftyCount++;} + hostVolume[i] = data[i + 0*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH]; + if (data[i + 0*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH] + epsilon >= infty) {hostVolume[i] = -infty; inftyCount++;} } std::cout << "inftyCount: " << inftyCount << std::endl; + // Reverse the order of hostVolume + for (int i = 0; i < VOLUME_WIDTH; i++) { + for (int j = 0; j < VOLUME_HEIGHT; j++) { + for (int k = 0; k < VOLUME_DEPTH/2; k++) { + float temp = hostVolume[i + j*VOLUME_WIDTH + k*VOLUME_WIDTH*VOLUME_HEIGHT]; + hostVolume[i + j*VOLUME_WIDTH + k*VOLUME_WIDTH*VOLUME_HEIGHT] = hostVolume[i + j*VOLUME_WIDTH + (VOLUME_DEPTH - 1 - k)*VOLUME_WIDTH*VOLUME_HEIGHT]; + hostVolume[i + j*VOLUME_WIDTH + (VOLUME_DEPTH - 1 - k)*VOLUME_WIDTH*VOLUME_HEIGHT] = temp; + } + } + } + + + // Store the half-way up slice data into a file + std::ofstream myfile; + myfile.open("halfwayup.txt"); + for (int i = 0; i < VOLUME_WIDTH; i++) { + for (int j = 0; j < VOLUME_HEIGHT; j++) { + myfile << hostVolume[i + j*VOLUME_WIDTH + VOLUME_DEPTH/2*VOLUME_WIDTH*VOLUME_HEIGHT] << " "; + } + myfile << std::endl; + } + myfile.close(); + 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; From 417913af2b52aed901394f7368469348cf82cedd Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 10:59:08 +0100 Subject: [PATCH 15/25] Determined volume coordinates direction --- src/illumination/Raycaster.cu | 59 ----------------------------------- 1 file changed, 59 deletions(-) diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index 8376a5b..33d0344 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -14,30 +14,6 @@ // TODO: Probbably move this transfer function business into a different file // Samples the voxel nearest to the given coordinates. TODO: Can be re-used in other places so move __device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz) { - // TODO: I actually do not know which screen coordinates corresponds to which volume coordinates ... so let's try all options - - // x <-> width, y <-> height, z <-> depth - // if (vx < 0) vx = 0; - // if (vy < 0) vy = 0; - // if (vz < 0) vz = 0; - // if (vx >= volW) vx = volW - 1; - // if (vy >= volH) vy = volH - 1; - // if (vz >= volD) vz = volD - 1; - - // int idx = vz * volW * volH + vy * volD + vx; - // return volumeData[idx]; - - // x <-> width, y <-> depth, z <-> height - // if (vx < 0) vx = 0; - // if (vy < 0) vy = 0; - // if (vz < 0) vz = 0; - // if (vx >= volW) vx = volW - 1; - // if (vy >= volH) vy = volH - 1; - // if (vz >= volD) vz = volD - 1; - - // int idx = vy * volW * volD + vz * volW + vx; - // return volumeData[idx]; - // x <-> height, y <-> width, z <-> depth <--- So far this is the best one if (vx < 0) vx = 0; if (vy < 0) vy = 0; @@ -48,41 +24,6 @@ __device__ float sampleVolumeNearest(float* volumeData, const int volW, const in int idx = vz * volW * volH + vx * volW + vy; return volumeData[idx]; - - // x <-> height, y <-> depth, z <-> width - // if (vx < 0) vx = 0; - // if (vy < 0) vy = 0; - // if (vz < 0) vz = 0; - // if (vx >= volH) vx = volH - 1; - // if (vy >= volD) vy = volD - 1; - // if (vz >= volW) vz = volW - 1; - - // int idx = vy * volW * volH + vz * volW + vx; - // return volumeData[idx]; - - - // x <-> depth, y <-> width, z <-> height - // if (vx < 0) vx = 0; - // if (vy < 0) vy = 0; - // if (vz < 0) vz = 0; - // if (vx >= volD) vx = volD - 1; - // if (vy >= volW) vy = volW - 1; - // if (vz >= volH) vz = volH - 1; - - // int idx = vz * volD * volW + vx * volW + vy; - // return volumeData[idx]; - - // x <-> depth, y <-> height, z <-> width - // if (vx < 0) vx = 0; - // if (vy < 0) vy = 0; - // if (vz < 0) vz = 0; - // if (vx >= volD) vx = volD - 1; - // if (vy >= volH) vy = volH - 1; - // if (vz >= volW) vz = volW - 1; - - // int idx = vy * volD * volW + vz * volD + vx; - // return volumeData[idx]; - } // tri-linear interpolation - ready if necessary (but no visible improvement for full volume) From 1711bcf7566f6a7608919171e1b988c9018d0015 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 11:50:19 +0100 Subject: [PATCH 16/25] Moved transfer function functions to their own file --- src/consts.cu | 35 +++++ src/consts.h | 14 ++ src/illumination/Raycaster.cu | 206 +-------------------------- src/illumination/transferFunction.cu | 129 +++++++++++++++++ src/illumination/transferFunction.h | 20 +++ 5 files changed, 200 insertions(+), 204 deletions(-) create mode 100644 src/illumination/transferFunction.cu create mode 100644 src/illumination/transferFunction.h diff --git a/src/consts.cu b/src/consts.cu index c55069f..8022c53 100644 --- a/src/consts.cu +++ b/src/consts.cu @@ -1,5 +1,31 @@ #include "consts.h" +// ----------------------- Colour mapping ----------------------- +__constant__ ColorStop d_stopsPythonLike[5]; +__constant__ ColorStop d_stopsGrayscale[2]; +__constant__ ColorStop d_stopsBluePurleRed[3]; + +const ColorStop h_stopsPythonLike[] = { + { 0.0f, Color3::init(0.2298057f, 0.29871797f, 0.75368315f) }, // Dark Blue + { 0.25f, Color3::init(0.23437708f, 0.30554173f, 0.75967953f) }, // Mid Blue + { 0.5f, Color3::init(0.27582712f, 0.36671692f, 0.81255294f) }, // White + { 0.75f, Color3::init(0.79606387f, 0.84869321f, 0.93347147f) }, // Light Orange + { 1.0f, Color3::init(0.70567316f, 0.01555616f, 0.15023281f) } // Red +}; + +const ColorStop h_stopsGrayscale[] = { + { 0.0f, Color3::init(0.0f, 0.0f, 0.0f) }, // No colour + { 1.0f, Color3::init(1.0f, 1.0f, 1.0f) } // White +}; + +const ColorStop h_stopsBluePurleRed[] = { + { 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 +}; + +// ----------------------- Camera and Light ----------------------- + __device__ Point3 d_cameraPos; __device__ Vec3 d_cameraDir; __device__ Vec3 d_cameraUp; @@ -12,7 +38,16 @@ Vec3 h_cameraDir = (center - h_cameraPos).normalize(); Vec3 h_cameraUp = Vec3::init(0.0, 0.0, 1.0).normalize(); Point3 h_lightPos = Point3::init(1.5, 2.0, -1.0); + +// Copy the above values to the device void copyConstantsToDevice() { + // ----------------------- Colour mapping ----------------------- + cudaMemcpyToSymbol(d_stopsPythonLike, h_stopsPythonLike, sizeof(h_stopsPythonLike)); + cudaMemcpyToSymbol(d_stopsGrayscale, h_stopsGrayscale, sizeof(h_stopsGrayscale)); + cudaMemcpyToSymbol(d_stopsBluePurleRed, h_stopsBluePurleRed, sizeof(h_stopsBluePurleRed)); + + + // ----------------------- Camera and Light ----------------------- cudaMemcpyToSymbol(d_cameraPos, &h_cameraPos, sizeof(Point3)); cudaMemcpyToSymbol(d_cameraDir, &h_cameraDir, sizeof(Vec3)); cudaMemcpyToSymbol(d_cameraUp, &h_cameraUp, sizeof(Vec3)); diff --git a/src/consts.h b/src/consts.h index dd2e4f6..efd4131 100644 --- a/src/consts.h +++ b/src/consts.h @@ -48,6 +48,20 @@ extern __device__ Vec3 d_cameraDir; extern __device__ Vec3 d_cameraUp; extern __device__ Point3 d_lightPos; + +// --------------------------- Transfer Function Constants --------------------------- +struct ColorStop { + float pos; // in [0,1] + Color3 color; +}; + +const int lenStopsPythonLike = 5; +const int lenStopsGrayscale = 2; +const int lenStopsBluePurpleRed = 3; +extern __constant__ ColorStop d_stopsPythonLike[5]; +extern __constant__ ColorStop d_stopsGrayscale[2]; +extern __constant__ ColorStop d_stopsBluePurleRed[3]; + // --------------------------- Functions for handling external constants --------------------------- void copyConstantsToDevice(); diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index 33d0344..ca16136 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -5,214 +5,12 @@ #include "linalg/linalg.h" #include "consts.h" +#include "transferFunction.h" #include "cuda_error.h" -#include "shading.h" + #include -#include "objs/sphere.h" #include -// TODO: Probbably move this transfer function business into a different file -// Samples the voxel nearest to the given coordinates. TODO: Can be re-used in other places so move -__device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz) { - // x <-> height, y <-> width, z <-> depth <--- So far this is the best one - if (vx < 0) vx = 0; - if (vy < 0) vy = 0; - if (vz < 0) vz = 0; - if (vx >= volH) vx = volH - 1; - if (vy >= volW) vy = volW - 1; - if (vz >= volD) vz = volD - 1; - - int idx = vz * volW * volH + vx * volW + vy; - return volumeData[idx]; -} - -// tri-linear interpolation - ready if necessary (but no visible improvement for full volume) -__device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz) { - int ix = (int)floorf(fx); - int iy = (int)floorf(fy); - int iz = (int)floorf(fz); - - // Clamp indices to valid range - int ix1 = min(ix + 1, volW - 1); - int iy1 = min(iy + 1, volH - 1); - int iz1 = min(iz + 1, volD - 1); - ix = max(ix, 0); - iy = max(iy, 0); - iz = max(iz, 0); - - // Compute weights - float dx = fx - ix; - float dy = fy - iy; - float dz = fz - iz; - - // Sample values - float c00 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz) * (1.0f - dx) + - sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy, iz) * dx; - float c10 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy1, iz) * (1.0f - dx) + - sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy1, iz) * dx; - float c01 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz1) * (1.0f - dx) + - sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy, iz1) * dx; - float c11 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy1, iz1) * (1.0f - dx) + - sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy1, iz1) * dx; - - float c0 = c00 * (1.0f - dy) + c10 * dy; - float c1 = c01 * (1.0f - dy) + c11 * dy; - - return c0 * (1.0f - dz) + c1 * dz; -} - - -__device__ float opacityFromGradient(const Vec3 &grad) { - float gradMag = grad.length(); // magnitude - float k = 1e-4f; // 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; -} - -// Monochromatic colormap for speed -__device__ Color3 colorMapMonochrome(float normalizedSpeed) { - // Define the color stops: black (0.0) to white (1.0) - ColorStop speedStops[] = { - { 0.0f, Color3::init(0.0f, 0.0f, 0.0f) }, // No colour - { 1.0f, Color3::init(1.0f, 1.0f, 1.0f) } // White - }; - - // Clamp to [0,1] - if (normalizedSpeed < 0.0f) normalizedSpeed = 0.0f; - if (normalizedSpeed > 1.0f) normalizedSpeed = 1.0f; - - // Single interval (N = 2 for black to white) - const int N = 2; - for (int i = 0; i < N - 1; ++i) - { - float start = speedStops[i].pos; - float end = speedStops[i + 1].pos; - - if (normalizedSpeed >= start && normalizedSpeed <= end) - { - float localT = (normalizedSpeed - start) / (end - start); - return interpolate(speedStops[i].color, speedStops[i + 1].color, localT); - } - } - - // Fallback - return speedStops[N - 1].color; -} - -// Colormap function for gradient from blue to white to red -__device__ Color3 colorMapPythonLike(float normalizedValue) { - // Define control points for the colormap (approximation of coolwarm) - ColorStop coolwarmStops[] = { - { 0.0f, Color3::init(0.2298057f, 0.29871797f, 0.75368315f) }, // Dark Blue - { 0.25f, Color3::init(0.23437708f, 0.30554173f, 0.75967953f) }, // Mid Blue - { 0.5f, Color3::init(0.27582712f, 0.36671692f, 0.81255294f) }, // White - { 0.75f, Color3::init(0.79606387f, 0.84869321f, 0.93347147f) }, // Light Orange - { 1.0f, Color3::init(0.70567316f, 0.01555616f, 0.15023281f) } // Red - }; - - // Clamp the scalar value to [0, 1] - normalizedValue = fminf(fmaxf(normalizedValue, 0.0f), 1.0f); - - // Interpolate between the defined color stops - const int N = 5; // Number of control points - for (int i = 0; i < N - 1; ++i) { - float start = coolwarmStops[i].pos; - float end = coolwarmStops[i + 1].pos; - - if (normalizedValue >= start && normalizedValue <= end) { - float localT = (normalizedValue - start) / (end - start); - return interpolate(coolwarmStops[i].color, coolwarmStops[i + 1].color, localT); - } - } - - // Fallback (shouldn't reach here due to clamping) - return coolwarmStops[N - 1].color; -} - - -// Transfer function -__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 - - // Color3 baseColor = Color3::init(density, 0.1f*density, 1.f - density); // TODO: Implement a proper transfer function - // Color3 baseColor = temperatureToRGB(density); - - float normDensity = (density - MIN_TEMP) / (MAX_TEMP - MIN_TEMP); - // float normDensity = (density - MIN_SPEED) / (MAX_SPEED - MIN_SPEED); - - normDensity = clamp(normDensity, 0.0f, 1.0f); - - // Color3 baseColor = colorMapViridis(normDensity); - // Color3 baseColor = colorMapMonochrome(normDensity); - Color3 baseColor = colorMapPythonLike(normDensity); - - - float alpha = opacityFromGradient(grad); - // alpha = 0.1f; - alpha = 1.0f / (1.0f + expf(-250.f * (normDensity - 0.5f))); // This is also quite nice, but the exponent should be parameterized - float alphaSample = density * alpha; // TODO: Decide whether to keep alpha here or not - - Vec3 normal = -grad.normalize(); - - Vec3 lightDir = (d_lightPos - pos).normalize(); - Vec3 viewDir = -rayDir.normalize(); - - // Apply Phong - Vec3 shadedColor = phongShading(normal, lightDir, viewDir, baseColor); - - // Compose - 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: 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; - result.w = 1.0f; - } - - return result; -} diff --git a/src/illumination/transferFunction.cu b/src/illumination/transferFunction.cu new file mode 100644 index 0000000..9490078 --- /dev/null +++ b/src/illumination/transferFunction.cu @@ -0,0 +1,129 @@ +#include "transferFunction.h" + + +// Samples the voxel nearest to the given coordinates. +__device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz) { + // x <-> height, y <-> width, z <-> depth <--- So far this is the best one + if (vx < 0) vx = 0; + if (vy < 0) vy = 0; + if (vz < 0) vz = 0; + if (vx >= volH) vx = volH - 1; + if (vy >= volW) vy = volW - 1; + if (vz >= volD) vz = volD - 1; + + int idx = vz * volW * volH + vx * volW + vy; + return volumeData[idx]; +} + +// tri-linear interpolation - ready if necessary (but no visible improvement for full volume) +__device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz) { + int ix = (int)floorf(fx); + int iy = (int)floorf(fy); + int iz = (int)floorf(fz); + + // Clamp indices to valid range + int ix1 = min(ix + 1, volW - 1); + int iy1 = min(iy + 1, volH - 1); + int iz1 = min(iz + 1, volD - 1); + ix = max(ix, 0); + iy = max(iy, 0); + iz = max(iz, 0); + + // Compute weights + float dx = fx - ix; + float dy = fy - iy; + float dz = fz - iz; + + // Sample values + float c00 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz) * (1.0f - dx) + + sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy, iz) * dx; + float c10 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy1, iz) * (1.0f - dx) + + sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy1, iz) * dx; + float c01 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy, iz1) * (1.0f - dx) + + sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy, iz1) * dx; + float c11 = sampleVolumeNearest(volumeData, volW, volH, volD, ix, iy1, iz1) * (1.0f - dx) + + sampleVolumeNearest(volumeData, volW, volH, volD, ix1, iy1, iz1) * dx; + + float c0 = c00 * (1.0f - dy) + c10 * dy; + float c1 = c01 * (1.0f - dy) + c11 * dy; + + return c0 * (1.0f - dz) + c1 * dz; +} + + +__device__ float opacityFromGradient(const Vec3 &grad) { + float gradMag = grad.length(); + float k = 1e-4f; // tweak (the smaller the value, the less opacity) // TODO: Add a slider for this + float alpha = 1.0f - expf(-k * gradMag); + return alpha; +} + +__device__ float opacitySigmoid(float normDensity) { + return 1.0f / (1.0f + expf(-250.f * (normDensity - 0.5f))); // TODO: Parametrize and add sliders +} + +__device__ Color3 colorMap(float normalizedValues, const ColorStop stops[], int N) { + // clamp to [0,1] + normalizedValues = fminf(fmaxf(normalizedValues, 0.0f), 1.0f); + + // N stops => N-1 intervals + for (int i = 0; i < N - 1; ++i) { + float start = stops[i].pos; + float end = stops[i + 1].pos; + + if (normalizedValues >= start && normalizedValues <= end) { + float localT = (normalizedValues - start) / (end - start); + return interpolate(stops[i].color, stops[i + 1].color, localT); + } + } + + // fallback if something goes out of [0,1] or numerical issues + return stops[N - 1].color; +} + + +// Transfer function +__device__ float4 transferFunction(float density, const Vec3& grad, const Point3& pos, const Vec3& rayDir) { + + // --------------------------- Sample the volume --------------------------- + // TODO: Somehow pick if to use temp of speed normalization ... or pass extremas as params. + float normDensity = (density - MIN_TEMP) / (MAX_TEMP - MIN_TEMP); + // float normDensity = (density - MIN_SPEED) / (MAX_SPEED - MIN_SPEED); + + normDensity = clamp(normDensity, 0.0f, 1.0f); + + // --------------------------- Map density to color --------------------------- + // TODO: Add a way to pick stops here + Color3 baseColor = colorMap(normDensity, d_stopsPythonLike, lenStopsPythonLike); + + + float alpha = opacityFromGradient(grad); + // alpha = 0.1f; + alpha = opacitySigmoid(normDensity); + float alphaSample = density * alpha; // TODO: Decide whether to keep alpha here or not + + // --------------------------- Shading --------------------------- + // Apply Phong + Vec3 normal = -grad.normalize(); + Vec3 lightDir = (d_lightPos - pos).normalize(); + Vec3 viewDir = -rayDir.normalize(); + Vec3 shadedColor = phongShading(normal, lightDir, viewDir, baseColor); + + // Compose + 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 + + // --------------------------- Silhouettes --------------------------- + // 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 more pronounced edges + if (grad.length() > epsilon && fabs(grad.normalize().dot(viewDir)) < 0.2f) { + result.x = 0.0f; + result.y = 0.0f; + result.z = 0.0f; + result.w = 1.0f; + } + + return result; +} diff --git a/src/illumination/transferFunction.h b/src/illumination/transferFunction.h new file mode 100644 index 0000000..cd33929 --- /dev/null +++ b/src/illumination/transferFunction.h @@ -0,0 +1,20 @@ +#ifndef TRANSFER_FUNCTION_H +#define TRANSFER_FUNCTION_H + +#include "linalg/linalg.h" +#include "consts.h" +#include "shading.h" + +// --------------------------- Color mapping --------------------------- + + +// --------------------------- Volume sampling --------------------------- +__device__ float sampleVolumeNearest(float* volumeData, const int volW, const int volH, const int volD, int vx, int vy, int vz); +__device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz); + + +// --------------------------- Transfer function --------------------------- +__device__ float4 transferFunction(float density, const Vec3& grad, const Point3& pos, const Vec3& rayDir); + + +#endif // TRANSFER_FUNCTION_H \ No newline at end of file From a5c66fc6ca8bb4e29c691b2a23d417202040e36d Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 13:00:43 +0100 Subject: [PATCH 17/25] Removed some old comments --- src/consts.h | 3 +-- src/illumination/Raycaster.cu | 26 +++++--------------------- 2 files changed, 6 insertions(+), 23 deletions(-) diff --git a/src/consts.h b/src/consts.h index efd4131..6164aad 100644 --- a/src/consts.h +++ b/src/consts.h @@ -5,7 +5,6 @@ #include // --------------------------- Basic Constants --------------------------- -// 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 @@ -27,7 +26,7 @@ const float MAX_SPEED = 14.0f; // --------------------------- Raycasting Constants --------------------------- -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 int SAMPLES_PER_PIXEL = 1; const float alphaAcumLimit = 1.0f; // 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; diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index ca16136..317b788 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -12,8 +12,6 @@ #include - - // TODO: instead of IMAGEWIDTH and IMAGEHEIGHT this should reflect the windowSize; __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { int px = blockIdx.x * blockDim.x + threadIdx.x; @@ -28,7 +26,6 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { curandState randState; curand_init(1234, px + py * IMAGE_WIDTH, 0, &randState); - // Multiple samples per pixel // Multiple samples per pixel for (int s = 0; s < SAMPLES_PER_PIXEL; s++) { // Map to [-1, 1] @@ -37,7 +34,6 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { float u = ((px + 0.5f + jitterU) / IMAGE_WIDTH ) * 2.0f - 1.0f; float v = ((py + 0.5f + jitterV) / 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; @@ -86,7 +82,7 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { float colorR = 0.0f, colorG = 0.0f, colorB = 0.0f; float alphaAccum = 0.0f; - float t = tNear; + float t = tNear; // Front to back while (t < tFar && alphaAccum < alphaAcumLimit) { Point3 pos = d_cameraPos + rayDir * t; @@ -95,23 +91,16 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { int iy = (int)roundf(pos.y); int iz = (int)roundf(pos.z); - // Sample (pick appropriate method based on volume size) - float density = sampleVolumeNearest(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz); - // float density = sampleVolumeTrilinear(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, pos.x, pos.y, pos.z); + // Sample (pick appropriate method based on volume size) TODO: Add a way to pick this in GUI + // float density = sampleVolumeNearest(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz); + float density = sampleVolumeTrilinear(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, pos.x, pos.y, pos.z); // If density ~ 0, skip shading if (density > minAllowedDensity) { Vec3 grad = computeGradient(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, ix, iy, iz); float4 color = transferFunction(density, grad, pos, rayDir); // This already returns the alpha-weighted color - // // Accumulate color, respecting current transparency - // colorR += (1.0f - alphaAccum) * color.x * color.w; - // colorG += (1.0f - alphaAccum) * color.y * color.w; - // colorB += (1.0f - alphaAccum) * color.z * color.w; - - // // Update the total accumulated alpha - // alphaAccum += (1.0f - alphaAccum) * color.w; - + //Accumulate color, and alpha colorR = (1.0f - alphaAccum) * color.x + colorR; colorG = (1.0f - alphaAccum) * color.y + colorG; colorB = (1.0f - alphaAccum) * color.z + colorB; @@ -135,7 +124,6 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { accumR = accumR + leftover * 0.1f; accumG = accumG + leftover * 0.1f; accumB = accumB + leftover * 0.1f; - // accumA = accumA + leftover * 0.0f; // This is set to 0.9 instead of 1.0 to be able to see where the volume is a little bit } } @@ -148,10 +136,6 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { // Final colour framebuffer.writePixel(px, py, accumR, accumG, accumB, accumA); - // 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); } From 9afd242121d0cacc1a9973818a3ebf8163d4ada3 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 13:05:18 +0100 Subject: [PATCH 18/25] Added background colour to consts.h --- src/consts.cu | 3 +++ src/consts.h | 3 +++ src/illumination/Raycaster.cu | 14 +++++++------- 3 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/consts.cu b/src/consts.cu index 8022c53..160720c 100644 --- a/src/consts.cu +++ b/src/consts.cu @@ -30,6 +30,7 @@ __device__ Point3 d_cameraPos; __device__ Vec3 d_cameraDir; __device__ Vec3 d_cameraUp; __device__ Point3 d_lightPos; +__device__ Color3 d_backgroundColor; // 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) @@ -37,6 +38,7 @@ Vec3 center = Vec3::init((float)VOLUME_WIDTH/2.0f, (float)VOLUME_HEIGHT/2.0f, (f Vec3 h_cameraDir = (center - h_cameraPos).normalize(); Vec3 h_cameraUp = Vec3::init(0.0, 0.0, 1.0).normalize(); Point3 h_lightPos = Point3::init(1.5, 2.0, -1.0); +Color3 h_backgroundColor = Color3::init(0.1f, 0.1f, 0.1f); // Copy the above values to the device @@ -52,4 +54,5 @@ void copyConstantsToDevice() { cudaMemcpyToSymbol(d_cameraDir, &h_cameraDir, sizeof(Vec3)); cudaMemcpyToSymbol(d_cameraUp, &h_cameraUp, sizeof(Vec3)); cudaMemcpyToSymbol(d_lightPos, &h_lightPos, sizeof(Point3)); + cudaMemcpyToSymbol(d_backgroundColor, &h_backgroundColor, sizeof(Color3)); } diff --git a/src/consts.h b/src/consts.h index 6164aad..ef8cce1 100644 --- a/src/consts.h +++ b/src/consts.h @@ -35,6 +35,7 @@ const float stepSize = 0.02f; // --------------------------- Illumination Constants --------------------------- +// Shading consts const double ambientStrength = 0.3; const double diffuseStrength = 0.8; const double specularStrength = 0.5; @@ -47,6 +48,8 @@ extern __device__ Vec3 d_cameraDir; extern __device__ Vec3 d_cameraUp; extern __device__ Point3 d_lightPos; +// Background color +extern __device__ Color3 d_backgroundColor; // --------------------------- Transfer Function Constants --------------------------- struct ColorStop { diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index 317b788..763c4f2 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -72,9 +72,9 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { if (tNear > tFar) { // 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; + accumR = d_backgroundColor.x * (float)SAMPLES_PER_PIXEL; + accumG = d_backgroundColor.y * (float)SAMPLES_PER_PIXEL; + accumB = d_backgroundColor.z * (float)SAMPLES_PER_PIXEL; accumA = 1.0f * (float)SAMPLES_PER_PIXEL; } else { if (tNear < 0.0f) tNear = 0.0f; @@ -119,11 +119,11 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { accumB += colorB; accumA += alphaAccum; - // Blend with background (for transparency) // TODO: Put background colour in a constant + // Blend with background (for transparency) float leftover = 1.0 - alphaAccum; - accumR = accumR + leftover * 0.1f; - accumG = accumG + leftover * 0.1f; - accumB = accumB + leftover * 0.1f; + accumR = accumR + leftover * d_backgroundColor.x; + accumG = accumG + leftover * d_backgroundColor.y; + accumB = accumB + leftover * d_backgroundColor.z; } } From 4971d693fdfbd60a18dc7bf789330a70c92a2a29 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 13:25:24 +0100 Subject: [PATCH 19/25] Fixed some TODOs, removed old comments --- src/consts.h | 2 +- src/illumination/shading.cu | 1 - src/main.cu | 25 ++++++------------------- 3 files changed, 7 insertions(+), 21 deletions(-) diff --git a/src/consts.h b/src/consts.h index ef8cce1..d53abff 100644 --- a/src/consts.h +++ b/src/consts.h @@ -28,7 +28,7 @@ const float MAX_SPEED = 14.0f; // --------------------------- Raycasting Constants --------------------------- const int SAMPLES_PER_PIXEL = 1; -const float alphaAcumLimit = 1.0f; // TODO: Idk what a good accumulation value is <--- This finally does something when using alpha in both places at least +const float alphaAcumLimit = 1.0f; // TODO: Atm, this does not work very intuitively, other parameters control transparency const float minAllowedDensity = 0.001f; const float stepSize = 0.02f; diff --git a/src/illumination/shading.cu b/src/illumination/shading.cu index 41141f0..8fa8970 100644 --- a/src/illumination/shading.cu +++ b/src/illumination/shading.cu @@ -1,6 +1,5 @@ #include "shading.h" -// TODO: Consider wrapping this in a class (?) __device__ Vec3 phongShading(const Vec3& normal, const Vec3& lightDir, const Vec3& viewDir, const Vec3& baseColor) { Vec3 ambient = baseColor * ambientStrength; double diff = fmax(normal.dot(lightDir), 0.0); diff --git a/src/main.cu b/src/main.cu index 942e4ae..e3faac7 100644 --- a/src/main.cu +++ b/src/main.cu @@ -60,19 +60,15 @@ int main() { std::cout << "DATA size: " << data.size() << std::endl; - // TODO: Eveontually remove debug below (i.e., eliminate for-loop etc.) - // Generate debug volume data + // TODO: Eventually, we should not need to load the volume like this 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++) { - // Discard temperatures above a small star (supposedly, missing temperature values) hostVolume[i] = data[i + 0*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH]; - if (data[i + 0*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH] + epsilon >= infty) {hostVolume[i] = -infty; inftyCount++;} + // Discard missing values + if (data[i + 0*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH] + epsilon >= infty) hostVolume[i] = -infty; } - std::cout << "inftyCount: " << inftyCount << std::endl; - // Reverse the order of hostVolume + // Reverse the order of hostVolume TODO: Idk why the volume upside down, this can probably be deleted and substituted with a minus sign somehwere (camera position?) for (int i = 0; i < VOLUME_WIDTH; i++) { for (int j = 0; j < VOLUME_HEIGHT; j++) { for (int k = 0; k < VOLUME_DEPTH/2; k++) { @@ -83,8 +79,7 @@ int main() { } } - - // Store the half-way up slice data into a file + // Store the half-way up slice data into a file TODO: Remove this debug std::ofstream myfile; myfile.open("halfwayup.txt"); for (int i = 0; i < VOLUME_WIDTH; i++) { @@ -95,6 +90,7 @@ int main() { } myfile.close(); + // Print min, max, avg., and median values TODO: Remove this debug 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; @@ -102,15 +98,6 @@ int main() { }); 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); - // } - // // 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); // float sum = std::accumulate(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH, 0.0f); From b151e70167d9148f5831aa19c8eb68dd389b137e Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 13:26:05 +0100 Subject: [PATCH 20/25] Fixed mismatch between coordinates indices (rather min/max values thereof) --- src/illumination/Raycaster.cu | 3 ++- src/illumination/transferFunction.cu | 16 +++++++++------- src/linalg/mat.cu | 24 ++++++++++++++++-------- 3 files changed, 27 insertions(+), 16 deletions(-) diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index 763c4f2..c9015cf 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -22,7 +22,8 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { float accumG = 0.0f; float accumB = 0.0f; float accumA = 1.0f; - // Initialize random state + + // Initialize random state for ray scattering curandState randState; curand_init(1234, px + py * IMAGE_WIDTH, 0, &randState); diff --git a/src/illumination/transferFunction.cu b/src/illumination/transferFunction.cu index 9490078..1012555 100644 --- a/src/illumination/transferFunction.cu +++ b/src/illumination/transferFunction.cu @@ -22,8 +22,8 @@ __device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const int iz = (int)floorf(fz); // Clamp indices to valid range - int ix1 = min(ix + 1, volW - 1); - int iy1 = min(iy + 1, volH - 1); + int ix1 = min(ix + 1, volH - 1); + int iy1 = min(iy + 1, volW - 1); int iz1 = min(iz + 1, volD - 1); ix = max(ix, 0); iy = max(iy, 0); @@ -96,11 +96,12 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3 // TODO: Add a way to pick stops here Color3 baseColor = colorMap(normDensity, d_stopsPythonLike, lenStopsPythonLike); - + // TODO: Add a way to pick different function for alpha float alpha = opacityFromGradient(grad); - // alpha = 0.1f; + alpha = 0.1f; alpha = opacitySigmoid(normDensity); - float alphaSample = density * alpha; // TODO: Decide whether to keep alpha here or not + + float alphaSample = density * alpha; // --------------------------- Shading --------------------------- // Apply Phong @@ -114,11 +115,12 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3 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 + result.w = alpha; // --------------------------- Silhouettes --------------------------- // 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 more pronounced edges - if (grad.length() > epsilon && fabs(grad.normalize().dot(viewDir)) < 0.2f) { + // TODO: Add a way to adjust the treshold (0.2f atm) + if (grad.length() > epsilon && fabs(grad.normalize().dot(viewDir.normalize())) < 0.2f) { result.x = 0.0f; result.y = 0.0f; result.z = 0.0f; diff --git a/src/linalg/mat.cu b/src/linalg/mat.cu index c206cb1..2dd2564 100644 --- a/src/linalg/mat.cu +++ b/src/linalg/mat.cu @@ -10,19 +10,27 @@ __device__ Vec3 computeGradient(float* volumeData, const int volW, const int vol // Normal should point from higher to lower intensities int xm = max(x - 1, 0); - int xp = min(x + 1, volW - 1); + int xp = min(x + 1, volH - 1); int ym = max(y - 1, 0); - int yp = min(y + 1, volH - 1); + int yp = min(y + 1, volW - 1); int zm = max(z - 1, 0); int zp = min(z + 1, volD - 1); // Note: Assuming data is linearized (idx = z*w*h + y*w + x) TODO: Unlinearize if data not linear - float gx = volumeData[z * volW * volH + y * volW + xp] - - volumeData[z * volW * volH + y * volW + xm]; - float gy = volumeData[z * volW * volH + yp * volW + x ] - - volumeData[z * volW * volH + ym * volW + x ]; - float gz = volumeData[zp * volW * volH + y * volW + x ] - - volumeData[zm * volW * volH + y * volW + x ]; + // float gx = volumeData[z * volW * volH + y * volW + xp] + // - volumeData[z * volW * volH + y * volW + xm]; + // float gy = volumeData[z * volW * volH + yp * volW + x ] + // - volumeData[z * volW * volH + ym * volW + x ]; + // float gz = volumeData[zp * volW * volH + y * volW + x ] + // - volumeData[zm * volW * volH + y * volW + x ]; + + // Note: Assuming data is linearized (idx = z * W * H + x * W + y;) TODO: Unlinearize if data not linear + float gx = volumeData[z * volW * volH + x * volW + xp] + - volumeData[z * volW * volH + x * volW + xm]; + float gy = volumeData[z * volW * volH + y * volW + yp] + - volumeData[z * volW * volH + y * volW + ym]; + float gz = volumeData[zp * volW * volH + x * volW + y ] + - volumeData[zm * volW * volH + x * volW + y ]; return Vec3::init(gx, gy, gz); }; From 371438980371bb054e3d677fc0a93ace94cf00f0 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 14:39:52 +0100 Subject: [PATCH 21/25] Fixed gradient --- src/linalg/mat.cu | 17 ++++------------- src/linalg/mat.h | 1 + 2 files changed, 5 insertions(+), 13 deletions(-) diff --git a/src/linalg/mat.cu b/src/linalg/mat.cu index 2dd2564..2d1c8dc 100644 --- a/src/linalg/mat.cu +++ b/src/linalg/mat.cu @@ -7,7 +7,6 @@ using namespace std; __device__ Vec3 computeGradient(float* volumeData, const int volW, const int volH, const int volD, int x, int y, int z) { // Finite difference for partial derivatives. // For boundary voxels - clamp to the boundary. - // Normal should point from higher to lower intensities int xm = max(x - 1, 0); int xp = min(x + 1, volH - 1); @@ -16,19 +15,11 @@ __device__ Vec3 computeGradient(float* volumeData, const int volW, const int vol int zm = max(z - 1, 0); int zp = min(z + 1, volD - 1); - // Note: Assuming data is linearized (idx = z*w*h + y*w + x) TODO: Unlinearize if data not linear - // float gx = volumeData[z * volW * volH + y * volW + xp] - // - volumeData[z * volW * volH + y * volW + xm]; - // float gy = volumeData[z * volW * volH + yp * volW + x ] - // - volumeData[z * volW * volH + ym * volW + x ]; - // float gz = volumeData[zp * volW * volH + y * volW + x ] - // - volumeData[zm * volW * volH + y * volW + x ]; - // Note: Assuming data is linearized (idx = z * W * H + x * W + y;) TODO: Unlinearize if data not linear - float gx = volumeData[z * volW * volH + x * volW + xp] - - volumeData[z * volW * volH + x * volW + xm]; - float gy = volumeData[z * volW * volH + y * volW + yp] - - volumeData[z * volW * volH + y * volW + ym]; + float gx = volumeData[z * volW * volH + xp * volW + y] + - volumeData[z * volW * volH + xm * volW + y]; + float gy = volumeData[z * volW * volH + x * volW + yp] + - volumeData[z * volW * volH + x * volW + ym]; float gz = volumeData[zp * volW * volH + x * volW + y ] - volumeData[zm * volW * volH + x * volW + y ]; diff --git a/src/linalg/mat.h b/src/linalg/mat.h index 4d2a1d5..02234c0 100644 --- a/src/linalg/mat.h +++ b/src/linalg/mat.h @@ -2,6 +2,7 @@ #define MAT_H #include "vec.h" +#include "consts.h" __device__ Vec3 computeGradient(float* volumeData, const int volW, const int volH, const int volD, int x, int y, int z); From 1aa33e8e144a80894470a6da046be2ed8947188a Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 14:59:27 +0100 Subject: [PATCH 22/25] Fixed rays not interacting with missing volume incorrectly when number of samples > 1 --- src/consts.h | 4 ++-- src/illumination/Raycaster.cu | 3 ++- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/consts.h b/src/consts.h index d53abff..2686450 100644 --- a/src/consts.h +++ b/src/consts.h @@ -26,9 +26,9 @@ const float MAX_SPEED = 14.0f; // --------------------------- Raycasting Constants --------------------------- -const int SAMPLES_PER_PIXEL = 1; +const int SAMPLES_PER_PIXEL = 4; -const float alphaAcumLimit = 1.0f; // TODO: Atm, this does not work very intuitively, other parameters control transparency +const float alphaAcumLimit = 0.4f; // TODO: Atm, this only works with sigmoid const float minAllowedDensity = 0.001f; const float stepSize = 0.02f; diff --git a/src/illumination/Raycaster.cu b/src/illumination/Raycaster.cu index c9015cf..b2d1939 100644 --- a/src/illumination/Raycaster.cu +++ b/src/illumination/Raycaster.cu @@ -21,7 +21,7 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { float accumR = 0.0f; float accumG = 0.0f; float accumB = 0.0f; - float accumA = 1.0f; + float accumA = 1.0f * (float)SAMPLES_PER_PIXEL; // Initialize random state for ray scattering curandState randState; @@ -77,6 +77,7 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer) { accumG = d_backgroundColor.y * (float)SAMPLES_PER_PIXEL; accumB = d_backgroundColor.z * (float)SAMPLES_PER_PIXEL; accumA = 1.0f * (float)SAMPLES_PER_PIXEL; + } else { if (tNear < 0.0f) tNear = 0.0f; From 3313e8cf84c7cf34c52488838e8f5b37e3a510d0 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 15:09:42 +0100 Subject: [PATCH 23/25] Played around with transfer function -> needs GUI --- src/illumination/transferFunction.cu | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/illumination/transferFunction.cu b/src/illumination/transferFunction.cu index 1012555..fd8db0c 100644 --- a/src/illumination/transferFunction.cu +++ b/src/illumination/transferFunction.cu @@ -53,13 +53,13 @@ __device__ float sampleVolumeTrilinear(float* volumeData, const int volW, const __device__ float opacityFromGradient(const Vec3 &grad) { float gradMag = grad.length(); - float k = 1e-4f; // tweak (the smaller the value, the less opacity) // TODO: Add a slider for this + float k = 1e-6f; // tweak (the smaller the value, the less opacity) // TODO: Add a slider for this float alpha = 1.0f - expf(-k * gradMag); return alpha; } -__device__ float opacitySigmoid(float normDensity) { - return 1.0f / (1.0f + expf(-250.f * (normDensity - 0.5f))); // TODO: Parametrize and add sliders +__device__ float opacitySigmoid(float val) { + return 1.0f / (1.0f + expf(-250.f * (val - 0.5f))); // TODO: Parametrize and add sliders } __device__ Color3 colorMap(float normalizedValues, const ColorStop stops[], int N) { @@ -98,10 +98,11 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3 // TODO: Add a way to pick different function for alpha float alpha = opacityFromGradient(grad); - alpha = 0.1f; + // alpha = 0.1f; alpha = opacitySigmoid(normDensity); + // alpha = (1.0f - fabs(grad.normalize().dot(rayDir.normalize()))) * 0.8f + 0.2f; - float alphaSample = density * alpha; + float alphaSample = density * alpha * 0.1; // --------------------------- Shading --------------------------- // Apply Phong @@ -120,12 +121,13 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3 // --------------------------- Silhouettes --------------------------- // 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 more pronounced edges // TODO: Add a way to adjust the treshold (0.2f atm) - if (grad.length() > epsilon && fabs(grad.normalize().dot(viewDir.normalize())) < 0.2f) { - result.x = 0.0f; - result.y = 0.0f; - result.z = 0.0f; - result.w = 1.0f; - } + // TODO: I don't think we should literally be doing this => use gradient based opacity => delete the below if-statement + // if (fabs(grad.normalize().dot(rayDir.normalize())) < 0.2f) { + // result.x = 0.0f; + // result.y = 0.0f; + // result.z = 0.0f; + // result.w = 1.0f; + // } return result; } From 079e0797e9d3588c586c1dcb6e5b000baf8dd7a2 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 15:11:27 +0100 Subject: [PATCH 24/25] Removed debug picture exportanig --- src/gui/MainWindow.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/gui/MainWindow.cpp b/src/gui/MainWindow.cpp index 3ab2dc6..f7700a9 100644 --- a/src/gui/MainWindow.cpp +++ b/src/gui/MainWindow.cpp @@ -64,17 +64,17 @@ int Window::init(float* data) { if (init_quad(data)) return -1; this->last_frame = std::chrono::steady_clock::now(); - // 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; + while (!glfwWindowShouldClose(window)) { + Window::tick(); + } + // TODO: Remove this, this was just for ray-casting debug + // 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; From 340efd641da4d6e0e71eaf0c5590ce515ef6caf9 Mon Sep 17 00:00:00 2001 From: Martin Opat Date: Sat, 11 Jan 2025 15:13:23 +0100 Subject: [PATCH 25/25] Made camera not upside down --- src/consts.cu | 4 ++-- src/main.cu | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/consts.cu b/src/consts.cu index 160720c..36a2a59 100644 --- a/src/consts.cu +++ b/src/consts.cu @@ -33,10 +33,10 @@ __device__ Point3 d_lightPos; __device__ Color3 d_backgroundColor; // 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) +Point3 h_cameraPos = Point3::init(50.0f, -50.0f, -75.0f); // Camera for partially trimmed data set 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, 0.0, 1.0).normalize(); +Vec3 h_cameraUp = Vec3::init(0.0, 1.0, 0.0).normalize(); Point3 h_lightPos = Point3::init(1.5, 2.0, -1.0); Color3 h_backgroundColor = Color3::init(0.1f, 0.1f, 0.1f); diff --git a/src/main.cu b/src/main.cu index e3faac7..0ff59bb 100644 --- a/src/main.cu +++ b/src/main.cu @@ -68,7 +68,7 @@ int main() { if (data[i + 0*VOLUME_DEPTH*VOLUME_HEIGHT*VOLUME_WIDTH] + epsilon >= infty) hostVolume[i] = -infty; } - // Reverse the order of hostVolume TODO: Idk why the volume upside down, this can probably be deleted and substituted with a minus sign somehwere (camera position?) + // Reverse the order of hostVolume - why is it upside down anyway? for (int i = 0; i < VOLUME_WIDTH; i++) { for (int j = 0; j < VOLUME_HEIGHT; j++) { for (int k = 0; k < VOLUME_DEPTH/2; k++) {