5
$\begingroup$

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?

Improve this question
$\endgroup$
4
  • $\begingroup$ RandomPoint[Cone[{{0, 0, 0}, {0, 0, h}}, r], 100] $\endgroup$
    – cvgmt
    Commented May 11, 2022 at 10:31
  • 2
    $\begingroup$ @cvgmt, I wrote the same comment, then realized that (s)he is not looking for random points, but for a grid! :-) $\endgroup$
    – Domen
    Commented May 11, 2022 at 10:33
  • $\begingroup$ Thanks for your comments, but yes: I do not want random points, but points with a specific spacing. $\endgroup$
    – rowsi
    Commented May 11, 2022 at 10:55
  • $\begingroup$ 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. $\endgroup$
    – MikeY
    Commented May 11, 2022 at 13:51

4 Answers 4

Reset to default
4
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you for your answer, this looks promising! Upvoted. $\endgroup$
    – rowsi
    Commented May 11, 2022 at 12:13
  • $\begingroup$ 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. $\endgroup$
    – rowsi
    Commented May 12, 2022 at 5:54
5
$\begingroup$

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}]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – rowsi
    Commented May 11, 2022 at 12:41
  • $\begingroup$ While I realy like this, I ultimately accepted the answer by @I.M., as it is better suited for upscaling. Thanks again! $\endgroup$
    – rowsi
    Commented May 12, 2022 at 5:56
3
$\begingroup$

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]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks a lot for your answer. These do indeed work as intended, but are working slower on my machine than my attempt. Anyway: upvoted. $\endgroup$
    – rowsi
    Commented May 11, 2022 at 11:46
  • $\begingroup$ 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. $\endgroup$
    – rhermans
    Commented May 11, 2022 at 11:48
1
$\begingroup$

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 <= #[[1]] <= r/2 && -r/2 <= #[[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 *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – rowsi
    Commented May 13, 2022 at 6:38
  • $\begingroup$ @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. $\endgroup$
    – MikeY
    Commented May 13, 2022 at 14:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.