Is this mathematica: WignerSeitz3D[a_, b_, c_, plotpoints_, size_]:=In[9]:=Module[{points, R, pointplot, nearorigin, contourplot},(*define the lattice vectors*)R[u_, v_, w_]:=u a+v b+w c;points=Flatten[Table[R[u, v, w],{u,-2, 2},{v,-2, 2},{w,-2, 2}], 2];pointplot=ListPointPlot3D[points,PlotStylePointSize[0.02], AxesFalse, AspectRatioAutomatic];(*define the function that returns one if the point isnearer to the origin than any other*)nearorigin[p1_?NumericQ, p2_?NumericQ, p3_?NumericQ]:=If[Nearest[points,{p1, p2, p3}][[1]] {0, 0, 0}, 1, 0];(*now plot the contour for the value nearorgin1/2*)(*contourplot=ContourPlot3D[nearorigin[{x,y,z}]1/2,{x,-2,2},{y,-2,2},{z,-2,2},MeshFalse];*)RegionPlot3D[nearorigin[x, y, z] 1,{x,-size, size},{y,-size, size},{z,-size, size}, MeshFalse, PlotPointsplotpoints,AxesFalse, BoxedFalse, BoxRatiosAutomatic