Unitize time, ray direction. Add CI tests for when ray is within the sphere.

[cgyurgyik]

Description

The yt version of walk_volume does not have a t_begin, and unitizes time and the ray direction. Our functionality should be similar. Doing so will require some refactoring.

Changes

• Change t_end to max_t.
• Set cythonized version max_t default value to 1.0
• Remove t_begin. This is not used in yt's walk_volume, and its removal allows for the removal of unnecessary calculations. If one wants the ray to start at a different t_begin, the origin can be changed.
• Unitize the ray direction.
• Unitize the ray time. This is done by using max_t. The documentation changes in the following way:

  // max_t is the unitized time for which the ray may travel. It is used in the following
  // linear function: ray_travel_duration = time_to_entrance + sphere.diameter() * max_t
  // If the ray origin is within the sphere, then time_to_entrance is 0. Its
  // expected values are within bounds [0.0, 1.0]. For example, if max_t <= 0.0,
  // then no voxels will be traversed. If max_t >= 1.0, then the entire sphere
  // will be traversed.

• Added CI tests with the ray placed randomly inside the sphere. Conducts similar basic property checks, except for those that require the ray to traverse the entire sphere.
• Removed RadialHitMetadata. It is simpler and easier to read to just provide a boolean value named is_radial_hit_transition. This is then passed by reference into radialHit().
• Reduced number of calls to svr::isEqual() in minimumIntersection(). Now, this is only called a maximum of 3 times rather than (worst case) 8 times. I also save previous floating point comparison in variables to avoid calling the same comparison twice.
• Reduce number of calls to svr::isEqual() in angularHit(). Since svr::lessThan() checks to see if !svr::isEqual(), it is better to just save the value of svr::isEqual() in a variable and re-use it later.
• Remove .within_bounds from HitParameter. We can simply exit early if none of the values are within bounds.
• Add a constexpr value DOUBLE_MAX which represents the maximum value a double can be. This is the value tMax will be set to if no intersection occurs.

Type of change

• [X] New feature (non-breaking change which adds functionality)
• [X] This change requires a documentation update

How Has This Been Tested?

Below are a list of tests for a different max_t values. We can also get simple First Quadrant tests passing simply by setting max_t = 0.5. This does not mean that sectored sphere traversal works, just a use case of the changes in this PR. All tests listed below are also added for the cythonized code.

TEST(SphericalCoordinateTraversal, RayOutsideSphereAndMaxTGreaterThanOne) {
  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 = 4;
  const svr::SphereBound max_bound = {
      .radial = sphere_max_radius, .polar = TAU, .azimuthal = TAU};
  const svr::SphericalVoxelGrid grid(MIN_BOUND, max_bound, num_radial_sections,
                                     num_polar_sections, num_azimuthal_sections,
                                     sphere_center);
  const BoundVec3 ray_origin(-13.0, -13.0, -13.0);
  const UnitVec3 ray_direction(1.0, 1.0, 1.0);
  const Ray ray(ray_origin, ray_direction);
  const auto actual_voxels = walkSphericalVolume(ray, grid, /*max_t=*/10.0);
  const std::vector<int> expected_radial_voxels = {1, 2, 3, 4, 4, 3, 2, 1};
  const std::vector<int> expected_theta_voxels = {2, 2, 2, 2, 0, 0, 0, 0};
  const std::vector<int> expected_phi_voxels = {2, 2, 2, 2, 0, 0, 0, 0};
  verifyEqualVoxels(actual_voxels, expected_radial_voxels,
                    expected_theta_voxels, expected_phi_voxels);
}

TEST(SphericalCoordinateTraversal, RayInsideSphereAndMaxTGreaterThanOne) {
  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 = 4;
  const svr::SphereBound max_bound = {
      .radial = sphere_max_radius, .polar = TAU, .azimuthal = TAU};
  const svr::SphericalVoxelGrid grid(MIN_BOUND, max_bound, num_radial_sections,
                                     num_polar_sections, num_azimuthal_sections,
                                     sphere_center);
  const BoundVec3 ray_origin(0.0, 0.0, 0.0);
  const UnitVec3 ray_direction(1.0, 1.0, 1.0);
  const Ray ray(ray_origin, ray_direction);
  const auto actual_voxels = walkSphericalVolume(ray, grid, /*max_t=*/10.0);
  const std::vector<int> expected_radial_voxels = {4, 3, 2, 1};
  const std::vector<int> expected_theta_voxels = {0, 0, 0, 0};
  const std::vector<int> expected_phi_voxels = {0, 0, 0, 0};
  verifyEqualVoxels(actual_voxels, expected_radial_voxels,
                    expected_theta_voxels, expected_phi_voxels);
}

