spherical-volume-rendering / svr-algorithm

A spherical volume rendering algorithm that performs ray casting through a spherical voxel grid.
Other
6 stars 7 forks source link

Radial hit discriminant calculation optimization #30

Closed cgyurgyik closed 4 years ago

cgyurgyik commented 4 years ago
...
else
    r_a = r_a + delta_radius;
    discr = r_a^2 - (dot(ray_circle_vector,ray_circle_vector) - v^2);
    d = sqrt(discr);
    ta = (v-d);
    tb = (v+d);
    pa = ray_origin + ta.*ray_unit_vector;
    pb = ray_origin + tb.*ray_unit_vector;
    t1 = (pa(1) - ray_origin(1))/ray_direction(1);
    t2 = (pb(1) - ray_origin(1))/ray_direction(1);
    time_array_a(1) = t1;
    time_array_a(2) = t2;
end

discr = r_b^2 - (dot(ray_circle_vector,ray_circle_vector) - v^2);
if (discr >= 0)        
    d = sqrt(discr);
    ta = (v-d);
    tb = (v+d);
    pa = ray_origin + ta.*ray_unit_vector;
    pb = ray_origin + tb.*ray_unit_vector;
    t1 = (pa(1) - ray_origin(1))/ray_direction(1);
    t2 = (pb(1) - ray_origin(1))/ray_direction(1);
    time_array_b(1) = t1;
    time_array_b(2) = t2;
end

Just sharing some thoughts as I look through the code for this next progress report. This isn't something we need to worry about now, but it is an optimization we can consider in the future. We know if discr <= 0 there is no hit. If it is not, then we have an intersection.

We can then calculate the first discriminant for r_a, and do an early exit if either t < t1 < t_end or t < t2 < t_end. This allows for early exit if r_a time fulfills our requirement. In the current code, we always calculate for both r_a and r_b

This pseudocode would look something like this:

...
else
 r_a = r_a + delta_radius;
    discr = r_a^2 - (dot(ray_circle_vector,ray_circle_vector) - v^2);
    d = sqrt(discr);
    ta = (v-d);
    tb = (v+d);
    pa = ray_origin + ta.*ray_unit_vector;
    pb = ray_origin + tb.*ray_unit_vector;
    t1 = (pa(1) - ray_origin(1))/ray_direction(1);
    t2 = (pb(1) - ray_origin(1))/ray_direction(1);
    time_array_a(1) = t1;
    time_array_a(2) = t2;
    if t1 > t || t2 > t
        tMaxR = time_array_a(time_array_a > t);
    end
end

time_array_b = [0, 0];
if (tMaxR has not been calculated):
    discr = r_b^2 - (dot(ray_circle_vector, ray_circle_vector) - v^2);
    if (discr >= 0)        
        d = sqrt(discr);
        ta = (v-d);
        tb = (v+d);
        pa = ray_origin + ta.*ray_unit_vector;
        pb = ray_origin + tb.*ray_unit_vector;
        t1 = (pa(1) - ray_origin(1))/ray_direction(1);
        t2 = (pb(1) - ray_origin(1))/ray_direction(1);
        time_array_b(1) = t1;
        time_array_b(2) = t2;
    end

In summary, this would allow us to totally skip the calculation of the discr and other necessary equations for time_array_b.

ak-2485 commented 4 years ago

There is a really important point to make here. The statement "We know if discr <= 0 there is no hit." is not completely correct. Consider the first test case in our test suite: the value of the discriminant that we want to consider is positive; the ray does intersect the circle, just not in the time frame that we consider. In particular, the Graphics Gems algorithm is a base algorithm that we've developed our initialization phase on; the GG algorithm doesn't make any mention of time.

ak-2485 commented 4 years ago

I'm not sure if the remaining comments are relevant after that point.

cgyurgyik commented 4 years ago

Ok I posted this a while back so I'll have to look into this again