Skip to content
Closed
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 44 additions & 1 deletion matlab/spherical-traversal/sphericalCoordinateTraversal.m
Original file line number Diff line number Diff line change
Expand Up @@ -194,6 +194,7 @@
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
Expand Down Expand Up @@ -251,7 +252,16 @@
else
t_start = t_begin;
end


% The chord length test. We first find the exit time and entrance
% time for the ray. We also need to check that if the ray
% begins/end within the grid, not the outside. Acc is used
% to keep track of total time spent on traversal
if verbose
fprintf('\nExpected chord length is: %d', min(t2,t_end)-max(t1,t_begin))
acc = 0.0
end

% TRAVERSAL PHASE
t = t_start;
t_end = min(t_grid, t_end);
Expand Down Expand Up @@ -281,6 +291,9 @@
% 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
if verbose
acc = acc + tMaxTheta - t
end
t = tMaxTheta;
current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections);
change_r = 0;
Expand All @@ -290,6 +303,9 @@
strictlyLess(tMaxR,t_end, 1e-12,1e-8) ...
&& ~rStepViolation
% Case tMaxR is minimum
if verbose
acc = acc + tMaxR - t
end
t = tMaxR;
current_voxel_ID_r = current_voxel_ID_r + tStepR;
if tStepR ~= 0; change_r = 1; else; change_r = 0; end
Expand All @@ -298,6 +314,9 @@
strictlyLess(t,tMaxPhi,1e-12,1e-8) && ...
strictlyLess(tMaxPhi,t_end,1e-12,1e-8)
% Case tMaxPhi is minimum
if verbose
acc = acc + tMaxPhi - t
end
t = tMaxPhi;
current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections);
change_r = 0;
Expand All @@ -306,6 +325,9 @@
strictlyLess(t,tMaxR,1e-12,1e-8) && ...
strictlyLess(tMaxR,t_end,1e-12,1e-8) && ~rStepViolation
% Triple boundary intersection
if verbose
acc = acc + tMaxPhi - t
end
t = tMaxPhi;
current_voxel_ID_r = current_voxel_ID_r + tStepR;
current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections);
Expand All @@ -318,10 +340,16 @@
if strictlyLess(tMaxR, tMaxPhi, 1e-12,1e-8) && ...
strictlyLess(t,tMaxR,1e-12,1e-8) && ~rStepViolation
% R min
if verbose
acc = acc + tMaxR - t
end
t = tMaxR;
current_voxel_ID_r = current_voxel_ID_r + tStepR;
if tStepR ~= 0; change_r = 1; else; change_r = 0; end
else
if verbose
acc = acc + tMaxPhi - t
end
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);
Expand All @@ -335,10 +363,16 @@
if strictlyLess(tMaxPhi,tMaxTheta,1e-12,1e-8) && ...
strictlyLess(t,tMaxPhi,1e-12,1e-8)
% Phi min
if verbose
acc = acc + tMaxPhi - t
end
t = tMaxPhi;
current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections);
change_r = 0;
else
if verbose
acc = acc + tMaxTheta - t
end
t = tMaxTheta;
current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections);
current_voxel_ID_r = current_voxel_ID_r + tStepR;
Expand All @@ -351,10 +385,16 @@
if strictlyLess(tMaxTheta,tMaxR,1e-12,1e-8) && ...
strictlyLess(t,tMaxTheta, 1e-12,1e-8)
% Theta min
if verbose
acc = acc + tMaxTheta - t
end
t = tMaxTheta;
current_voxel_ID_theta = mod(current_voxel_ID_theta + tStepTheta, num_angular_sections);
change_r = 0;
else
if verbose
acc = acc + tMaxR - t
end
t = tMaxR;
current_voxel_ID_phi = mod(current_voxel_ID_phi + tStepPhi, num_azimuthal_sections);
current_voxel_ID_r = current_voxel_ID_r + tStepR;
Expand All @@ -363,6 +403,9 @@
else
return;
end
if verbose:
fprintf('\nCurrently, we have spent %d time on traversal\n',acc)
end
% If the radial voxel changes, update the angular voxel boundary
% segments for intersection checks in angular hit.
if change_r == 1
Expand Down