# Calculate list of points inside a 3D region with specific spacing

I'm trying to get a list of points inside a 3D region. However, I want these points to have a certain spacing in x-, y- and z-directions.

As of now I generate a list of points with the appropriate spacing and use Select with an appropriate RegionMemberFunction to filter for the desired points. This generates quite a bit of useless data, as millions of points get discarded after the use of Select.

An example:

h = 100;
r = 100;
rm = RegionMember[Cone[{{0, 0, 0}, {0, 0, h}}, r]];
pts = ParallelTable[{x, y, z}, {x, -r/2, r/2, r/100}, {y, -r/2, r/2,
r/100}, {z, 0, h, h/100}];
selpts = Select[rm]@Flatten[pts, 2]


Around 1 million points are generated by pts, but around 50% are discarded for the desired points selpts. Is there a way to generate only points inside the desired 3D region without the overhead?

• RandomPoint[Cone[{{0, 0, 0}, {0, 0, h}}, r], 100] May 11 at 10:31
• @cvgmt, I wrote the same comment, then realized that (s)he is not looking for random points, but for a grid! :-) May 11 at 10:33
• Thanks for your comments, but yes: I do not want random points, but points with a specific spacing. May 11 at 10:55
• What you'd like to do is to create a set of points in a cylinder (or cube, etc.), with the points distributed in such a way that when you transform from the cylinder to the cone, the points are evenly spaced. Haven't thought through the details of it yet. May 11 at 13:51

You can get an explicit constraint:

Simplify[RegionMember[Cone[{{0, 0, 0}, {0, 0, h}}, r], {x, y, z}], Assumptions-> {r > 0, h >0, 0<=z/h <= 1, Element[_, Reals]}]


Then:

h = 10;
r = 10;
ps = Flatten[Table[If[x^2 + y^2 <= r^2*(h -z)^2/h^2,{x, y, z}, Nothing], {x, -r/2, r/2, r/10}, {y, -r/2, r/2,r/10}, {z, 0, h, h/10}] , 2];
ListPointPlot3D[ps,PlotRange->All]


Perhabs it is possible to move constraint to ranges.

• Thank you for your answer, this looks promising! Upvoted. May 11 at 12:13
• Accepted this answer, as it is best suited for scaling up to several millions of points with ParallelTable and the list structure of the result works in my favor for further computation. May 12 at 5:54

Pick seems faster.

h = 100;
r = 100;
reg = Cone[{{0, 0, 0}, {0, 0, h}}, r];
rm = RegionMember[reg];
n = 30;
pts = Tuples[Subdivide[##, n] & @@@ RegionBounds[reg]];
selpts = Pick[pts, rm[pts]];
Graphics3D[{Blue, Point@selpts, Opacity[.2], reg}] • Thank you very much for your answer. Indeed, Pick works faster here than Select, but fundamentally it does the same as my code: generating a bunch of points and discarding the unwanted ones afterwards. Anyway: upvoted, as it does work as intended. May 11 at 12:41
• While I realy like this, I ultimately accepted the answer by @I.M., as it is better suited for upscaling. Thanks again! May 12 at 5:56

This works but is not fast. I really hope to see a variety of answers to this question, and hopefully, somebody will benchmark them all.

pnts = Block[{h=20,r=10,rmf},
rmf = Reduce@RegionMember[Cone[{{0, 0, 0}, {0, 0, h}}, r]][{x,y,z}];
ReplaceAll[
{x,y,z},
Solve[rmf,{x,y,z}∈Integers]
]
]


Or a bit faster for version 12.3 and after

pnts = Block[{h=20,r=10,rmf},
rmf = Reduce@ RegionMember[Cone[{{0, 0, 0}, {0, 0, h}}, r]][{x,y,z}];
NSolveValues[rmf,{x,y,z}∈Integers]
];


Display

ListPointPlot3D[pnts,PlotRange->All] • Thanks a lot for your answer. These do indeed work as intended, but are working slower on my machine than my attempt. Anyway: upvoted. May 11 at 11:46
• Thanks for the vote, I do realize this solution is slow. Probably a cubic grid Mesh could be faster, but that's beyond my knowledge at the moment. I honestly hope you will get a variety of answers. May 11 at 11:48

Not an answer, just an analysis. Some details left out, just considering an alternate approach.

If you can accept a pretty close approximation, you can kind of use RandomPoint[ ]. As an example, I generate 10,000,000 points in the cone, taking just the integer part of the float number to get a sampling of the grid. This can be fast. I wonder about an additional transformation to make it almost surely cover the grid with fewer samples.

The slow part is then chopping off the edges of the cone to get the rectangular base.

r = h = 100;
cone = Cone[{{0, 0, 0}, {0, 0, h}}, r];
prePts = Union@(IntegerPart /@ RandomPoint[cone, 10000000]);) // Timing
(*  {3.26563, Null} *)

Length@prePts
(* 1045777 *)

ptsCube = Select[prePts, (-r/2 <= #[] <= r/2 && -r/2 <= #[] <= r/2 ) &]; // Timing
(*  {3.92188, Null} *)

Length@ptsCube
(* 625310 *)


The length of the true grid is

Length@selpts
(* 625310 *)


Also tried this approach, directly forming a region that is an intersection of the spaces. The sampling algorithm obviously has to wrestle with the odd space shape.

cube = Cuboid[{-r/2, -r/2, 0}, {r/2, r/2, h}];
ri = RegionIntersection[cube, cone];
coneCubePts = Union@(IntegerPart /@ RandomPoint[ri, 2000000]); // Timing
(* {9.79688, Null} *)

Length@coneCubePts
(* 582744 *)

• Thank you very much for your answer. Indeed, its a clever approach to use RandomPoint for fast generation of points and only take the integer parts of the generated points. However, this wouldn't work for my application, when I want a spacing smaller than integers. But to develop this further: Another idea is to use Round on the list of random points. This works around 30% faster on my machine than IntegerPart, altough some points will be outside the defined region. Ultimately, a threedimensional matrix works better for my further applications than a one-dimensional list of points. May 13 at 6:38
• @rowsi, you could come up with an approach for variable spacing. It was convenient that your problem was nicely spaced on integers. I played with Round[ ], but as you noted, it requires another constraint check, which adds a bunch of time to it. Ultimately, just an interesting experiment. May 13 at 14:28