From bda1f448c37285d3d371fcc7f440890a0d3e2007 Mon Sep 17 00:00:00 2001 From: Chris Gyurgyik Date: Sat, 23 May 2020 12:48:13 -0400 Subject: [PATCH 1/5] Add sampler function to cpp version. --- cpp/spherical_volume_rendering_util.cpp | 25 ++++++++++++++++--------- cpp/spherical_volume_rendering_util.h | 12 ++++++++++-- 2 files changed, 26 insertions(+), 11 deletions(-) diff --git a/cpp/spherical_volume_rendering_util.cpp b/cpp/spherical_volume_rendering_util.cpp index ddbb37d..0fc457f 100644 --- a/cpp/spherical_volume_rendering_util.cpp +++ b/cpp/spherical_volume_rendering_util.cpp @@ -424,8 +424,11 @@ inline void initializeVoxelBoundarySegments( } std::vector walkSphericalVolume( - const Ray &ray, const svr::SphericalVoxelGrid &grid, - double max_t) noexcept { + const Ray &ray, const svr::SphericalVoxelGrid &grid, double max_t, + const std::function &sampler, + void *data) noexcept { if (max_t <= 0.0) return {}; const FreeVec3 rsv = grid.sphereCenter() - ray.pointAtParameter(0.0); // Ray Sphere Vector. @@ -499,6 +502,7 @@ std::vector walkSphericalVolume( 0.0, ray.timeOfIntersectionAt(grid.sphereCenter())}; RaySegment ray_segment(max_t, ray); + double previous_intersection_t; bool radial_step_has_transitioned = false; while (true) { const auto radial = @@ -516,6 +520,7 @@ std::vector walkSphericalVolume( } const auto voxel_intersection = minimumIntersection(radial, polar, azimuthal); + previous_intersection_t = t; switch (voxel_intersection) { case Radial: { t = radial.tMax; @@ -566,14 +571,16 @@ std::vector walkSphericalVolume( break; } } - if (voxels.back().radial == current_radial_voxel && - voxels.back().polar == current_polar_voxel && - voxels.back().azimuthal == current_azimuthal_voxel) { - continue; + const svr::SphericalVoxel intersected_voxel = { + .radial = current_radial_voxel, + .polar = current_polar_voxel, + .azimuthal = current_azimuthal_voxel}; + if (voxels.back() == intersected_voxel) continue; + if (sampler != nullptr) { + sampler(grid, ray, /*enter_t=*/previous_intersection_t, /*exit_t=*/t, + intersected_voxel, data); } - voxels.push_back({.radial = current_radial_voxel, - .polar = current_polar_voxel, - .azimuthal = current_azimuthal_voxel}); + voxels.push_back(intersected_voxel); } } diff --git a/cpp/spherical_volume_rendering_util.h b/cpp/spherical_volume_rendering_util.h index 0e64711..cd94f3f 100644 --- a/cpp/spherical_volume_rendering_util.h +++ b/cpp/spherical_volume_rendering_util.h @@ -8,12 +8,16 @@ #include "vec3.h" namespace svr { - // Represents a spherical voxel coordinate. struct SphericalVoxel { int radial; int polar; int azimuthal; + + inline bool operator==(SphericalVoxel v) { + return this->radial == v.radial && this->polar == v.polar && + this->azimuthal == v.azimuthal; + } }; // A spherical coordinate voxel traversal algorithm. The algorithm traces the @@ -26,7 +30,11 @@ struct SphericalVoxel { // then no voxels will be traversed. If max_t >= 1.0, then the entire sphere // will be traversed. std::vector walkSphericalVolume( - const Ray &ray, const svr::SphericalVoxelGrid &grid, double max_t) noexcept; + const Ray &ray, const svr::SphericalVoxelGrid &grid, double max_t, + const std::function& + sampler = nullptr, void* data = nullptr) noexcept; // Simplified parameters to Cythonize the function; implementation remains the // same as above. From 4aaf54d5a83ee3bf7af46419de4a0e88cd9f3445 Mon Sep 17 00:00:00 2001 From: Chris Gyurgyik Date: Sat, 23 May 2020 12:48:48 -0400 Subject: [PATCH 2/5] Add sampler function to cpp version. --- cpp/spherical_volume_rendering_util.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/spherical_volume_rendering_util.cpp b/cpp/spherical_volume_rendering_util.cpp index 0fc457f..c122707 100644 --- a/cpp/spherical_volume_rendering_util.cpp +++ b/cpp/spherical_volume_rendering_util.cpp @@ -576,7 +576,7 @@ std::vector walkSphericalVolume( .polar = current_polar_voxel, .azimuthal = current_azimuthal_voxel}; if (voxels.back() == intersected_voxel) continue; - if (sampler != nullptr) { + if (sampler != nullptr && data != nullptr) { sampler(grid, ray, /*enter_t=*/previous_intersection_t, /*exit_t=*/t, intersected_voxel, data); } From 52f4b4805db78e93079a0bece8230ad5830a2f7d Mon Sep 17 00:00:00 2001 From: Chris Gyurgyik Date: Sat, 23 May 2020 12:52:01 -0400 Subject: [PATCH 3/5] const noexcept function. --- cpp/spherical_volume_rendering_util.h | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/cpp/spherical_volume_rendering_util.h b/cpp/spherical_volume_rendering_util.h index cd94f3f..cd2202c 100644 --- a/cpp/spherical_volume_rendering_util.h +++ b/cpp/spherical_volume_rendering_util.h @@ -14,7 +14,7 @@ struct SphericalVoxel { int polar; int azimuthal; - inline bool operator==(SphericalVoxel v) { + inline bool operator==(SphericalVoxel v) const noexcept { return this->radial == v.radial && this->polar == v.polar && this->azimuthal == v.azimuthal; } @@ -31,10 +31,11 @@ struct SphericalVoxel { // will be traversed. std::vector walkSphericalVolume( const Ray &ray, const svr::SphericalVoxelGrid &grid, double max_t, - const std::function& - sampler = nullptr, void* data = nullptr) noexcept; + const std::function + &sampler = nullptr, + void *data = nullptr) noexcept; // Simplified parameters to Cythonize the function; implementation remains the // same as above. From c48bac2568ca50913f99ea20791dd968386f7684 Mon Sep 17 00:00:00 2001 From: Chris Gyurgyik Date: Sat, 23 May 2020 12:52:53 -0400 Subject: [PATCH 4/5] Add todo. --- cpp/spherical_volume_rendering_util.h | 1 + 1 file changed, 1 insertion(+) diff --git a/cpp/spherical_volume_rendering_util.h b/cpp/spherical_volume_rendering_util.h index cd2202c..944628e 100644 --- a/cpp/spherical_volume_rendering_util.h +++ b/cpp/spherical_volume_rendering_util.h @@ -29,6 +29,7 @@ struct SphericalVoxel { // 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. +// todo(cgyurgyik): Add description of sampler, data. std::vector walkSphericalVolume( const Ray &ray, const svr::SphericalVoxelGrid &grid, double max_t, const std::function Date: Fri, 7 Aug 2020 14:25:10 -0400 Subject: [PATCH 5/5] Const correctness for spherical_voxel_grid. Makes some slight changes to the algorithm. --- cpp/spherical_volume_rendering_util.cpp | 175 +++++++++-------- cpp/spherical_volume_rendering_util.h | 2 + cpp/spherical_voxel_grid.h | 244 ++++++++++++++---------- 3 files changed, 233 insertions(+), 188 deletions(-) diff --git a/cpp/spherical_volume_rendering_util.cpp b/cpp/spherical_volume_rendering_util.cpp index 01bc6f4..86c481f 100644 --- a/cpp/spherical_volume_rendering_util.cpp +++ b/cpp/spherical_volume_rendering_util.cpp @@ -9,6 +9,7 @@ namespace svr { +namespace { constexpr double DOUBLE_MAX = std::numeric_limits::max(); // The type corresponding to the voxel(s) with the minimum tMax value for a @@ -56,8 +57,8 @@ struct RaySegment { inline double intersectionTimeAt(double intersect_parameter, const Ray &ray) const noexcept { return (P1_[NZDI_] + ray_segment_[NZDI_] * intersect_parameter - - ray.origin()[NZDI_]) * - ray.invDirection()[NZDI_]; + ray.origin()[NZDI_]) * + ray.invDirection()[NZDI_]; } inline const BoundVec3 &P1() const noexcept { return P1_; } @@ -98,7 +99,7 @@ inline int calculateAngularVoxelIDFromPoints( const double Y_p1_diff = angular_max[j].P1 - p1; const double Y_p2_diff = angular_max[j].P2 - p2; const double d1d2 = (X_p1_diff * X_p1_diff) + (X_p2_diff * X_p2_diff) + - (Y_p1_diff * Y_p1_diff) + (Y_p2_diff * Y_p2_diff); + (Y_p1_diff * Y_p1_diff) + (Y_p2_diff * Y_p2_diff); const double d3 = (X_diff * X_diff) + (Y_diff * Y_diff); if (d1d2 < d3 || svr::isEqual(d1d2, d3)) return i; } @@ -107,22 +108,22 @@ inline int calculateAngularVoxelIDFromPoints( // Returns true if the "step" taken from the current voxel ID remains in // the grid bounds. -inline bool inBoundsAzimuthal(const SphericalVoxelGrid &grid, - const int step, +inline bool inBoundsAzimuthal(const SphericalVoxelGrid &grid, const int step, const int azi_voxel) noexcept { - const double radian = (azi_voxel + 1)*grid.deltaPhi(); + const double radian = (azi_voxel + 1) * grid.deltaPhi(); const double angval = radian - std::abs(step * grid.deltaPhi()); - return angval <= grid.sphereMaxBoundAzi() && angval >= grid.sphereMinBoundAzi() ; + return angval <= grid.sphereMaxBoundAzi() && + angval >= grid.sphereMinBoundAzi(); } // Returns true if the "step" taken from the current voxel ID remains in // the grid bounds. -inline bool inBoundsPolar(const SphericalVoxelGrid &grid, - const int step, +inline bool inBoundsPolar(const SphericalVoxelGrid &grid, const int step, const int pol_voxel) noexcept { - const double radian = (pol_voxel + 1)*grid.deltaTheta(); + const double radian = (pol_voxel + 1) * grid.deltaTheta(); const double angval = radian - std::abs(step * grid.deltaTheta()); - return angval <= grid.sphereMaxBoundPolar() && angval >= grid.sphereMinBoundPolar() ; + return angval <= grid.sphereMaxBoundPolar() && + angval >= grid.sphereMinBoundPolar(); } // Initializes an angular voxel ID. For polar initialization, *_2 represents @@ -163,7 +164,7 @@ inline HitParameters radialHit(const Ray &ray, if (radial_step_has_transitioned) { const double d_b = std::sqrt(grid.deltaRadiiSquared(current_radial_voxel - 1) - - rsvd_minus_v_squared); + rsvd_minus_v_squared); const double intersection_t = ray.timeOfIntersectionAt(v + d_b); if (intersection_t < max_t) return {.tMax = intersection_t, .tStep = -1}; } else { @@ -172,7 +173,7 @@ inline HitParameters radialHit(const Ray &ray, grid.numRadialSections() - 1); const double r_a = grid.deltaRadiiSquared( previous_idx - - (grid.deltaRadiiSquared(previous_idx) < rsvd_minus_v_squared)); + (grid.deltaRadiiSquared(previous_idx) < rsvd_minus_v_squared)); const double d_a = std::sqrt(r_a - rsvd_minus_v_squared); const double t_entrance = ray.timeOfIntersectionAt(v - d_a); const double t_exit = ray.timeOfIntersectionAt(v + d_a); @@ -213,12 +214,12 @@ HitParameters angularHit( const std::vector &P_max, int current_voxel) noexcept { const bool is_parallel_min = svr::isEqual(perp_uv_min, 0.0); const bool is_collinear_min = is_parallel_min && - svr::isEqual(perp_uw_min, 0.0) && - svr::isEqual(perp_vw_min, 0.0); + svr::isEqual(perp_uw_min, 0.0) && + svr::isEqual(perp_vw_min, 0.0); const bool is_parallel_max = svr::isEqual(perp_uv_max, 0.0); const bool is_collinear_max = is_parallel_max && - svr::isEqual(perp_uw_max, 0.0) && - svr::isEqual(perp_vw_max, 0.0); + svr::isEqual(perp_uw_max, 0.0) && + svr::isEqual(perp_vw_max, 0.0); double a, b; double t_min = collinear_times[is_collinear_min]; bool is_intersect_min = false; @@ -227,7 +228,7 @@ HitParameters angularHit( a = perp_vw_min * inv_perp_uv_min; b = perp_uw_min * inv_perp_uv_min; if (!((svr::lessThan(a, 0.0) || svr::lessThan(1.0, a)) || - svr::lessThan(b, 0.0) || svr::lessThan(1.0, b))) { + svr::lessThan(b, 0.0) || svr::lessThan(1.0, b))) { is_intersect_min = true; t_min = ray_segment.intersectionTimeAt(b, ray); } @@ -239,12 +240,11 @@ HitParameters angularHit( a = perp_vw_max * inv_perp_uv_max; b = perp_uw_max * inv_perp_uv_max; if (!((svr::lessThan(a, 0.0) || svr::lessThan(1.0, a)) || - svr::lessThan(b, 0.0) || svr::lessThan(1.0, b))) { + svr::lessThan(b, 0.0) || svr::lessThan(1.0, b))) { is_intersect_max = true; t_max = ray_segment.intersectionTimeAt(b, ray); } } - const bool t_t_max_eq = svr::isEqual(t, t_max); const bool t_max_within_bounds = t < t_max && !t_t_max_eq && t_max < max_t; const bool t_t_min_eq = svr::isEqual(t, t_min); @@ -276,9 +276,9 @@ HitParameters angularHit( const int next_step = std::abs( current_voxel - calculateAngularVoxelIDFromPoints(P_max, p1, p2)); return {.tMax = t_max, - .tStep = ray.direction().x() < 0.0 || ray_direction_2 < 0.0 - ? next_step - : -next_step}; + .tStep = ray.direction().x() < 0.0 || ray_direction_2 < 0.0 + ? next_step + : -next_step}; } if (t_min_within_bounds && ((t_min < t_max && !min_max_eq) || t_t_max_eq)) { return {.tMax = t_min, .tStep = -1}; @@ -309,15 +309,15 @@ inline HitParameters polarHit(const Ray &ray, const FreeVec3 w_min = p_one - ray_segment.P1(); const FreeVec3 w_max = p_two - ray_segment.P1(); const double perp_uv_min = u_min->x() * ray_segment.vector().y() - - u_min->y() * ray_segment.vector().x(); + u_min->y() * ray_segment.vector().x(); const double perp_uv_max = u_max->x() * ray_segment.vector().y() - - u_max->y() * ray_segment.vector().x(); + u_max->y() * ray_segment.vector().x(); const double perp_uw_min = u_min->x() * w_min.y() - u_min->y() * w_min.x(); const double perp_uw_max = u_max->x() * w_max.y() - u_max->y() * w_max.x(); const double perp_vw_min = ray_segment.vector().x() * w_min.y() - - ray_segment.vector().y() * w_min.x(); + ray_segment.vector().y() * w_min.x(); const double perp_vw_max = ray_segment.vector().x() * w_max.y() - - ray_segment.vector().y() * w_max.x(); + ray_segment.vector().y() * w_max.x(); return angularHit(grid, ray, perp_uv_min, perp_uv_max, perp_uw_min, perp_uw_max, perp_vw_min, perp_vw_max, ray_segment, collinear_times, t, max_t, ray.direction().y(), @@ -346,15 +346,15 @@ inline HitParameters azimuthalHit(const Ray &ray, const FreeVec3 w_min = p_one - ray_segment.P1(); const FreeVec3 w_max = p_two - ray_segment.P1(); const double perp_uv_min = u_min->x() * ray_segment.vector().z() - - u_min->z() * ray_segment.vector().x(); + u_min->z() * ray_segment.vector().x(); const double perp_uv_max = u_max->x() * ray_segment.vector().z() - - u_max->z() * ray_segment.vector().x(); + u_max->z() * ray_segment.vector().x(); const double perp_uw_min = u_min->x() * w_min.z() - u_min->z() * w_min.x(); const double perp_uw_max = u_max->x() * w_max.z() - u_max->z() * w_max.x(); const double perp_vw_min = ray_segment.vector().x() * w_min.z() - - ray_segment.vector().z() * w_min.x(); + ray_segment.vector().z() * w_min.x(); const double perp_vw_max = ray_segment.vector().x() * w_max.z() - - ray_segment.vector().z() * w_max.x(); + ray_segment.vector().z() * w_max.x(); return angularHit(grid, ray, perp_uv_min, perp_uv_max, perp_uw_min, perp_uw_max, perp_vw_min, perp_vw_max, ray_segment, collinear_times, t, max_t, ray.direction().z(), @@ -373,7 +373,7 @@ inline HitParameters azimuthalHit(const Ray &ray, // 6. tMaxR, tMaxPhi equal intersection. // 7. tMaxTheta, tMaxPhi equal intersection. // For each case, the following must hold: t < tMax < max_t -// For reference in shorthand naming: +// For reference, uses the following shortform naming: // RP = Radial - Polar // RA = Radial - Azimuthal // PA = Polar - Azimuthal @@ -398,10 +398,12 @@ inline VoxelIntersectionType minimumIntersection( // Initialize an array of values representing the points of intersection between // the lines corresponding to voxel boundaries and a given radial voxel in the -// XY plane and XZ plane. Here, P_* represents these points with a given radius -// 'current_radius'. The case where the number of polar voxels is equal to the -// number of azimuthal voxels is also checked to reduce the number of -// trigonometric and floating point calculations. +// XY plane and XZ plane. Here, P_* represents these points with a given radius. +// +// The calculations used for P_polar are: +// P1 = current_radius * trig_value.cosine + sphere_center.x() +// P2 = current_radius * trig_value.sine + sphere_center.y() +// Similar for P_azimuthal, but uses Z-axis instead of Y-axis. inline void initializeVoxelBoundarySegments( std::vector &P_polar, std::vector &P_azimuthal, bool ray_origin_is_outside_grid, @@ -411,38 +413,24 @@ inline void initializeVoxelBoundarySegments( P_azimuthal = grid.pMaxAzimuthal(); return; } - if (grid.numPolarSections() == grid.numAzimuthalSections()) { - std::transform( - grid.polarTrigValues().cbegin(), grid.polarTrigValues().cend(), - P_polar.begin(), P_azimuthal.begin(), - [current_radius, &grid](const TrigonometricValues &tv, - LineSegment &polar_LS) -> LineSegment { - const double px_value = - current_radius * tv.cosine + grid.sphereCenter().x(); - const double current_radius_times_sin = current_radius * tv.sine; - polar_LS = {.P1 = px_value, - .P2 = current_radius_times_sin + grid.sphereCenter().y()}; - return {.P1 = px_value, - .P2 = current_radius_times_sin + grid.sphereCenter().z()}; - }); - return; - } std::transform( grid.polarTrigValues().cbegin(), grid.polarTrigValues().cend(), P_polar.begin(), [current_radius, &grid](const TrigonometricValues &tv) -> LineSegment { return {.P1 = current_radius * tv.cosine + grid.sphereCenter().x(), - .P2 = current_radius * tv.sine + grid.sphereCenter().y()}; + .P2 = current_radius * tv.sine + grid.sphereCenter().y()}; }); std::transform( grid.azimuthalTrigValues().cbegin(), grid.azimuthalTrigValues().cend(), P_azimuthal.begin(), [current_radius, &grid](const TrigonometricValues &tv) -> LineSegment { return {.P1 = current_radius * tv.cosine + grid.sphereCenter().x(), - .P2 = current_radius * tv.sine + grid.sphereCenter().z()}; + .P2 = current_radius * tv.sine + grid.sphereCenter().z()}; }); } +} // namespace + std::vector walkSphericalVolume( const Ray &ray, const svr::SphericalVoxelGrid &grid, double max_t) noexcept { @@ -461,7 +449,7 @@ std::vector walkSphericalVolume( const double entry_radius_squared = grid.deltaRadiiSquared(vector_index); const double entry_radius = grid.deltaRadius() * - static_cast(grid.numRadialSections() - vector_index); + static_cast(grid.numRadialSections() - vector_index); const double rsvd = rsv.dot(rsv); const double v = rsv.dot(ray.direction().to_free()); const double rsvd_minus_v_squared = rsvd - v * v; @@ -480,8 +468,8 @@ std::vector walkSphericalVolume( const FreeVec3 ray_sphere = ray_origin_is_outside_grid - ? grid.sphereCenter() - ray.pointAtParameter(t_ray_entrance) - : SED_from_center == 0.0 ? rsv - ray.direction().to_free() : rsv; + ? grid.sphereCenter() - ray.pointAtParameter(t_ray_entrance) + : SED_from_center == 0.0 ? rsv - ray.direction().to_free() : rsv; int current_polar_voxel = initializeAngularVoxelID( grid, grid.numPolarSections(), ray_sphere, P_polar, ray_sphere.y(), @@ -501,14 +489,14 @@ std::vector walkSphericalVolume( std::vector voxels; voxels.reserve(grid.numRadialSections() + grid.numPolarSections() + - grid.numAzimuthalSections()); + grid.numAzimuthalSections()); voxels.push_back({.radial = current_radial_voxel, - .polar = current_polar_voxel, - .azimuthal = current_azimuthal_voxel}); + .polar = current_polar_voxel, + .azimuthal = current_azimuthal_voxel}); double t = t_ray_entrance * ray_origin_is_outside_grid; const double unitized_ray_time = max_t * grid.sphereMaxDiameter() + - t_ray_entrance * ray_origin_is_outside_grid; + 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; @@ -529,16 +517,13 @@ std::vector walkSphericalVolume( current_polar_voxel, t, max_t); const auto azimuthal = azimuthalHit(ray, grid, ray_segment, collinear_times, current_azimuthal_voxel, t, max_t); + if (current_radial_voxel + radial.tStep == 0 || (radial.tMax == DOUBLE_MAX && polar.tMax == DOUBLE_MAX && - azimuthal.tMax == DOUBLE_MAX)) { + azimuthal.tMax == DOUBLE_MAX)) { voxels.back().exit_t = t_ray_exit; return voxels; } - const bool in_azi_bounds = - inBoundsAzimuthal(grid,azimuthal.tStep,current_azimuthal_voxel); - const bool in_polar_bounds = - inBoundsPolar(grid,polar.tStep,current_polar_voxel); const auto voxel_intersection = minimumIntersection(radial, polar, azimuthal); switch (voxel_intersection) { @@ -549,21 +534,31 @@ std::vector walkSphericalVolume( } case Polar: { t = polar.tMax; - if (!in_polar_bounds) return voxels; + if (!inBoundsPolar(grid, polar.tStep, current_polar_voxel)) { + voxels.back().exit_t = t_ray_exit; + return voxels; + } current_polar_voxel = (current_polar_voxel + polar.tStep) % grid.numPolarSections(); break; } case Azimuthal: { - if (!in_azi_bounds) return voxels; + if (!inBoundsAzimuthal(grid, azimuthal.tStep, + current_azimuthal_voxel)) { + voxels.back().exit_t = t_ray_exit; + return voxels; + } t = azimuthal.tMax; current_azimuthal_voxel = (current_azimuthal_voxel + azimuthal.tStep) % - grid.numAzimuthalSections(); + grid.numAzimuthalSections(); break; } case RadialPolar: { t = radial.tMax; - if (!in_polar_bounds) return voxels; + if (!inBoundsPolar(grid, polar.tStep, current_polar_voxel)) { + voxels.back().exit_t = t_ray_exit; + return voxels; + } current_radial_voxel += radial.tStep; current_polar_voxel = (current_polar_voxel + polar.tStep) % grid.numPolarSections(); @@ -571,29 +566,43 @@ std::vector walkSphericalVolume( } case RadialAzimuthal: { t = radial.tMax; - if (!in_azi_bounds) return voxels; + if (!inBoundsAzimuthal(grid, azimuthal.tStep, + current_azimuthal_voxel)) { + voxels.back().exit_t = t_ray_exit; + return voxels; + } current_radial_voxel += radial.tStep; current_azimuthal_voxel = (current_azimuthal_voxel + azimuthal.tStep) % - grid.numAzimuthalSections(); + grid.numAzimuthalSections(); break; } case PolarAzimuthal: { t = polar.tMax; - if (!in_azi_bounds) return voxels; + if (!inBoundsAzimuthal(grid, azimuthal.tStep, + current_azimuthal_voxel) || + !(inBoundsPolar(grid, polar.tStep, current_polar_voxel))) { + voxels.back().exit_t = t_ray_exit; + return voxels; + } current_polar_voxel = (current_polar_voxel + polar.tStep) % grid.numPolarSections(); current_azimuthal_voxel = (current_azimuthal_voxel + azimuthal.tStep) % - grid.numAzimuthalSections(); + grid.numAzimuthalSections(); break; } case RadialPolarAzimuthal: { t = radial.tMax; - if (!in_azi_bounds) return voxels; + if (!inBoundsAzimuthal(grid, azimuthal.tStep, + current_azimuthal_voxel) || + !(inBoundsPolar(grid, polar.tStep, current_polar_voxel))) { + voxels.back().exit_t = t_ray_exit; + return voxels; + } current_radial_voxel += radial.tStep; current_polar_voxel = (current_polar_voxel + polar.tStep) % grid.numPolarSections(); current_azimuthal_voxel = (current_azimuthal_voxel + azimuthal.tStep) % - grid.numAzimuthalSections(); + grid.numAzimuthalSections(); break; } } @@ -604,9 +613,9 @@ std::vector walkSphericalVolume( } voxels.back().exit_t = t; voxels.push_back({.radial = current_radial_voxel, - .polar = current_polar_voxel, - .azimuthal = current_azimuthal_voxel, - .enter_t = t}); + .polar = current_polar_voxel, + .azimuthal = current_azimuthal_voxel, + .enter_t = t}); } } @@ -621,11 +630,11 @@ std::vector walkSphericalVolume( UnitVec3(ray_direction[0], ray_direction[1], ray_direction[2])), svr::SphericalVoxelGrid( svr::SphereBound{.radial = min_bound[0], - .polar = min_bound[1], - .azimuthal = min_bound[2]}, + .polar = min_bound[1], + .azimuthal = min_bound[2]}, svr::SphereBound{.radial = max_bound[0], - .polar = max_bound[1], - .azimuthal = max_bound[2]}, + .polar = max_bound[1], + .azimuthal = max_bound[2]}, num_radial_voxels, num_polar_voxels, num_azimuthal_voxels, BoundVec3(sphere_center[0], sphere_center[1], sphere_center[2])), max_t); diff --git a/cpp/spherical_volume_rendering_util.h b/cpp/spherical_volume_rendering_util.h index d3984f9..b146808 100644 --- a/cpp/spherical_volume_rendering_util.h +++ b/cpp/spherical_volume_rendering_util.h @@ -14,6 +14,8 @@ struct SphericalVoxel { int radial; int polar; int azimuthal; + + // Entrance and exit time into the given voxel. double enter_t; double exit_t; }; diff --git a/cpp/spherical_voxel_grid.h b/cpp/spherical_voxel_grid.h index db573cf..453255f 100644 --- a/cpp/spherical_voxel_grid.h +++ b/cpp/spherical_voxel_grid.h @@ -6,7 +6,6 @@ #include "vec3.h" namespace svr { -constexpr double TAU = 2 * M_PI; // Represents the boundary for the sphere. This is used to determine the minimum // and maximum boundaries for a sectored traversal. @@ -29,11 +28,117 @@ struct TrigonometricValues { double sine; }; +namespace { + +constexpr double TAU = 2 * M_PI; + +// Initializes the delta radii squared. These values are used for radial hit +// calculations in the main spherical volume algorithm. This calculates +// delta_radius^2 for num_radial_sections + 1 iterations. The delta radius value +// begins at max_radius, and subtracts delta_radius with each index. +// For example, +// +// Given: num_radial_voxels = 3, max_radius = 6, delta_radius = 2 +// Returns: { 6*6, 4*4, 2*2, 0*0 } +std::vector initializeDeltaRadiiSquared( + const std::size_t num_radial_voxels, const double max_radius, + const double delta_radius) { + std::vector delta_radii_squared(num_radial_voxels + 1); + + double current_delta_radius = max_radius; + std::generate(delta_radii_squared.begin(), delta_radii_squared.end(), + [&]() -> double { + const double old_delta_radius = current_delta_radius; + current_delta_radius -= delta_radius; + return old_delta_radius * old_delta_radius; + }); + return delta_radii_squared; +} + +// Returns a vector of TrigonometricValues for the given number of voxels. +// This begins with min_bound, and increments by a value of delta +// for num_voxels + 1 iterations. For example, +// +// Given: num_voxels = 2, min_bound = 0.0, delta = pi/2 +// Returns: { {.cosine=1.0, .sine=0.0}, +// {.cosine=0.0, .sine=1.0}, +// {.cosine=1.0, .sine=0.0} } +std::vector initializeTrigonometricValues( + const std::size_t num_voxels, const double min_bound, const double delta) { + std::vector trig_values(num_voxels + 1); + + double radians = min_bound; + std::generate(trig_values.begin(), trig_values.end(), + [&]() -> TrigonometricValues { + const double cos = std::cos(radians); + const double sin = std::sin(radians); + radians += delta; + return {.cosine = cos, .sine = sin}; + }); + return trig_values; +} + +// Returns a vector of maximum radius line segments for the given trigonometric +// values. The following predicate should be true: +// num_voxels + 1 == trig_values.size() + 1. +// +// The LineSegment points P1 and P2 are calculated with the following equation: +// .P1 = max_radius * trig_value.cosine + center.x(). +// .P2 = max_radius * trig_value.sine + center.y(). +std::vector initializeMaxRadiusLineSegments( + const std::size_t num_voxels, const BoundVec3 ¢er, + const double max_radius, + const std::vector &trig_values) { + std::vector line_segments(num_voxels + 1); + std::transform(trig_values.cbegin(), trig_values.cend(), + line_segments.begin(), + [&](const TrigonometricValues &trig_value) -> LineSegment { + return {.P1 = max_radius * trig_value.cosine + center.x(), + .P2 = max_radius * trig_value.sine + center.y()}; + }); + return line_segments; +} + +// Initializes the vectors determined by the following calculation: +// sphere center - {X, Y, Z}, WHERE X, Y = P1, P2 for polar voxels. +std::vector initializeCenterToPolarPMaxVectors( + const std::vector &line_segments, const BoundVec3 ¢er) { + std::vector center_to_pmax_vectors; + center_to_pmax_vectors.reserve(line_segments.size()); + + for (const auto &points : line_segments) { + center_to_pmax_vectors.emplace_back(center - + FreeVec3(points.P1, points.P2, 0.0)); + } + return center_to_pmax_vectors; +} + +// Similar to above, but uses: +// sphere center - {X, Y, Z}, WHERE X, Z = P1, P2 for azimuthal voxels. +std::vector initializeCenterToAzimuthalPMaxVectors( + const std::vector &line_segments, const BoundVec3 ¢er) { + std::vector center_to_pmax_vectors; + center_to_pmax_vectors.reserve(line_segments.size()); + + for (const auto &points : line_segments) { + center_to_pmax_vectors.emplace_back(center - + FreeVec3(points.P1, 0.0, points.P2)); + } + return center_to_pmax_vectors; +} + +} // namespace + // Represents a spherical voxel grid used for ray casting. The bounds of the // grid are determined by min_bound and max_bound. The deltas are then // determined by (max_bound.X - min_bound.X) / num_X_sections. To minimize // calculation duplication, many calculations are completed once here and used // each time a ray traverses the spherical voxel grid. +// +// Note that the grid system currently does not align with one would expect +// from spherical coordinates. We represent both polar and azimuthal within +// bounds [0, 2pi]. +// TODO(cgyurgyik): Look into updating polar grid from [0, 2pi] -> [0, pi]. struct SphericalVoxelGrid { public: SphericalVoxelGrid(const SphereBound &min_bound, const SphereBound &max_bound, @@ -49,97 +154,35 @@ struct SphericalVoxelGrid { sphere_min_bound_polar_(min_bound.polar), sphere_max_bound_azimuthal_(max_bound.azimuthal), sphere_min_bound_azimuthal_(min_bound.azimuthal), + // TODO(cgyurgyik): Verify we want the sphere_max_radius to simply be + // max_bound.radial. sphere_max_radius_(max_bound.radial), sphere_max_diameter_(sphere_max_radius_ * 2.0), delta_radius_((max_bound.radial - min_bound.radial) / num_radial_sections), delta_theta_((max_bound.polar - min_bound.polar) / num_polar_sections), delta_phi_((max_bound.azimuthal - min_bound.azimuthal) / - num_azimuthal_sections) { - double current_delta_radius = max_bound.radial - min_bound.radial; - delta_radii_sq_.resize(num_radial_sections + 1); - std::generate(delta_radii_sq_.begin(), delta_radii_sq_.end(), - [&]() -> double { - const double old_delta_radius = current_delta_radius; - current_delta_radius -= delta_radius_; - return old_delta_radius * old_delta_radius; - }); - P_max_polar_.resize(num_polar_sections + 1); - P_max_azimuthal_.resize(num_azimuthal_sections + 1); - center_to_polar_bound_vectors_.reserve(num_polar_sections + 1); - center_to_azimuthal_bound_vectors_.reserve(num_azimuthal_sections + 1); - if (num_polar_sections == num_azimuthal_sections) { - double radians = 0.0; - polar_trig_values_.resize(num_polar_sections + 1); - std::generate(polar_trig_values_.begin(), polar_trig_values_.end(), - [&]() -> TrigonometricValues { - const double cos = std::cos(radians); - const double sin = std::sin(radians); - radians += delta_theta_; - return {.cosine = cos, .sine = sin}; - }); - std::transform(polar_trig_values_.cbegin(), polar_trig_values_.cend(), - P_max_polar_.begin(), P_max_azimuthal_.begin(), - [&](const TrigonometricValues &tv, - LineSegment &ang_LS) -> LineSegment { - const double px_max_value = - sphere_max_radius_ * tv.cosine + sphere_center.x(); - const double max_radius_times_s = - sphere_max_radius_ * tv.sine; - ang_LS = {.P1 = px_max_value, - .P2 = max_radius_times_s + sphere_center.y()}; - return {.P1 = px_max_value, - .P2 = max_radius_times_s + sphere_center.z()}; - }); - for (const auto &points : P_max_polar_) { - const BoundVec3 v = - sphere_center - FreeVec3(points.P1, points.P2, points.P2); - center_to_polar_bound_vectors_.push_back(v); - center_to_azimuthal_bound_vectors_.push_back(v); - } - return; - } - double radians = min_bound.polar; - polar_trig_values_.resize(num_polar_sections + 1); - std::generate(polar_trig_values_.begin(), polar_trig_values_.end(), - [&]() -> TrigonometricValues { - const double cos = std::cos(radians); - const double sin = std::sin(radians); - radians += delta_theta_; - return {.cosine = cos, .sine = sin}; - }); - radians = min_bound.azimuthal; - azimuthal_trig_values_.resize(num_azimuthal_sections + 1); - std::generate(azimuthal_trig_values_.begin(), azimuthal_trig_values_.end(), - [&]() -> TrigonometricValues { - const double cos = std::cos(radians); - const double sin = std::sin(radians); - radians += delta_phi_; - return {.cosine = cos, .sine = sin}; - }); - std::transform( - polar_trig_values_.cbegin(), polar_trig_values_.cend(), - P_max_polar_.begin(), - [&](const TrigonometricValues &ang_tv) -> LineSegment { - return {.P1 = sphere_max_radius_ * ang_tv.cosine + sphere_center.x(), - .P2 = sphere_max_radius_ * ang_tv.sine + sphere_center.y()}; - }); - std::transform( - azimuthal_trig_values_.cbegin(), azimuthal_trig_values_.cend(), - P_max_azimuthal_.begin(), - [&](const TrigonometricValues &azi_tv) -> LineSegment { - return {.P1 = sphere_max_radius_ * azi_tv.cosine + sphere_center.x(), - .P2 = sphere_max_radius_ * azi_tv.sine + sphere_center.z()}; - }); - for (const auto &points : P_max_polar_) { - center_to_polar_bound_vectors_.emplace_back( - sphere_center - FreeVec3(points.P1, points.P2, 0.0)); - } - for (const auto &points : P_max_azimuthal_) { - center_to_azimuthal_bound_vectors_.emplace_back( - sphere_center - FreeVec3(points.P1, 0.0, points.P2)); - } - } + num_azimuthal_sections), + // TODO(cgyurgyik): Verify this is actually what we want for + // 'max_radius'. The other option is simply using max_bound.radial + delta_radii_sq_(initializeDeltaRadiiSquared( + num_radial_sections, + /*max_radius=*/max_bound.radial - min_bound.radial, delta_radius_)), + polar_trig_values_(initializeTrigonometricValues( + num_polar_sections, min_bound.polar, delta_theta_)), + azimuthal_trig_values_(initializeTrigonometricValues( + num_azimuthal_sections, min_bound.azimuthal, delta_phi_)), + P_max_polar_(initializeMaxRadiusLineSegments( + num_polar_sections, sphere_center, sphere_max_radius_, + polar_trig_values_)), + P_max_azimuthal_(initializeMaxRadiusLineSegments( + num_azimuthal_sections, sphere_center, sphere_max_radius_, + azimuthal_trig_values_)), + center_to_polar_bound_vectors_( + initializeCenterToPolarPMaxVectors(P_max_polar_, sphere_center)), + center_to_azimuthal_bound_vectors_( + initializeCenterToAzimuthalPMaxVectors(P_max_azimuthal_, + sphere_center)) {} inline std::size_t numRadialSections() const noexcept { return this->num_radial_sections_; @@ -258,28 +301,19 @@ struct SphericalVoxelGrid { // sections respectively. const double delta_theta_, delta_phi_; - // The delta radii squared ranging from 0...num_radial_voxels. - std::vector delta_radii_sq_; - - // The maximum radius line segments for polar voxels. - std::vector P_max_polar_; - - // The maximum radius line segments for azimuthal voxels. - std::vector P_max_azimuthal_; - - // The trigonometric values for each delta theta. - std::vector polar_trig_values_; + // The delta radii squared calculated for use in radial hit calculations. + const std::vector delta_radii_sq_; - // The trigonometric values for each delta phi. In the case where delta theta - // is equal to delta phi, this is left uninitialized and polar_trig_values_ is - // used. - std::vector azimuthal_trig_values_; + // The trigonometric values calculated for the azimuthal and polar voxels. + const std::vector azimuthal_trig_values_, + polar_trig_values_; - // The vectors represented by sphere_center_ - P_max_polar_[i]. - std::vector center_to_polar_bound_vectors_; + // The maximum radius line segments for polar and azimuthal voxels. + const std::vector P_max_polar_, P_max_azimuthal_; - // The vectors represented by sphere_center_ - P_max_azimuthal_[i]. - std::vector center_to_azimuthal_bound_vectors_; + // The vectors represented by the vector sphere center - P_max[i]. + const std::vector center_to_polar_bound_vectors_, + center_to_azimuthal_bound_vectors_; }; } // namespace svr