How to generate random points in a region?

Question

The Mathematica 10 documentation was updated for FindInstance adding support for regions.

In my use case, I'm trying to sample points in a set of disks:

region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}]
FindInstance[{x, y} ∈ region, {x, y}, Reals, 2] // N

However the above code fails and generates the following error:

FindInstance::elemc: "Unable to resolve the domain or region membership condition {x,y}∈"

What's going wrong here?

Answer 1

There are already good answers, but I'm going to improve the performance, generalize to any region in any dimensions and make the function more convenient. The main idea is to use DirichletDistribution (the uniform distribution on a simplex, e.g. triangle or tetrahedron). This idea was implemented by PlatoManiac and me in the related question obtaining random element of a set given by multiple inequalities (there is also Metropolis algorithm, but it is not suitable here).

The code is relatively short:

RegionDistribution /: 
 Random`DistributionVector[RegionDistribution[reg_MeshRegion], n_Integer, prec_?Positive] :=
  Module[{d = RegionDimension@reg, cells, measures, s, m},
   cells = Developer`ToPackedArray@MeshPrimitives[reg, d][[All, 1]];
   s = RandomVariate[DirichletDistribution@ConstantArray[1, d + 1], n];
   measures = PropertyValue[{reg, d}, MeshCellMeasure];
   m = RandomVariate[#, n] &@EmpiricalDistribution[measures -> Range@Length@cells];
   #[[All, 1]] (1 - Total[s, {2}]) + Total[#[[All, 2 ;;]] s, {2}] &@
   cells[[m]]]

Examples

Random disks (2D in 2D)

SeedRandom[0];
region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}];
pts = RandomVariate[RegionDistribution[region], 10000]; // AbsoluteTiming
ListPlot[pts, AspectRatio -> Automatic]

{0.004473, Null}

Precise test

pts = RandomVariate[RegionDistribution[region], 200000000]; // AbsoluteTiming

{85.835022, Null}

Histogram3D[pts, 50, "PDF", BoxRatios -> {Automatic, Automatic, 1.5}]

It is fast for $2\cdot10^8$ points and the distribution is really flat!

Intervals (1D in 1D)

region = DiscretizeRegion[Interval[{0, 1}, {2, 4}]];
pts = RandomVariate[RegionDistribution[region], 100000]; // AbsoluteTiming
Histogram[Flatten@pts]

{0.062430, Null}

Random circles (1D in 2D)

region = DiscretizeRegion@RegionUnion[Circle /@ RandomReal[10, {100, 2}]];
pts = RandomVariate[RegionDistribution[region], 10000]; // AbsoluteTiming
ListPlot[pts, AspectRatio -> Automatic]

{0.006216, Null}

Balls (3D in 3D)

region = DiscretizeRegion@RegionUnion[Ball[{0, 0, 0}], Ball[{1.5, 0, 0}], Ball[{3, 0, 0}]];
pts = RandomVariate[RegionDistribution[region], 10000]; // AbsoluteTiming
ListPointPlot3D[pts, BoxRatios -> Automatic]

{0.082202, Null}

Surface cow disctribution (2D in 3D)

region = DiscretizeGraphics@ExampleData[{"Geometry3D", "Cow"}];
pts = RandomVariate[RegionDistribution[region], 2000]; // AbsoluteTiming
ListPointPlot3D[pts, BoxRatios -> Automatic]

{0.026357, Null}

Line in space (1D in 3D)

region = DiscretizeGraphics@ParametricPlot3D[{Sin[2 t], Cos[3 t], Cos[5 t]}, {t, 0, 2 π}];
pts = RandomVariate[RegionDistribution[region], 1000]; // AbsoluteTiming
ListPointPlot3D[pts, BoxRatios -> Automatic]

{0.005056, Null}

Answer 2

It would be nice if UniformDistribution worked on arbitrary regions, then we could simply do RandomVariate[UniformDistribution[region]]. Someone at Wolfram should get on that.

