hello, I’ve just started learning OpenGL.
- According to an OpenGL tutorial: [GPU-accelerated single-pass volumetric raycasting in Qt and OpenGL], I implemented 3D visualization of weather radar data using volume raycasting technique. However, the rendered effect consists of regularly arranged pillars instead of continuous shapes representing meteorological echoes.
-
I don’t know if it’s a data issue or if there’s a problem with OpenGL settings.I will provide the relevant OpenGL code along with how I generated the data.
-
The following code is mostly sourced from open-source projects: volume-raycasting, with only some modifications made by me according to the requirements.
The .frag GLSL code:
#version 130
out vec4 a_colour;
uniform mat4 ViewMatrix;
uniform mat3 NormalMatrix;
uniform float focal_length;
uniform float aspect_ratio;
uniform vec2 viewport_size;
uniform vec3 ray_origin;
uniform vec3 top;
uniform vec3 bottom;
uniform vec3 background_colour;
uniform vec3 material_colour;
uniform vec3 light_position;
uniform float step_length;
uniform float threshold;
uniform sampler3D volume;
uniform sampler2D jitter;
uniform float gamma;
uniform vec3 colors[255];// Rainbow_r, new added
// Ray
struct Ray {
vec3 origin;
vec3 direction;
};
// Axis-aligned bounding box
struct AABB {
vec3 top;
vec3 bottom;
};
// Estimate normal from a finite difference approximation of the gradient
vec3 normal(vec3 position, float intensity)
{
float d = step_length;
float dx = texture(volume, position + vec3(d,0,0)).r - intensity;
float dy = texture(volume, position + vec3(0,d,0)).r - intensity;
float dz = texture(volume, position + vec3(0,0,d)).r - intensity;
return -normalize(NormalMatrix * vec3(dx, dy, dz));
}
// Slab method for ray-box intersection
void ray_box_intersection(Ray ray, AABB box, out float t_0, out float t_1)
{
vec3 direction_inv = 1.0 / ray.direction;
vec3 t_top = direction_inv * (box.top - ray.origin);
vec3 t_bottom = direction_inv * (box.bottom - ray.origin);
vec3 t_min = min(t_top, t_bottom);
vec2 t = max(t_min.xx, t_min.yz);
t_0 = max(0.0, max(t.x, t.y));
vec3 t_max = max(t_top, t_bottom);
t = min(t_max.xx, t_max.yz);
t_1 = min(t.x, t.y);
}
// A very simple colour transfer function
vec4 colour_transfer(float intensity)
{
vec3 high = vec3(1.0, 0.0, 0.0);
vec3 low = vec3(0.0, 0.0, 0.0);
vec3 color_new = vec3(0.0, 0.0, 0.0);
float index = intensity * 255;
if(index < 255.){
color_new = colors[int(index)];
}
float alpha = intensity;
if (intensity < 0.1){
alpha = 0;
}
if (intensity > 0.8){
alpha = 1;
}
return vec4(color_new, alpha);
}
void main()
{
vec3 ray_direction;
ray_direction.xy = 2.0 * gl_FragCoord.xy / viewport_size - 1.0;
ray_direction.x *= aspect_ratio;
ray_direction.z = -focal_length;
ray_direction = (vec4(ray_direction, 0) * ViewMatrix).xyz;
float t_0, t_1;
Ray casting_ray = Ray(ray_origin, ray_direction);
AABB bounding_box = AABB(top, bottom);
ray_box_intersection(casting_ray, bounding_box, t_0, t_1);
vec3 ray_start = (ray_origin + ray_direction * t_0 - bottom) / (top - bottom);
vec3 ray_stop = (ray_origin + ray_direction * t_1 - bottom) / (top - bottom);
vec3 ray = ray_stop - ray_start;
float ray_length = length(ray);
vec3 step_vector = step_length * ray / ray_length;
// Random jitter
ray_start += step_vector * texture(jitter, gl_FragCoord.xy / viewport_size).r;
vec3 position = ray_start;
vec4 colour = vec4(0.0);
// Ray march until reaching the end of the volume, or colour saturation
while (ray_length > 0 && colour.a < 1.0) {
float intensity = texture(volume, position).r;
vec4 c = colour_transfer(intensity);
// Alpha-blending
colour.rgb = c.a * c.rgb + (1 - c.a) * colour.a * colour.rgb;
colour.a = c.a + (1 - c.a) * colour.a;
ray_length -= step_length;
position += step_vector;
}
// Blend background
colour.rgb = colour.a * colour.rgb + (1 - colour.a) * pow(background_colour, vec3(gamma)).rgb;
// Gamma correction
a_colour.rgb = pow(colour.rgb, vec3(1.0 / gamma));
a_colour.a = colour.a;
}
The texture binding code:
glDeleteTextures(1, &m_volume_texture);
glGenTextures(1, &m_volume_texture);
glBindTexture(GL_TEXTURE_3D, m_volume_texture);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_BASE_LEVEL, 0);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
glPixelStorei(GL_UNPACK_ALIGNMENT, 1); // The array on the host has 1 byte alignment
glTexImage3D(GL_TEXTURE_3D, 0, GL_R8, m_size.x(), m_size.y(), m_size.z(), 0, GL_RED, GL_UNSIGNED_BYTE, data.data());
glBindTexture(GL_TEXTURE_3D, 0);
The following is the raycasting para:
m_shaders[shader]->setUniformValue("ViewMatrix", m_viewMatrix);
m_shaders[shader]->setUniformValue("ModelViewProjectionMatrix", m_modelViewProjectionMatrix);
m_shaders[shader]->setUniformValue("NormalMatrix", m_normalMatrix);
m_shaders[shader]->setUniformValue("aspect_ratio", m_aspectRatio);
m_shaders[shader]->setUniformValue("focal_length", m_focalLength);
m_shaders[shader]->setUniformValue("viewport_size", m_viewportSize);
m_shaders[shader]->setUniformValue("ray_origin", m_rayOrigin);
m_shaders[shader]->setUniformValue("top", m_raycasting_volume->top());
m_shaders[shader]->setUniformValue("bottom", m_raycasting_volume->bottom());
m_shaders[shader]->setUniformValue("background_colour", to_vector3d(m_background));
m_shaders[shader]->setUniformValue("light_position", m_lightPosition);
m_shaders[shader]->setUniformValue("material_colour", m_diffuseMaterial);
m_shaders[shader]->setUniformValue("step_length", m_stepLength);
m_shaders[shader]->setUniformValue("threshold", m_threshold);
m_shaders[shader]->setUniformValue("gamma", m_gamma);
m_shaders[shader]->setUniformValue("volume", 0);
m_shaders[shader]->setUniformValue("jitter", 1);
m_shaders[shader]->setUniformValueArray("colors", m_rainbow_colormap.data(), m_rainbow_colormap.size()); //new added
glClearColor(m_background.redF(), m_background.greenF(), m_background.blueF(), m_background.alphaF());
glClear(GL_COLOR_BUFFER_BIT);
m_raycasting_volume->paint();
the mesh code:
// Generates and populates a VBO for the vertices
glDeleteBuffers(1, &(m_vertex_VBO));
glGenBuffers(1, &(m_vertex_VBO));
glBindBuffer(GL_ARRAY_BUFFER, m_vertex_VBO);
glBufferData(GL_ARRAY_BUFFER, vertices.size() * sizeof(vertices[0]), vertices.data(), GL_STATIC_DRAW);
// Generates and populates a VBO for the element indices
glDeleteBuffers(1, &(m_index_VBO));
glGenBuffers(1, &(m_index_VBO));
glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, m_index_VBO);
glBufferData(GL_ELEMENT_ARRAY_BUFFER, indices.size() * sizeof(indices[0]), indices.data(), GL_STATIC_DRAW);
// Creates a vertex array object (VAO) for drawing the mesh
glDeleteVertexArrays(1, &(m_vao));
glGenVertexArrays(1, &(m_vao));
glBindVertexArray(m_vao);
glBindBuffer(GL_ARRAY_BUFFER, m_vertex_VBO);
glEnableVertexAttribArray(POSITION);
glVertexAttribPointer(POSITION, 3, GL_FLOAT, GL_FALSE, 0, nullptr);
glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, m_index_VBO);
glBindVertexArray(0);
The header of the VTK data file I created:
# vtk DataFile Version 5.1
vtk output
BINARY
DATASET STRUCTURED_POINTS
DIMENSIONS 65 501 501
SPACING 1 1 1
ORIGIN 0 0 0
POINT_DATA 16315065
SCALARS scalars float
LOOKUP_TABLE default
And here is the Python code I used to create the VTK file：
import vtk
import netCDF4
import numpy as np
nc_file = r"D:\pyCode\xiongan\code\2023072108_merge.nc"
nc_ds = netCDF4.Dataset(nc_file)
ref = np.array(nc_ds['ref'][:])
data_matrix = np.zeros(ref.shape, dtype=np.float64, order='F')
data_matrix = ref
Nx, Ny, Nz = data_matrix.shape
dataImporter = vtk.vtkImageImport()
data_string = data_matrix.tostring()
dataImporter.CopyImportVoidPointer(data_string, data_matrix.nbytes)
dataImporter.SetDataScalarTypeToFloat()
dataImporter.SetNumberOfScalarComponents(1)
dataImporter.SetDataExtent(0, Nx-1, 0, Ny-1, 0, Nz-1)
dataImporter.SetWholeExtent(0, Nx-1, 0, Ny-1, 0, Nz-1)
spacing = (0, 0, 0)
dataImporter.SetDataSpacing(spacing)
writer = vtk.vtkDataSetWriter()
writer.SetInputConnection(dataImporter.GetOutputPort())
writer.SetFileName('ref_2023072108.vtk')
writer.SetFileTypeToBinary()
writer.Update()
writer.Write()
Could you give me any advice to help me towards my goal? Thank you very much！