• Advertisement
Sign in to follow this  

Volume Raycasting algorithm question

This topic is 3629 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

I am implementing a Volume raycaster, amd I am stuck when trying to create an intensity per pixel. I am trying to implement the model illustrated here: http://en.wikipedia.org/wiki/Volume_ray_casting. For each point in my volume's scalar field, I have a function that obtains an intensity and an opacity for a given density by mapping the density points on two 1d textures: 1 texture for opacity, and one for intensity. For each ray that intersects the volume, I calculate the entry and exit points of the ray. The purpose of the loop is to accumulate intensity along the ray from back to front, but I don't understand how to do this. As of now, the pseudocode I have is
volume v;
ray r;
color accumulatedColor = black;
double tDelta;

for (double t = exit; t > entrance; t-= tDelta)
{
    point p = r.getPoint(t)
    double d = v.getDensity(p)
    color c  = v.getColor(d)
    double o = v.getOpacity(d)
    ??? -- something with accumulated color
}
return accumulatedColor;
Does anyone have any suggestions as to how to perform this computation?

Share this post


Link to post
Share on other sites
Advertisement
It helps if your volume's boundary is simple, like an axis-aligned bounding box.

In this scenario I would detect collision between the eye ray and the bounding box. If there's a hit, then I do a second collision test to determine the exit point of the eye ray. Once you have the line that travels through the bounding box, you can set up a bunch of points along it to do your opacity/shadowing etc.

Here is the GLSL code that I use to obtain the entry and exit points for an axis-aligned bounding box.

The combined_aabb_min and max variables describe the dimensions of the bounding box.

The inf constant is used to denote a non-hit.

The *_normal constants are used to help with detecting collision. It also means that we are using an AABB, since the normals of the box are axis-aligned.


uniform vec3 combined_aabb_min;
uniform vec3 combined_aabb_max;

const float inf = 1e30;

const vec3 x_plus_normal = vec3(1, 0, 0);
const vec3 x_minus_normal = vec3(-1, 0, 0);
const vec3 y_plus_normal = vec3(0, 1, 0);
const vec3 y_minus_normal = vec3(0, -1, 0);
const vec3 z_plus_normal = vec3(0, 0, 1);
const vec3 z_minus_normal = vec3(0, 0, -1);

bool within(const in float a, const in float min, const in float max);
void ray_aabb_entry(const in vec3 start, const in vec3 normalized_dir, const in vec3 aabb_min, const in vec3 aabb_max, inout vec3 entry_point);
void ray_aabb_exit(const in vec3 start, const in vec3 normalized_dir, const in vec3 aabb_min, const in vec3 aabb_max, inout vec3 exit_point);
void ray_aabb_intersection(const in vec3 start, const in vec3 normalized_dir, const in vec3 aabb_min, const in vec3 aabb_max, inout vec3 entry_point, inout vec3 exit_point);


bool within(const in float a, const in float min, const in float max)
{
if(a >= min && a <= max)
return true;

return false;
}

void ray_aabb_entry(const in vec3 start, const in vec3 normalized_dir, const in vec3 aabb_min, const in vec3 aabb_max, inout vec3 entry_point)
{
// If within the box, the entry point is the ray origin
if(within(start.x, aabb_min.x, aabb_max.x) && within(start.y, aabb_min.y, aabb_max.y) && within(start.z, aabb_min.z, aabb_max.z))
{
entry_point = start;
return;
}

// Find which plane normal points opposite to normalized_dir

// Test against +x plane
if(start.x > aabb_max.x && 0.0 > dot(normalized_dir, x_plus_normal))
{
// Extend the unit line segment to whatever length that meets the +x plane
vec3 final_dir = normalized_dir * distance(start.x, aabb_max.x) / distance(start.x, start.x + normalized_dir.x);
vec3 potential_point = start + final_dir;

// Is the end of the extended line segment within the volume defined by the infinitely extended yz plane?
if(within(potential_point.y, aabb_min.y, aabb_max.y) && within(potential_point.z, aabb_min.z, aabb_max.z))
{
entry_point = potential_point;
return;
}
}

// Test again -x plane
else if(start.x < aabb_min.x && 0.0 > dot(normalized_dir, x_minus_normal))
{
vec3 final_dir = normalized_dir * distance(start.x, aabb_min.x) / distance(start.x, start.x + normalized_dir.x);
vec3 potential_point = start + final_dir;

if(within(potential_point.y, aabb_min.y, aabb_max.y) && within(potential_point.z, aabb_min.z, aabb_max.z))
{
entry_point = potential_point;
return;
}
}

// Test against +y plane
if(start.y > aabb_max.y && 0.0 > dot(normalized_dir, y_plus_normal))
{
vec3 final_dir = normalized_dir * distance(start.y, aabb_max.y) / distance(start.y, start.y + normalized_dir.y);
vec3 potential_point = start + final_dir;

if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.z, aabb_min.z, aabb_max.z))
{
entry_point = potential_point;
return;
}
}