TEST(SphericalCoordinateTraversal, MaxTHalvedAndRayOutsideSphere) {
  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 = 4;
  const svr::SphereBound max_bound = {
      .radial = sphere_max_radius, .polar = TAU, .azimuthal = TAU};
  const svr::SphericalVoxelGrid grid(MIN_BOUND, max_bound, num_radial_sections,
                                     num_polar_sections, num_azimuthal_sections,
                                     sphere_center);
  const BoundVec3 ray_origin(-13.0, -13.0, -13.0);
  const UnitVec3 ray_direction(1.0, 1.0, 1.0);
  const Ray ray(ray_origin, ray_direction);
  const auto actual_voxels = walkSphericalVolume(ray, grid, /*max_t=*/0.5);
  const std::vector<int> expected_radial_voxels = {1, 2, 3, 4, 4};
  const std::vector<int> expected_theta_voxels = {2, 2, 2, 2, 0};
  const std::vector<int> expected_phi_voxels = {2, 2, 2, 2, 0};
  verifyEqualVoxels(actual_voxels, expected_radial_voxels,
                    expected_theta_voxels, expected_phi_voxels);
}

TEST(SphericalCoordinateTraversal, MaxTHalvedAndRayInsideSphere) {
  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 = 4;
  const svr::SphereBound max_bound = {
      .radial = sphere_max_radius, .polar = TAU, .azimuthal = TAU};
  const svr::SphericalVoxelGrid grid(MIN_BOUND, max_bound, num_radial_sections,
                                     num_polar_sections, num_azimuthal_sections,
                                     sphere_center);
  const BoundVec3 ray_origin(0.0, 0.0, 0.0);
  const UnitVec3 ray_direction(1.0, 1.0, 1.0);
  const Ray ray(ray_origin, ray_direction);
  const auto actual_voxels = walkSphericalVolume(ray, grid, /*max_t=*/0.5);
  const std::vector<int> expected_radial_voxels = {4, 3, 2, 1};
  const std::vector<int> expected_theta_voxels = {0, 0, 0, 0};
  const std::vector<int> expected_phi_voxels = {0, 0, 0, 0};
  verifyEqualVoxels(actual_voxels, expected_radial_voxels,
                    expected_theta_voxels, expected_phi_voxels);
}

TEST(SphericalCoordinateTraversal, MaxTAtOrLessThanZero) {
  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 = 4;
  const svr::SphereBound max_bound = {
      .radial = sphere_max_radius, .polar = TAU, .azimuthal = TAU};
  const svr::SphericalVoxelGrid grid(MIN_BOUND, max_bound, num_radial_sections,
                                     num_polar_sections, num_azimuthal_sections,
                                     sphere_center);
  const BoundVec3 ray_origin(0.0, 0.0, 0.0);
  const UnitVec3 ray_direction(1.0, 1.0, 1.0);
  const Ray ray(ray_origin, ray_direction);
  const auto v1 = walkSphericalVolume(ray, grid, /*max_t=*/0.0);
  const auto v2 = walkSphericalVolume(ray, grid, /*max_t=*/-0.1);
  EXPECT_EQ(0, v1.size());
  EXPECT_EQ(0, v2.size());
}

Checklist:

• [X] I have commented my code, particularly in hard-to-understand areas
• [X] I have made corresponding changes to the documentation
• [X] New and existing unit tests pass locally with my changes

[cgyurgyik]

@matthewturk As discussed on Slack, here's the meat of changing to a unitized time where max_t can be (0.0, 1.0]. This means no matter where the ray is, inside or outside the sphere, if max_t = 1.0, it will traverse the entire sphere. If you want to traverse half the sphere, set max_t = 0.5.

// Quick notes: ray_origin_is_outside_grid is a boolean so its value is 0 or 1.           

// Set current t.
  double t = t_ray_entrance * ray_origin_is_outside_grid;

// Set max_t.
  const double unitized_ray_time = max_t * grid.sphereMaxRadius() * 2.0 +
                                   t_ray_entrance * ray_origin_is_outside_grid;
  max_t = ray_origin_is_outside_grid ? std::min(t_ray_exit, unitized_ray_time)
                                     : unitized_ray_time;