spherical-volume-rendering / svr-algorithm

A spherical volume rendering algorithm that performs ray casting through a spherical voxel grid.

Fast Voxel Traversal Algorithm Over Spherical Grids

About

This project extends the yt open-source data analysis and visualization package, providing an enhanced, integrated user interface for data exploration and enabling the visualization of physical data that comes from non-cartesian grids. Currently, yt implements a fast voxel traversal over a cartesian coordinate grid. The objective is to develop a fast voxel traversal over a spherical coordinate grid, based on ideas from Amanatides and Woo’s seminal paper on fast voxel traversal for ray tracing.

Authors

• Chris Gyurgyik (cpg49 at cornell.edu)
• Ariel Kellison (ak2485 at cornell.edu)
• Youhan Yuan (yy435 at cornell.edu)

Initial Benchmarks*

# Rays	# Voxels	CPU Mean (ms)	CPU Median (ms)	CPU Std Dev (ms)
128^2	64^3	103	103	1.1
256^2	64^3	407	406	1.3
512^2	64^3	1616	1613	2.3
128^2	128^3	201	201	0.8
256^2	128^3	796	794	2.4
512^2	128^3	3187	3180	8.3

^{*Run on (8 X 1400 MHz CPUs).
CPU Caches:
L1 Data 32 KiB (x4), L1 Instruction 32 KiB (x4), L2 Unified 256 KiB (x4), L3 Unified 6144 KiB (x1)}

C++ Build Requirements

• CMake 3.7 or later
• C++11-standard-compliant compiler

To run the benchmarks:

