# Select the surface points of a body on a cubic lattice

I have a conceptually simple problem. Given an arbitrary shaped body on a cubic lattice, I have to select the points corresponding to the surface. An arbitrary point may have up to 6 neighbors and one just has to select the ones, which have less than 6 neighbors. I'm going to denote the lattice spacing with $$\delta L$$. My solution is

SurfaceQ[list_, volume_] := Block[{tmp, tlist},

tlist =
ConstantArray[list,
6] + {{0. + \[Delta]L, 0., 0.}, {0., 0. + \[Delta]L, 0.}, {0.,
0., 0. + \[Delta]L}, {0. - \[Delta]L, 0., 0.}, {0.,
0. - \[Delta]L, 0.}, {0., 0., 0. - \[Delta]L}};
tmp = Complement[tlist, volume, SameTest -> (#1 == #2 &)];
If[Length[tmp] > 0, list, {}]
]


One inputs the point in the variable 'list' and the data corresponding to the total body 'volume'. Then the function SurfaceQ computes the complement set of neighboring points of 'list' and 'volume'. If it is empty, i.e. all the neighboring points can be found in 'volume', then we detected a bulk point. If the complement set is not empty, then we detected a surface point. My main problem is, that Complement and related functions does not seem to do their jobs properly!

In order to check whether this function is correct I tested it in a simple example:

\[Delta]L = 0.01;
nPts = 2;
cubify[point_?ListQ] :=
Tuples[{Range[point[] - nPts*\[Delta]L,
point[] + nPts*\[Delta]L, \[Delta]L],
Range[point[] - nPts*\[Delta]L,
point[] + nPts*\[Delta]L, \[Delta]L],
Range[point[] - nPts*\[Delta]L,
point[] + nPts*\[Delta]L, \[Delta]L]}]

testdat = cubify[{0, 0, 0}];
result = DeleteCases[
SurfaceQ[testdat[[#]], testdat] & /@ Range[Length@testdat], {}];


Now, the problem is that this result is far from correct.

facePts[x_] := 2 x^2 + 2 x (x - 2) + 2 (x - 2) (x - 2)
ListPointPlot3D[{testdat, result}]
facePts[2 nPts + 1] == Length@result


facePts gives the number of surface points of a cube on a cubic lattice, which should obviously be equal to the number of points given by my function. The incorrect behavior can also be seen from the plot.

My question is: why does Complement not work in this scenario? (It only works for nPts = 1, i.e. a 3x3x3 cube). Or is there a better way to perform this task? I feel like there exists a 200 IQ one-liner solution... :D

• Does cubify do almost the same thing as CoordinateBoundingBoxArray? Oct 17 at 19:19
• @MichaelE2 Yes, I was not aware of CoordinateBoundingBoxArray until now. I really want to avoid using mesh, because I need the surface points with the precision of the lattice spacing, and in general, the 3D body on the lattice will not be a convex shape. A solution using a meshing function is OK so far as it recovers the surface points with an error $\pm\delta L /10$ and i do not know how to control the precision for meshes. Oct 17 at 19:31

Update: You may want to use MorphologicalPerimeter as EdgeDetect can sometimes lead to stray voxels. You can use the CornerNeighbors -> False option of MorphologicalPerimeter to specify you don't want to include the corners, i.e 6-connected, instead of 26-connected.

It would make a lot more sense to me to use Mathematica's EdgeDetect on an Image3D.

So suppose our cubic lattice is this:

lattice = CoordinateBoundsArray[{{-2, 2}, {-2, 2}, {-2, 2}}, 0.1];


Then we can define a solid ball of points in this cubic lattice. This is just for show - in practice you would just have some arbitrary rank-3 tensor full of 1s and 0s:

rmf = RegionMember[Ball[]];
ball = Image3D@Map[Boole@*rmf, lattice, {3}]


To get the surface elements you just use EdgeDetect like this:

surf = EdgeDetect@ball


... and you can see from the slices that the shape is now hollowed out:

Image3DSlices@surf


Now to get the lattice points back you can just find the positions of 1s and extract from the lattice:

ListPointPlot3D[
Extract[lattice, Position[ImageData@surf, 1]]
, BoxRatios -> 1]

• Nice solution, thank you! Do you have any idea, why the complement function does not work properly? If no other answers appear, I'll accept your answer in the evening. Oct 18 at 10:25
• @dzsoga Your Complement appears to work fine, it's your expectations of the number of points that's wrong. e.g with dL=1 volume = Flatten[ CoordinateBoundsArray[{{-2, 2}, {-2, 2}, {-2, 2}}, 1], 2]; Total[(SurfaceQ[#, volume] & /@ volume) /. {List[] -> 0, List[_, _, _] -> 1}]  returns 98 which is a 5x5x5 cube minus the 3x3x3 interior. You don't need the 2 nPts + 1 I think. But in any case, it's less efficient to use set operations like Complement, and far quicker to use optimized built-ins like I've done above. Oct 18 at 11:52
• @flinfty I have one last question: how do I make a Region/RegionFunction out of a list of numbers/a rank 3 tensor containing the coordinates of the elements of the body? I obtain a set of data( a closed volume) from a parameter space scan and the best I can tell is the max and min coordinates in each direction and of course the resolution (dL). The scan takes some seed values and checks a cubic volume around the seed for some conditions. The final paramer space is a closed volume, but it is not embedded in a nxnxn cube. Oct 19 at 11:19
• @dzsoga I'm having trouble understanding that. Is your input data a list of points i.e a list of 3d coordinates, or a rank-3 tensor, i.e a list-of-lists-of-lists of 1s and 0s? You don't need a region function if your data is already in a rank-3 tensor form to achieve the results in my answer. Oct 19 at 11:43
• @dzsoga if you need to take a list of points and voxelize them into a rank-3 tensor, then you can do this with BinCounts + Unitize. Oct 19 at 11:49