I am now trying to calculate the poisson sphere distribution(a 3D version of the poisson disk) using python and then plug in the result to POV-RAY so that I can generate some random distributed packing rocks. I am following these two links:
0.Create an n-dimensional grid array and cell size = r/sqrt(n) where r is the minimum distance between each sphere. All arrays are set to be default -1 which stands for 'without point'
1.Create an initial sample. (it should be placed randomly but I choose to put it in the middle). Put it in the grid array. Also, intialize an active array. Put the initial sample in the active array.
2.While the active list is not empty, pick a random index. Generate points near it and make sure the points are not overlapping with nearby points(only test with the nearby arrays). If no sample can be created near the 'random index', kick the 'random index' out. Loop the process.
And here is my code:
import math
import numpy
from random import uniform
import random
from math import floor
r = 1
k = 30
grid = []
w = r / math.sqrt(2)
active = []
width = 100
height = 100
depth = 100
cols = floor(width / w)
rows = floor(height / w)
deps = floor(depth / w)
default = numpy.array((-1,-1,-1))
for i in range(cols * rows * deps):
grid.append(default)
x = width / 2
y = height / 2
z = depth / 2
i = floor(x / w)
j = floor(y / w)
k = floor(z / w)
pos = numpy.array((x,y,z))
grid[i + cols * (j + rows * k)] = pos
active.append(pos)
while (len(active) > 0) and (len(grid[grid == -1]) > 0):
randIndex = floor(uniform(0, len(active)))
pos = active[randIndex]
found = False
for n in range(k):
m1 = uniform(-2 * r, 2 * r)
m2 = uniform(-2 * r, 2 * r)
m3 = uniform(-2 * r, 2 * r)
m = numpy.array((m1,m2,m3))
sample = numpy.add(pos, m)
col = floor(sample[0] / w)
row = floor(sample[1] / w)
dep = floor(sample[2] / w)
if (col > -1 and row > -1 and dep > -1 and col < cols and row < rows and dep < deps and numpy.all([grid[col + cols * (row + rows * dep)],default])==True):
ok = True
for i in range(-1,2):
for j in range(-1, 2):
for k in range(-1, 2):
index = (col + i) + cols * ((row + j) + rows * (dep + k))
if col + i > -1 and row + j > -1 and dep + k > -1 and col + i < cols and row + j < rows and dep + k < deps:
neighbor = grid[index]
if numpy.all([neighbor, default]) == False:
d = numpy.linalg.norm(sample - neighbor)
if (d < r):
ok = False
if ok == True:
found = True
grid[col + cols * (row + rows * dep)] = sample
active.append(sample)
if found == False:
del active[randIndex]
print(len(active))
for printout in range(len(grid)):
print("<" + str(active[printout][0]) + "," + str(active[printout][1]) + "," + str(active[printout][2]) + ">")
print(len(grid))
My code seems to run forever and do not obey my condition(distance of two spheres must be larger than 2 * radius) as shown in the visualization by POV-RAY.(picture in comment)
Therefore I tried to add a print(len(active)) in the last of the while loop.
Surprisingly, I think I discovered the bug as the length of the active list just keep increasing! (It is supposed to be the same length as the grid) I think the problem is caused by the active.append(), but I can't figure out where is the problem as the code is literally the 90% the same as the one made by Mr.Shiffman.
I don't want to free ride this but I have already checked again and again while correcting again and again for this code :(. Still, I don't know where the bug is. (why do the active[] keep appending!?)
Thank you for the precious time.
