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