Raycast sphere sampling mockup

output.txt

scanradius = 50000.0
set point distance = 256.0
num_pts = 610352
scantime = 2441798s = 40696.6m = 678.3h
too many points to run k nearest neighbors

--

scanradius = 25000.0
set point distance = 256.0
num_pts = 152588
scantime = 305224s = 5087.1m = 84.8h
too many points to run k nearest neighbors

---

scanradius = 5000.0
set point distance = 256.0
num_pts = 6104
scantime = 2441s = 40.7m = 0.7h
average point distance = 226.89593558252824

raycast.py

from numpy import pi, cos, sin, arccos, arange, round, column_stack, sort, mean
import mpl_toolkits.mplot3d
import matplotlib.pyplot as pp

scanradius = 50000.0
pointdistance = 256.0
print("scanradius =", scanradius)
print("set point distance =", pointdistance)

num_pts = int(round((4.0 / ((pointdistance/2)**2)) * scanradius**2))
scantime = (num_pts * scanradius) / (2083 * 6)
print("num_pts =", num_pts)
print("scantime = %ds = %.1fm = %.1fh" % (scantime, scantime / 60.0, scantime / 3600.0))

indices = arange(0, num_pts, dtype=float) + 0.5

phi = arccos(1 - 2*indices/num_pts)
theta = pi * (1 + 5**0.5) * indices

x = scanradius * cos(theta) * sin(phi)
y = scanradius * sin(theta) * sin(phi)
z = scanradius * cos(phi);

##pp.figure().add_subplot(111, projection='3d').scatter(x, y, z);
##pp.show()

from scipy.spatial import distance
if num_pts < 10000:
    points = column_stack((x,y,z))
    all_distances = distance.squareform(distance.pdist(points))
    closest_distances = sort(all_distances, axis=1)[:, 1:4]
    average_distance = mean(mean(closest_distances, axis=1))
    print("average point distance =", average_distance)
else:
    print("too many points to run k nearest neighbors")