// Test against -y plane
else if(start.y < aabb_min.y && 0.0 > dot(normalized_dir, y_minus_normal))
{
vec3 final_dir = normalized_dir * distance(start.y, aabb_min.y) / distance(start.y, start.y + normalized_dir.y);
vec3 potential_point = start + final_dir;

if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.z, aabb_min.z, aabb_max.z))
{
entry_point = potential_point;
return;
}
}

// Test against +z plane
if(start.z > aabb_max.z && 0.0 > dot(normalized_dir, z_plus_normal))
{
vec3 final_dir = normalized_dir * distance(start.z, aabb_max.z) / distance(start.z, start.z + normalized_dir.z);
vec3 potential_point = start + final_dir;

if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.y, aabb_min.y, aabb_max.y))
{
entry_point = potential_point;
return;
}
}

// Test against -z plane
else if(start.z < aabb_min.z && 0.0 > dot(normalized_dir, z_minus_normal))
{
vec3 final_dir = normalized_dir * distance(start.z, aabb_min.z) / distance(start.z, start.z + normalized_dir.z);
vec3 potential_point = start + final_dir;

if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.y, aabb_min.y, aabb_max.y))
{
entry_point = potential_point;
return;
}
}

// If here, no match found, and entry_point is invalidated
entry_point.x = inf;
}

void ray_aabb_exit(const in vec3 start, const in vec3 normalized_dir, const in vec3 aabb_min, const in vec3 aabb_max, inout vec3 exit_point)
{
// Find which plane normal points in the same direction as normalized_dir

// Test against +x plane
if(start.x < aabb_max.x && 0.0 < dot(normalized_dir, x_plus_normal))
{
// Extend the unit line segment to whatever length that meets the +x plane
vec3 final_dir = normalized_dir * distance(start.x, aabb_max.x) / distance(start.x, start.x + normalized_dir.x);
vec3 potential_point = start + final_dir;

// Is the end of the extended line segment within the volume defined by the infinitely extended yz plane?
if(within(potential_point.y, aabb_min.y, aabb_max.y) && within(potential_point.z, aabb_min.z, aabb_max.z))
{
exit_point = potential_point;
return;
}
}

// Test against -x plane
else if(start.x > aabb_min.x && 0.0 < dot(normalized_dir, x_minus_normal))
{
vec3 final_dir = normalized_dir * distance(start.x, aabb_min.x) / distance(start.x, start.x + normalized_dir.x);
vec3 potential_point = start + final_dir;

if(within(potential_point.y, aabb_min.y, aabb_max.y) && within(potential_point.z, aabb_min.z, aabb_max.z))
{
exit_point = potential_point;
return;
}
}

// Test against +y plane
if(start.y < aabb_max.y && 0.0 < dot(normalized_dir, y_plus_normal))
{
vec3 final_dir = normalized_dir * distance(start.y, aabb_max.y) / distance(start.y, start.y + normalized_dir.y);
vec3 potential_point = start + final_dir;

if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.z, aabb_min.z, aabb_max.z))
{
exit_point = potential_point;
return;
}
}

