From dcc478740ba965c41d518c23dfa52fc465c9a753 Mon Sep 17 00:00:00 2001 From: ak2485 Date: Fri, 17 Apr 2020 13:27:27 -0400 Subject: [PATCH 1/3] Add collinear cases to angular/azimuthal hit; update tests; remove pointArray update from main traversal. --- .../SphericalCoordinateTraversalExample.m | 10 +- matlab/spherical-traversal/angular_hit.m | 121 ++++-- matlab/spherical-traversal/angular_hit.m~ | 212 ++++++++++ matlab/spherical-traversal/azimuthal_hit.m | 133 ++++--- .../sphericalCoordinateTraversal.m | 66 +--- .../sphericalCoordinateTraversal.m~ | 371 ++++++++++++++++++ .../sphericalCoordinateTraversalTest.m | 61 +-- 7 files changed, 815 insertions(+), 159 deletions(-) create mode 100644 matlab/spherical-traversal/angular_hit.m~ create mode 100644 matlab/spherical-traversal/sphericalCoordinateTraversal.m~ diff --git a/matlab/spherical-traversal/SphericalCoordinateTraversalExample.m b/matlab/spherical-traversal/SphericalCoordinateTraversalExample.m index 70ef7ab..38a4b85 100644 --- a/matlab/spherical-traversal/SphericalCoordinateTraversalExample.m +++ b/matlab/spherical-traversal/SphericalCoordinateTraversalExample.m @@ -1,7 +1,7 @@ - min_bound = [-20, -20.0, -20.0]; +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.0, 1.0]; + ray_origin = [13.0, -15.0, -15.0]; + ray_direction = [-1.0, 1.0, 1.0]; sphere_center = [0.0, 0.0, 0.0]; sphere_max_radius = 10.0; @@ -10,7 +10,7 @@ num_azimuthal_sections = 4; t_begin = 0.0; t_end = 30.0; - verbose = true; + verbose = false; [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) + sphere_center, sphere_max_radius, num_radial_sections, num_angular_sections, num_azimuthal_sections, t_begin, t_end, verbose) \ No newline at end of file diff --git a/matlab/spherical-traversal/angular_hit.m b/matlab/spherical-traversal/angular_hit.m index 3ecef71..c4cca8f 100644 --- a/matlab/spherical-traversal/angular_hit.m +++ b/matlab/spherical-traversal/angular_hit.m @@ -1,5 +1,5 @@ function [tMaxTheta, tStepTheta] = angular_hit(ray_origin, ray_direction, current_voxel_ID_theta,... - num_angular_sections, sphere_center, pointArray, t, t_end, verbose) + sphere_center, P_max, sphere_max_radius, t, t_end, verbose) % Determines whether an angular hit occurs for the given ray. % Input: % ray_origin: vector of the origin of the ray in cartesian coordinate. @@ -22,7 +22,7 @@ if verbose fprintf("\nCurrent Voxel ID Theta: %d\n", current_voxel_ID_theta) end - +pointArray = [P_max(current_voxel_ID_theta+1,:) ; P_max(current_voxel_ID_theta+2,:)]; % Get current ray location to use as start point p = ray_origin + ray_direction * t; % Get ray exit point to use as end point @@ -43,12 +43,34 @@ perpvw_max = v(1) * w_max(2) - v(2) * w_max(1); % check to see if ray is parallel to each angular voxel boundary -% if the ray is parallel to or collinear with a given voxel boundary, -% then intersection does not occur between the given voxel boundary and -% the ray. -if approximatelyEqual(perpuv_min,0.0,1e-12,1e-8) ; parallel_min = 1; else; parallel_min = 0; end -if approximatelyEqual(perpuv_max,0.0,1e-12,1e-8) ; parallel_max = 1; else; parallel_max = 0; end - +if approximatelyEqual(perpuv_min,0.0,1e-12,1e-8) + parallel_min = 1; + % check to see if the ray is colliner with the angular voxel boundary + if approximatelyEqual(perpuw_min,0.0,1e-12,1e-8) && ... + approximatelyEqual(perpvw_min,0.0,1e-12,1e-8) + collinear_min = 1; + p_min = sphere_center; + else + collinear_min = 0; + end +else + parallel_min = 0; + collinear_min = 0; +end +if approximatelyEqual(perpuv_max,0.0,1e-12,1e-8) + parallel_max = 1; + if approximatelyEqual(perpuw_max,0.0,1e-12,1e-8) && ... + approximatelyEqual(perpvw_max,0.0,1e-12,1e-8) + collinear_max = 1; + p_max = sphere_center; + else + collinear_max = 0; + end +else + parallel_max = 0; + collinear_max = 0; +end +% Determine intersection points, p_min and p_max if parallel_min == 0 a = perpvw_min / perpuv_min; b = perpuw_min / perpuv_min; @@ -62,7 +84,6 @@ else intersect_min = 0; end - if parallel_max == 0 a = perpvw_max / perpuv_max; b = perpuw_max / perpuv_max; @@ -76,71 +97,95 @@ else intersect_max = 0; end - +% Compute inverses of slopes if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) inv_dir = 1/ray_direction(1); else inv_dir = 1/ray_direction(2); end - -if intersect_max == 1 +% Compute intersection times +if intersect_max == 1 || collinear_max == 1 if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) - inv_dir = 1/ray_direction(1); t_max = (p_max(1) - ray_origin(1)) * inv_dir; else - inv_dir = 1/ray_direction(2); t_max = (p_max(2) - ray_origin(2)) * inv_dir; end end - -if intersect_min == 1 +if intersect_min == 1 || collinear_min == 1 if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) - inv_dir = 1/ray_direction(1); t_min = (p_min(1) - ray_origin(1)) * inv_dir; else - inv_dir = 1/ray_direction(2); t_min = (p_min(2) - ray_origin(2)) * inv_dir; end end - -if intersect_max == 1 && intersect_min == 0 && ... - t < t_max && t_max < t_end && ... - ~approximatelyEqual(t,t_max,1e-12,1e-8) +% Determine tStepTheta, tMaxTheta +if intersect_max == 1 && intersect_min == 0 && collinear_min == 0 &&... + strictlyLess(t,t_max,1e-12,1e-8) && strictlyLess(t_max,t_end,1e-12,1e-8) tStepTheta = 1; tMaxTheta = t_max; -elseif intersect_min == 1 && intersect_max == 0 && ... - t < t_min && t_min < t_end && ... - ~approximatelyEqual(t,t_min,1e-12,1e-8) +elseif intersect_min == 1 && intersect_max == 0 && collinear_max == 0 && ... + strictlyLess(t,t_min,1e-12,1e-8) && ... + strictlyLess(t_min,t_end,1e-12,1e-8) tStepTheta = -1; tMaxTheta = t_min; -elseif intersect_min == 1 && intersect_max == 1 +elseif (intersect_min == 1 && intersect_max == 1) || ... + (intersect_min == 1 && collinear_max == 1) || ... + (collinear_min == 1 && intersect_max == 1) % hit the min & max boundary simultaneously. if approximatelyEqual(t_min,t_max,1e-12,1e-8) && ... - t < t_min && t_min < t_end && ... - ~approximatelyEqual(t,t_min,1e-12,1e-8) + strictlyLess(t,t_min,1e-12,1e-8) && ... + strictlyLess(t_min,t_end,1e-12,1e-8) tMaxTheta = t_max; if verbose fprintf("hit origin t\n") + end + pert_t = 0.1; + pert_x = sphere_center(1) + ray_direction(1) * pert_t; + pert_y = sphere_center(2) + ray_direction(2) * pert_t; + a = sphere_center(1) - pert_x; + b = sphere_center(2) - pert_y; + l = sqrt(a^2 + b^2); + sphere_center_xy = [sphere_center(1), sphere_center(2)]; + p1 = sphere_center_xy - (sphere_max_radius/l) .* [a b]; + i = 1; + while i < length(P_max) + d1 = (P_max(i,1)-p1(1))^2 + (P_max(i,2)-p1(2))^2; + d2 = (P_max(i+1,1)-p1(1))^2 + (P_max(i+1,2)-p1(2))^2; + d3 = (P_max(i,1)-P_max(i+1,1))^2 + (P_max(i,2)-P_max(i+1,2))^2; + if strictlyLess(d1+d2,d3,1e-12,1e-8) || ... + approximatelyEqual(d1+d2,d3,1e-12,1e-8) + next_vox = i - 1 ; + i = length(P_max); + end + i = i + 1; end - if ~approximatelyEqual(ray_direction(2),0.0,1e-12,1e-8) - % change in theta is dependent on the slope of the line - tStepTheta = - num_angular_sections/2; + step = abs(current_voxel_ID_theta - next_vox); + if strictlyLess(ray_direction(2),0.0,1e-12,1e-8) + tStepTheta = step; + elseif strictlyLess(0.0,ray_direction(2),1e-12,1e-8) + tStepTheta = -step; else - tStepTheta = num_angular_sections/2; + if strictlyLess(ray_direction(1),0.0,1e-12,1e-8) + tStepTheta = step; + else + tStepTheta = -step; + end end % hit min bound - elseif t < t_min && t_min < t_end && ... - (t_min < t_max || approximatelyEqual(t,t_max,1e-12,1e-8)) && ... - ~approximatelyEqual(t,t_min,1e-12,1e-8) + elseif strictlyLess(t,t_min,1e-12,1e-8) && ... + strictlyLess(t_min,t_end,1e-12,1e-8) && ... + (strictlyLess(t_min,t_max,1e-12,1e-8) || ... + approximatelyEqual(t,t_max,1e-12,1e-8)) tStepTheta = -1; tMaxTheta = t_min; if verbose fprintf("hit min bound\n") end % hit max bound - elseif t < t_max && t_max < t_end && ... - (t_max < t_min || approximatelyEqual(t,t_min,1e-12,1e-8)) && ... - ~approximatelyEqual(t,t_max,1e-12,1e-8) + elseif strictlyLess(t,t_max,1e-12,1e-8) && ... + strictlyLess(t_max,t_end,1e-12,1e-8) && ... + (strictlyLess(t_max,t_min,1e-12,1e-8) || ... + approximatelyEqual(t,t_min,1e-12,1e-8)) tStepTheta = 1; tMaxTheta = t_max; if verbose diff --git a/matlab/spherical-traversal/angular_hit.m~ b/matlab/spherical-traversal/angular_hit.m~ new file mode 100644 index 0000000..af5777b --- /dev/null +++ b/matlab/spherical-traversal/angular_hit.m~ @@ -0,0 +1,212 @@ +function [tMaxTheta, tStepTheta] = angular_hit(ray_origin, ray_direction, current_voxel_ID_theta,... + sphere_center, P_max, sphere_max_radius, t, t_end, verbose) +% Determines whether an angular hit occurs for the given ray. +% Input: +% ray_origin: vector of the origin of the ray in cartesian coordinate. +% ray_direction: vector of the direction of the ray in cartesian coordinate. +% current_voxel_ID_theta: the (angular) ID of current voxel. +% num_angular_sections: number of total angular sections on the grid. +% sphere_center: The center of the circle. +% pointArray: +% t: The current time parameter of the ray. +% t_end: +% verbose: Determines whether debug print statements are enabled. +% +% Returns: +% tMaxTheta: is the time at which a hit occurs for the ray at the next point of intersection. +% tStepTheta: The direction the theta voxel steps. +1, -1, or 0. +if verbose + fprintf("\n-- angular_hit --") +end + +if verbose + fprintf("\nCurrent Voxel ID Theta: %d\n", current_voxel_ID_theta) +end +pointArray = [P_max(current_voxel_ID_theta+1,:) ; P_max(current_voxel_ID_theta+2,:)]; +% Get current ray location to use as start point +p = ray_origin + ray_direction * t; +% Get ray exit point to use as end point +p_f = ray_origin + ray_direction * t_end; +% Caculate the ray segment vector +v = p_f - p; +% Calculate the voxel boundary vectors +u_min = [sphere_center(1) sphere_center(2)] - pointArray(1,:); +u_max = [sphere_center(1) sphere_center(2)] - pointArray(2,:); +% http://geomalgorithms.com/a05-_intersect-1.html#intersect2D_2Segments() +w_min = pointArray(1,:) - [p(1) p(2)]; +w_max = pointArray(2,:) - [p(1) p(2)]; +perpuv_min = u_min(1) * v(2) - u_min(2) * v(1); +perpuv_max = u_max(1) * v(2) - u_max(2) * v(1); +perpuw_min = u_min(1) * w_min(2) - u_min(2) * w_min(1); +perpuw_max = u_max(1) * w_max(2) - u_max(2) * w_max(1); +perpvw_min = v(1) * w_min(2) - v(2) * w_min(1); +perpvw_max = v(1) * w_max(2) - v(2) * w_max(1); + +% check to see if ray is parallel to each angular voxel boundary +if approximatelyEqual(perpuv_min,0.0,1e-12,1e-8) + parallel_min = 1; + % check to see if the ray is colliner with the angular voxel boundary + if approximatelyEqual(perpuw_min,0.0,1e-12,1e-8) && ... + approximatelyEqual(perpvw_min,0.0,1e-12,1e-8) + collinear_min = 1; + p_min = sphere_center; + else + collinear_min = 0; + end +else + parallel_min = 0; + collinear_min = 0; +end +if approximatelyEqual(perpuv_max,0.0,1e-12,1e-8) + parallel_max = 1; + if approximatelyEqual(perpuw_max,0.0,1e-12,1e-8) && ... + approximatelyEqual(perpvw_max,0.0,1e-12,1e-8) + collinear_max = 1; + p_max = sphere_center; + else + collinear_max = 0; + end +else + parallel_max = 0; + collinear_max = 0; +end + +if parallel_min == 0 + a = perpvw_min / perpuv_min; + b = perpuw_min / perpuv_min; + if (strictlyLess(a,0.0,1e-12,1e-8) || strictlyLess(1.0,a,1e-12,1e-8)) || ... + (strictlyLess(b,0.0,1e-12,1e-8) || strictlyLess(1.0,b,1e-12,1e-8)) + intersect_min = 0; + else + intersect_min = 1; + p_min = p + b * v; + end +else + intersect_min = 0; +end + +if parallel_max == 0 + a = perpvw_max / perpuv_max; + b = perpuw_max / perpuv_max; + if (strictlyLess(a,0.0,1e-12,1e-8) || strictlyLess(1.0,a,1e-12,1e-8)) || ... + (strictlyLess(b,0.0,1e-12,1e-8) || strictlyLess(1.0,b,1e-12,1e-8)) + intersect_max = 0; + else + intersect_max = 1; + p_max = p + b * v; + end +else + intersect_max = 0; +end + +if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) + inv_dir = 1/ray_direction(1); +else + inv_dir = 1/ray_direction(2); +end + +if intersect_max == 1 || collinear_max == 1 + if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) + t_max = (p_max(1) - ray_origin(1)) * inv_dir; + else + t_max = (p_max(2) - ray_origin(2)) * inv_dir; + end +end + +if intersect_min == 1 || collinear_min == 1 + if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) + t_min = (p_min(1) - ray_origin(1)) * inv_dir; + else + t_min = (p_min(2) - ray_origin(2)) * inv_dir; + end +end + +if intersect_max == 1 && intersect_min == 0 && collinear_min == 0 &&... + strictlyLess(t,t_max,1e-12,1e-8) && strictlyLess(t_max,t_end,1e-12,1e-8) + tStepTheta = 1; + tMaxTheta = t_max; +elseif intersect_min == 1 && intersect_max == 0 && collinear_max == 0 && ... + strictlyLess(t,t_min,1e-12,1e-8) && ... + strictlyLess(t_min,t_end,1e-12,1e-8) + tStepTheta = -1; + tMaxTheta = t_min; +elseif (intersect_min == 1 && intersect_max == 1) || ... + (intersect_min == 1 && collinear_max == 1) || ... + (collinear_min == 1 && intersect_max == 1) + % hit the min & max boundary simultaneously. + if approximatelyEqual(t_min,t_max,1e-12,1e-8) && ... + strictlyLess(t,t_min,1e-12,1e-8) && ... + strictlyLess(t_min,t_end,1e-12,1e-8) + tMaxTheta = t_max; + if verbose + fprintf("hit origin t\n") + end + pert_t = 0.1; + pert_x = sphere_center(1) + ray_direction(1) * pert_t; + pert_y = sphere_center(2) + ray_direction(2) * pert_t; + a = sphere_center(1) - pert_x; + b = sphere_center(2) - pert_y; + l = sqrt(a^2 + b^2); + p1 = sphere_center - (sphere_max_radius/l) .* [a b]; + i = 1; + while i < length(P_max) + d1 = (P_max(i,1)-p1(1))^2 + (P_max(i,2)-p1(2))^2; + d2 = (P_max(i+1,1)-p1(1))^2 + (P_max(i+1,2)-p1(2))^2; + d3 = (P_max(i,1)-P_max(i+1,1))^2 + (P_max(i,2)-P_max(i+1,2))^2; + if strictlyLess(d1+d2,d3,1e-12,1e-8) || ... + approximatelyEqual(d1+d2,d3,1e-12,1e-8) + next_vox = i - 1 ; + i = length(P_max); + end + i = i + 1; + end + step = abs(current_voxel_ID_theta - next_vox); + if strictlyLess(ray_direction(2),0.0,1e-12,1e-8) + tStepTheta = step; + elseif strictlyLess(0.0,ray_direction(2),1e-12,1e-8) + tStepTheta = -step; + else + if strictlyLess(ray_direction(1),0.0,1e-12,1e-8) + tStepTheta = step; + else + tStepTheta = -step; + end + end + elseif t < t_min && t_min < t_end && ... + (t_min < t_max || approximatelyEqual(t,t_max,1e-12,1e-8)) && ... + ~approximatelyEqual(t,t_min,1e-12,1e-8) + tStepTheta = -1; + tMaxTheta = t_min; + if verbose + fprintf("hit min bound\n") + end + % hit max bound + elseif t < t_max && t_max < t_end && ... + (t_max < t_min || approximatelyEqual(t,t_min,1e-12,1e-8)) && ... + ~approximatelyEqual(t,t_max,1e-12,1e-8) + tStepTheta = 1; + tMaxTheta = t_max; + if verbose + fprintf("hit max bound\n") + end + else + tStepTheta = 0; + tMaxTheta = t_end; + if verbose + fprintf("no hit within time bound\n") + end + end + else + tStepTheta = 0; + tMaxTheta = t_end; + if verbose + fprintf("no hit within voxel bound\n") + end + end + + +if verbose + fprintf(['\ntMaxTheta: %d \n' ... + 'tStepTheta: %d \n'], tMaxTheta, tStepTheta); +end +end diff --git a/matlab/spherical-traversal/azimuthal_hit.m b/matlab/spherical-traversal/azimuthal_hit.m index f9cc89a..53735a4 100644 --- a/matlab/spherical-traversal/azimuthal_hit.m +++ b/matlab/spherical-traversal/azimuthal_hit.m @@ -1,5 +1,5 @@ function [tMaxPhi, tStepPhi] = azimuthal_hit(ray_origin, ray_direction, current_voxel_ID_phi,... - num_azimuthal_sections, sphere_center, pointArray, t, t_end, verbose) + sphere_center, P_max, sphere_max_radius, t, t_end, verbose) % Determines whether an azimuthal hit occurs for the given ray. % Input: % ray_origin: vector of the origin of the ray in cartesian coordinate. @@ -22,7 +22,7 @@ if verbose fprintf("\nCurrent Voxel ID Phi: %d\n", current_voxel_ID_phi) end -pointArray +pointArray = [P_max(current_voxel_ID_phi+1,:) ; P_max(current_voxel_ID_phi+2,:)]; % Get current ray location to use as start point p = ray_origin + ray_direction * t; % Get ray exit point to use as end point @@ -35,20 +35,42 @@ % http://geomalgorithms.com/a05-_intersect-1.html#intersect2D_2Segments() w_min = pointArray(1,:) - [p(1) p(3)]; w_max = pointArray(2,:) - [p(1) p(3)]; -perpuv_min = u_min(1) * v(2) - u_min(2) * v(1); -perpuv_max = u_max(1) * v(2) - u_max(2) * v(1); +perpuv_min = u_min(1) * v(3) - u_min(2) * v(1); +perpuv_max = u_max(1) * v(3) - u_max(2) * v(1); perpuw_min = u_min(1) * w_min(2) - u_min(2) * w_min(1); perpuw_max = u_max(1) * w_max(2) - u_max(2) * w_max(1); -perpvw_min = v(1) * w_min(2) - v(2) * w_min(1); -perpvw_max = v(1) * w_max(2) - v(2) * w_max(1); - -% check to see if ray is parallel to each azimuthal voxel boundary -% if the ray is parallel to or collinear with a given voxel boundary, -% then intersection does not occur between the given voxel boundary and -% the ray. -if approximatelyEqual(perpuv_min,0.0,1e-12,1e-8) ; parallel_min = 1; else; parallel_min = 0; end -if approximatelyEqual(perpuv_max,0.0,1e-12,1e-8) ; parallel_max = 1; else; parallel_max = 0; end +perpvw_min = v(1) * w_min(2) - v(3) * w_min(1); +perpvw_max = v(1) * w_max(2) - v(3) * w_max(1); +% check to see if ray is parallel to each angular voxel boundary +if approximatelyEqual(perpuv_min,0.0,1e-12,1e-8) + parallel_min = 1; + % check to see if the ray is colliner with the angular voxel boundary + if approximatelyEqual(perpuw_min,0.0,1e-12,1e-8) && ... + approximatelyEqual(perpvw_min,0.0,1e-12,1e-8) + collinear_min = 1; + p_min = sphere_center; + else + collinear_min = 0; + end +else + parallel_min = 0; + collinear_min = 0; +end +if approximatelyEqual(perpuv_max,0.0,1e-12,1e-8) + parallel_max = 1; + if approximatelyEqual(perpuw_max,0.0,1e-12,1e-8) && ... + approximatelyEqual(perpvw_max,0.0,1e-12,1e-8) + collinear_max = 1; + p_max = sphere_center; + else + collinear_max = 0; + end +else + parallel_max = 0; + collinear_max = 0; +end +% Determine intersection points, p_min and p_max if parallel_min == 0 a = perpvw_min / perpuv_min; b = perpuw_min / perpuv_min; @@ -62,7 +84,6 @@ else intersect_min = 0; end - if parallel_max == 0 a = perpvw_max / perpuv_max; b = perpuw_max / perpuv_max; @@ -76,71 +97,95 @@ else intersect_max = 0; end - +% Compute inverses of slopes if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) inv_dir = 1/ray_direction(1); else inv_dir = 1/ray_direction(3); end - -if intersect_max == 1 +% Compute intersection times +if intersect_max == 1 || collinear_max == 1 if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) - inv_dir = 1/ray_direction(1); t_max = (p_max(1) - ray_origin(1)) * inv_dir; else - inv_dir = 1/ray_direction(3); t_max = (p_max(2) - ray_origin(3)) * inv_dir; end end - -if intersect_min == 1 +if intersect_min == 1 || collinear_min == 1 if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) - inv_dir = 1/ray_direction(1); t_min = (p_min(1) - ray_origin(1)) * inv_dir; else - inv_dir = 1/ray_direction(3); t_min = (p_min(2) - ray_origin(3)) * inv_dir; end end - -if intersect_max == 1 && intersect_min == 0 && ... - t < t_max && t_max < t_end && ... - ~approximatelyEqual(t,t_max,1e-12,1e-8) +% Determine tStepPhi, tMaxPhi +if intersect_max == 1 && intersect_min == 0 && collinear_min == 0 &&... + strictlyLess(t,t_max,1e-12,1e-8) && strictlyLess(t_max,t_end,1e-12,1e-8) tStepPhi = 1; tMaxPhi = t_max; -elseif intersect_min == 1 && intersect_max == 0 && ... - t < t_min && t_min < t_end && ... - ~approximatelyEqual(t,t_min,1e-12,1e-8) +elseif intersect_min == 1 && intersect_max == 0 && collinear_max == 0 && ... + strictlyLess(t,t_min,1e-12,1e-8) && ... + strictlyLess(t_min,t_end,1e-12,1e-8) tStepPhi = -1; tMaxPhi = t_min; -elseif intersect_min == 1 && intersect_max == 1 +elseif (intersect_min == 1 && intersect_max == 1) || ... + (intersect_min == 1 && collinear_max == 1) || ... + (collinear_min == 1 && intersect_max == 1) % hit the min & max boundary simultaneously. - if approximatelyEqual(t_min,t_max,1e-12,1e-8) && ... - t < t_min && t_min < t_end && ... - ~approximatelyEqual(t,t_min,1e-12,1e-8) + if approximatelyEqual(t_min,t_max,1e-12,1e-8) && ... + strictlyLess(t,t_min,1e-12,1e-8) && ... + strictlyLess(t_min,t_end,1e-12,1e-8) tMaxPhi = t_max; if verbose fprintf("hit origin t\n") + end + pert_t = 0.1; + pert_x = sphere_center(1) + ray_direction(1) * pert_t; + pert_z = sphere_center(3) + ray_direction(3) * pert_t; + a = sphere_center(1) - pert_x; + b = sphere_center(2) - pert_z; + l = sqrt(a^2 + b^2); + sphere_center_xz = [sphere_center(1), sphere_center(3)]; + p1 = sphere_center_xz - (sphere_max_radius/l) .* [a b]; + i = 1; + while i < length(P_max) + d1 = (P_max(i,1)-p1(1))^2 + (P_max(i,2)-p1(2))^2; + d2 = (P_max(i+1,1)-p1(1))^2 + (P_max(i+1,2)-p1(2))^2; + d3 = (P_max(i,1)-P_max(i+1,1))^2 + (P_max(i,2)-P_max(i+1,2))^2; + if strictlyLess(d1+d2,d3,1e-12,1e-8) || ... + approximatelyEqual(d1+d2,d3,1e-12,1e-8) + next_vox = i - 1 ; + i = length(P_max); + end + i = i + 1; end - if ~approximatelyEqual(ray_direction(3),0.0,1e-12,1e-8) - % change in phi is dependent on the slope of the line - tStepPhi = - num_azimuthal_sections/2; + step = abs(current_voxel_ID_phi - next_vox); + if strictlyLess(ray_direction(3),0.0,1e-12,1e-8) + tStepPhi = step; + elseif strictlyLess(0.0,ray_direction(3),1e-12,1e-8) + tStepPhi = -step; else - tStepPhi = num_azimuthal_sections/2; + if strictlyLess(ray_direction(1),0.0,1e-12,1e-8) + tStepPhi = step; + else + tStepPhi = -step; + end end % hit min bound - elseif t < t_min && t_min < t_end && ... - (t_min < t_max || approximatelyEqual(t,t_max,1e-12,1e-8)) && ... - ~approximatelyEqual(t,t_min,1e-12,1e-8) + elseif strictlyLess(t,t_min,1e-12,1e-8) && ... + strictlyLess(t_min,t_end,1e-12,1e-8) && ... + (strictlyLess(t_min,t_max,1e-12,1e-8) || ... + approximatelyEqual(t,t_max,1e-12,1e-8)) tStepPhi = -1; tMaxPhi = t_min; if verbose fprintf("hit min bound\n") end % hit max bound - elseif t < t_max && t_max < t_end && ... - (t_max < t_min || approximatelyEqual(t,t_min,1e-12,1e-8)) && ... - ~approximatelyEqual(t,t_max,1e-12,1e-8) + elseif strictlyLess(t,t_max,1e-12,1e-8) && ... + strictlyLess(t_max,t_end,1e-12,1e-8) && ... + (strictlyLess(t_max,t_min,1e-12,1e-8) || ... + approximatelyEqual(t,t_min,1e-12,1e-8)) tStepPhi = 1; tMaxPhi = t_max; if verbose diff --git a/matlab/spherical-traversal/sphericalCoordinateTraversal.m b/matlab/spherical-traversal/sphericalCoordinateTraversal.m index 21abe9c..0d7a848 100644 --- a/matlab/spherical-traversal/sphericalCoordinateTraversal.m +++ b/matlab/spherical-traversal/sphericalCoordinateTraversal.m @@ -106,21 +106,23 @@ % the lines corresponding to angular voxels boundaries and the initial % radial voxel of the ray. Note that spherical coordinates are unnecessary, % this is just marking the voxel boundaries in the plane. -P_ang = []; -k_ang = 0; -while k_ang <= 2*pi - pt = [r * cos(k_ang) + sphere_center(1), r * sin(k_ang) + sphere_center(2)]; - P_ang = [P_ang ; pt]; - k_ang = k_ang + delta_theta; +k = 0; +trig_ang = []; +while k <= 2*pi + trig_ang = [trig_ang; cos(k), sin(k)]; + k = k + delta_theta; end +P_max_ang = sphere_max_radius .* trig_ang + [sphere_center(1) sphere_center(2)]; +P_ang = r .* trig_ang + [sphere_center(1) sphere_center(2)]; -P_azi = []; -k_azi = 0; -while k_azi <= 2*pi - pt = [r * cos(k_azi) + sphere_center(1), r * sin(k_azi) + sphere_center(3)]; - P_azi = [P_azi ; pt]; - k_azi = k_azi + delta_phi; +j=0; +trig_azi = []; +while j <= 2*pi + trig_azi = [trig_azi; cos(j), sin(j)]; + j = j + delta_phi; end +P_max_azi = sphere_max_radius .* trig_azi + [sphere_center(1) sphere_center(3)]; +P_azi = r .* trig_azi + [sphere_center(1) sphere_center(3)]; % Find the point of intersection between the vector created by the % ray intersection with the initial sphere radius and the sphere center. @@ -160,7 +162,8 @@ d1 = (P_ang(i,1)-p1_xy(1))^2 + (P_ang(i,2)-p1_xy(2))^2; d2 = (P_ang(i+1,1)-p1_xy(1))^2 + (P_ang(i+1,2)-p1_xy(2))^2; d3 = (P_ang(i,1)-P_ang(i+1,1))^2 + (P_ang(i,2)-P_ang(i+1,2))^2; - if d1 + d2 <= d3 + if strictlyLess(d1+d2,d3,1e-12,1e-8) || ... + approximatelyEqual(d1+d2,d3,1e-12,1e-8) current_voxel_ID_theta = i - 1; i = length(P_ang); end @@ -202,8 +205,9 @@ while i < length(P_azi) d1 = (P_azi(i,1)-p1_xz(1))^2 + (P_azi(i,2)-p1_xz(2))^2; d2 = (P_azi(i+1,1)-p1_xz(1))^2 + (P_azi(i+1,2)-p1_xz(2))^2; - d3 = (P_azi(i,1)-P_ang(i+1,1))^2 + (P_azi(i,2)-P_ang(i+1,2))^2; - if d1 + d2 <= d3 + d3 = (P_azi(i,1)-P_azi(i+1,1))^2 + (P_azi(i,2)-P_azi(i+1,2))^2; + if strictlyLess(d1+d2,d3,1e-12,1e-8) || ... + approximatelyEqual(d1+d2,d3,1e-12,1e-8) current_voxel_ID_phi = i - 1; i = length(P_azi); end @@ -256,20 +260,15 @@ t = t_start; t_end = min(t_grid, t_end); previous_transition_flag = false; -change_r = 0; while t < t_end % 1. Calculate tMaxR, tMaxTheta, tMaxPhi [tMaxR, tStepR, previous_transition_flag] = radial_hit(ray_origin, ray_direction, ... current_voxel_ID_r, sphere_center, sphere_max_radius, delta_radius, t, ray_unit_vector, ... ray_sphere_vector, v, previous_transition_flag, verbose); [tMaxTheta, tStepTheta] = angular_hit(ray_origin, ray_direction, current_voxel_ID_theta,... - num_angular_sections, sphere_center, ... - [P_ang(current_voxel_ID_theta+1,:) ; P_ang(current_voxel_ID_theta+2,:)], ... - t, t_end, verbose); + sphere_center, P_max_ang, sphere_max_radius, t, t_end, verbose); [tMaxPhi, tStepPhi] = azimuthal_hit(ray_origin, ray_direction, current_voxel_ID_phi,... - num_azimuthal_sections, sphere_center, ... - [P_azi(current_voxel_ID_phi+1,:) ; P_azi(current_voxel_ID_phi+2,:)], ... - t, t_end, verbose); + sphere_center, P_max_azi, sphere_max_radius, t, t_end, verbose); rStepViolation = current_voxel_ID_r + tStepR == 0; @@ -283,7 +282,6 @@ % angular boundary, we need the second half of disjunction t = tMaxTheta; current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections); - change_r = 0; elseif strictlyLess(tMaxR, tMaxTheta, 1e-12,1e-8) && ... strictlyLess(tMaxR, tMaxPhi, 1e-12,1e-8) && ... strictlyLess(t,tMaxR,1e-12,1e-8) && ... @@ -292,7 +290,6 @@ % Case tMaxR is minimum t = tMaxR; current_voxel_ID_r = current_voxel_ID_r + tStepR; - if tStepR ~= 0; change_r = 1; else; change_r = 0; end elseif strictlyLess(tMaxPhi,tMaxTheta,1e-12,1e-8) && ... strictlyLess(tMaxPhi,tMaxR,1e-12,1e-8) && ... strictlyLess(t,tMaxPhi,1e-12,1e-8) && ... @@ -300,7 +297,6 @@ % Case tMaxPhi is minimum t = tMaxPhi; current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections); - change_r = 0; elseif approximatelyEqual(tMaxPhi,tMaxTheta,1e-12,1e-8) && ... approximatelyEqual(tMaxPhi,tMaxR,1e-12,1e-8) && ... strictlyLess(t,tMaxR,1e-12,1e-8) && ... @@ -310,7 +306,6 @@ current_voxel_ID_r = current_voxel_ID_r + tStepR; 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); - if tStepR ~= 0; change_r = 1; else; change_r = 0; end elseif approximatelyEqual(tMaxPhi, tMaxTheta, 1e-12,1e-8) && ... strictlyLess(t,tMaxPhi,1e-12,1e-8) && ... strictlyLess(tMaxPhi,t_end, 1e-12,1e-8) @@ -320,12 +315,10 @@ % R min t = tMaxR; current_voxel_ID_r = current_voxel_ID_r + tStepR; - if tStepR ~= 0; change_r = 1; else; change_r = 0; end 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); - change_r = 0; end elseif approximatelyEqual(tMaxTheta,tMaxR,1e-12,1e-8) && ... strictlyLess(t,tMaxR,1e-12,1e-8) && ... @@ -337,12 +330,10 @@ % Phi min t = tMaxPhi; current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections); - change_r = 0; 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; - if tStepR ~= 0; change_r = 1; else; change_r = 0; end end elseif approximatelyEqual(tMaxR,tMaxPhi,1e-12,1e-8) && ... strictlyLess(t,tMaxR,1e-12,1e-8) && ... @@ -353,29 +344,14 @@ % Theta min t = tMaxTheta; current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections); - change_r = 0; 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; - if tStepR ~= 0; change_r = 1; else; change_r = 0; end end else return; end - % If the radial voxel changes, update the angular voxel boundary - % segments for intersection checks in angular hit. - if change_r == 1 - r = sphere_max_radius - delta_radius * (current_voxel_ID_r - 1); - a = [sphere_center(1) sphere_center(2)] - P_ang; - l1 = (a(:,1).^2 + a(:,2).^2).^-0.5; - c = [sphere_center(1) sphere_center(3)] - P_azi; - l2 = (c(:,1).^2 + c(:,2).^2).^-0.5; - P_ang(:,1) = sphere_center(1) - (r*l1) .* (sphere_center(1) - P_ang(:,1)); - P_ang(:,2) = sphere_center(2) - (r*l1) .* (sphere_center(2) - P_ang(:,2)); - P_azi(:,1) = sphere_center(1) - (r*l2) .* (sphere_center(1) - P_azi(:,1)); - P_azi(:,2) = sphere_center(3) - (r*l2) .* (sphere_center(3) - P_azi(:,2)); - 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]; diff --git a/matlab/spherical-traversal/sphericalCoordinateTraversal.m~ b/matlab/spherical-traversal/sphericalCoordinateTraversal.m~ new file mode 100644 index 0000000..577365e --- /dev/null +++ b/matlab/spherical-traversal/sphericalCoordinateTraversal.m~ @@ -0,0 +1,371 @@ +function [radial_voxels, angular_voxels, azimuthal_voxels] = 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) +% Input: +% min_bound: The lower left corner of the bounding box. +% max_bound: The upper right corner of the bounding box. +% ray origin: The origin of the ray in (x, y, z) coordinates. +% ray direction: The direction of the ray in (x, y, z) coordinates. +% sphere_center: The x, y location of the center of the sphere. +% sphere_max_radius: The largest that encompasses the sphere. +% num_radial_sections: The number of radial sections in the sphere. +% num_angular_sections: The number of angular sections in the sphere. +% The angular sections lie on the x-y plane. +% num_azimuthal_sections: The number of azimuthal sections in the sphere. +% The azimuthal sections lie on the x-z plane. +% t_begin: The beginning time of the ray. +% t_end: The end time of the ray. +% +% Requires: +% max_bound > min_bound +% The entire sphere is within min_bound and max_bound. +% t_end > t_begin >= 0.0 +% sphere_max_radius > 0 +% num_radial_sections > 0 +% num_angular_sections > 0 +% num_azimuthal_sections > 0 +% +% Returns: +% radial_voxels: A list of the radial voxels that were hit by the ray. +% angular_voxels: A list of the angular voxels that were hit by the ray. +% azimuthal_voxels: A list of the phi voxels that were hit by the ray. +% +% Note: These lists, used in conjunction, will produce the path of the ray +% through the voxels using each voxel transition. For example, +% [radial_voxels(1), angular_voxels(1), azimuthal_voxels(1)] +% is the first voxel the ray travels through. If the next traversal is a radial hit, +% angular_voxels(2) and azimuthal_voxels(2) will remain the same. + +angular_voxels = []; +radial_voxels = []; +azimuthal_voxels = []; + +% INITIALIZATION PHASE +delta_radius = sphere_max_radius / num_radial_sections; + +% Determine ray location at t_begin. +p = ray_origin + t_begin .* ray_direction; +ray_sphere_vector = [sphere_center(1) - p(1); sphere_center(2) - p(2); sphere_center(3) - p(3)]'; +% Find the radial shell containing the ray at t_begin. +r = delta_radius; +while (ray_sphere_vector(1)^2 + ray_sphere_vector(2)^2 + ray_sphere_vector(3)^2 > r^2) && r < sphere_max_radius + r = r + delta_radius; +end + +% Find the intersection times for the ray and the radial shell containing +% the ray at t_begin; in order to determine if the ray intersects the grid. +ray_unit_vector = 1 / sqrt(ray_direction(1)^2 + ray_direction(2)^2 + ray_direction(3)^2)... + .* [ray_direction(1); ray_direction(2); ray_direction(3)]'; +v = dot(ray_sphere_vector, ray_unit_vector); +discr = r^2 - (dot(ray_sphere_vector,ray_sphere_vector) - v^2); +d = sqrt(discr); +pa = ray_origin + (v-d) .* ray_unit_vector; +pb = ray_origin + (v+d) .* ray_unit_vector; +tol = 10^-16; +% Calculate the time of entrance and exit of the ray. +if ~approximatelyEqual(ray_direction(2),0.0,1e-12,1e-8) + % Use the y-direction if it is non-zero. + t1 = (pa(2) - ray_origin(2)) / ray_direction(2); + t2 = (pb(2) - ray_origin(2)) / ray_direction(2); +elseif ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) + % Use the x-direction if it is non-zero. + t1 = (pa(1) - ray_origin(1)) / ray_direction(1); + t2 = (pb(1) - ray_origin(1)) / ray_direction(1); +else + % Use the z-direction if it is non-zero. + t1 = (pa(3) - ray_origin(3)) / ray_direction(3); + t2 = (pb(3) - ray_origin(3)) / ray_direction(3); +end + +% The ray may not intersect the grid at all. +% In particular, if the ray is outside the grid at t_begin. +if t1 < t_begin && t2 < t_begin + if verbose + fprintf("\nRay does not intersect spherical grid for t_begin.") + end + return; +end + +% It may be a tangent hit. + if approximatelyEqual(t1,t2,1e-12,1e-8) + if verbose + fprintf("\nTangent hit.") + end + return; +end + +% If there is a ray/shell intersection, then set the radial voxel ID. +current_voxel_ID_r = 1 + (sphere_max_radius - r) / delta_radius; +if verbose + fprintf('RADIAL HIT INITIALIZED.\n') +end + +% Calculate Voxel ID Theta. +delta_theta = 2 * pi/ num_angular_sections; +delta_phi = 2 * pi/ num_azimuthal_sections; +% Create an array of values representing the points of intersection between +% the lines corresponding to angular voxels boundaries and the initial +% radial voxel of the ray. Note that spherical coordinates are unnecessary, +% this is just marking the voxel boundaries in the plane. +k = 0; +while k <= 2*pi + trig_ang = [trig_ang; cos(k), sin(k)]; + k = k + delta_theta; +end +P_max_ang = circle_max_radius .* trig_ang + circle_center; +P_ang = r .* trig_ang + circle_center; + +j=0; +while j <= 2*pi + trig_azi = [trig_azi; cos(j), sin(j)]; + j = j + delta_theta; +end +P_max_azi = circle_max_radius .* trig_azi + circle_center; +P_azi = r .* trig_azi + circle_center; + +% Find the point of intersection between the vector created by the +% ray intersection with the initial sphere radius and the sphere center. +ray_origin_xy = [ray_origin(1), ray_origin(2)]; +sphere_center_xy = [sphere_center(1), sphere_center(2)]; +if approximatelyEqual(ray_origin_xy,sphere_center_xy,1e-12,1e-8) + % 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(1) * pert_t; + pert_y = ray_origin(2) + ray_direction(2) * pert_t; + a = sphere_center(1) - pert_x; + b = sphere_center(2) - pert_y; +elseif approximatelyEqual(r,sphere_max_radius,1e-12,1e-8) +% If the ray origin is outside the grid, this will snap to the grid and use the +% resulting intersection point as the ray origin for this calculation. + a = sphere_center(1) - pa(1); + b = sphere_center(2) - pa(2); +else + a = sphere_center(1) - ray_origin(1); + b = sphere_center(2) - ray_origin(2); +end +l = sqrt(a^2 + b^2); +if approximatelyEqual(l,0.0,1e-12,1e-8) + % if l is 0.0 then only the z-dir of the ray is non-zero; set the + % initial point s.t. voxel_ID_theta init is 0. + p1_xy = [sphere_center(1) + r, sphere_center(2)]; +else + p1_xy = sphere_center_xy - (r/l) .* [a b]; +end +% This point will lie between two angular voxel boundaries iff the angle between +% it and the angular boundary intersection points along the circle of +% max radius is obtuse. Equality represents the case when the +% point lies on an angular boundary. +i = 1; +while i < length(P_ang) + d1 = (P_ang(i,1)-p1_xy(1))^2 + (P_ang(i,2)-p1_xy(2))^2; + d2 = (P_ang(i+1,1)-p1_xy(1))^2 + (P_ang(i+1,2)-p1_xy(2))^2; + d3 = (P_ang(i,1)-P_ang(i+1,1))^2 + (P_ang(i,2)-P_ang(i+1,2))^2; + if d1 + d2 <= d3 + current_voxel_ID_theta = i - 1; + i = length(P_ang); + end + i = i + 1; +end + +ray_origin_xz = [ray_origin(1), ray_origin(3)]; +sphere_center_xz = [sphere_center(1), sphere_center(3)]; +if approximatelyEqual(ray_origin_xz,sphere_center_xz,1e-12,1e-8) + % 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(1) * pert_t; + pert_z = ray_origin(3) + ray_direction(3) * pert_t; + a = sphere_center(1) - pert_x; + c = sphere_center(3) - pert_z; +elseif approximatelyEqual(r,sphere_max_radius,1e-12,1e-8) +% If the ray origin is outside the grid, this will snap to the grid and use the +% resulting intersection point as the ray origin for this calculation. + a = sphere_center(1) - pa(1); + c = sphere_center(3) - pa(3); +else + a = sphere_center(1) - ray_origin(1); + c = sphere_center(3) - ray_origin(3); +end +l = sqrt(a^2 + c^2); +if approximatelyEqual(l,0.0,1e-12,1e-8) + % if l is 0.0 then only the z-dir of the ray is non-zero; set the + % initial point s.t. voxel_ID_theta init is 0. + p1_xz = [sphere_center(1) + r, sphere_center(3)]; +else + p1_xz = sphere_center_xz - (r/l) .* [a c]; +end +% This point will lie between two angular voxel boundaries iff the angle between +% it and the angular boundary intersection points along the circle of +% max radius is obtuse. Equality represents the case when the +% point lies on an angular boundary. +i = 1; +while i < length(P_azi) + d1 = (P_azi(i,1)-p1_xz(1))^2 + (P_azi(i,2)-p1_xz(2))^2; + d2 = (P_azi(i+1,1)-p1_xz(1))^2 + (P_azi(i+1,2)-p1_xz(2))^2; + d3 = (P_azi(i,1)-P_azi(i+1,1))^2 + (P_azi(i,2)-P_azi(i+1,2))^2; + if d1 + d2 <= d3 + current_voxel_ID_phi = i - 1; + i = length(P_azi); + end + i = i + 1; +end + +azimuthal_voxels = [current_voxel_ID_phi]; +angular_voxels = [current_voxel_ID_theta]; +radial_voxels = [current_voxel_ID_r]; + +if verbose + fprintf('\nInitial phi voxel: %d', current_voxel_ID_phi) + fprintf('\nInitial theta voxel: %d', current_voxel_ID_theta) + fprintf('\nInitial radial voxel: %d \n', current_voxel_ID_r) +end + +% Find the maximum time the ray will be in the grid. +discr = sphere_max_radius^2 - (dot(ray_sphere_vector,ray_sphere_vector) - v^2); +d = sqrt(discr); +pa_max = ray_origin + (v-d) .* ray_unit_vector; +pb_max = ray_origin + (v+d) .* ray_unit_vector; +if ~approximatelyEqual(ray_direction(2),0.0,1e-12,1e-8) + t1 = (pa_max(2) - ray_origin(2)) / ray_direction(2); + t2 = (pb_max(2) - ray_origin(2)) / ray_direction(2); +elseif ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) + t1 = (pa_max(1) - ray_origin(1)) / ray_direction(1); + t2 = (pb_max(1) - ray_origin(1)) / ray_direction(1); +else + t1 = (pa_max(3) - ray_origin(3)) / ray_direction(3); + t2 = (pb_max(3) - ray_origin(3)) / ray_direction(3); +end +t_grid = max(t1,t2); + +% Determine the correct time to begin the traversal phase. If the ray +% starts outside the grid at t_begin, snap to the grid and find the +% corresopnding start time. +if approximatelyEqual(r,sphere_max_radius,1e-12,1e-8) + if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) + t_start = (p1_xy(1) - ray_origin(1))/ray_direction(1); + elseif ~approximatelyEqual(ray_direction(2),0.0,1e-12,1e-8) + t_start = (p1_xy(2) - ray_origin(2))/ray_direction(2); + else + t_start = (p1_xz(2) - ray_origin(3))/ray_direction(3); + end +else + t_start = t_begin; +end + +% TRAVERSAL PHASE +t = t_start; +t_end = min(t_grid, t_end); +previous_transition_flag = false; +while t < t_end + % 1. Calculate tMaxR, tMaxTheta, tMaxPhi + [tMaxR, tStepR, previous_transition_flag] = radial_hit(ray_origin, ray_direction, ... + current_voxel_ID_r, sphere_center, sphere_max_radius, delta_radius, t, ray_unit_vector, ... + ray_sphere_vector, v, previous_transition_flag, verbose); + [tMaxTheta, tStepTheta] = angular_hit(ray_origin, ray_direction, current_voxel_ID_theta,... + sphere_center, P_max_ang, circle_max_radius, t, t_end, verbose); + [tMaxPhi, tStepPhi] = azimuthal_hit(ray_origin, ray_direction, current_voxel_ID_phi,... + sphere_center, P_max_azi, circle_max_radius, t, t_end, verbose); + + rStepViolation = current_voxel_ID_r + tStepR == 0; + + % 2. Compare tMaxR, tMaxTheta, tMaxPhi + if ((strictlyLess(tMaxTheta,tMaxR,1e-12,1e-8) && ... + strictlyLess(tMaxR, tMaxPhi,1e-12,1e-8)) || rStepViolation) ... + && strictlyLess(t,tMaxTheta,1e-12,1e-8) && ... + strictlyLess(tMaxTheta,t_end,1e-12,1e-8) + % Note 1: Case tMaxTheta is a minimum + % Note 2: When the ray only intersects one radial shell but crosses an + % angular boundary, we need the second half of disjunction + t = tMaxTheta; + current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections); + elseif strictlyLess(tMaxR, tMaxTheta, 1e-12,1e-8) && ... + strictlyLess(tMaxR, tMaxPhi, 1e-12,1e-8) && ... + strictlyLess(t,tMaxR,1e-12,1e-8) && ... + strictlyLess(tMaxR,t_end, 1e-12,1e-8) ... + && ~rStepViolation + % Case tMaxR is minimum + t = tMaxR; + current_voxel_ID_r = current_voxel_ID_r + tStepR; + if tStepR ~= 0; change_r = 1; else; change_r = 0; end + elseif strictlyLess(tMaxPhi,tMaxTheta,1e-12,1e-8) && ... + strictlyLess(tMaxPhi,tMaxR,1e-12,1e-8) && ... + strictlyLess(t,tMaxPhi,1e-12,1e-8) && ... + strictlyLess(tMaxPhi,t_end,1e-12,1e-8) + % Case tMaxPhi is minimum + t = tMaxPhi; + current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections); + elseif approximatelyEqual(tMaxPhi,tMaxTheta,1e-12,1e-8) && ... + approximatelyEqual(tMaxPhi,tMaxR,1e-12,1e-8) && ... + strictlyLess(t,tMaxR,1e-12,1e-8) && ... + strictlyLess(tMaxR,t_end,1e-12,1e-8) && ~rStepViolation + % Triple boundary intersection + t = tMaxPhi; + current_voxel_ID_r = current_voxel_ID_r + tStepR; + 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); + if tStepR ~= 0; change_r = 1; else; change_r = 0; end + elseif approximatelyEqual(tMaxPhi, tMaxTheta, 1e-12,1e-8) && ... + strictlyLess(t,tMaxPhi,1e-12,1e-8) && ... + strictlyLess(tMaxPhi,t_end, 1e-12,1e-8) + % Phi, Theta equal + if strictlyLess(tMaxR, tMaxPhi, 1e-12,1e-8) && ... + strictlyLess(t,tMaxR,1e-12,1e-8) && ~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 approximatelyEqual(tMaxTheta,tMaxR,1e-12,1e-8) && ... + strictlyLess(t,tMaxR,1e-12,1e-8) && ... + strictlyLess(tMaxR, t_end,1e-12,1e-8) && ... + ~rStepViolation + % R, Theta equal + if strictlyLess(tMaxPhi,tMaxTheta,1e-12,1e-8) && ... + strictlyLess(t,tMaxPhi,1e-12,1e-8) + % 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 approximatelyEqual(tMaxR,tMaxPhi,1e-12,1e-8) && ... + strictlyLess(t,tMaxR,1e-12,1e-8) && ... + strictlyLess(tMaxR, t_end, 1e-12,1e-8) && ~rStepViolation + % R, Phi equal + if strictlyLess(tMaxTheta,tMaxR,1e-12,1e-8) && ... + strictlyLess(t,tMaxTheta, 1e-12,1e-8) + % 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 + % If the radial voxel changes, update the angular voxel boundary + % segments for intersection checks in angular hit. + if change_r == 1 + r = sphere_max_radius - delta_radius * (current_voxel_ID_r - 1); + a = [sphere_center(1) sphere_center(2)] - P_ang; + l1 = (a(:,1).^2 + a(:,2).^2).^-0.5; + c = [sphere_center(1) sphere_center(3)] - P_azi; + l2 = (c(:,1).^2 + c(:,2).^2).^-0.5; + P_ang(:,1) = sphere_center(1) - (r*l1) .* (sphere_center(1) - P_ang(:,1)); + P_ang(:,2) = sphere_center(2) - (r*l1) .* (sphere_center(2) - P_ang(:,2)); + P_azi(:,1) = sphere_center(1) - (r*l2) .* (sphere_center(1) - P_azi(:,1)); + P_azi(:,2) = sphere_center(3) - (r*l2) .* (sphere_center(3) - P_azi(:,2)); + 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 + +end diff --git a/matlab/spherical-traversal/sphericalCoordinateTraversalTest.m b/matlab/spherical-traversal/sphericalCoordinateTraversalTest.m index 1a1d314..23c7671 100644 --- a/matlab/spherical-traversal/sphericalCoordinateTraversalTest.m +++ b/matlab/spherical-traversal/sphericalCoordinateTraversalTest.m @@ -66,11 +66,16 @@ function testRayOffsetInXYPlane(testCase) num_radial_sections = 4; num_angular_sections = 4; - num_azimuthal_sections = 4 + num_azimuthal_sections = 4; + t_begin = 0.0; + t_end = 30.0; + verbose = false; + + [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); - verifyEqual(testCase, rVoxels, [1,2,2,3,2,2,2,1]); - verifyEqual(testCase, thetaVoxels, [2,2,1,1,1,1,0,0]); - verifyEqual(testCase, phiVoxels, [2,2,2,2,2,0,0,0]); + verifyEqual(testCase, rVoxels, [1,2,2,3,2,2,1]); + verifyEqual(testCase, thetaVoxels, [2,2,1,1,1,0,0]); + verifyEqual(testCase, phiVoxels, [2,2,2,2,2,0,0]); end % Ray travels up z-axis @@ -91,9 +96,11 @@ function testRayTravelsAlongZAxis(testCase) [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); - verifyEqual(testCase, rVoxels, [1,2,3,4,3,2,1]); - verifyEqual(testCase, thetaVoxels, [0,0,0,0,0,0,0]); - verifyEqual(testCase, phiVoxels, [3,3,3,3,3,3,3]); + verifyEqual(testCase, rVoxels, [1,2,3,4,4,3,2,1]); + verifyEqual(testCase, thetaVoxels, [0,0,0,0,0,0,0,0]); + % Note that phi_voxels transition is somewhat ambiguous, as long as + % voxels both lie along the z-axis this should be counted as correct + verifyEqual(testCase, phiVoxels, [2,2,2,2,0,0,0,0]); end % Ray travels along x-axis @@ -114,9 +121,9 @@ function testRayTravelsAlongXAxis(testCase) [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); - verifyEqual(testCase, rVoxels, [1,2,3,4,3,2,1]); - verifyEqual(testCase, thetaVoxels, [1,1,1,1,0,0,0]); - verifyEqual(testCase, phiVoxels, [1,1,1,1,0,0,0]); + verifyEqual(testCase, rVoxels, [1,2,3,4,4,3,2,1]); + verifyEqual(testCase, thetaVoxels, [3,3,3,3,0,0,0,0]); + verifyEqual(testCase, phiVoxels, [1,1,1,1,0,0,0,0]); end % Ray travels along y-axis @@ -137,9 +144,9 @@ function testRayTravelsAlongYAxis(testCase) [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); - verifyEqual(testCase, rVoxels, [1,2,3,4,3,2,1]); - verifyEqual(testCase, thetaVoxels, [1,1,1,1,0,0,0]); - verifyEqual(testCase, phiVoxels, [0,0,0,0,0,0,0]); + verifyEqual(testCase, rVoxels, [1,2,3,4,4,3,2,1]); + verifyEqual(testCase, thetaVoxels, [5,5,5,5,1,1,1,1]); + verifyEqual(testCase, phiVoxels, [0,0,0,0,0,0,0,0]); end % Ray parallel to the XY Plane. @@ -162,7 +169,7 @@ function testRayParallelToXYPlane(testCase) 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,4,4,3,2,1]); verifyEqual(testCase, thetaVoxels, [2,2,2,2,0,0,0,0]); - verifyEqual(testCase, phiVoxels, [2,2,2,2,0,0,0,0]); + verifyEqual(testCase, phiVoxels, [1,1,1,1,0,0,0,0]); end % Ray parallel to the XZ Plane. @@ -184,7 +191,7 @@ function testRayParallelToXZPlane(testCase) [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); verifyEqual(testCase, rVoxels, [1,2,3,4,4,3,2,1]); - verifyEqual(testCase, thetaVoxels, [2,2,2,2,0,0,0,0]); + verifyEqual(testCase, thetaVoxels, [1,1,1,1,0,0,0,0]); verifyEqual(testCase, phiVoxels, [2,2,2,2,0,0,0,0]); end @@ -207,8 +214,8 @@ function testRayParallelToYZPlane(testCase) [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); verifyEqual(testCase, rVoxels, [1,2,3,4,4,3,2,1]); - verifyEqual(testCase, thetaVoxels, [3,3,3,3,0,0,0,0]); - verifyEqual(testCase, phiVoxels, [3,3,3,3,0,0,0,0]); + verifyEqual(testCase, thetaVoxels, [2,2,2,2,0,0,0,0]); + verifyEqual(testCase, phiVoxels, [2,2,2,2,0,0,0,0]); end % Ray with negative X positive YZ direction @@ -229,13 +236,13 @@ function testRayNegXPosYZ(testCase) [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); - verifyEqual(testCase, rVoxels, [1,2,3,4,4,4,3,2,1]); - verifyEqual(testCase, thetaVoxels, [3,3,3,3,3,2,1,1,1]); - verifyEqual(testCase, phiVoxels, [3,3,3,3,3,2,1,1,1]); + verifyEqual(testCase, rVoxels, [1,2,3,3,4,4,3,2,1]); + verifyEqual(testCase, thetaVoxels, [3,3,3,2,2,1,1,1,1]); + verifyEqual(testCase, phiVoxels, [3,3,3,2,2,1,1,1,1]); end % Ray with negative Y positive XZ direction -function testRayNegXPosYZ(testCase) +function testRayNegYPosXZ(testCase) min_bound = [-20, -20.0, -20.0]; max_bound = [20.0, 20.0, 20.0]; ray_origin = [-13.0, 17.0, -15.0]; @@ -252,9 +259,9 @@ function testRayNegXPosYZ(testCase) [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); - verifyEqual(testCase, rVoxels, [1,2,3,4,4,4,3,2,1]); - verifyEqual(testCase, thetaVoxels, [1,1,1,1,0,3,3,3,3]); - verifyEqual(testCase, phiVoxels, [2,2,2,2,1,0,0,0,0]); + verifyEqual(testCase, rVoxels, [1,2,3,3,4,4,3,3,2,1]); + verifyEqual(testCase, thetaVoxels, [1,1,1,1,1,0,0,3,3,3]); + verifyEqual(testCase, phiVoxels, [2,2,2,1,1,0,0,0,0,0]); end % Ray with negative Z positive XY direction, NOT COMPLETED @@ -275,9 +282,9 @@ function testRayNegZPosXY(testCase) [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); - verifyEqual(testCase, rVoxels, [1,2,3,3,3,2,1]); - verifyEqual(testCase, thetaVoxels, [2,2,2,1,1,0,0]); - verifyEqual(testCase, phiVoxels, [2,2,2,2,1,0,0,0,0]); + verifyEqual(testCase, rVoxels, [1,1,2,2,1]); + verifyEqual(testCase, thetaVoxels, [2,1,1,0,0]); + verifyEqual(testCase, phiVoxels, [1,1,1,0,0]); end % Ray with negative XY positive Z direction From 34ce1ed8822f9cef73acff1268de4b73efede8e2 Mon Sep 17 00:00:00 2001 From: Ariel Kellison Date: Fri, 17 Apr 2020 14:41:39 -0400 Subject: [PATCH 2/3] Delete angular_hit.m~ --- matlab/spherical-traversal/angular_hit.m~ | 212 ---------------------- 1 file changed, 212 deletions(-) delete mode 100644 matlab/spherical-traversal/angular_hit.m~ diff --git a/matlab/spherical-traversal/angular_hit.m~ b/matlab/spherical-traversal/angular_hit.m~ deleted file mode 100644 index af5777b..0000000 --- a/matlab/spherical-traversal/angular_hit.m~ +++ /dev/null @@ -1,212 +0,0 @@ -function [tMaxTheta, tStepTheta] = angular_hit(ray_origin, ray_direction, current_voxel_ID_theta,... - sphere_center, P_max, sphere_max_radius, t, t_end, verbose) -% Determines whether an angular hit occurs for the given ray. -% Input: -% ray_origin: vector of the origin of the ray in cartesian coordinate. -% ray_direction: vector of the direction of the ray in cartesian coordinate. -% current_voxel_ID_theta: the (angular) ID of current voxel. -% num_angular_sections: number of total angular sections on the grid. -% sphere_center: The center of the circle. -% pointArray: -% t: The current time parameter of the ray. -% t_end: -% verbose: Determines whether debug print statements are enabled. -% -% Returns: -% tMaxTheta: is the time at which a hit occurs for the ray at the next point of intersection. -% tStepTheta: The direction the theta voxel steps. +1, -1, or 0. -if verbose - fprintf("\n-- angular_hit --") -end - -if verbose - fprintf("\nCurrent Voxel ID Theta: %d\n", current_voxel_ID_theta) -end -pointArray = [P_max(current_voxel_ID_theta+1,:) ; P_max(current_voxel_ID_theta+2,:)]; -% Get current ray location to use as start point -p = ray_origin + ray_direction * t; -% Get ray exit point to use as end point -p_f = ray_origin + ray_direction * t_end; -% Caculate the ray segment vector -v = p_f - p; -% Calculate the voxel boundary vectors -u_min = [sphere_center(1) sphere_center(2)] - pointArray(1,:); -u_max = [sphere_center(1) sphere_center(2)] - pointArray(2,:); -% http://geomalgorithms.com/a05-_intersect-1.html#intersect2D_2Segments() -w_min = pointArray(1,:) - [p(1) p(2)]; -w_max = pointArray(2,:) - [p(1) p(2)]; -perpuv_min = u_min(1) * v(2) - u_min(2) * v(1); -perpuv_max = u_max(1) * v(2) - u_max(2) * v(1); -perpuw_min = u_min(1) * w_min(2) - u_min(2) * w_min(1); -perpuw_max = u_max(1) * w_max(2) - u_max(2) * w_max(1); -perpvw_min = v(1) * w_min(2) - v(2) * w_min(1); -perpvw_max = v(1) * w_max(2) - v(2) * w_max(1); - -% check to see if ray is parallel to each angular voxel boundary -if approximatelyEqual(perpuv_min,0.0,1e-12,1e-8) - parallel_min = 1; - % check to see if the ray is colliner with the angular voxel boundary - if approximatelyEqual(perpuw_min,0.0,1e-12,1e-8) && ... - approximatelyEqual(perpvw_min,0.0,1e-12,1e-8) - collinear_min = 1; - p_min = sphere_center; - else - collinear_min = 0; - end -else - parallel_min = 0; - collinear_min = 0; -end -if approximatelyEqual(perpuv_max,0.0,1e-12,1e-8) - parallel_max = 1; - if approximatelyEqual(perpuw_max,0.0,1e-12,1e-8) && ... - approximatelyEqual(perpvw_max,0.0,1e-12,1e-8) - collinear_max = 1; - p_max = sphere_center; - else - collinear_max = 0; - end -else - parallel_max = 0; - collinear_max = 0; -end - -if parallel_min == 0 - a = perpvw_min / perpuv_min; - b = perpuw_min / perpuv_min; - if (strictlyLess(a,0.0,1e-12,1e-8) || strictlyLess(1.0,a,1e-12,1e-8)) || ... - (strictlyLess(b,0.0,1e-12,1e-8) || strictlyLess(1.0,b,1e-12,1e-8)) - intersect_min = 0; - else - intersect_min = 1; - p_min = p + b * v; - end -else - intersect_min = 0; -end - -if parallel_max == 0 - a = perpvw_max / perpuv_max; - b = perpuw_max / perpuv_max; - if (strictlyLess(a,0.0,1e-12,1e-8) || strictlyLess(1.0,a,1e-12,1e-8)) || ... - (strictlyLess(b,0.0,1e-12,1e-8) || strictlyLess(1.0,b,1e-12,1e-8)) - intersect_max = 0; - else - intersect_max = 1; - p_max = p + b * v; - end -else - intersect_max = 0; -end - -if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) - inv_dir = 1/ray_direction(1); -else - inv_dir = 1/ray_direction(2); -end - -if intersect_max == 1 || collinear_max == 1 - if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) - t_max = (p_max(1) - ray_origin(1)) * inv_dir; - else - t_max = (p_max(2) - ray_origin(2)) * inv_dir; - end -end - -if intersect_min == 1 || collinear_min == 1 - if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) - t_min = (p_min(1) - ray_origin(1)) * inv_dir; - else - t_min = (p_min(2) - ray_origin(2)) * inv_dir; - end -end - -if intersect_max == 1 && intersect_min == 0 && collinear_min == 0 &&... - strictlyLess(t,t_max,1e-12,1e-8) && strictlyLess(t_max,t_end,1e-12,1e-8) - tStepTheta = 1; - tMaxTheta = t_max; -elseif intersect_min == 1 && intersect_max == 0 && collinear_max == 0 && ... - strictlyLess(t,t_min,1e-12,1e-8) && ... - strictlyLess(t_min,t_end,1e-12,1e-8) - tStepTheta = -1; - tMaxTheta = t_min; -elseif (intersect_min == 1 && intersect_max == 1) || ... - (intersect_min == 1 && collinear_max == 1) || ... - (collinear_min == 1 && intersect_max == 1) - % hit the min & max boundary simultaneously. - if approximatelyEqual(t_min,t_max,1e-12,1e-8) && ... - strictlyLess(t,t_min,1e-12,1e-8) && ... - strictlyLess(t_min,t_end,1e-12,1e-8) - tMaxTheta = t_max; - if verbose - fprintf("hit origin t\n") - end - pert_t = 0.1; - pert_x = sphere_center(1) + ray_direction(1) * pert_t; - pert_y = sphere_center(2) + ray_direction(2) * pert_t; - a = sphere_center(1) - pert_x; - b = sphere_center(2) - pert_y; - l = sqrt(a^2 + b^2); - p1 = sphere_center - (sphere_max_radius/l) .* [a b]; - i = 1; - while i < length(P_max) - d1 = (P_max(i,1)-p1(1))^2 + (P_max(i,2)-p1(2))^2; - d2 = (P_max(i+1,1)-p1(1))^2 + (P_max(i+1,2)-p1(2))^2; - d3 = (P_max(i,1)-P_max(i+1,1))^2 + (P_max(i,2)-P_max(i+1,2))^2; - if strictlyLess(d1+d2,d3,1e-12,1e-8) || ... - approximatelyEqual(d1+d2,d3,1e-12,1e-8) - next_vox = i - 1 ; - i = length(P_max); - end - i = i + 1; - end - step = abs(current_voxel_ID_theta - next_vox); - if strictlyLess(ray_direction(2),0.0,1e-12,1e-8) - tStepTheta = step; - elseif strictlyLess(0.0,ray_direction(2),1e-12,1e-8) - tStepTheta = -step; - else - if strictlyLess(ray_direction(1),0.0,1e-12,1e-8) - tStepTheta = step; - else - tStepTheta = -step; - end - end - elseif t < t_min && t_min < t_end && ... - (t_min < t_max || approximatelyEqual(t,t_max,1e-12,1e-8)) && ... - ~approximatelyEqual(t,t_min,1e-12,1e-8) - tStepTheta = -1; - tMaxTheta = t_min; - if verbose - fprintf("hit min bound\n") - end - % hit max bound - elseif t < t_max && t_max < t_end && ... - (t_max < t_min || approximatelyEqual(t,t_min,1e-12,1e-8)) && ... - ~approximatelyEqual(t,t_max,1e-12,1e-8) - tStepTheta = 1; - tMaxTheta = t_max; - if verbose - fprintf("hit max bound\n") - end - else - tStepTheta = 0; - tMaxTheta = t_end; - if verbose - fprintf("no hit within time bound\n") - end - end - else - tStepTheta = 0; - tMaxTheta = t_end; - if verbose - fprintf("no hit within voxel bound\n") - end - end - - -if verbose - fprintf(['\ntMaxTheta: %d \n' ... - 'tStepTheta: %d \n'], tMaxTheta, tStepTheta); -end -end From fc3f01069000de39d8e8dff13ccbc22e4d9a7073 Mon Sep 17 00:00:00 2001 From: Ariel Kellison Date: Fri, 17 Apr 2020 14:41:50 -0400 Subject: [PATCH 3/3] Delete sphericalCoordinateTraversal.m~ --- .../sphericalCoordinateTraversal.m~ | 371 ------------------ 1 file changed, 371 deletions(-) delete mode 100644 matlab/spherical-traversal/sphericalCoordinateTraversal.m~ diff --git a/matlab/spherical-traversal/sphericalCoordinateTraversal.m~ b/matlab/spherical-traversal/sphericalCoordinateTraversal.m~ deleted file mode 100644 index 577365e..0000000 --- a/matlab/spherical-traversal/sphericalCoordinateTraversal.m~ +++ /dev/null @@ -1,371 +0,0 @@ -function [radial_voxels, angular_voxels, azimuthal_voxels] = 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) -% Input: -% min_bound: The lower left corner of the bounding box. -% max_bound: The upper right corner of the bounding box. -% ray origin: The origin of the ray in (x, y, z) coordinates. -% ray direction: The direction of the ray in (x, y, z) coordinates. -% sphere_center: The x, y location of the center of the sphere. -% sphere_max_radius: The largest that encompasses the sphere. -% num_radial_sections: The number of radial sections in the sphere. -% num_angular_sections: The number of angular sections in the sphere. -% The angular sections lie on the x-y plane. -% num_azimuthal_sections: The number of azimuthal sections in the sphere. -% The azimuthal sections lie on the x-z plane. -% t_begin: The beginning time of the ray. -% t_end: The end time of the ray. -% -% Requires: -% max_bound > min_bound -% The entire sphere is within min_bound and max_bound. -% t_end > t_begin >= 0.0 -% sphere_max_radius > 0 -% num_radial_sections > 0 -% num_angular_sections > 0 -% num_azimuthal_sections > 0 -% -% Returns: -% radial_voxels: A list of the radial voxels that were hit by the ray. -% angular_voxels: A list of the angular voxels that were hit by the ray. -% azimuthal_voxels: A list of the phi voxels that were hit by the ray. -% -% Note: These lists, used in conjunction, will produce the path of the ray -% through the voxels using each voxel transition. For example, -% [radial_voxels(1), angular_voxels(1), azimuthal_voxels(1)] -% is the first voxel the ray travels through. If the next traversal is a radial hit, -% angular_voxels(2) and azimuthal_voxels(2) will remain the same. - -angular_voxels = []; -radial_voxels = []; -azimuthal_voxels = []; - -% INITIALIZATION PHASE -delta_radius = sphere_max_radius / num_radial_sections; - -% Determine ray location at t_begin. -p = ray_origin + t_begin .* ray_direction; -ray_sphere_vector = [sphere_center(1) - p(1); sphere_center(2) - p(2); sphere_center(3) - p(3)]'; -% Find the radial shell containing the ray at t_begin. -r = delta_radius; -while (ray_sphere_vector(1)^2 + ray_sphere_vector(2)^2 + ray_sphere_vector(3)^2 > r^2) && r < sphere_max_radius - r = r + delta_radius; -end - -% Find the intersection times for the ray and the radial shell containing -% the ray at t_begin; in order to determine if the ray intersects the grid. -ray_unit_vector = 1 / sqrt(ray_direction(1)^2 + ray_direction(2)^2 + ray_direction(3)^2)... - .* [ray_direction(1); ray_direction(2); ray_direction(3)]'; -v = dot(ray_sphere_vector, ray_unit_vector); -discr = r^2 - (dot(ray_sphere_vector,ray_sphere_vector) - v^2); -d = sqrt(discr); -pa = ray_origin + (v-d) .* ray_unit_vector; -pb = ray_origin + (v+d) .* ray_unit_vector; -tol = 10^-16; -% Calculate the time of entrance and exit of the ray. -if ~approximatelyEqual(ray_direction(2),0.0,1e-12,1e-8) - % Use the y-direction if it is non-zero. - t1 = (pa(2) - ray_origin(2)) / ray_direction(2); - t2 = (pb(2) - ray_origin(2)) / ray_direction(2); -elseif ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) - % Use the x-direction if it is non-zero. - t1 = (pa(1) - ray_origin(1)) / ray_direction(1); - t2 = (pb(1) - ray_origin(1)) / ray_direction(1); -else - % Use the z-direction if it is non-zero. - t1 = (pa(3) - ray_origin(3)) / ray_direction(3); - t2 = (pb(3) - ray_origin(3)) / ray_direction(3); -end - -% The ray may not intersect the grid at all. -% In particular, if the ray is outside the grid at t_begin. -if t1 < t_begin && t2 < t_begin - if verbose - fprintf("\nRay does not intersect spherical grid for t_begin.") - end - return; -end - -% It may be a tangent hit. - if approximatelyEqual(t1,t2,1e-12,1e-8) - if verbose - fprintf("\nTangent hit.") - end - return; -end - -% If there is a ray/shell intersection, then set the radial voxel ID. -current_voxel_ID_r = 1 + (sphere_max_radius - r) / delta_radius; -if verbose - fprintf('RADIAL HIT INITIALIZED.\n') -end - -% Calculate Voxel ID Theta. -delta_theta = 2 * pi/ num_angular_sections; -delta_phi = 2 * pi/ num_azimuthal_sections; -% Create an array of values representing the points of intersection between -% the lines corresponding to angular voxels boundaries and the initial -% radial voxel of the ray. Note that spherical coordinates are unnecessary, -% this is just marking the voxel boundaries in the plane. -k = 0; -while k <= 2*pi - trig_ang = [trig_ang; cos(k), sin(k)]; - k = k + delta_theta; -end -P_max_ang = circle_max_radius .* trig_ang + circle_center; -P_ang = r .* trig_ang + circle_center; - -j=0; -while j <= 2*pi - trig_azi = [trig_azi; cos(j), sin(j)]; - j = j + delta_theta; -end -P_max_azi = circle_max_radius .* trig_azi + circle_center; -P_azi = r .* trig_azi + circle_center; - -% Find the point of intersection between the vector created by the -% ray intersection with the initial sphere radius and the sphere center. -ray_origin_xy = [ray_origin(1), ray_origin(2)]; -sphere_center_xy = [sphere_center(1), sphere_center(2)]; -if approximatelyEqual(ray_origin_xy,sphere_center_xy,1e-12,1e-8) - % 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(1) * pert_t; - pert_y = ray_origin(2) + ray_direction(2) * pert_t; - a = sphere_center(1) - pert_x; - b = sphere_center(2) - pert_y; -elseif approximatelyEqual(r,sphere_max_radius,1e-12,1e-8) -% If the ray origin is outside the grid, this will snap to the grid and use the -% resulting intersection point as the ray origin for this calculation. - a = sphere_center(1) - pa(1); - b = sphere_center(2) - pa(2); -else - a = sphere_center(1) - ray_origin(1); - b = sphere_center(2) - ray_origin(2); -end -l = sqrt(a^2 + b^2); -if approximatelyEqual(l,0.0,1e-12,1e-8) - % if l is 0.0 then only the z-dir of the ray is non-zero; set the - % initial point s.t. voxel_ID_theta init is 0. - p1_xy = [sphere_center(1) + r, sphere_center(2)]; -else - p1_xy = sphere_center_xy - (r/l) .* [a b]; -end -% This point will lie between two angular voxel boundaries iff the angle between -% it and the angular boundary intersection points along the circle of -% max radius is obtuse. Equality represents the case when the -% point lies on an angular boundary. -i = 1; -while i < length(P_ang) - d1 = (P_ang(i,1)-p1_xy(1))^2 + (P_ang(i,2)-p1_xy(2))^2; - d2 = (P_ang(i+1,1)-p1_xy(1))^2 + (P_ang(i+1,2)-p1_xy(2))^2; - d3 = (P_ang(i,1)-P_ang(i+1,1))^2 + (P_ang(i,2)-P_ang(i+1,2))^2; - if d1 + d2 <= d3 - current_voxel_ID_theta = i - 1; - i = length(P_ang); - end - i = i + 1; -end - -ray_origin_xz = [ray_origin(1), ray_origin(3)]; -sphere_center_xz = [sphere_center(1), sphere_center(3)]; -if approximatelyEqual(ray_origin_xz,sphere_center_xz,1e-12,1e-8) - % 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(1) * pert_t; - pert_z = ray_origin(3) + ray_direction(3) * pert_t; - a = sphere_center(1) - pert_x; - c = sphere_center(3) - pert_z; -elseif approximatelyEqual(r,sphere_max_radius,1e-12,1e-8) -% If the ray origin is outside the grid, this will snap to the grid and use the -% resulting intersection point as the ray origin for this calculation. - a = sphere_center(1) - pa(1); - c = sphere_center(3) - pa(3); -else - a = sphere_center(1) - ray_origin(1); - c = sphere_center(3) - ray_origin(3); -end -l = sqrt(a^2 + c^2); -if approximatelyEqual(l,0.0,1e-12,1e-8) - % if l is 0.0 then only the z-dir of the ray is non-zero; set the - % initial point s.t. voxel_ID_theta init is 0. - p1_xz = [sphere_center(1) + r, sphere_center(3)]; -else - p1_xz = sphere_center_xz - (r/l) .* [a c]; -end -% This point will lie between two angular voxel boundaries iff the angle between -% it and the angular boundary intersection points along the circle of -% max radius is obtuse. Equality represents the case when the -% point lies on an angular boundary. -i = 1; -while i < length(P_azi) - d1 = (P_azi(i,1)-p1_xz(1))^2 + (P_azi(i,2)-p1_xz(2))^2; - d2 = (P_azi(i+1,1)-p1_xz(1))^2 + (P_azi(i+1,2)-p1_xz(2))^2; - d3 = (P_azi(i,1)-P_azi(i+1,1))^2 + (P_azi(i,2)-P_azi(i+1,2))^2; - if d1 + d2 <= d3 - current_voxel_ID_phi = i - 1; - i = length(P_azi); - end - i = i + 1; -end - -azimuthal_voxels = [current_voxel_ID_phi]; -angular_voxels = [current_voxel_ID_theta]; -radial_voxels = [current_voxel_ID_r]; - -if verbose - fprintf('\nInitial phi voxel: %d', current_voxel_ID_phi) - fprintf('\nInitial theta voxel: %d', current_voxel_ID_theta) - fprintf('\nInitial radial voxel: %d \n', current_voxel_ID_r) -end - -% Find the maximum time the ray will be in the grid. -discr = sphere_max_radius^2 - (dot(ray_sphere_vector,ray_sphere_vector) - v^2); -d = sqrt(discr); -pa_max = ray_origin + (v-d) .* ray_unit_vector; -pb_max = ray_origin + (v+d) .* ray_unit_vector; -if ~approximatelyEqual(ray_direction(2),0.0,1e-12,1e-8) - t1 = (pa_max(2) - ray_origin(2)) / ray_direction(2); - t2 = (pb_max(2) - ray_origin(2)) / ray_direction(2); -elseif ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) - t1 = (pa_max(1) - ray_origin(1)) / ray_direction(1); - t2 = (pb_max(1) - ray_origin(1)) / ray_direction(1); -else - t1 = (pa_max(3) - ray_origin(3)) / ray_direction(3); - t2 = (pb_max(3) - ray_origin(3)) / ray_direction(3); -end -t_grid = max(t1,t2); - -% Determine the correct time to begin the traversal phase. If the ray -% starts outside the grid at t_begin, snap to the grid and find the -% corresopnding start time. -if approximatelyEqual(r,sphere_max_radius,1e-12,1e-8) - if ~approximatelyEqual(ray_direction(1),0.0,1e-12,1e-8) - t_start = (p1_xy(1) - ray_origin(1))/ray_direction(1); - elseif ~approximatelyEqual(ray_direction(2),0.0,1e-12,1e-8) - t_start = (p1_xy(2) - ray_origin(2))/ray_direction(2); - else - t_start = (p1_xz(2) - ray_origin(3))/ray_direction(3); - end -else - t_start = t_begin; -end - -% TRAVERSAL PHASE -t = t_start; -t_end = min(t_grid, t_end); -previous_transition_flag = false; -while t < t_end - % 1. Calculate tMaxR, tMaxTheta, tMaxPhi - [tMaxR, tStepR, previous_transition_flag] = radial_hit(ray_origin, ray_direction, ... - current_voxel_ID_r, sphere_center, sphere_max_radius, delta_radius, t, ray_unit_vector, ... - ray_sphere_vector, v, previous_transition_flag, verbose); - [tMaxTheta, tStepTheta] = angular_hit(ray_origin, ray_direction, current_voxel_ID_theta,... - sphere_center, P_max_ang, circle_max_radius, t, t_end, verbose); - [tMaxPhi, tStepPhi] = azimuthal_hit(ray_origin, ray_direction, current_voxel_ID_phi,... - sphere_center, P_max_azi, circle_max_radius, t, t_end, verbose); - - rStepViolation = current_voxel_ID_r + tStepR == 0; - - % 2. Compare tMaxR, tMaxTheta, tMaxPhi - if ((strictlyLess(tMaxTheta,tMaxR,1e-12,1e-8) && ... - strictlyLess(tMaxR, tMaxPhi,1e-12,1e-8)) || rStepViolation) ... - && strictlyLess(t,tMaxTheta,1e-12,1e-8) && ... - strictlyLess(tMaxTheta,t_end,1e-12,1e-8) - % Note 1: Case tMaxTheta is a minimum - % Note 2: When the ray only intersects one radial shell but crosses an - % angular boundary, we need the second half of disjunction - t = tMaxTheta; - current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections); - elseif strictlyLess(tMaxR, tMaxTheta, 1e-12,1e-8) && ... - strictlyLess(tMaxR, tMaxPhi, 1e-12,1e-8) && ... - strictlyLess(t,tMaxR,1e-12,1e-8) && ... - strictlyLess(tMaxR,t_end, 1e-12,1e-8) ... - && ~rStepViolation - % Case tMaxR is minimum - t = tMaxR; - current_voxel_ID_r = current_voxel_ID_r + tStepR; - if tStepR ~= 0; change_r = 1; else; change_r = 0; end - elseif strictlyLess(tMaxPhi,tMaxTheta,1e-12,1e-8) && ... - strictlyLess(tMaxPhi,tMaxR,1e-12,1e-8) && ... - strictlyLess(t,tMaxPhi,1e-12,1e-8) && ... - strictlyLess(tMaxPhi,t_end,1e-12,1e-8) - % Case tMaxPhi is minimum - t = tMaxPhi; - current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections); - elseif approximatelyEqual(tMaxPhi,tMaxTheta,1e-12,1e-8) && ... - approximatelyEqual(tMaxPhi,tMaxR,1e-12,1e-8) && ... - strictlyLess(t,tMaxR,1e-12,1e-8) && ... - strictlyLess(tMaxR,t_end,1e-12,1e-8) && ~rStepViolation - % Triple boundary intersection - t = tMaxPhi; - current_voxel_ID_r = current_voxel_ID_r + tStepR; - 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); - if tStepR ~= 0; change_r = 1; else; change_r = 0; end - elseif approximatelyEqual(tMaxPhi, tMaxTheta, 1e-12,1e-8) && ... - strictlyLess(t,tMaxPhi,1e-12,1e-8) && ... - strictlyLess(tMaxPhi,t_end, 1e-12,1e-8) - % Phi, Theta equal - if strictlyLess(tMaxR, tMaxPhi, 1e-12,1e-8) && ... - strictlyLess(t,tMaxR,1e-12,1e-8) && ~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 approximatelyEqual(tMaxTheta,tMaxR,1e-12,1e-8) && ... - strictlyLess(t,tMaxR,1e-12,1e-8) && ... - strictlyLess(tMaxR, t_end,1e-12,1e-8) && ... - ~rStepViolation - % R, Theta equal - if strictlyLess(tMaxPhi,tMaxTheta,1e-12,1e-8) && ... - strictlyLess(t,tMaxPhi,1e-12,1e-8) - % 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 approximatelyEqual(tMaxR,tMaxPhi,1e-12,1e-8) && ... - strictlyLess(t,tMaxR,1e-12,1e-8) && ... - strictlyLess(tMaxR, t_end, 1e-12,1e-8) && ~rStepViolation - % R, Phi equal - if strictlyLess(tMaxTheta,tMaxR,1e-12,1e-8) && ... - strictlyLess(t,tMaxTheta, 1e-12,1e-8) - % 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 - % If the radial voxel changes, update the angular voxel boundary - % segments for intersection checks in angular hit. - if change_r == 1 - r = sphere_max_radius - delta_radius * (current_voxel_ID_r - 1); - a = [sphere_center(1) sphere_center(2)] - P_ang; - l1 = (a(:,1).^2 + a(:,2).^2).^-0.5; - c = [sphere_center(1) sphere_center(3)] - P_azi; - l2 = (c(:,1).^2 + c(:,2).^2).^-0.5; - P_ang(:,1) = sphere_center(1) - (r*l1) .* (sphere_center(1) - P_ang(:,1)); - P_ang(:,2) = sphere_center(2) - (r*l1) .* (sphere_center(2) - P_ang(:,2)); - P_azi(:,1) = sphere_center(1) - (r*l2) .* (sphere_center(1) - P_azi(:,1)); - P_azi(:,2) = sphere_center(3) - (r*l2) .* (sphere_center(3) - P_azi(:,2)); - 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 - -end