See links at the end for documentation links of Mathematica functions used in this answer.
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.
The more advanced approach
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
The more advanced method :
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}
]
The more advanced method:
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 ?
Consider the following possible answers/heuristics:
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
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
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 .
{RegionWithin,MeshConnectivityGraph,AbsoluteTiming,Table,Null,Graphics,Opacity,Magenta,BoundaryDiscretizeRegion,RGBColor,PointSize,Point,All,With,DiscretizeRegion,Clear,Module,TransformedRegion,TranslationTransform,RandomPoint,While,Not,If,Print,Default,Abort,Options,OptionsPattern,ConnectedMeshComponents,Map,GraphCenter,First,Association,OptionValue,GraphDistance,VertexList,Length,RandomInteger,MeshCoordinates,Last,RandomChoice,Values,Keys,Do,MaximalBy,RegionDistance,RegionUnion,AppendTo,Reap,Nest,RegionDifference,Sow,In,For,HighlightMesh,WeaklyConnectedGraphComponents,Flatten,Extract,List,Join,Take,Links}