1. Install CMake version 3.7 or higher (https://cmake.org/)
2. Clone the repository and build the benchmarks:

git clone https://github.com/spherical-volume-rendering/svr-algorithm.git && 
cd svr-algorithm/cpp/benchmarks && mkdir build && cd build && cmake .. && make

3. Run the benchmarks:

cd .. && ./bin/benchmark_svr

C++ Example

#include "spherical_volume_rendering_util.h"

const BoundVec3 sphere_center(0.0, 0.0, 0.0);
const double sphere_max_radius = 10.0;
const svr::SphereBound min_bound = { .radial=0.0, .polar=0.0, .azimuthal=0.0 };
const svr::SphereBound max_bound = { .radial=sphere_max_radius, .polar=2*M_PI, .azimuthal=2*M_PI };
const svr::SphericalVoxelGrid grid(min_bound, max_bound, 
                                   /*num_radial_sections=*/4, 
                                   /*num_polar_sections=*/4,
                                   /*num_azimuthal_sections=*/4, 
                                   sphere_center);
const BoundVec3 ray_origin(-13.0, -13.0, -13.0);
const FreeVec3 ray_direction(1.0, 1.0, 1.0);
const Ray ray(ray_origin, ray_direction);
const auto voxels = svr::walkSphericalVolume(ray, grid, /*t_begin=*/0.0, /*t_end=*/30.0);

Cython Build Requirements

• Python3
• Cython
• Numpy
• distutils

Cython Example

#   Compile code before use:
#   python3 cython_SVR_setup.py build_ext --inplace

import cython_SVR
import numpy as np

ray_origin    = np.array([-13.0, -13.0, -13.0])
ray_direction = np.array([1.0, 1.0, 1.0])
sphere_center = np.array([0.0, 0.0, 0.0])
sphere_max_radius      = 10.0
num_radial_sections    =  4
num_polar_sections     =  4
num_azimuthal_sections =  4
min_bound = np.array([0.0, 0.0, 0.0])
max_bound = np.array([sphere_max_radius, 2 * np.pi, 2 * np.pi])
t_begin   = 0.0
t_end     = 30.0
voxels = cython_SVR.walk_spherical_volume(ray_origin, ray_direction, min_bound, max_bound, 
                                          num_radial_sections, num_polar_sections, 
                                          num_azimuthal_sections, sphere_center, t_begin, t_end)

# Expected voxels: [ [1, 2, 2], [2, 2, 2], [3, 2, 2], [4, 2, 2],
#                    [4, 0, 0], [3, 0, 0], [2, 0, 0], [1, 0, 0] ]

Project Links

• Initial Proposal
• Feasibility Study
• Progress Report
• Final Report
• Final Handover Package

References

• John Amanatides and Andrew Woo. A fast voxel traversal algorithm for ray tracing. In Eurographics ’87, pages 3–10, 1987.
• James Foley, Andries van Dam, Steven Feiner & John Hughes, "Clipping Lines" in Computer Graphics (3rd Edition) (2013)
• Paul S. Heckbert, editor. Graphics Gems IV. Academic Press Professional, Inc., USA, 1994.
• Donald. E. Knuth, 1998, Addison-Wesley Longman, Inc., ISBN 0-201-89684-2, Addison-Wesley Professional; 3rd edition.
• Joseph O'Rourke, "Search and Intersection" in Computational Geometry in C (2nd Edition) (1998)

svr-algorithm's Issues

Chris' list of fixes

Here are my goals for Thursday:

• Early exit from function if ray doesn’t go into any voxel. (#12)
• Early exit from loop if we’ve hit the maximum radius (#11)
• A few test cases to hit both general and corner cases. This will make our lives easier for implementation in other languages as well. (#13 for initial test suite, but still WIP)
• Re-name / understand jenkyR. Related to #14, #15

This will require me to edit polarCoordinateTraversal and radial_hit. I will send each of these updates as small PRs and wait for approval from one of you so if you see a new PR please take a look ASAP.

XZ plane for Azimuthal Hit.

https://github.com/spherical-volume-rendering/algorithm-team-collaboration/blob/4d4774b694a9f84e27060b31de3fc1a46adb917e/matlab/spherical-traversal/azimuthal_hit.m#L74

Should this be checking the z-direction() since we're now living in the XZ plane for azimuthal hit?

Ray collinear with voxel boundary.

TEST(SphericalCoordinateTraversal, RayTravelsAlongXAxis) {
        const BoundVec3 min_bound(0.0, 0.0, 0.0);
        const BoundVec3 max_bound(30.0, 30.0, 30.0);
        const BoundVec3 sphere_center(0.0, 0.0, 0.0);
        const double sphere_max_radius = 10.0;
        const std::size_t num_radial_sections = 4;
        const std::size_t num_angular_sections = 8;
        const std::size_t num_azimuthal_sections = 4;
        const SphericalVoxelGrid grid(min_bound, max_bound, num_radial_sections,
                                      num_angular_sections,
                                      num_azimuthal_sections, sphere_center, sphere_max_radius);
        const BoundVec3 ray_origin(-15.0, 0.0, 0.0);
        const FreeVec3 ray_direction(1.0, 0.0, 0.0);
        const Ray ray(ray_origin, ray_direction);
        const double t_begin = 0.0;
        const double t_end = 30.0;

        const auto actual_voxels = sphericalCoordinateVoxelTraversal(ray, grid, t_begin, t_end);
        const std::vector<int> expected_radial_voxels = {1,2,3,4,4,3,2,1};
        const std::vector<int> expected_theta_voxels = {4,4,4,4,5,5,5,5};
        const std::vector<int> expected_phi_voxels = {2,2,2,2,3,3,3,3};
        expectEqualVoxels(actual_voxels, expected_radial_voxels, expected_theta_voxels, expected_phi_voxels);
    }

This is the correct theta voxels:
expected_theta_voxels = {4,4,4,4,7,7,7,7};

Parallel lines are not intersections in our current implementation. In this case though, they are parallel and collinear.

We may want to consider that collinear lines are intersecting. The issue with this is that intersection will occur at EACH step of the traversal.

Calculating tStepR

https://github.com/spherical-volume-rendering/algorithm-team-collaboration/blob/8db18b6e4eff0b2cf82d3267dc3688a6a6c1c43a/polarCoordinateTraversal.m#L173

Is there a reason tStepR is calculated in the while loop rather than in radial_hit?

More tests for polar case

We need additional test cases for the polar algorithm:

1. We need test cases where the center of the gird occurs in all four quadrants, and at the origin. Currently, I think all (or most) of the test cases occur for the same grid!
2. We need test cases for very large grids - say radius is 100 and the number of radial voxels is 100; also a large number of angular voxels so that the voxel sizes are small.

Renaming files

Let's rename files; something better than the _old append. I'll PR for a radial_hit_dotprod and polar_coordinate_traversal_dotprod. For the originals let's get rid of the _old and otherwise remain the same. Thoughts?

Math for sphere intersection.

A while back, I did a primer on ray tracing. Part of this dealt with ray intersection with spheres. They use similar equations to the one we're using to calculate intersection with spheres. More information can be found here:
https://raytracing.github.io/books/RayTracingInOneWeekend.html#addingasphere

And here's the code I wrote (it was originally C++, but turned it into pseudocode since a lot of the C++ keywords will cause more confusion than anything.

    // Determines whether a ray has hit a sphere in the boundaries (minimum, maximum)
    // with the given center and radius. The formula is generated as follows:
    // -> (x - Cx)^2 + (y - Cy)^2 + (z - Cz)^2 = R^2  (sphere centered at Cx, Cy, Cz)
    // -> Since a vector from center C to a point P = (x, y, z) is (P - C)
    // -> Then, dot((p - C), (p - C)) = R^2
    // -> = dot((p(t) - C), p(t) - C)) = R^2
    // -> = dot((A + t * B - C), A + t * B - C)) = R^2
    // -> = t^2 * dot(B, B) + 2t * dot(B, A - C) + dot (A - C, A - C) - R^2 = 0
    bool hit(ray, t_min, t_max) {
        oc = ray.origin() - circle_center;
        direction = ray.direction()
        a = direction.dot(direction);
        b = direction.dot(oc);
        c = oc.dot(oc) - (radius_ * radius_);
        discriminant = (b * b) - (a * c);
        if (discriminant <= 0) return false;

        hit_point_one = (-b - sqrt(discriminant)) / a;
        if (hit_point_one > t_min && hit_point_one < t_max) {
            // Given the hit point and ray, we can also calculate where along the ray it intersects.
            return true;
        }
        hit_point_two = (-b + sqrt(discriminant)) / a;
        if (hit_point_two > t_min && hit_point_two < t_max) {
            return true;
        }
        return false;
    }

[bug] radial_hit transition to other side of the origin.

Let's say, num_radial_sections = 4;

So currently, when the ray first travels into the origin of the circle (i.e. current_radial_voxel = 4, it will update the current_radial_voxel to 4. When it then travels again out of the same voxel (when it enters the center voxel, it has to exit it as well), it updates current_radial_voxel to 5. This is incorrect. Instead, this should be 3.

So our current_voxel_r, if it were to travel left to right through the center, the current_radial_voxel
should do the following: 1 -> 2 -> 3 -> 4 -> 3 -> 2 -> 1

Currently, ours does: 1 -> 2 -> 3 -> 4 -> 5 and then stops. It stops likely because we determine the current radius boundaries by

current_radius = circle_max_radius - (delta_radius * (current_radial_voxel - 1));

This is nonsensical when the current_radial_voxel is greater than the number of voxels.

Radial hit discriminant calculation optimization

...
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.

Fix theta voxel initialization phase.

Issue found with small voxel count:

thetaVoxels =  1     2     2     2     2     2     2     2     2     2     3     
               3     3     3     3     4     4     5     6     7     7     8     
               8     8     9     9     9     9     9     9     9     9     10    
               10    10    10

The algorithm produces thetaVoxel as being 1 during initialization phase (using both old version with atan2 and new version with no trigonometric function), yet this image clearly depicts it being 2.
If I double the num_angular_sections from 20 to 40, it is off by one again.
My guess is that our current calculations are not taking into account the position along the ray at first entry. So, we need:

if abs(ray_origin - circle_center) < tol 
    % If the ray starts at the origin, we need to perturb slightly along its 
    % path to find the correct angular voxel.
    pert_t = 0.1;
    pert_x = ray_origin(1) + ray_direction_x * pert_t;
    pert_y = ray_origin(2) + ray_direction_y * pert_t;
    a = circle_center(1) - pert_x;
    b = circle_center(2) - pert_y;
else
    % Otherwise, we want the theta voxel where the ray intersects with the first radial boundary. 
    ray_circle_intersection_x = ray_origin(1) + ray_direction_x * t1; % t1 calculated in initialization phase
    ray_circle_intersection_y = ray_origin(2) + ray_direction_y * t1;
    a = circle_center(1) - ray_circle_intersection_x;
    b = circle_center(2) - ray_circle_intersection_y;
end

We then need to update t to be t = max(t_begin, t1);, so the first angular_hit does not go back.
And it works! At least for this case. We're assuming here though that it'll hit a radial boundary before it traverses a theta boundary. So what we may want to do is have the following:

if abs(ray_origin - circle_center) < tol 
    % If the ray starts at the origin, we need to perturb slightly along its 
    % path to find the correct angular voxel.
    pert_t = 0.1;
    pert_x = ray_origin(1) + ray_direction_x * pert_t;
    pert_y = ray_origin(2) + ray_direction_y * pert_t;
    a = circle_center(1) - pert_x;
    b = circle_center(2) - pert_y;
elseif the ray begins inside the circle, do what we were originally doing.
else
    % Otherwise, we want the theta voxel where the ray intersects with the first radial boundary. 
    ray_circle_intersection_x = ray_origin(1) + ray_direction_x * t1; % t1 calculated in initialization phase
    ray_circle_intersection_y = ray_origin(2) + ray_direction_y * t1;
    a = circle_center(1) - ray_circle_intersection_x;
    b = circle_center(2) - ray_circle_intersection_y;
end

Its important to note that all of our tests pass with the first code segment, which means we need to add a test that hits this case.

New PR to follow soon.

Originally posted by @cgyurgyik in #58 (comment)

3-d traversal versus visualization

Test

function SphericalCoordinateTraversalExample
   min_bound = [-20, -20.0, -20.0];
    max_bound = [20.0, 20.0, 20.0];
    ray_origin = [-13.0, -13.0, -13.0];
    ray_direction = [1.0, 1.5, 1.0];
    sphere_center = [0.0, 0.0, 0.0];
    sphere_max_radius = 10.0;
    
    num_radial_sections = 4;
    num_angular_sections = 4;
    num_azimuthal_sections = 4;
    t_begin = 0.0;
    t_end = 30.0;
    verbose = true;
    
    [rVoxels, thetaVoxels, phiVoxels] = sphericalCoordinateTraversal(min_bound, max_bound, ray_origin, ray_direction, ...
    sphere_center, sphere_max_radius, num_radial_sections, num_angular_sections, num_azimuthal_sections, t_begin, t_end, verbose);
    
    rVoxels
    thetaVoxels
    phiVoxels
end

spherePlot.py code

# max radius of the sphere
max_radius = 10
# number of radial sections
num_rad = 4
# number of angular sections
num_ang = 4
# number of azimuthal sections
num_azi = 4
# sphere center
origin_sphere = np.array([0, 0, 0])
# ray Start
origin_ray = np.array([-13, -13, -13])
# ray direction
ray_dir = np.array([1, 2, 1])
# t_begin
t_begin = 0.0
# t_end
t_end = 30.0

Ok so using the visualization, I've changed the ray direction from what our only test has, which is:

ray_direction = [1.0, 1.0, 1.0];

to

ray_direction = [1.0, 1.5, 1.0];

This produces the following voxels:

The XY, XZ Plane are pictured as follows:

The discrepancy:

• According to the pictures, the ray never traverses into the 0th theta voxel, yet according to our algorithm, it does this twice at the end.
• This means either our algorithm or the picture is incorrect.

Clarification on azimuthal and polar angle bounds

Hi! I was working on some of the previous work when I ran into something that I was unsure of. In the test_svr.cpp file, I see lots of areas where the sphere bounds are set to 0 .. 2*pi for both azimuthal and polar angles. My understanding is that we should usually be setting them to 0..pi for polar and 0..2pi for azimuth.

Is the code set up to address this somewhere internally? Or does it need to be modified somehow to take this into account?

Performance bottleneck with equality.

    inline bool isKnEqual(double a, double b) noexcept {
        const double diff = std::abs(a - b);
        if (diff <= ABS_EPSILON) { return true; }
        return diff <= std::max(std::abs(a), std::abs(b)) * REL_EPSILON;
    }

This is currently the performance bottleneck. One consideration would be to only use return std::abs(a - b) <= ABS_EPSILON, with degradation in accuracy.
FWIW, all our current tests pass using only ABS_EPSILON

[Code reduction enhancement] Removing flags

When we fix our current issues, one optimization that can occur is just setting the tMaxR or tMaxTheta to inf when no hit occurs. This guarantees that it will always be greater than its counterpart when checking:
if tMaxTheta < tMaxR: ...

It also removes a lot of redundant checks with the current flags is_radial_hit and is_angular_hit.

Our code, in practice, should guarantee that either hit will always occur with each iteration. Otherwise, one of two cases is true:

1. We have a hole in our circle
2. The ray has exited the circle and we're still iterating

Fix ray origin inside sphere, other corner cases.

Edit: This has been fixed. This will remain open until the rest of the tests are implemented.

When writing the spherical algorithm, we didn't include any tests for when the ray origin is inside the sphere. This is something that needs to be accounted for. Here is a list of tests I've compiled to verify the correctness of some corner cases in our algorithm:

1. ray origin is within the sphere in radial voxel 1.
2. ray origin is within sphere in radial voxel N, where N = num_radial_voxels (NOT directly centered at origin (0.0, 0.0, 0.0); this is covered in case 6).
3. ray origin is within sphere and within bounds [1, N] (i.e. exclusive). For example, if N = 4, the ray origin should be 2 or 3.
4. ray origin is within sphere and t_begin != 0.0
5. ray origin is within sphere and t_begin = 0.0
6. ray origin is (0.0, 0.0, 0.0)
7. ray origin begins on an angular boundary.
8. ray origin begins on an azimuthal boundary.
9. ray origin begins on a radial_boundary.
10. ray origin begins inside sphere and ends inside sphere.
11. ray origin begins outside sphere and ends inside sphere.
12. ray origin begins inside sphere and ends outside sphere.
13. ray origin is (0.0, 0.0, 0.0) and t_begin != 0.0

Tangential Hit Bug - No Progress

I believe I may have some issues with the tangential hit.

        const BoundVec3 min_bound(-200000.0, -200000.0, -200000.0);
        const BoundVec3 max_bound(200000.0, 200000.0, 200000.0);
        const BoundVec3 sphere_center(0.0, 0.0, 0.0);
        const double sphere_max_radius = 10e3;
        const std::size_t num_radial_sections = 128;
        const std::size_t num_angular_sections = 1;
        const std::size_t num_azimuthal_sections = 1;
        const svr::SphericalVoxelGrid grid(min_bound, max_bound, num_radial_sections, num_angular_sections,
                                           num_azimuthal_sections, sphere_center, sphere_max_radius);
        const double t_begin = 0.0;
        const double t_end = sphere_max_radius * 3;
        const BoundVec3 ray_origin(-421.875, -562.5, -(sphere_max_radius + 1.0));
        const FreeVec3 ray_direction(0.0, 0.0, 1.0);
        const Ray ray(ray_origin, ray_direction);
        const auto actual_voxels = sphericalCoordinateVoxelTraversal(ray, grid, t_begin, t_end);

The expected behavior for radial_voxels is 1, 2, 3, ..., X, ..., 3, 2, 1. This works fine when num_radial_sections = 64, num_radial_sections = 127, and num_radial_sections = 129. When it is changed to 128, we get the following behavior:

...
Radial Hit Entrance: < Current Voxel: 113, Current t: 8967.503394 > - {it[0]: 9063.500000, it[1]: 10938.500000, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9063.500000, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 114, Current t: 9063.500000 > - {it[0]: 9163.202711, it[1]: 10838.797289, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9163.202711, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 115, Current t: 9163.202711 > - {it[0]: 9268.122538, it[1]: 10733.877462, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9268.122538, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 116, Current t: 9268.122538 > - {it[0]: 9380.902036, it[1]: 10621.097964, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9380.902036, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 117, Current t: 9380.902036 > - {it[0]: 9506.894116, it[1]: 10495.105884, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9506.894116, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 118, Current t: 9506.894116 > - {it[0]: 9660.461020, it[1]: 10341.538980, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9660.461020, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 119, Current t: 9660.461020 > - {it[0]: 10001.000000, it[1]: 10001.000000, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 10001.000000, is_radial_transition: False, is_tangential_hit: True
Radial Hit Entrance: < Current Voxel: 119, Current t: 10001.000000 > - {it[0]: 10001.000000, it[1]: 10001.000000, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 19976.250134, is_radial_transition: False, is_tangential_hit: True
// algorithm ends

Here, I've printed out the data for the last few radial hits. The issue arises with the second to last radial hit:

Radial Hit Entrance: < Current Voxel: 119, Current t: 9660.461020 > - {it[0]: 10001.000000, it[1]: 10001.000000, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 10001.000000, is_radial_transition: False, is_tangential_hit: True
Radial Hit Entrance: < Current Voxel: 119, Current t: 10001.000000 > - {it[0]: 10001.000000, it[1]: 10001.000000, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 19976.250134, is_radial_transition: False, is_tangential_hit: True

The algorithm just ends after this. It should continue here; since the ray passes entirely through the sphere, the first and last radial voxel should be 1.

At first I thought it was because the sphere maximum radius was small and the number of radial sections large, but I don't think this is the case. Let's set num_radial_sections to 129:

...
Radial Hit Entrance: < Current Voxel: 113, Current t: 8886.419490 > - {it[0]: 8979.244490, it[1]: 11022.755510, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 8979.244490, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 114, Current t: 8979.244490 > - {it[0]: 9074.880439, it[1]: 10927.119561, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9074.880439, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 115, Current t: 9074.880439 > - {it[0]: 9174.302312, it[1]: 10827.697688, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9174.302312, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 116, Current t: 9174.302312 > - {it[0]: 9279.072578, it[1]: 10722.927422, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9279.072578, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 117, Current t: 9279.072578 > - {it[0]: 9391.945120, it[1]: 10610.054880, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9391.945120, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 118, Current t: 9391.945120 > - {it[0]: 9518.573841, it[1]: 10483.426159, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9518.573841, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 119, Current t: 9518.573841 > - {it[0]: 9674.594333, it[1]: 10327.405667, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 9674.594333, is_radial_transition: False, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 120, Current t: 9674.594333 > - {it[0]: 9674.594333, it[1]: 10327.405667, it[2]: 25.749866, it[3]: 19976.250134} - intersection_time: 10327.405667, is_radial_transition: True, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 119, Current t: 10327.405667 > - {it[0]: 9674.594333, it[1]: 10327.405667, it[2]: 9518.573841, it[3]: 10483.426159} - intersection_time: 10483.426159, is_radial_transition: True, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 118, Current t: 10483.426159 > - {it[0]: 9518.573841, it[1]: 10483.426159, it[2]: 9391.945120, it[3]: 10610.054880} - intersection_time: 10610.054880, is_radial_transition: True, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 117, Current t: 10610.054880 > - {it[0]: 9391.945120, it[1]: 10610.054880, it[2]: 9279.072578, it[3]: 10722.927422} - intersection_time: 10722.927422, is_radial_transition: True, is_tangential_hit: False
Radial Hit Entrance: < Current Voxel: 116, Current t: 10722.927422 > - {it[0]: 9279.072578, it[1]: 10722.927422, it[2]: 9174.302312, it[3]: 10827.697688} - intersection_time: 10827.697688, is_radial_transition: True, is_tangential_hit: False
...
// algorithm continues to radial voxel 1.

No issues at all.

Difference in behavior when num_polar_sections != num_azimuthal_sections

I'm working through things related to #203 and #204 and I believe I have identified the principle problem, which is that there is some behavioral difference in the initialization of center_to_polar_bound_vectors_ and center_to_azimuthal_bound_vectors_ when num_polar_sections == num_azimuthal_sections and when the converse. The more correct (appearing) behavior is when they are equal, but I'm not sure I can completely convince myself of the correctness of that. The thing that gives me the most pause is that in the case that the count is equal, I don't see why it would get the correct answer for the values for both polar and azimuthal, since they span different values.

Anyway, I'm working on it, and have a partially completed pull request in the works.

Sectored sphere: quarter hemisphere traversal.

Describe the feature request
Allow the ray to traverse only in a quarter of the hemisphere. For example, the first quadrant of a sphere. This will eventually lead to the generalized of sectored sphere traversal.

Test(s) that the feature request must pass

    TEST(SphericalCoordinateTraversal, FirstQuadrantHit) {
        const BoundVec3 sphere_center(0.0, 0.0, 0.0);
        const double sphere_max_radius = 10.0;
        const std::size_t num_radial_sections = 4;
        const std::size_t num_polar_sections = 4;
        const std::size_t num_azimuthal_sections = 8;
        const svr::SphereBound max_bound = {.radial=sphere_max_radius, .polar=M_PI/2.0, .azimuthal=M_PI/2.0};
        const svr::SphericalVoxelGrid grid(MIN_BOUND, max_bound, num_radial_sections, num_polar_sections,
                                           num_azimuthal_sections, sphere_center);
        const double t_begin = 0.0;
        const double t_end = 30.0;
        const auto actual_voxels = walkSphericalVolume(Ray(BoundVec3(13.0, 13.0, 13.0), FreeVec3(-1.0, -1.0, -1.0)),
                                                       grid, t_begin, t_end);
        const std::vector<int> expected_radial_voxels = {1, 2, 3, 4};
        const std::vector<int> expected_theta_voxels = {0, 0, 0, 0};
        const std::vector<int> expected_phi_voxels = {0, 0, 0, 0};
        expectEqualVoxels(actual_voxels, expected_radial_voxels, expected_theta_voxels, expected_phi_voxels);
    }

    TEST(SphericalCoordinateTraversal, FirstQuadrantMiss) {
        const BoundVec3 sphere_center(0.0, 0.0, 0.0);
        const double sphere_max_radius = 10.0;
        const std::size_t num_radial_sections = 4;
        const std::size_t num_polar_sections = 4;
        const std::size_t num_azimuthal_sections = 8;
        const svr::SphereBound max_bound = {.radial=sphere_max_radius, .polar=M_PI/2.0, .azimuthal=M_PI/2.0};
        const svr::SphericalVoxelGrid grid(MIN_BOUND, max_bound, num_radial_sections, num_polar_sections,
                                           num_azimuthal_sections, sphere_center);
        const double t_begin = 0.0;
        const double t_end = 30.0;
        const auto actual_voxels = walkSphericalVolume(Ray(BoundVec3(13.0, -13.0, 13.0), FreeVec3(-1.0, 1.0, -1.0)),
                                                       grid, t_begin, t_end);
        EXPECT_EQ(actual_voxels.size(), 0);
    }

Expected behavior
The ray should only traverse voxels within the quarter hemisphere. Once the ray exits the quarter hemisphere, the traversal algorithm ends and returns said voxels.

Screenshots
N/A

Additional context
This is in relation to the sectored sphere traversal discussed in #102.

Potential solutions (optional)
I implemented this very simply using radial voxels as the limitation. For example, to traverse the first quadrant, I'd set the [min radial, max radial] bounds to [1, N] respectively. Then, the ray would traverse until radial voxel == N. A limitation of this is that we want the ray to traverse even when it is within the radial voxel N; it should only stop when it reaches the planes made by the X, Y, Z axes. We could calculate the time it takes for the ray to reach such a plane, but this would need to be generalized for any sectored sphere.

Angular hit increment/decrement

Check out theta voxel values for example
ray_origin = [25, 17];
ray_direction = [-1.0, -1.0];

Optimize initializeVoxelBoundarySegments() for orthographic case.

Currently, we call initializeVoxelBoundarySegments() once per ray. If we're following an orthographic projection, we know the ray will always enter through the outermost radius, and thus the boundary segments will be the same for each ray. Then, we can call initializeVoxelBoundarySegments() once in SphericalVoxelGrid.

This is pretty simple to tackle. We always need to calculate pMaxAngular, pMaxAzimuthal; this is done in SphericalVoxelGrid. In the algorithm, we just need to check if the ray origin is outside the grid. If it is, we know pAngular == pMaxAngular (similar for azimuthal).

Corner case tests

Here is a list of corner case examples that I believe should be covered in our test suite. These tests help us verify the correctness of our code, and eliminate any discrepancies I may have produced in the cpp version.

1. ray origin is within sphere and within bounds [1, N] (i.e. exclusive). For example, if N = 4, the ray origin should be within radial voxel 2 or 3.
2. ray origin begins on an angular boundary.
3. ray origin begins on an azimuthal boundary.
4. ray origin begins on a radial_boundary within the sphere.
5. Double intersection with radial/polar voxel

Edit: One more. We want to hit the following case in angular_hit:
(is_intersect_max && is_collinear_min)) . This might be more up Ariel's alley.

updateVoxelBoundarySegments() profiling

This is called during each iteration when a radial step occurs. It updates the boundary segments for the azimuthal and angular boundaries.

Location in the traversal phase:

The function:

Profiling Call Tree:

As shown, it is taking a large chunk of time.
Some potential fixes:

• Replace grid.sphereCenter().x() with const double sphere_center_x (and similar for y, z) this MIGHT avoid extra lookups.
• Raw for loops require the use of vector::size() and [] operator No single algorithm will be able to conduct what we're doing here, but I may be able to write an algorithm that does the following:
  Given two vectors A,B transform A, B while also using components from both itself and the other.

These are just micro-optimizations at the end of the day. Thinking about changes to the implementation may be necessary to truly reduce the footprint.

Maybe some kind of table look up to reduce calculations, or since voxel sizes are the same for each radial section other than orientation, we can use that information to minimize calculations.

3-dimensional spherical coordinate case.

Currently, we've been using theta and r in the 2-dimensional case as pictured above.

With the 3-dimensional case, we'll be adding azimuthal_voxels (phi).
This will begin in a new directory named spherical_coordinate_traversal.

Floating Point Limitations

We need to discuss floating point limitations of our algorithm. For example, when a large number of radial sections exist, but the sphere max radius is not "large enough," this can cause the algorithm to break down. For example, I send a ray orthographically near the center of the sphere. Let the sphere max radius be 10e4 divided into 128 radial sections. On radial section 120, it will say the next radial hit is a tangential hit; this simply isn't true.

Transferring the SVR function to yt

The intent of this project was to use it in yt. I think now is a good time to open some dialogue on what else still needs to occur for this to happen.

Tangential hit, intersection of the same voxel twice.

Describe the bug
There should not be a double hit on the voxel in which the tangential hit occurs. This is happening since tStepR = 0 for the tangential hit, and thus the voxel intersected remains the same if the next hit after a tangential hit is a radial hit. Just changing tStepR to -1 in the case that it is a tangential hit breaks tests.

Test that reproduces the bug

% Ray with a tangential hit
function testRayTangentialHit2(testCase)
    min_bound = [-20, -20.0, -20.0];
    max_bound = [20.0, 20.0, 20.0];
    ray_origin = [-2.5, 0.0, 10.0];
    ray_direction = [0.0,0.0,-1.0];
    sphere_center = [0.0, 0.0, 0.0];
    sphere_max_radius = 10.0;

    num_radial_sections = 4;
    num_angular_sections = 1;
    num_azimuthal_sections = 1;
    t_begin = 0.0;
    t_end = 30.0;
    verbose = false;

    [rVoxels, thetaVoxels, phiVoxels, tTest, tTraversal] = sphericalCoordinateTraversal(min_bound, max_bound, ray_origin, ray_direction, ...
    sphere_center, sphere_max_radius, num_radial_sections, num_angular_sections, num_azimuthal_sections, t_begin, t_end, verbose);
    verifyEqual(testCase, rVoxels,     [1,2,3,2,1]);
    verifyEqual(testCase, thetaVoxels, [1,1,1,1,1]);
    verifyEqual(testCase, phiVoxels,   [1,1,1,1,1]);
    tRelError = (tTest-tTraversal)/(tTest^2)
end

Expected behavior
First, theta and phi voxels should be 0 as Ariel has already fixed.
Secondly, the test was initially made by Youhan as: verifyEqual(testCase, rVoxels, [1,2,3,3,2,1]);
The correct version should be:
verifyEqual(testCase, rVoxels, [1,2,3,2,1]);. A tangential hit should not cause voxel [3,1,1] to be hit twice.

Screenshots
N/A

Additional context
N/A

Potential solution(s)
A not-so-elegant fix is just to simply check to see if the next voxel being added is the same as the last one. If this is true, then don't append it.

/*before the while loop*/ 
 int previous_radial_voxel_ID, previous_theta_voxel_ID, previous_phi_voxel_ID;
 previous_radial_voxel_ID = previous_theta_voxel_ID = previous_phi_voxel_ID = -1;
...

/* in the while loop */
 if (previous_radial_voxel_ID == current_voxel_ID_r &&
     previous_theta_voxel_ID == current_voxel_ID_theta &&
     previous_phi_voxel_ID == current_voxel_ID_phi) { continue;  }
            previous_radial_voxel_ID = current_voxel_ID_r;
            previous_theta_voxel_ID = current_voxel_ID_theta;
            previous_phi_voxel_ID = current_voxel_ID_phi;
            voxels.push_back({.radial_voxel=current_voxel_ID_r,
                              .polar_voxel=current_voxel_ID_theta,
                              .azimuthal_voxel=current_voxel_ID_phi});

Furthermore, since this only occurs when the radial voxel remains the same and its solely a radial voxel hit, we need only check that the previous and current voxel ID are the same in the case of a Radial intersection.

switch (voxel_intersection) {
                case Radial: {
                    t = radial_params.tMaxR;
                    current_voxel_ID_r += radial_params.tStepR;
                    if (previous_radial_voxel_ID == current_voxel_ID_r) continue;
                    previous_radial_voxel_ID = current_voxel_ID_r;
                    break;
                }
                // All other cases.
                ...

Change angular -> polar in cpp code.

Refactor the cpp code to switch from angular -> polar. Example, angularHit -> polarHit. Angular was initially used since we used polar to describe the 2D case.

Issues when ray does not go through innermost radial voxel

With this example:

    min_bound = [0.0, 0.0];
    max_bound = [30.0, 30.0];
    
    ray_origin = [3.0, 5.0];
    ray_direction = [2.3, 1.0];
    
    circle_center = [15.0, 15.0];
    circle_max_radius = 10.0;
    num_radial_sections = 3;
    num_angular_sections = 4;
    
    t_begin = 0.0;
    t_end = 10.0;
    
    verbose = true;
    
    [rVoxels, thetaVoxels] = polarCoordinateTraversal_old(min_bound, max_bound, ray_origin, ...
        ray_direction, circle_center, circle_max_radius, ...
        num_radial_sections, num_angular_sections, t_begin, t_end, verbose);

The ray does not go through the innermost radial voxel, the function ends early before traversal ends. This is because our hot fix jenkyR only goes into action when current_radius == 0, a.k.a. when we go through the inner most radius.

In the example above,
The current radius cr is: 11.1111, and the intersection sect is

   5.9777 - 1.1366i
   5.9777 + 1.1366i

Uh-oh. We need to generalize this so that our step switches even when we don't go through this voxel.

Inside sphere corner case tests

Here is a list of corner case examples that I believe should be covered in our test suite. These tests help us verify the correctness of our code, and eliminate any discrepancies I may have produced in the cpp version.

1. ray origin is within sphere and within bounds [1, N] (i.e. exclusive). For example, if N = 4, the ray origin should be 2 or 3.
2. ray origin begins on an angular boundary.
3. ray origin begins on an azimuthal boundary.
4. ray origin begins on a radial_boundary within the sphere.
5. ray origin is (0.0, 0.0, 0.0) and t_begin != 0.0

Edit: One more. We want to hit the following case in angular_hit:
(is_intersect_max && is_collinear_min)) . This might be more up Ariel's alley.

Edge case with ray beginning at a theta boundary test.

% Ray begins at a theta boundary (and goes in orthogonal directions).
function testBeginAtThetaBoundaryInOrthogonal(testCase)
    fprintf("Ray begins at a theta boundary (and goes in orthogonal directions\n")
    min_bound = [0.0, 0.0];
    max_bound = [30.0, 30.0];
    ray_origin = [10.0, 15.0];
    ray_direction = [0.0, 7.0];
    circle_center = [15.0, 15.0];
    circle_max_radius = 10.0;
    num_radial_sections = 3;
    num_angular_sections = 4;
    t_begin = 0.0;
    t_end = 20.0;
    verbose = false;
    
    [rVoxels, thetaVoxels] = polarCoordinateTraversal(min_bound, max_bound, ray_origin, ray_direction, ...
        circle_center, circle_max_radius, num_radial_sections, num_angular_sections, t_begin, t_end, verbose);
        expected_rVoxels     = [2,1];
        expected_thetaVoxels = [1,1];
        
        verifyEqual(testCase, rVoxels, expected_rVoxels);
        verifyEqual(testCase, thetaVoxels, expected_thetaVoxels);
end

Here, we get radial_voxels: [2, 2, 1]
when we should be expecting
radial_voxels: [2, 1]

This means the theta voxel begins in theta voxel 2, and traverses into 1, which causes an extra hit.

We seem to agree that as long as we aren't intersecting with objects that aren't in the ray's path, and we aren't causing any unnecessary extension (or deletion) in the chord line, we should be ok.

AMR rendering issues

I'm running the system with an AMR dataset, and I'm getting some results that confuse me somewhat. I have a suspicion why they are doing what they are doing, but I cannot yet verify this.

To set the stage, I've loaded up a "fake" AMR dataset using yt.testing.fake_amr_ds(geometry = "spherical") and I'm rendering the field "ones".

With a random normal, [0.18617054, 0.15260178, 0.83324573], I get an image that looks like the attached. The min/max are 0 and 3.971; this is just under what I'd expect the max to be (4.0, for through-the-center).

The AMR grids are well-defined, and the total enclosed volume is correct. It looks to me, however, like it's not fully traversing some of the grids. I am not sure if this is because they need to be entered by rays twice (which is definitely possible), but I am working on figuring that out. So I will get back to you. I just wanted to give an update.

I'm starting to suspect that I'll need to use the walkSphericalVolume function to compute the order of traversal for the grids.

Early exit for initialization phase

We want to exit as early as possible if there is no ray intersection.
This is as simple as calculating the discriminant early, and determining whether
discriminant <= 0. This can be done before any square rooting occurs as well.

Redundancy in traversal phase

 elseif abs(tMaxPhi - tMaxTheta) < tol && t < tMaxPhi && tMaxPhi < t_end
        % Phi, Theta equal
        if tMaxR < tMaxPhi && t < tMaxR && ~rStepViolation
            % R min
            t = tMaxR;
            current_voxel_ID_r = current_voxel_ID_r + tStepR;
        else
            t = tMaxPhi;
            current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections);
            current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections);
        end
    elseif abs(tMaxTheta - tMaxR) < tol && t < tMaxR && tMaxR < t_end && ...
            ~rStepViolation
        % R, Theta equal
        if tMaxPhi < tMaxTheta && t < tMaxPhi
            % Phi min
            t = tMaxPhi;
            current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections);
        else
            t = tMaxTheta;
            current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections);
            current_voxel_ID_r = current_voxel_ID_r + tStepR;
        end
    elseif abs(tMaxR - tMaxPhi) < tol && t < tMaxR && tMaxR < t_end && ...
            ~rStepViolation
        % R, Phi equal
        if tMaxTheta < tMaxR && t < tMaxTheta
            % Theta min
            t = tMaxTheta;
            current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections);
        else
            t = tMaxR;
            current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections);
            current_voxel_ID_r = current_voxel_ID_r + tStepR;
        end
    else
       return;
    end    
    radial_voxels = [radial_voxels, current_voxel_ID_r];
    angular_voxels = [angular_voxels, current_voxel_ID_theta];
    azimuthal_voxels = [azimuthal_voxels, current_voxel_ID_phi];
end

This is the last part of the traversal phase, checking for equal intersection between two tMax's.
Let's look at Phi, Theta equal. In the beginning of this case you have:

       if tMaxR < tMaxPhi && t < tMaxR && ~rStepViolation
            % R min
            t = tMaxR;
            current_voxel_ID_r = current_voxel_ID_r + tStepR;
        else
            t = tMaxPhi;
            current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections);
            current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections);
        end

Isn't this redundant since we check for minimum in the beginning anyways?
I removed these in each case, and the one test we have still works. Just curious if there's something I'm missing here.

Edge case with along theta boundary test

% Ray goes along a theta boundary.
function testRayGoesAlongThetaBoundary(testCase)
    fprintf("Ray goes along a theta boundary\n")
    min_bound = [0.0, 0.0];
    max_bound = [30.0, 30.0];
    ray_origin = [25.0, 15.0];
    ray_direction = [-1.0, 0.0];
    circle_center = [15.0, 15.0];
    circle_max_radius = 10.0;
    num_radial_sections = 3;
    num_angular_sections = 4;
    t_begin = 0.0;
    t_end = 20.0;
    verbose = false;
    
    [rVoxels, thetaVoxels] = polarCoordinateTraversal(min_bound, max_bound, ray_origin, ray_direction, ...
        circle_center, circle_max_radius, num_radial_sections, num_angular_sections, t_begin, t_end, verbose);
        expected_rVoxels     = [1,2,3,3,2,1];
        expected_thetaVoxels = [0,0,0,1,1,1];
        
        verifyEqual(testCase, rVoxels, expected_rVoxels);
        verifyEqual(testCase, thetaVoxels, expected_thetaVoxels);
end

The actual_thetaVoxels is instead [0,0,0,2,2,2]. This is making a jump to 2 instead of 1 since its along the theta boundary. This issue is just to ensure we are aware in case it comes up again in the future.

Reduce std::cos, std::sin calls in SphericalVoxelGrid

Since cosine and sine are symmetric in both the x- and y- direction, we can reduce the number of calls even further. This shouldn't be too difficult, I just need some time to think about how to do so elegantly.

There is no restrictions on num_angular_voxels, num_azimuthal_voxels other than its > 0, so odd numbers also need to be considered.

Currently, this is still far from being the performance bottleneck, so I'm keeping in mind Amdahl's law for now.

Divide by zero in vec3

Current division operators allow division by zero.

Reduce trigonometric calculations in angular_hit.

During our presentation, Matt was asking about calls to expensive functions, and other optimizations we can conduct. Below I present an optimization in the angular_hit function. I confess I am hesitant to implement any optimizations yet, at least not until we have a working 3-dimensional implementation.

Right now in angular_hit, we have the following:

delta_theta = 2 * pi / num_angular_sections;
interval_theta = [current_voxel_ID_theta * delta_theta, (current_voxel_ID_theta + 1) * delta_theta];

xmin = cos(min(interval_theta));
xmax = cos(max(interval_theta));
ymin = sin(min(interval_theta));
ymax = sin(max(interval_theta));

if (single(tan(min(interval_theta))) == ray_direction(2)/ray_direction(1) ...
        && single(tan(max(interval_theta))) == ray_direction(2)/ray_direction(1))
        fprintf("parallel");
        ...
end

In a single call to radial_hit, this ends up being the following :
function | #
cos | 2 calls
sin | 2 calls
tan | 2 calls
We also calculate delta_theta (float multiplication, float division) and interval_theta (2 float multiplications).

Multiply this by N, where N is the number of angular section traversals.

We can reduce this to 1 by adding the following calculations to the initialization phase in polarCoordinateTraversal.m:

delta_theta = 2 * pi / num_angular_sections;
interval_thetas = zeros(num_angular_sections);
interval_theta_cos = zeros(num_angular_sections);
interval_theta_sin = zeros(num_angular_sections);
interval_theta_tan = zeros(num_angular_sections);

% Calculate cosine min/max, sine min/max for angular sections.
for k = 1 : num_angular_sections
    i = (k - 1) * delta_theta;  % Recall: theta goes from 0 to N - 1, but MATLAB begins at 1.
    interval_thetas(k) = i;
    interval_theta_cos(k) = cos(i);
    interval_theta_sin(k) = sin(i);
    interval_theta_tan(k) = sin(i);
end

And then, in angular_hit, we need only to reference these pre-calculated numbers:

% Note: must be careful about off-by-one error.
voxel_one = current_voxel_ID_theta + 1;
voxel_two = current_voxel_ID_theta + 2;
interval_theta = [interval_thetas(voxel_one), interval_thetas(voxel_two)];

xmin = min(interval_theta_cos(voxel_one),  interval_theta_cos(voxel_two));
xmax = max(interval_theta_cos(voxel_one),  interval_theta_cos(voxel_two));
ymin = min(interval_theta_sin(voxel_one),  interval_theta_sin(voxel_two));
ymax = max(interval_theta_sin(voxel_one),  interval_theta_sin(voxel_two));

y_over_x = ray_direction(2) / ray_direction(1);
is_parallel_min = single(min(interval_theta_tan(voxel_one), interval_theta_tan(voxel_two) )) == y_over_x;
is_parallel_max = single(max(interval_theta_tan[voxel_one], interval_theta_tan(voxel_two) )) == y_over_x;
if (is_parallel_min && is_parallel_max)
   fprintf("parallel\n");
   ...
end

This cuts out a lot of the meat in our hot loop (see exact numbers above), at the cost of a computationally more expensive start-up time. Given many rays, this will be worth it though.

With regards to atan2 (asked about in the presentation), we call this once per ray in the initialization phase as well.

Traversal over a segment of the sphere.

One thing we talked about in our meeting this week was only traversing through a "slice" of the sphere. We need to discuss how this can be implemented, and what needs to change in our current algorithm.

Spherical traversal phase

In the spherical traversal phase, when angular and azimuthal boundaries are hit at the same time, the code currently only increments/decrements the angular voxels. An example case is given in the spherical test file. According to the logic of the traversal phase, we only enter increment/decrement phi if it is a strict minimum. Furthermore, the angular/azimuthal hit functions set tMaxTheta and tMaxPhi to inf in the case when the calculated intersection time is equal to the global time. You will see this if you add relevant print statements to the above mentioned example.

Review of cpp code.

@matthewturk I pushed the current PR since we'll be working on it today.
As mentioned on Slack, any early failure will pay dividends in the future. Right now, we're still trying to translate our MATLAB version into C++, work out some kinks of the spherical traversal, and verify with testing.

Please take a look at the cpp code, and the Cython code as well, where my knowledge dips. I've been relying solely on Cython documentation for now.

Optimizing the angular hits

Currently, the angular hit functions are the computational bottlenecks with many, many floating point operations. So the question is then, how can we reduce the number of calculations, branch predictors, and complexity without affecting the soundness of the algorithm.

What I've done so far:

• Tried to reduce as many calculations as possible from per ray to per grid. See all of these in SphericalVoxelGrid.h

Singular matrix when ray is parallel to a theta boundary

Currently, we are getting a singular matrix when our ray is parallel to a theta boundary (i.e. infinite number of solutions).

We're also having issues when cos(theta) = sin(theta) at pi/4 and 3pi/4 since det = 0.
@ak-2485, Let me know if this is correct.

Potential solution:
https://www.mathworks.com/matlabcentral/answers/412513-number-of-solutions-of-a-system-of-linear-equations

If this works, then the next question is how to translate this to other langs.

current_radius floating point error

So in the following example,

  min_bound = [0.0, 0.0];
    max_bound = [30.0, 30.0];
    
    ray_origin = [3.0, 10.0];
    ray_direction = [2.3, 1.0];
    
    circle_center = [15.0, 15.0];
    circle_max_radius = 10.0;
    num_radial_sections = 6;
    num_angular_sections = 8;
    
    t_begin = 0.0;
    t_end = 10.0;
    
    verbose = true;
    
    [rVoxels, thetaVoxels] = polarCoordinateTraversal_old(min_bound, max_bound, ray_origin, ray_direction, circle_center, circle_max_radius, num_radial_sections, num_angular_sections, t_begin, t_end, verbose);

If you try to run it with our current version of polarCoordinateTraversal_old,
it will stop at the origin.

Printing the variable cr (current radius) will yield a number similar to 3.32e-34.
I presume this is some floating point error. A simple fix for this might be:

Replace
if (cr == 0), cr = current_radius^2; end
with
if (abs(cr) < eps), cr = current_radius^2; end

Let me know your thoughts either way.