-1

I am trying to create a Voronoi diagram using a 3d array representing voxels where the coordinates are represented by the indices of the array. For example:

points = np.random.random([10,3])
vox_len = 0.1
lx = ly = lz = 11
hull_space=np.zeros([lx,ly,lz],dtype=np.int16)
hull_space.fill(-1)
for i in range(lx):
    for j in range(ly):
        for k in range(lz):
            coord = np.array([i,j,k]).astype(float)*vox_len
            diff = points - coord
            dist = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)
            closest = np.argmin(dist)
            hull_space[i][j][k]=closest

The problem is that with larger arrays the code is very slow. I've also tried:

grid = np.indices((lx,ly,lz)).astype(float)*vox_len
hull_space = np.ones([lx,ly,lz],dtype=np.int16)
hull_space.fill(-1)
closest_dist = np.ones(hull_space.shape)
closest_dist.fill(999)
for index,point in enumerate(points):
    dist2 = (grid[0]-point[0])**2 + (grid[1]-point[1])**2 + (grid[2]-point[2])**2
    hull_space[dist2 < closest_dist]=index
    closest_dist[dist2 < closest_dist]=dist2[dist2 < closest_dist]

But results were even worse. I am aware of scipy's voronoi diagram methods but I need to build a voxel image in order to do distance transforms of the edges of the voronoi regions. My goal is to create an image of an entangled fibrous structure to use in a pore network model and the voxel approach is useful for getting pore volumes.

Any help making the algorithm run faster would be appreciated. For a 300x300x300 diagram with 500 points it takes about 50 minutes on my machine

TomTranter
  • 3
  • 3
  • Hi, if you just embed your process with a njit decorator and change the types of the zeros and coord variables to numpy specific types you will see a noticeable speedup. – Oier Arcelus Jan 05 '22 at 09:06

1 Answers1

0

You can try to sort the points. Treat the co-ordinate as a binary and interleave it. Then sort it both ways. In fact Cgal uses the same algorithm.

Micromega
  • 12,486
  • 7
  • 35
  • 72
  • Thanks for the reply. Could you add a code snippet as an example of how to do this please. I'm not too familiar with manipulating binary data. – TomTranter Feb 03 '15 at 14:56
  • @TomTranter:Can you elaborate? Python isn't my language either also asking for code isn't the right place here?!? – Micromega Feb 03 '15 at 17:01
  • I don't really understand why or how your method works and the best way to see would have been to get an example but if that breaks some sort of protocol I'm sorry. If you could point me towards some theory to help me understand that would be great. – TomTranter Feb 04 '15 at 00:45
  • @TomTranter:The idea is to reduce the complexity,i.e to transform a n-dimension to a 1 dimension. In mathematics there are some very useful fractals, or patterns. One is I described above. You can also look for space filling curve, morton curve and hilbert curve. It works on any dimension and fills it completely. When you traverse the space in 1 dimension you can use it to sort nearby points. – Micromega Feb 04 '15 at 12:32
  • Thanks for the extra info. I found this article which has helped me understand why the method works for anyone reading this without much knowledge of computer science. http://www.forceflow.be/2013/10/07/morton-encodingdecoding-through-bit-interleaving-implementations/ – TomTranter Feb 04 '15 at 14:34
  • @TomTranter:If my answer is helpful please consider to accept it! – Micromega Feb 04 '15 at 15:35
  • I'm actually still finding it quite difficult to implement this suggestion as I need a method for sorting the points once in morton order. I have managed to speed up the process by invoking a subdivision process so instead of evaluating the closest neighbor point for every voxel I can evaluate it at the corners of a box containing a subset of voxels. If all corners are nearest the same point then all voxels within the box are nearest that point. Starting with a box covering the whole domain the process can be run recursively dividing the box into smaller ones. I will still try your method – TomTranter Feb 05 '15 at 16:35