# How to randomly insert a region inside another region like RandomPoint?

With RandomPoint one can find a random point inside of a region. I would like to instead obtain a region randomly placed inside another region. One possibility would be to use:

(Note …=\[Ellipsis])

GeometricTransformation[smaller…region,
TranslationTransform[RandomPoint[bigger…region]]]


And to keep trying until RegionWithin

Is there a better way ?

Possible Test regions:

The choices below are optional you can use whatever regions you want while making sure the method is not restricted to that particular case

Annulus[]


Edit:

Random bubbles :

RegionUnion[Disk[RandomReal[1, {10, 2}], RandomReal[]] // Thread]


The example below is taken from the neat example section of CylindricalDecomposition.

disk[{x0_, y0_}, r_] := (x - x0)^2 + (y - y0)^2 < r^2;

smiley =
ImplicitRegion[
disk[{-7/4, 9/4}, 1/3] ||
disk[{9/4, 9/4},
1/3] || (disk[{0, 0}, 5] && ! disk[{2, 2}, 1] && !
disk[{-2, 2}, 1] && (! disk[{0, -1}, Sqrt[5]] ||
disk[{0, 1}, 3])), {x, y}];

scaled…smiley =
GeometricTransformation[smiley,
RescalingTransform[
smiley // BoundingRegion // Apply[List] //
Transpose, {{0, 0.1}, {0, 0.1}}]];

cat = (ImageCases@
Entity["Species", "Species:FelisCatus"]["Image"])[[1, 1]]

cat…outline =
ColorDetect[cat, RGBColor[
0.5177182138426386, 0.5201882005011664, 0.2585894129044725]] //
ColorNegate // SkeletonTransform // InverseDistanceTransform //
MorphologicalTransform[#, {Min, Max}] & //
SelectComponents[#, Large] & // ImageMesh

scaled…cat =
GeometricTransformation[cat…outline,
RescalingTransform[
cat…outline // BoundingRegion // Apply[List] //
Transpose, {{0, 0.1}, {0, 0.1}}]];


Examples:

shifted =
GeometricTransformation[scaled…smiley,
TranslationTransform[RandomPoint[Annulus[]]]]

Show[BoundaryDiscretizeRegion@shifted,
Graphics[{Magenta, Opacity[0.2], Annulus[]}]]


In the image above the smiley does not fit entirely within the annulus

But in the image below:

The cat fits nicely inside the annulus

EDIT :

The issue I have with a direct brute force algorithm that keeps trying until RegionWithin is that for complex geometries the loop might try lots of times in the same region where it can't fit. For example a medium size circle trying to fit in the eye of the smiley. Ideally I am hoping for some method that maybe exploits meshing in an optimal way that avoids checking the same area multiple times or guesses beforehand if the region would fit. Maybe some kind of iterative or heuristic method that exploits information about properties of the region like the size or the meshes not sure.

Try RegionWithin in a while-loop.

Examplary I try to find a disk which is within Annulus[]

reg = Annulus[]
p = RandomPoint[reg]
While[! RegionWithin[reg, disk = Disk[p, 1/8]], p = RandomPoint[reg]]

Graphics[{FaceForm[Gray], EdgeForm[Blue], reg,FaceForm[Darker[Green]], disk}]


• Hi. thanks, a while loop is kind of what I meant when I wrote "to keep trying until RegionWithin " but it's my fault for not being explicit enough. I made an edit to the question in an attempt to be more explicit. Commented Nov 19, 2022 at 11:57

Note 1 : I noticed that the code runs a lot faster if regions are discretized before using RegionWithin.

Note 2 : The more advanced code below uses MeshConnectivityGraph in version 13.1. You might need to replicate that function in the simpler case of a graph of mesh points if you wish to use the code below in previous versions.

Note 3 : The code below does not check whether it is actually possible to fit one region in another.

### Examples:

(functions defined in next section with code)

The direct approach

The example below shows the output of 8 calls to a function that implements the direct method suggested in the question "just keep trying":

AbsoluteTiming[
points…and…regions =
Table[random…region…position•direct[
bubbles, scaled…annulus], 8];]


{0.449325, Null}

Visualizing the regions :

Graphics[{Opacity[0.9], Magenta, BoundaryDiscretizeRegion@bubbles,
Opacity[0.8], RGBColor[0.148, 0.33, 0.54], PointSize[0.02],
Point@points…and…regions[[All, 1]], RGBColor[
0.386403, 0.743518, 0.934237],
points…and…regions[[All, 2]]}]


With 5 times more bubbles and 40 small region positions instead of 8:

timing: 10.55 s

Notice that regions may accumulate in the same zone of the bigger region and can also intersect. There are two codes in the code section that try to tackle that and examples are given below.

8 calls :

AbsoluteTiming[
points…and…regions =
randomRegionPosition[bubbles, scaled…annulus, 8];]


{0.2574, Null}

Visualization:

With 5 times more bubbles and 40 small region positions instead of 8:

timing: 1.18432 s (* instead of the 10 seconds from the direct approach *)

Favoring configurations with less overlaps

Note : the method in the next subsection can be significantly faster and depending on the case more reliable

We can try to maximize the distance between regions by taking 10 region positions, for example, and taking the position that is furthest away from previous positions. The function that implements that is : weakly…self…avoiding…regions. The syntax is :

weakly…self…avoiding…regions[generator,number…of…points]


where generator is a function that takes as input an integer representing the number of region positions to output and outputs the the region positions and shifted regions.

With the direct approach :

AbsoluteTiming[
points…and…regions =
weakly…self…avoiding…regions[
Table[random…region…position•direct[
bubbles, scaled…annulus], #] &, 8];]


{4.32627, Null}

Note: timing can be decreased by decreasing the default integer valued avoidance parameter (representing the sample size over which distance is maximized)

Visualization:

Graphics[{Opacity[0.9], Magenta, BoundaryDiscretizeRegion@bubbles,
Opacity[0.8], RGBColor[0.148, 0.33, 0.54], PointSize[0.02],
Point@points…and…regions[[All, 1]], RGBColor[
0.386403, 0.743518, 0.934237],
points…and…regions[[All, 2]]}]


With 5 times more bubbles and 40 small region positions instead of 8:

timing: 1 min 45 s

Going directly to the case of 5 times more bubbles and 40 small region positions instead of 8 :

AbsoluteTiming[
points…and…regions =
weakly…self…avoiding…regions[
randomRegionPosition[bubbles…more,
scaled…annulus, #] &, 40];]


{35.0745, Null}

Notice the code is roughly 3 times faster in this case. Results may differ depending on the geometry of the regions.

Imposing no overlaps

The method below imposes no overlaps between regions. The syntax is :

strongly…self…avoiding…regions[
bigger…region, smaller…region,
generator, number…of…points]


To improve efficiency, all regions should be discretized. The output regions are also discretized. generator expects as input a bigger…region and outputs the smaller region with position shifted to a new location in the bigger region (as usual). The positions of the regions are not outputted just the regions themselves.

As an example:

AbsoluteTiming[
points…and…regions =
strongly…self…avoiding…regions[
DiscretizeRegion@bubbles,
DiscretizeRegion@scaled…annulus,
random…region…position•direct[#,
DiscretizeRegion@scaled…annulus][[2]] &, 8];]


{0.27135, Null}

5 times more bubbles and 40 small region positions instead of 8 :

timing : 3.4 s

The timing is significantly faster and the regions avoid one another

### Code

(discussion and explanation of some of the codes in next sections)

The direct method for finding random points (overlaps can happen)

Clear[random…region…position•direct];
random…region…position•direct[
big…region_, small…region_,
max…iterations_ : 50] :=

Module[{small…region•centered,
small…region…center,
small…region•shifted, counter,
discretized…region},

small…region•centered =
TransformedRegion[small…region,
TranslationTransform[-RegionCentroid@small…region]
];

small…region…center =
RandomPoint[big…region];

small…region•shifted =
TransformedRegion[small…region•centered,
TranslationTransform[
small…region…center]];
counter = 0;

discretized…region = DiscretizeRegion@big…region;

While[Not@RegionWithin[discretized…region,
DiscretizeRegion@small…region•shifted],

If[counter > max…iterations,

Print["Too many attempts. \n
The region might not be able to fit. \n
Consider also increasing option max iterations. \n
Default value for max iterations is 50. "]; Abort[]
];

small…region…center =
RandomPoint[big…region];

small…region•shifted =
TransformedRegion[small…region•centered,

TranslationTransform[ small…region…center]];

counter++;

];
{small…region…center,
small…region•shifted}

]


Options[randomRegionPosition]=
{"distribution"->(1/(1+#^2) &),"max iterations"-> 50};

Clear[randomRegionPosition];
randomRegionPosition[big…region_, small…region_,
how…many_, OptionsPattern[]]:=

Module[{centers,graph…components,weight…table,
region•discretized•list,mesh…graphs,small…region•centered,
small…region…center,small…region•shifted,
counter,random…component,discretized…region},

region•discretized•list= big…region //BoundaryDiscretizeRegion //ConnectedMeshComponents//Map[DiscretizeRegion] ;

mesh…graphs= region•discretized•list // Map[MeshConnectivityGraph];

centers=mesh…graphs //Map[GraphCenter] //Map[First];

weight…table=
Table[
Association@Table[
v-> OptionValue["distribution"]@GraphDistance[
mesh…graphs[[component]],
centers[[component]],
v],

{v,VertexList[mesh…graphs[[component]]]}
],

{component,Length@mesh…graphs}
]
;

small…region•centered=
TransformedRegion[small…region,
TranslationTransform[-RegionCentroid@small…region]
]
;
discretized…region=DiscretizeRegion@big…region;

Table[

random…component=RandomInteger[{1,Length@mesh…graphs}];

small…region…center=
MeshCoordinates[region•discretized•list[[random…component]]][[
Last@RandomChoice[
(Values@weight…table[[random…component]])->
(Keys@weight…table[[random…component]])
]
]]
;

small…region•shifted=
TransformedRegion[small…region•centered,
TranslationTransform[small…region…center]];

counter=0;

While[Not@RegionWithin[discretized…region,
DiscretizeRegion@small…region•shifted],

If[counter> OptionValue["max iterations"],

Print["Too many attempts. \n
The region might not be able to fit. \n
Consider also increasing option max iterations. \n
Default value for max iterations is 50. "];
Abort[];
];

random…component=RandomInteger[{1,Length@mesh…graphs}];

small…region…center=
MeshCoordinates[region•discretized•list[[random…component]]][[
Last@RandomChoice[
(Values@weight…table[[random…component]])->
(Keys@weight…table[[random…component]])
]
]]
;

small…region•shifted=
TransformedRegion[small…region•centered,
TranslationTransform[small…region…center]];

counter++;

];

{small…region…center, small…region•shifted}

, how…many]

];

randomRegionPosition[big…region_,
small…region_,
opt: OptionsPattern[]]:= randomRegionPosition[big…region, small…region,1, opt]


Weakly avoiding overlaps:

Clear[weakly…self…avoiding…regions];

weakly…self…avoiding…regions[random…position…generator_,how…many_,avoidance_:10]:=

Module[{random…pos…list,possible…new…points,
chosen…position,shifted,small…region•centered},

random…pos…list=random…position…generator[1];

Do[
possible…new…points = random…position…generator[avoidance];

chosen…position=
First@MaximalBy[possible…new…points,
RegionDistance[RegionUnion[random…pos…list[[All,2]]], #[[1]]]&];

AppendTo[random…pos…list,chosen…position];

,
how…many-1
];

random…pos…list

]


Strongly avoiding overlaps :

strongly…self…avoiding…regions[
big…region_, small…region_,
random…position…generator_,
how…many_] :=
Last@Last@
Reap@Nest[
RegionDifference[#,
Sow[random…position…generator[#,
small…region]]] &, big…region,
how…many]


### Discussion

I spent some time trying to find a way to position a smaller region in a larger region efficiently.

Two issues that might come to mind are:

• How to maximize the chances that on the first try the smaller region fits into the larger region

• If that fails, what is the next best step to take and the next after that ? In particular, how to avoid trying over and over again in parts of the region that were already explored ?

• Taking points too close to the boundary of the bigger region will likely cause the smaller region not to fit inside. Prefer points "in the bulk".

• Store all previous guesses and avoid regions near them.

The first point implies that the random distribution can not be uniform and should decrease with a so far undefined bulk center. For point 2 we can just take n samples of the distribution without the self avoiding walk like constraint and take the one that maximizes distance from previous steps.

The code below chooses the points in a mesh of the larger region obtained from DiscretizeRegion

### Coding the idea above

• First question

What is the "bulk center" and how to find it ?

The method proposed here consists in making a mesh of the larger region and finding the GraphCenter of the MeshConnectivityGraph (version 13.1):

The cat (in the question above not in the hat):

cat…outline•discretized =
cat…outline // DiscretizeRegion;

centers =
cat…outline•discretized // MeshConnectivityGraph //
GraphCenter


The output is a list of mesh cells. We can view their positions using

HighlightMesh[cat…outline•discretized, centers]


In the case of the smiley there are disconnected regions. We will separate them and find the GraphCenter for each component :

smiley•discretized = smiley // DiscretizeRegion

centers =
smiley•discretized // MeshConnectivityGraph //
WeaklyConnectedGraphComponents // Map[GraphCenter] //
Flatten[#, 1] &

HighlightMesh[smiley•discretized, centers]


The vertices of the graph are the cell labels in the mesh and have the form {0,m}. To find the coordinate of the vertex I used

pts=smiley•discretized // MeshCoordinates //
Extract[List /@ Last /@ centers] // Point

• Second question

How to define the random distribution ?

The distribution will be defined on the mesh points and uses the weight option of RandomChoice. The weights below are defined to be (decreasing) functions of the GraphDistance from the GraphCenter. Only one GraphCenter is chosen per region component. The function below allows the user to specify how they would like to distribute the weights and has a default Lorentzian/Cauchy distribution.

( Note The code below has been changed in the code section )

center…distance…weighted…choice[mesh…graph_,distribution_:(1/(1+#^2) &)]:=

Module[{centers,graph…components,weight…table},

graph…components=mesh…graph//WeaklyConnectedGraphComponents;

centers=graph…components //Map[GraphCenter] //Map[First];

weight…table=
Table[
Association@Table[
v-> distribution@GraphDistance[
graph…components[[component]],
centers[[component]],
v],

{v,VertexList[graph…components[[component]]]}
],

{component,Length@graph…components}
];

RandomChoice[
(Join@@Values/@weight…table)->(Join@@Keys/@weight…table)
]

]


As an example;

smiley•discretized•graph =
smiley•discretized // MeshConnectivityGraph;

HighlightMesh[smiley•discretized,
center…distance…weighted…choice[
smiley•discretized•graph]]


The mesh coordinate can be obtained using:

MeshCoordinates[
smiley•discretized][[Last@
center…distance…weighted…choice[
smiley•discretized•graph]]]


To prevent picking the same points too often we may store previous points and choose points that are furthest from previous choices. Consider "previously" chosen points:

pts = Table[
MeshCoordinates[
smiley•discretized][[Last@
center…distance…weighted…choice[
smiley•discretized•graph]]], 10]


Then some new candidates to add to the list:

possible…new…points =
Table[MeshCoordinates[
smiley•discretized][[Last@
center…distance…weighted…choice[
smiley•discretized•graph]]], 10]


Take then the point furthest away from previous points:

MaximalBy[possible…new…points,
RegionDistance[Point@pts]]


Links below are generated automatically using Mathematica on the text of this answer. May contain errors .