you can see the hurricane a bit

This commit is contained in:
Robin 2025-01-19 00:26:15 +01:00
commit 9352a30bdf
16 changed files with 324 additions and 192 deletions

View File

@ -4,8 +4,7 @@
First, initialize the imGui submodule:
```bash
git submodule init imgui
git submodule update imgui
git submodule update --init
```
Then, compile using cmake:

View File

@ -4,7 +4,7 @@ stdenv.mkDerivation {
name = "cuda-raytracer";
src = ./.;
nativeBuildInputs = with cudaPackages; [cmake autoAddDriverRunpath autoPatchelfHook ];
buildInputs = with pkgs; with cudaPackages; [ libGL glfw netcdf netcdfcxx4 cuda_nvcc cuda_cudart cuda_cccl libcurand];
buildInputs = with pkgs; with cudaPackages; [ libGL glfw netcdf netcdfcxx4 cuda_nvcc cuda_cudart cuda_cccl libcurand libdecor];
postConfigure = ''
export netCDFCxx_DIR=${netcdfcxx4}/lib/cmake/netCDFCxx

2
notes.md Normal file
View File

@ -0,0 +1,2 @@
# Notes to self (probably consider for report TODO: Delete this later probably)
- Temp. "zero" is at celsius 0, but the interval of values of temps. is not symmetric so technically the color distrubution is skewed -> Maybe discuss this in the report.

View File

@ -2,16 +2,40 @@
#include "cuda_error.h"
// ----------------------- Colour mapping -----------------------
__device__ __constant__ ColorStop d_stopsPythonLike[5];
// __device__ __constant__ ColorStop d_stopsPythonLike[5];
__device__ __constant__ ColorStop d_stopsPythonLike[11];
__device__ __constant__ ColorStop d_stopsGrayscale[2];
__device__ __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_stopsPythonLike[] = {
// { 0.0f, Color3::init(0.2298057f, 0.29871797f, 0.75368315f) }, // Dark Blue
// { 0.85f, Color3::init(0.23437708f, 0.30554173f, 0.75967953f) }, // Mid Blue
// { 0.90f, Color3::init(0.27582712f, 0.36671692f, 0.81255294f) }, // White
// { 0.95f, Color3::init(0.79606387f, 0.84869321f, 0.93347147f) }, // Light Orange
// { 1.0f, Color3::init(0.70567316f, 0.01555616f, 0.15023281f) } // Red
// };
// Python "jet" colour scheme
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
{ 0.00f, Color3::init(0.00000000f, 0.00000000f, 0.50000000f) },
{ 0.82f, Color3::init(0.00000000f, 0.00000000f, 0.94563280f) },
{ 0.84f, Color3::init(0.00000000f, 0.30000000f, 1.00000000f) },
{ 0.86f, Color3::init(0.00000000f, 0.69215686f, 1.00000000f) },
{ 0.88f, Color3::init(0.16129032f, 1.00000000f, 0.80645161f) },
{ 0.90f, Color3::init(0.49019608f, 1.00000000f, 0.47754586f) },
{ 0.92f, Color3::init(0.80645161f, 1.00000000f, 0.16129032f) },
{ 0.94f, Color3::init(1.00000000f, 0.77051561f, 0.00000000f) },
{ 0.96f, Color3::init(1.00000000f, 0.40740741f, 0.00000000f) },
{ 0.98f, Color3::init(0.94563280f, 0.02977487f, 0.00000000f) },
{ 1.00f, Color3::init(0.50000000f, 0.00000000f, 0.00000000f) },
};
const ColorStop h_stopsGrayscale[] = {
@ -52,5 +76,14 @@ void copyConstantsToDevice() {
// ----------------------- TransferFunction -----------------------
__device__ float d_opacityK;
__device__ float d_sigmoidOne;
__device__ float d_sigmoidTwo;
__device__ float d_sigmoidShift;
__device__ float d_sigmoidExp;
__device__ float d_alphaAcumLimit;
__device__ int d_tfComboSelected;
__device__ int d_tfComboSelectedColor;
__device__ float d_opacityConst;
__device__ bool d_showSilhouettes;
__device__ float d_silhouettesThreshold;
// ----------------------- Raycasting -----------------------
__device__ int d_samplesPerPixel;

View File

@ -5,12 +5,6 @@
#include <cmath>
// --------------------------- Basic Constants ---------------------------
// 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 INITIAL_WINDOW_WIDTH = 1200;
const int INITIAL_WINDOW_HEIGHT = 900;
@ -18,17 +12,26 @@ const double epsilon = 1e-10f;
const double infty = 1e15f; // This value is used to represent missing values in data
// --------------------------- Dataset Constants ---------------------------
// const int VOLUME_WIDTH = 576; // lon
// const int VOLUME_WIDTH = 97; // lon
const int VOLUME_WIDTH = 57; // lon
// const int VOLUME_HEIGHT = 361; // lat
// const int VOLUME_HEIGHT = 71; // lat
const int VOLUME_HEIGHT = 121; // lat
const int VOLUME_DEPTH = 42; // lev
const float DLON = 35.0f / VOLUME_WIDTH; // 35 for current trimmed data set range
const float DLAT = 60.0f / VOLUME_HEIGHT; // 60 for current trimmed data set range
const float DLEV = 1000.0f / VOLUME_DEPTH; // 1000 from max pressure (hPa) but not sure here
const float MIN_TEMP = 210.0f;
const float MAX_TEMP = 240.0f;
const float MAX_TEMP = 293.0f;
const float MIN_SPEED = 0.0F;
const float MAX_SPEED = 14.0f;
// --------------------------- Raycasting Constants ---------------------------
const int SAMPLES_PER_PIXEL = 4;
const float alphaAcumLimit = 0.4f; // TODO: Atm, this only works with sigmoid
const float minAllowedDensity = 0.001f;
const float stepSize = 0.02f;
@ -58,16 +61,28 @@ struct ColorStop {
Color3 color;
};
// factor for the opacity function
// factor for the gradient opacity function
extern __device__ float d_opacityK;
// sigmoid function variables
extern __device__ float d_sigmoidOne;
extern __device__ float d_sigmoidTwo;
extern __device__ float d_sigmoidShift;
extern __device__ float d_sigmoidExp;
// alpha accumulation limit
extern __device__ float d_alphaAcumLimit;
// combo box index
extern __device__ int d_tfComboSelected;
extern __device__ int d_tfComboSelectedColor;
// constant opacity option
extern __device__ float d_opacityConst;
// samples per pixel
extern __device__ int d_samplesPerPixel;
// Silhouettes
extern __device__ bool d_showSilhouettes;
extern __device__ float d_silhouettesThreshold;
const int lenStopsPythonLike = 5;
const int lenStopsPythonLike = 11;
const int lenStopsGrayscale = 2;
const int lenStopsBluePurpleRed = 3;
extern __constant__ ColorStop d_stopsPythonLike[5];
extern __constant__ ColorStop d_stopsPythonLike[11];
extern __constant__ ColorStop d_stopsGrayscale[2];
extern __constant__ ColorStop d_stopsBluePurleRed[3];

View File

@ -9,16 +9,6 @@
#include "cuda_error.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;
@ -70,14 +60,6 @@ int Window::init(float* data) {
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;

View File

@ -13,11 +13,15 @@ Widget::Widget(GLFWwindow* window) {
ImGui_ImplGlfw_InitForOpenGL(window, true);
ImGui_ImplOpenGL3_Init();
this->cameraPos = Point3::init(50.0f, -50.0f, -75.0f); // Camera for partially trimmed data set
this->cameraPos = Point3::init(64.0f, 42.5f, 250.0f); // Camera for partially trimmed data set
// this->cameraPos = Point3::init(300.0f, 200.0f, -700.0f); // Camera for full data set
this->cameraDir = Point3::init(0.074f, -0.301f, -2.810f);
Vec3 h_center = Vec3::init((float)VOLUME_WIDTH/2.0f, (float)VOLUME_HEIGHT/2.0f, (float)VOLUME_DEPTH/2.0f);
this->cameraDir = (h_center - this->cameraPos).normalize();
// Vec3 h_center = Vec3::init((float)VOLUME_WIDTH/2.0f, (float)VOLUME_HEIGHT/2.0f, (float)VOLUME_DEPTH/2.0f);
// this->cameraDir = (h_center - this->cameraPos).normalize();
// Vec3 h_center = Vec3::init((float)VOLUME_WIDTH, (float)VOLUME_HEIGHT, (float)VOLUME_DEPTH);
// this->cameraDir = (h_center - this->cameraPos).normalize();
this->bgColor = Color3::init(0.1f, 0.1f, 0.1f);
this->lightPos = Point3::init(1.5, 2.0, -1.0);
@ -25,14 +29,23 @@ Widget::Widget(GLFWwindow* window) {
this->fps = (char*)malloc(512*sizeof(char));
this->paused = true;
this->renderOnce = false;
this->samplesPerPixel = 1;
this->opacityK = 0;
this->sigmoidOne = 0.5f;
this->sigmoidTwo = -250.0f;
this->sigmoidShift = 0.5f;
this->sigmoidExp = -250.0f;
this->alphaAcumLimit = 0.4f;
this->tfComboSelected = 2;
this->opacityConst = 100;
this->showSilhouettes = false;
this->silhouettesThreshold = 0.02f;
};
// TODO: can be marginally improvement by only copying changed values to device - however we're dealing with individual floats here so i dont think the benefit would be all that obvious.
// REFACTOR: should probably not have all the logic in one function; something like a list of ImplementedWidgets with each a Render() function (a la interface) would be better.
// TODO: can be marginally improved by only copying changed values to device - however we're dealing with individual floats here so i dont think the benefit would be all that obvious.
// TODO: wrap basically all ImGui calls in if statements; better form + allows for checking return values / errors.
void Widget::tick(double fps) {
if (this->renderOnce) {
this->renderOnce = false;
this->paused = true;
@ -46,9 +59,42 @@ void Widget::tick(double fps) {
float min = -1, max = 1;
ImGui::Begin("Transfer Function Controls");
ImGui::DragInt("k (log [1e-10, 1])", &this->opacityK, 1, 0, 100, "%d%%", ImGuiSliderFlags_AlwaysClamp);
ImGui::DragFloat("sigmoidOne", &this->sigmoidOne, 0.01f, 0.0f, 1.0f, "%.2f");
ImGui::InputFloat("sigmoidTwo", &this->sigmoidTwo, 10.0f, 100.0f, "%.0f");
ImGui::DragInt("Grad. exp. (log [1e-10, 1])", &this->opacityK, 1, 0, 100, "%d%%", ImGuiSliderFlags_AlwaysClamp);
ImGui::DragFloat("Sig. shift", &this->sigmoidShift, 0.01f, 0.0f, 1.0f, "%.2f");
ImGui::InputFloat("Sig. sxp", &this->sigmoidExp, 10.0f, 100.0f, "%.0f");
ImGui::DragFloat("Alpha accum. limit", &this->alphaAcumLimit, 0.01f, 0.0f, 1.0f, "%.2f");
ImGui::DragInt("Opacity const. (log [1e-5, 1])", &this->opacityConst, 1, 0, 100, "%d%%", ImGuiSliderFlags_AlwaysClamp);
// the items[] contains the entries for the combobox. The selected index is stored as an int on this->tfComboSelected
// the default entry is set in the constructor, so if you want that to be a specific entry just change it
// whatever value is selected here is available on the gpu as d_tfComboSelected.
const char* items[] = {"Opacity - gradient", "Opacity - sigmoid", "Opacity - constant", "..."};
if (ImGui::BeginCombo("Transfer function", items[this->tfComboSelected])) {
for (int n = 0; n < IM_ARRAYSIZE(items); n++) {
const bool is_selected = (this->tfComboSelected == n);
if (ImGui::Selectable(items[n], is_selected))
this->tfComboSelected = n;
if (is_selected)
ImGui::SetItemDefaultFocus();
}
ImGui::EndCombo();
}
// Same comments as above apply
const char* items2[] = {"Python-like", "BPR", "Greyscale", "..."};
if (ImGui::BeginCombo("Color map", items2[this->tfComboSelectedColor])) {
for (int n = 0; n < IM_ARRAYSIZE(items2); n++) {
const bool is_selected = (this->tfComboSelectedColor == n);
if (ImGui::Selectable(items2[n], is_selected))
this->tfComboSelectedColor = n;
if (is_selected)
ImGui::SetItemDefaultFocus();
}
ImGui::EndCombo();
}
if (ImGui::Button(this->showSilhouettes ? "Hide Silhouettes" : "Show Silhouettes")) this->showSilhouettes = !this->showSilhouettes;
ImGui::DragFloat("Silhouettes threshold", &this->silhouettesThreshold, 0.001f, 0.0f, 0.5f, "%.3f");
ImGui::End();
ImGui::Begin("Light Controls");
@ -66,6 +112,7 @@ void Widget::tick(double fps) {
}
sprintf(this->fps, "%.3f fps\n", fps);
ImGui::Text("%s", this->fps);
ImGui::DragInt("Samples per pixel", &this->samplesPerPixel, 1, 1, 16, "%d", ImGuiSliderFlags_AlwaysClamp);
ImGui::End();
ImGui::Begin("Camera Controls");
@ -91,12 +138,23 @@ void Widget::copyToDevice() {
cudaMemcpyToSymbol(&d_lightPos, &this->lightPos, sizeof(Point3));
cudaMemcpyToSymbol(&d_backgroundColor, &this->bgColor, sizeof(Color3));
cudaMemcpyToSymbol(&d_samplesPerPixel, &this->samplesPerPixel, sizeof(int));
// cudaMemcpyToSymbol(&d_opacityK, &this->opacityK, sizeof(float));
this->opacityKReal = std::pow(10.0f, (-10 + 0.1 * this->opacityK));
cudaMemcpyToSymbol(&d_opacityK, &this->opacityKReal, sizeof(float));
cudaMemcpyToSymbol(&d_sigmoidOne, &this->sigmoidOne, sizeof(float));
cudaMemcpyToSymbol(&d_sigmoidTwo, &this->sigmoidTwo, sizeof(float));
cudaMemcpyToSymbol(&d_sigmoidShift, &this->sigmoidShift, sizeof(float));
cudaMemcpyToSymbol(&d_sigmoidExp, &this->sigmoidExp, sizeof(float));
cudaMemcpyToSymbol(&d_alphaAcumLimit, &this->alphaAcumLimit, sizeof(float));
cudaMemcpyToSymbol(&d_tfComboSelected, &this->tfComboSelected, sizeof(int));
cudaMemcpyToSymbol(&d_showSilhouettes, &this->showSilhouettes, sizeof(bool));
cudaMemcpyToSymbol(&d_silhouettesThreshold, &this->silhouettesThreshold, sizeof(float));
this->opacityConstReal = std::pow(10.0f, (-5 + 0.05 * this->opacityConst));
cudaMemcpyToSymbol(&d_opacityConst, &this->opacityConstReal, sizeof(float));
cudaMemcpyToSymbol(&d_tfComboSelectedColor, &this->tfComboSelectedColor, sizeof(int));
}
Widget::~Widget() {

View File

@ -17,11 +17,19 @@ public:
bool paused;
bool renderOnce;
char* fps;
int samplesPerPixel;
int tfComboSelected;
int tfComboSelectedColor;
int opacityK;
float opacityKReal;
float sigmoidOne;
float sigmoidTwo;
float sigmoidShift;
float sigmoidExp;
float alphaAcumLimit;
int opacityConst;
float opacityConstReal;
bool showSilhouettes;
float silhouettesThreshold;
ImGuiIO io;

View File

@ -20,14 +20,14 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer, const
float accumR = 0.0f;
float accumG = 0.0f;
float accumB = 0.0f;
float accumA = 1.0f * (float)SAMPLES_PER_PIXEL;
float accumA = 1.0f * (float)d_samplesPerPixel;
// Initialize random state for ray scattering
curandState randState;
curand_init(1234, px + py * width, 0, &randState);
// Multiple samples per pixel
for (int s = 0; s < SAMPLES_PER_PIXEL; s++) {
for (int s = 0; s < d_samplesPerPixel; s++) {
// Map to [-1, 1]
float jitterU = (curand_uniform(&randState) - 0.5f) / width;
float jitterV = (curand_uniform(&randState) - 0.5f) / height;
@ -71,11 +71,11 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer, const
intersectAxis(d_cameraPos.z, rayDir.z, 0.0f, (float)VOLUME_DEPTH);
if (tNear > tFar) {
// No intersection -> Set to brackground color (multiply by SAMPLES_PER_PIXEL because we divide by it later)
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;
// No intersection -> Set to brackground color (multiply by d_samplesPerPixel because we divide by it later)
accumR = d_backgroundColor.x * (float)d_samplesPerPixel;
accumG = d_backgroundColor.y * (float)d_samplesPerPixel;
accumB = d_backgroundColor.z * (float)d_samplesPerPixel;
accumA = 1.0f * (float)d_samplesPerPixel;
} else {
if (tNear < 0.0f) tNear = 0.0f;
@ -84,7 +84,7 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer, const
float alphaAccum = 0.0f;
float t = tNear; // Front to back
while (t < tFar && alphaAccum < alphaAcumLimit) {
while (t < tFar && alphaAccum < d_alphaAcumLimit) {
Point3 pos = d_cameraPos + rayDir * t;
// Convert to volume indices
@ -92,13 +92,13 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer, const
int iy = (int)roundf(pos.y);
int iz = (int)roundf(pos.z);
// Sample (pick appropriate method based on volume size) TODO: Add a way to pick this in GUI
// Sample (pick appropriate method based on volume size) TODO: Consider adding 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);
Vec3 grad = computeGradient(volumeData, VOLUME_WIDTH, VOLUME_HEIGHT, VOLUME_DEPTH, pos.x, pos.y, pos.z);
float4 color = transferFunction(density, grad, pos, rayDir); // This already returns the alpha-weighted color
//Accumulate color, and alpha
@ -130,10 +130,10 @@ __global__ void raycastKernel(float* volumeData, FrameBuffer framebuffer, const
// Average samples
accumR /= (float)SAMPLES_PER_PIXEL;
accumG /= (float)SAMPLES_PER_PIXEL;
accumB /= (float)SAMPLES_PER_PIXEL;
accumA /= (float)SAMPLES_PER_PIXEL;
accumR /= (float)d_samplesPerPixel;
accumG /= (float)d_samplesPerPixel;
accumB /= (float)d_samplesPerPixel;
accumA /= (float)d_samplesPerPixel;
// Final colour
framebuffer.writePixel(px, py, accumR, accumG, accumB, accumA);
@ -156,7 +156,6 @@ Raycaster::Raycaster(cudaGraphicsResource_t resources, int w, int h, float* data
void Raycaster::render() {
std::cout << "hello???\n";
check_cuda_errors(cudaGetLastError());
check_cuda_errors(cudaGraphicsMapResources(1, &this->resources));
check_cuda_errors(cudaGraphicsResourceGetMappedPointer((void**)&(this->fb->buffer), &(this->fb->buffer_size), resources));

View File

@ -4,64 +4,16 @@
#include <stdio.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, volH - 1);
int iy1 = min(iy + 1, volW - 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) {
__device__ float opacityFromGradient(const Vec3 &grad, const Vec3& rayDir) {
float gradMag = grad.length();
// float gradMag = grad.length()*(1-fabs(grad.normalize().dot(rayDir))); // Alternative, but not particularly better
float alpha = 1.0f - expf(-d_opacityK * gradMag);
return alpha;
}
__device__ float opacitySigmoid(float val) {
return 1.0f / (1.0f + expf(d_sigmoidTwo * (val - d_sigmoidOne)));
return 1.0f / (1.0f + expf(d_sigmoidExp * (val - d_sigmoidShift)));
}
__device__ Color3 colorMap(float normalizedValues, const ColorStop stops[], int N) {
@ -88,31 +40,63 @@ __device__ Color3 colorMap(float normalizedValues, const ColorStop stops[], int
__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);
// TODO: Somehow pick if to use temp of speed normalization ... or pass extremas as params. <-If we decide to visualize more than 1 type of data
// float normDensity = (density - MIN_TEMP) / (MAX_TEMP - MIN_TEMP);
float normDensity = (density - 273) / (MAX_TEMP - MIN_TEMP)+16.f/21.f; // Make zero match Celsius zero
// 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);
// Pick color map
Color3 baseColor;
switch (d_tfComboSelectedColor) {
case 0:
baseColor = colorMap(normDensity, d_stopsPythonLike, lenStopsPythonLike);
break;
// TODO: This is a Gui select element
// TODO: Add a way to pick different function for alpha
float alpha = opacityFromGradient(grad);
// alpha = 0.1f;
alpha = opacitySigmoid(normDensity);
// alpha = (1.0f - fabs(grad.normalize().dot(rayDir.normalize()))) * 0.8f + 0.2f;
case 1:
baseColor = colorMap(normDensity, d_stopsBluePurleRed, lenStopsBluePurpleRed);
break;
float alphaSample = density * alpha * 0.1;
case 2:
baseColor = colorMap(normDensity, d_stopsGrayscale, lenStopsGrayscale);
break;
default:
baseColor = colorMap(normDensity, d_stopsPythonLike, lenStopsPythonLike);
break;
}
// Pick opacity function
float alpha;
switch (d_tfComboSelected) {
case 0:
alpha = opacityFromGradient(grad, rayDir);
break;
case 1:
alpha = opacitySigmoid(normDensity);
break;
case 2:
alpha = d_opacityConst;
break;
default:
alpha = 1.0f; // This should not be reached anyway
break;
}
float alphaSample = density * alpha * 0.1; // TODO: Why is this still 0.1?
// --------------------------- Shading ---------------------------
// Apply Phong
Vec3 normal = -grad.normalize();
Vec3 lightDir = (d_lightPos - pos).normalize();
Vec3 viewDir = -rayDir.normalize();
Vec3 shadedColor = phongShading(normal, lightDir, viewDir, baseColor);
Vec3 shadedColor = phongShading(normal, lightDir, viewDir, baseColor); // TODO: Check if still pixelated
// Compose
float4 result;
@ -122,15 +106,13 @@ __device__ float4 transferFunction(float density, const Vec3& grad, const Point3
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
// TODO: Add a way to adjust the treshold (0.2f atm)
// 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;
// }
Vec3 N = grad.normalize();
if (d_showSilhouettes && grad.length() > 0.2f && fabs(N.dot(viewDir)) < d_silhouettesThreshold) {
result.x = 0.0f;
result.y = 0.0f;
result.z = 0.0f;
result.w = alpha;
}
return result;
}

View File

@ -7,12 +7,6 @@
// --------------------------- 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);

View File

@ -3,7 +3,7 @@
#include <fstream>
void saveImage(const char* filename, unsigned char* framebuffer, int width, int height) { // TODO: Figure out a better way to do this
void saveImage(const char* filename, unsigned char* framebuffer, int width, int height) {
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++) {

View File

@ -1,29 +1,86 @@
#include "mat.h"
#include "consts.h"
#include <vector>
#include <algorithm>
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.
// 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) {
// For boundary voxels - clamp to the boundary.
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 xm = max(x - 1, 0);
int xp = min(x + 1, volH - 1);
int ym = max(y - 1, 0);
int yp = min(y + 1, volW - 1);
int zm = max(z - 1, 0);
int zp = min(z + 1, volD - 1);
// x <-> height, y <-> width, z <-> depth
int idx = vz * volW * volH + vx * volW + vy;
return volumeData[idx];
}
// Note: Assuming data is linearized (idx = z * W * H + x * W + y;) TODO: Unlinearize if data not linear
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 ];
// 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);
return Vec3::init(gx, gy, gz);
// Clamp indices to valid range
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);
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__ Vec3 computeGradient(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz) {
// Compute gradient using central differencing with trilinear interpolation
float hx = DLAT; // x => height => lat
float hy = DLON; // y => width => lon
float hz = DLEV; // z => depth => alt
// Default
float dfdx = (sampleVolumeTrilinear(volumeData, volW, volH, volD, fx + hx, fy, fz) -
sampleVolumeTrilinear(volumeData, volW, volH, volD, fx - hx, fy, fz)) / (2.0f * hx);
float dfdy = (sampleVolumeTrilinear(volumeData, volW, volH, volD, fx, fy + hy, fz) -
sampleVolumeTrilinear(volumeData, volW, volH, volD, fx, fy - hy, fz)) / (2.0f * hy);
float dfdz = (sampleVolumeTrilinear(volumeData, volW, volH, volD, fx, fy, fz + hz) -
sampleVolumeTrilinear(volumeData, volW, volH, volD, fx, fy, fz - hz)) / (2.0f * hz);
// // DEBUG (TODO: Delete) - Back to nearest
// float dfdx = (sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx + 1), (int)roundf(fy), (int)roundf(fz)) -
// sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx - 1), (int)roundf(fy), (int)roundf(fz))) / (2.0f * hx);
// float dfdy = (sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx), (int)roundf(fy + 1), (int)roundf(fz)) -
// sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx), (int)roundf(fy - 1), (int)roundf(fz))) / (2.0f * hy);
// float dfdz = (sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx), (int)roundf(fy), (int)roundf(fz + 1)) -
// sampleVolumeNearest(volumeData, volW, volH, volD, (int)roundf(fx), (int)roundf(fy), (int)roundf(fz - 1))) / (2.0f * hz);
return Vec3::init(dfdx, dfdy, dfdz);
};
// TESTING: haven't tested this function at all tbh
@ -35,7 +92,7 @@ __device__ unsigned int packUnorm4x8(float r, float g, float b, float a) {
float len = sqrtf(r*r + g*g + b*b + a*a);
// This is a Vec4 but i can't be bothered to make that its own struct/class; FIXME: maybe do that if we need to?
// This is a Vec4 but i can't be bothered to make that its own struct/class; FIXME: maybe do that if we need to? From Martin: We could use a Vec4 for rgba too, but I don't feel like it either
u.in[0] = round(r/len * 255.0f);
u.in[1] = round(g/len * 255.0f);
u.in[2] = round(b/len * 255.0f);

View File

@ -4,7 +4,10 @@
#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);
__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);
__device__ Vec3 computeGradient(float* volumeData, const int volW, const int volH, const int volD, float fx, float fy, float fz);
__device__ unsigned int packUnorm4x8(float r, float g, float b, float a);

View File

@ -3,7 +3,7 @@
#include <cuda_runtime.h>
#include <cmath>
struct Vec3 { // TODO: Maybe make this into a class ... maybe
struct Vec3 {
double x, y, z;
static __host__ __device__ Vec3 init(double x, double y, double z) {Vec3 v = {x, y, z}; return v;}

View File

@ -21,7 +21,7 @@ static float* d_volume = nullptr;
// * time controls - arbitrary skipping to specified point (would require some changes to gpubuffer) (could have)
void getTemperature(std::vector<float>& temperatureData, int idx = 0) {
std::string path = "../../../hurricane-sandy-data";
std::string path = "../../../../hurricane-sandy-data/trimmed";
// std::string path = "data";
std::string variable = "T";
DataReader dataReader(path, variable);
@ -31,7 +31,7 @@ void getTemperature(std::vector<float>& temperatureData, int idx = 0) {
}
void getSpeed(std::vector<float>& speedData, int idx = 0) {
std::string path = "../../../hurricane-sandy-data";
std::string path = "../../../../hurricane-sandy-data/trimmed";
// std::string path = "data";
std::string varU = "U";
std::string varV = "V";
@ -55,7 +55,7 @@ void getSpeed(std::vector<float>& speedData, int idx = 0) {
int main() {
std::vector<float> data;
getTemperature(data, 0);
getTemperature(data, 7); // 20121028
// getSpeed(data, 294);
std::cout << "DATA size: " << data.size() << std::endl;
@ -79,25 +79,25 @@ int main() {
}
}
// 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++) {
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();
// // 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++) {
// 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();
// 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;
return a < b;
});
float maxVal = *std::max_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH);
std::cout << "minVal: " << minVal << " maxVal: " << maxVal << std::endl;
// // 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;
// return a < b;
// });
// float maxVal = *std::max_element(hostVolume, hostVolume + VOLUME_WIDTH * VOLUME_HEIGHT * VOLUME_DEPTH);
// std::cout << "minVal: " << minVal << " maxVal: " << maxVal << std::endl;
// // 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);