How to quickly count the number of neighboring voxels?

Question

I have got a 3D grid (voxels), where some of the voxels are filled, and some are not. The 3D grid is sparsely filled, so I have got a set filledVoxels with coordinates (x, y, z) of the filled voxels. What I am trying to do is find out is for each filled voxel, how many neighboring voxels are filled too.

Here is an example:

• filledVoxels contains the voxels (1, 1, 1), (1, 2, 1), and (1, 3, 1).
• Therefore, the neighbor counts are:
  • (1,1,1) has 1 neighbor
  • (1,2,1) has 2 neighbors
  • (1,3,1) has 1 neighbor.

Right now I have this algorithm:

voxelCount = new Map<Voxel, Integer>();

for (voxel v in filledVoxels)
  count = checkAllNeighbors(v, filledVoxels);
  voxelCount[v] = count;
end

checkAllNeighbors() looks up all 26 surrounding voxels. So in total I am doing 26*filledVoxels.size() lookups, which is quite slow.

Is there any way to cut down the number of required lookups? When you look at the above example you can see that I am checking the same voxels several times, so it might be possible to get rid of lookups with some clever caching.

If this helps in any way, the voxels represent a voxelized 3D surface (but there might be holes in it). I usually want to get a list of all voxels that have 5 or 6 neighbors.

Answer 1

If you're marching along the voxels one at a time, you can keep a lookup table corresponding to the grid, so that after you've checked it once using IsFullVoxel() you put the value in this grid. For each voxel you're marching in you can check if its lookup table value is valid, and only call IsFullVoxel() it it isn't.

OTOH it seems like you can't avoid iterating over all neighboring voxels, either using IsFullVoxel() or the LUT. If you had some more a priori information it could help. For instance, if you knew that there were at most x neighboring filled voxels, or you knew that there were at most y neighboring filled voxels in each direction. For instance, if you know you're looking for voxels with 5 to 6 neighbors, you can stop after you've found 7 full neighbors or 22 empty neighbors.

I'm assuming that a function IsFullVoxel() exists that returns true if a voxel is full.

Answer 2

Um, your question is not very clear. I'm assuming you just have a list of the filled points. In that case, this is going to be very slow, because you have to iterate through it (or use some kind of tree structure such as a kd-tree, but this will still be O(log n)).

If you can (i.e. the grid is not too big), just make a 3d array of bools. 26 lookups in a 3d array shouldn't really take that long (and there really is no way to cut down on the number of lookups).

Actually, now that I think of it, you could make it a 3d array of longs (64 bits). Each 64 bit block would hold 64 (4 x 4 x 4) voxels. When you are checking the neighbors of a voxel in the middle of the block, you could do a single 64 bit read (which would be much faster).

Answer 3

You can transform your voxel space into a octree in which every node contains a flag that specifies whether it contains filled voxels at all.

When a node does not contain filled voxels, you don't need to check any of its descendants.

Answer 4

If most of the moves in your iteration were to neighbors, you could reduce your checking by around 25% by not looking back at the ones you just checked before you made the step.

Answer 5

I'd say if each of your lookups is slow (O(size)), you should optimize it by binary search in an ordered list (O(log(size))).

The constant 26, I wouldn't worry much. If you iterate smarter, you could cache something and have 26 -> 10 or something, I think, but unless you have profiled the whole application and found out decisively that it is the bottleneck I would concentrate on something else.

Answer 6

As ilya states, there's not much you can do to get around the 26 neighbor look-ups. You have to make your biggest gains in efficiently identifying whether a given neighbor is filled or not. Given that the brute force solution is essentially O(N^2), you have a lot of possible ground to gain in that area. Since you have to iterate over all filled voxels at least once, I would take an approach similar to the following:

voxelCount = new Map<Voxel, Integer>();
visitedVoxels = new EfficientSpatialDataType();

for (voxel v in filledVoxels)
  for (voxel n in neighbors(v))
    if (visitedVoxels.contains(n))
      voxelCount[v]++;
      voxelCount[n]++;
    end
  next
  visitedVoxels.add(v);
next

For your efficient spatial data type, a kd-tree, as Zifre suggested, might be a good idea. In any case, you're going to want to reduce your search space by binning visited voxels.

Answer 7

You may find a Z-order curve a useful concept here. It lets you (with certain provisos) keep a sliding window of data around the point you're currently querying, so that when you move to the next point, you don't have to throw away many of the queries you've already performed.