From c8515e5006ae329b914ae01eacb4290ad9f42545 Mon Sep 17 00:00:00 2001 From: Youhan Date: Wed, 15 Apr 2020 12:14:23 -0400 Subject: [PATCH 1/2] Add Chord Length test --- .../sphericalCoordinateTraversal.m | 42 ++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/matlab/spherical-traversal/sphericalCoordinateTraversal.m b/matlab/spherical-traversal/sphericalCoordinateTraversal.m index 21abe9c..0efae66 100644 --- a/matlab/spherical-traversal/sphericalCoordinateTraversal.m +++ b/matlab/spherical-traversal/sphericalCoordinateTraversal.m @@ -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 @@ -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); @@ -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; @@ -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 @@ -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; @@ -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); @@ -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); @@ -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; @@ -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; From 2718a633cb9478db0c19eb84000fa9e163df57d6 Mon Sep 17 00:00:00 2001 From: Youhan Date: Wed, 15 Apr 2020 12:19:59 -0400 Subject: [PATCH 2/2] Add display statement for acc --- matlab/spherical-traversal/sphericalCoordinateTraversal.m | 3 +++ 1 file changed, 3 insertions(+) diff --git a/matlab/spherical-traversal/sphericalCoordinateTraversal.m b/matlab/spherical-traversal/sphericalCoordinateTraversal.m index 0efae66..16fad64 100644 --- a/matlab/spherical-traversal/sphericalCoordinateTraversal.m +++ b/matlab/spherical-traversal/sphericalCoordinateTraversal.m @@ -403,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