In the meantime, it seems we have to write our own sampling routines. @m_goldberg's answer is very nice (vote it up!) and uses rejection sampling, which works for arbitrary regions. However it will become slow if the measure of the region is small compared to that of its bounding box, for example if the disks are small and far apart. On the other hand, since you have a MeshRegion, it is actually possible to generate uniformly distributed samples without any rejection. We will pick a mesh cell randomly with probability proportional to its measure, then pick a point uniformly within that cell.

Sorry, I haven't bothered to wrap the code up into a neat module.

SeedRandom[0];
region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}]

d = RegionDimension[region];
coords = MeshCoordinates[region];
cells = MeshCells[region, d];
cellMeasures = PropertyValue[{region, d}, MeshCellMeasure];

randomCell[] := RandomChoice[cellMeasures -> cells]
randomBarycentric[] := Differences@Flatten@{0, Sort@RandomReal[1, d], 1}
randomRegionPoint[] := randomBarycentric[].coords[[First@randomCell[]]]

ListPlot[Table[randomRegionPoint[], {10000}], AspectRatio -> Automatic]

Answer 3

Good news! Version 10.2 of Mathematica has this built-in with the function RandomPoint[]. From the documentation:

• RandomPoint can generate random points for any RegionQ region that is also ConstantRegionQ.
• RandomPoint will generate points uniformly in the region reg.

The first example given is a simple disk, but there are a whole host of neat applications given in the documentation!

pts = RandomPoint[Disk[], 5000];
Graphics[{PointSize[Tiny], Point[pts]}]

Answer 4

Here is a work-around that will let you extract as many randoms points from a region as you wish.

SeedRandom[0]; 
region = 
  DiscretizeRegion @ RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}]

randomFromRegion[region_, trials_] :=
  Module[{bounds, randPts},
    bounds = RegionBounds@region;
    randPts = 
      Transpose @ {RandomReal[bounds[[1]], trials], RandomReal[bounds[[2]], trials]};
    Select[randPts, RegionMember[region]]]

Let's see what we get from 5000 trials.

pts = randomFromRegion[region, 5000];
Length @ pts

2550

We extracted 2550 points that lie within the region. Here is how they are distributed.

Graphics @ Point @ pts

Answer 5

In version 10.1 the undocumented function Random`RandomPointVector is useful:

region = DiscretizeRegion@RegionUnion@
 Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}];

Graphics@Point@
  Random`RandomPointVector[region, 1000, Automatic, Automatic]

The two Automatic arguments appear to be working precision and a method option - other allowed values are "Mesh" and "Rejection"

Answer 6

SeedRandom[0];

region = RegionUnion @ Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}];

DiscretizeRegion @ region

points = RandomReal[{-1, 5}, {10000, 2}];

circles = List @@ (region /. Disk -> Circle);

Graphics[{AbsolutePointSize[2],
  Transpose[{RegionMember[region, points] /.
     {False -> White, True -> Gray}, Point /@ points}], circles}, Frame -> True]

Answer 7

FindInstance have to apply to a expression like an equation or an inequality. DiscretizeRegion does not make an inequality. You should apply like the following example.

region = RegionUnion@Table[Disk[RandomReal[4, 2], RandomReal[1]],
    {10}
    ];
DiscretizeRegion[region]

RegionQ[region]

True

FindInstance[RegionMember[region, {x, y}], {x, y}, Reals, 2] // N

{{x -> 2.4597, y -> 4.09129}, {x -> 1.31552, y -> 2.00384}}

Answer 8

you can look at this question here

This method also may work fine:

p1 = {x, y} /. FindInstance[x^2 + y^2 < 1, {x, y}, Reals, 1000];
p2 = {y, x} /. FindInstance[x^2 + y^2 < 1, {x, y}, Reals, 1000];
point = RandomChoice[Join[p1, p2], 1000];
Graphics[{Circle[], {Red, Point[point]}}]