# Why are the 3D radar images I draw independent discrete pillars rather than continuous?

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);

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！