// Test against -y plane
else if(start.y > aabb_min.y && 0.0 < dot(normalized_dir, y_minus_normal))
{
vec3 final_dir = normalized_dir * distance(start.y, aabb_min.y) / distance(start.y, start.y + normalized_dir.y);
vec3 potential_point = start + final_dir;

if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.z, aabb_min.z, aabb_max.z))
{
exit_point = potential_point;
return;
}
}

// Test against +z plane
if(start.z < aabb_max.z && 0.0 < dot(normalized_dir, z_plus_normal))
{
vec3 final_dir = normalized_dir * distance(start.z, aabb_max.z) / distance(start.z, start.z + normalized_dir.z);
vec3 potential_point = start + final_dir;

if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.y, aabb_min.y, aabb_max.y))
{
exit_point = potential_point;
return;
}
}

// Test against -z plane
else if(start.z > aabb_min.z && 0.0 < dot(normalized_dir, z_minus_normal))
{
vec3 final_dir = normalized_dir * distance(start.z, aabb_min.z) / distance(start.z, start.z + normalized_dir.z);
vec3 potential_point = start + final_dir;

if(within(potential_point.x, aabb_min.x, aabb_max.x) && within(potential_point.y, aabb_min.y, aabb_max.y))
{
exit_point = potential_point;
return;
}
}

// If here, no match found, and exit_point is invalidated
exit_point.x = inf;
}

void ray_aabb_intersection(const in vec3 start, const in vec3 normalized_dir, const in vec3 aabb_min, const in vec3 aabb_max, inout vec3 entry_point, inout vec3 exit_point)
{
// Find potentially valid entry point
ray_aabb_entry(start, normalized_dir, aabb_min, aabb_max, entry_point);

// Is entry point valid?
if(inf != entry_point.x)
ray_aabb_exit(start, normalized_dir, aabb_min, aabb_max, exit_point);
}























To use this code, you call the function like this:


ray_aabb_intersection(camera_location, eye_ray_normal, combined_aabb_min, combined_aabb_max, entry_point, exit_point);























Basically the code takes the eye ray's starting point (the camera) and normalized direction, and then returns the entry and exit point where it runs through the AABB. If the entry_point.x member is equal to "inf", then no collision occurred (this only makes sense if the bounding box doesn't fill up the entire screen).

The eye ray's normalized direction is calculated beforehand using the image plane (this is standard raytracing stuff). Keenan Crane has some source code out on how to do quaternion Julia set raytracing using nvidia's CG -- including how to encode the entire set of eye ray directions into one texture so that the entire scene is raytraced using a single shader pass. The output can either be drawn directly to the screen, or to an off-screen texture, which allows for post-processing effects like constrast, levels, hue-brightness-saturation, etc, implemented as other shaders.

In practice, I found that using a single pass to do the entire scene can lead to unresponsive behaviour, and even shader calls that time out completely. So I cut the entire screen up into many smaller rendering regions. By running a larger number of shorter shaders, the mouse/keyboard/etc may be checked more frequently, and the risk of time outs is reduced.

I suppose I could have set up the primary function to return true/false if a collision is detected, but I was already setting entry_point.x to inf, so it seemed kind of redundant. Feel free to use/change the code in any way you wish.

Passing density data into the shader is a matter of stuffing an OpenGL 3D texture and passing it in as a uniform. The cool thing about today's shaders is that you can pass in a lot of textures at once, and return multiple textures too ("multiple render targets").

On the 6800 I believe you can return four RGBA textures, giving a total of 16 output scalars (when using floating point textures).

Some chipsets (6800 GS) do not perform linear interpolation for floating point input textures, so manual trilinear interpolation must be performed within the main shader when sampling from the density data. If you don't, your density data will be sampled using nearest neighbour lookup, and will end up looking like a giant Q*Bert level -- blocky.

[Edited by - taby on March 15, 2008 9:44:01 PM]

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement