# How to compute intersections of circles on a lattice

Take a set of points in the plane and draw a circle of a given radius at each point. The resulting pattern of circles may look something like this:

Graphics[Circle[#, 9/10] & /@ Tuples[Range, 2], PlotRange -> {{1, 6}, {1, 6}}] Notice all of the small regions that are formed by these intersecting circles. How can I get each of these regions as say, a list of Region objects?

EDIT: I've adapted the code by kglr to compute what I want. Here is the resulting code:

tuples = Tuples[Range, 2];
disks = Disk[#, 0.9] & /@ tuples;

intersecting[d_] :=
DeleteCases[Select[disks, ! RegionDisjoint[d, #] &], d]

funcs[a_, b_] :=
If[a == b, {BooleanCountingFunction[{a}, b]},
Or -> List]]

regs[n_] :=
DeleteDuplicatesBy[
Flatten[
Function[x,
Select[
RegionIntersection[
x,
BooleanRegion[#, intersecting[x]]
] & /@ funcs[n, Length[intersecting[x]]]
, Quiet[RegionDimension[#]] == 2 &
]
] /@ disks
],
RegionCentroid
]


The resulting regions look like this:

r = regs /@ Range;
colours = {Red, Green, Blue};
Show[Flatten[(Function[x, Region[x, BaseStyle -> FaceForm[colours[[#]]]]] /@ r[[#]]) & /@ Range]]; • Which ones specifically do you want? If you just make each Circle a Region using Disk it's pretty easy to do RegionIntersection over the pairs or triples or RegionDifference over the quadruples. Feb 8, 2018 at 23:42

Update: Disjoint regions corresponding to the intersection of exactly k disks for k = 1, 2, 3, 4.

Computation of 7-by-7 example is too large for free Wolfram Cloud, so I use a smaller example with 16 disks. Using Carl's method for identifying the neighbors of each disk

tuples = Tuples[Range @ 4, 2];
disks = Disk[#, 9/10] & /@ tuples;
circles = Circle[#, 9/10] & /@ tuples;

ClearAll[nF, boolReg]
nF[x_] := Module[{d = DeleteCases[disks, x]}, Pick[d, RegionDisjoint[#, x] & /@ d, False]]

boolReg[n_] := Module[{bCF = BooleanCountingFunction[{n}, Length @ nF @ #]},
DeleteCases[RegionIntersection[#, BooleanRegion[bCF, nF @ #]], _EmptyRegion]] &

r1 = Show[Region[#, BaseStyle -> Yellow]&/@boolReg /@ disks, Graphics[{Gray, circles}]] r2 = Show[Graphics[circles], Region[#, BaseStyle -> Blue] & /@ boolReg /@ disks] r3 = Show[Graphics[circles], Region[#, BaseStyle -> Red] & /@ boolReg /@ disks] r4 = Show[Graphics[circles], Region[#, BaseStyle -> Green] & /@ boolReg /@ disks] Show[r1, r2, r3, r4] intersections = DeleteCases[RegionIntersection @@@
Subsets[(Disk[#, 9/10] & /@ Tuples[Range, 2]), {2,4}], _EmptyRegion];

Show[Graphics[{Opacity[.5, Yellow], EdgeForm[{Gray,Thick}],
Disk[#, 9/10] & /@ Tuples[Range, 2]}],
RegionPlot[#, PlotStyle -> RandomColor[]]&/@intersections] Colorcoding points by the number of disks a point lies in:

ints = DeleteCases[RegionIntersection @@@ (Subsets[(Disk[#, 9/10] & /@
Tuples[Range, 2]), {#}]), _EmptyRegion]& /@ {2, 3, 4};
colors = {Red, Green, Blue};
ints2=Join @@ (Thread /@ Transpose[{ints, colors}]);

Legended[Show[Graphics[{Opacity[.5, Yellow], EdgeForm[{Gray, Thick}],
Disk[#, 9/10] & /@ Tuples[Range, 2]}],
RegionPlot[#, PlotStyle -> #2]& @@@ ints2, PlotRange -> {{1, 6}, {1, 6}}],
SwatchLegend[{Yellow, Red, Green, Blue}, {"1", "2","3","4"}]] • Same problem as one of the other answers - the regions aren't disjoint, they overlap each other.
– Ben
Feb 9, 2018 at 0:55
• @Ben, I see. I will update if I find a fix.
– kglr
Feb 9, 2018 at 1:20
• @Ben You should mention this in the body of the question (click "edit") or/and answer b2m3a1's comment. Feb 9, 2018 at 2:17
• @Ben, please see the update.
– kglr
Feb 9, 2018 at 13:37
• @kglr Thanks, this is almost perfect but there is still a small problem. Take the small red triangles in your image generated by Show[r1, r2, r3, r4]. Some of the regions which you are plotting actually contain several of these triangles together, and some of the triangles are contained in multiple regions, so they are still overlapping. I've managed to adapt your code to fix it, so I'll add the modified code to my original post and accept your answer.
– Ben
Feb 9, 2018 at 19:09

Update

Here is a revised version of my answer that is much faster. First, I define the function RegionPieces which takes a region and its neighbors, and creates the disjoint pieces when adding each neighbor and its complement:

RegionPieces[r_, neighbors_List] := Fold[iRegionPieces, {r}, neighbors]

iRegionPieces[r_, next_] := With[
{
new = Flatten[
{RegionIntersection[#, next], RegionDifference[#, next]}& /@ r
]
},

Pick[new, Unitize @ Map[Area] @ new, 1]
]


Then, I take the region difference between the current disk and already processed disks, and then find the pieces with the new, unprocessed neighbor disks:

pieces = Flatten @ Table[
cur = disks[[i]];
old = With[{d = disks[[;;i-1]]},
Pick[d, RegionDisjoint[#, cur]& /@ d, False]
];
new = With[{d = disks[[UpTo[i+1] ;;]]},
Pick[d, RegionDisjoint[#, cur]& /@ d, False]
];

RegionPieces[
If[Length[old] > 0,
RegionDifference[cur, RegionUnion @@ old],
cur
],
new
],
{i, Length @ disks}
]; //AbsoluteTiming

Length[pieces]


{33.7915, Null}

313

About 20 times faster than my previous answer, and the same number of pieces have been generated. This answer is much faster than the accepted answer.

Here are some visualizations. First a few of the 313 pieces:

Grid @ Table[Region[pieces[[7 i+j]], ImageSize->50], {i, 0, 6}, {j, 7}] Next, I put the pieces together with colors based on the area of the piece:

With[{area = N[Area /@ pieces, 10]},
colors = area /. Thread @ Rule[
DeleteDuplicates @ area,
{Red, Yellow, Orange, Blue, Green, Gray, Pink}
]
];

Show @ Table[Region[pieces[[i]], BaseStyle->FaceForm[colors[[i]]]], {i, 313}] regions = Flatten@Table[Disk[{i, j}, .9], {i, 1, 7}, {j, 1, 7}];

twos = (RegionIntersection @@@ Subsets[regions, {2}]);
threes = (RegionIntersection @@@ Subsets[regions, {3}]);
fours = (RegionIntersection @@@ Subsets[regions, {4}]);
twoInt = Show[
Table[RegionPlot[twos[[i]], PlotStyle -> Hue[RandomReal[]],
PlotRange -> {{0, 8}, {0, 8}}], {i, 1, Length[twos]}]];
threeInt =
Show[Table[
RegionPlot[threes[[i]], PlotStyle -> Hue[RandomReal[]],
PlotRange -> {{0, 8}, {0, 8}}], {i, 1, Length[threes]}]];
fourInt =
Show[Table[
RegionPlot[fours[[i]], PlotStyle -> Hue[RandomReal[]],
PlotRange -> {{0, 8}, {0, 8}}], {i, 1, Length[fours]}]];
Show[twoInt, threeInt, fourInt] • why the small squares are left out? Feb 9, 2018 at 0:01
• Because they are not the intersection of disks. If you want to include "uncovered" sections, then merely color the background. Alternatively, you could create regions of NOT region 1 AND region 2 (etc.). Feb 9, 2018 at 0:02
• This doesn't work correctly. The regions which are computed aren't disjoint, they all overlap, it just looks correct when plotted because the overlaps get covered up.
– Ben
Feb 9, 2018 at 0:53
• @Ben: Yeah, it looks like the list twoInt will contain all points that are within at least two circles, not within exactly two circles (which is what I think you're asking for.) I think you could fix this with RegionDifference, but it'd take some finagling in the code and I don't have Mathematica easily accessible on this machine. Feb 9, 2018 at 3:28

This is the same idea as @bill_s (with the same caveat for resolution). I am just making it a bit easier for the image processing function to work. Notice that if you plod opaque disks instead of circles, each region of 2, 3 or 4 intersecting disks will share the same opacity:

circ = Graphics[{Opacity[0.2], Blue, Disk[#, 9/10]} & /@ Tuples[Range, 2], PlotRange -> {{1, 6}, {1, 6}}] Then it's very easy and quick to get each region using ClusteringComponents:

im = Rasterize@circ;
clCom = ClusteringComponents[im, 4];


so here is the region where 2 circles intersect:

clCom /. (1 | 3 | 4) -> 0 // Colorize Another approach is to use the image processing functions. Using the original .png image gives a kind of stained glass version.

pat = Binarize[Import["https://i.stack.imgur.com/5yZTG.png"]];
MorphologicalComponents[Erosion[pat, 1]] // Colorize You can access all the different regions and their properties (area, location, size, best-fit-ellipse, etc) using ComponentMeasurements. The thick borders are easy to fix by using a higher resolution image.

img = Binarize[Image[Graphics[Circle[#, 9/10] & /@ Tuples[Range, 2],
PlotRange -> {{1, 6}, {1, 6}}], ImageSize -> 2000]];
MorphologicalComponents[Erosion[img, 1]] // Colorize Not very different from David G. Stork's answer but I had prepared it so... The following computes all the intersections of n circles (for n from 2 to 4):

disks = Disk[#, 9/10] & /@ Tuples[Range, 2];
intersect[n_] := RegionIntersection[Sequence @@ #] & /@ (disks[[#]] & /@
Subsets[Range[Length@disks], {n}]) // DeleteDuplicates
GraphicsRow[Table[Show[Region /@ intersect[i]], {i, 2, 4}]] Or, if you prefer:

circles = Circle[#, 9/10] & /@ Tuples[Range, 2];
GraphicsRow[
Table[Show[Region /@ intersect[i], Graphics[circles],
PlotRange -> {{1, 6}, {1, 6}}], {i, 2, 4}]] These can be computed via RegionIntersection and RegionDifference in general. Here's such an approach for some of them:

Table[
RegionDifference[
RegionDifference[
DiscretizeRegion@
RegionIntersection[circA[{i, j}], circA[{i - 1, j}]],
DiscretizeRegion@RegionUnion[
circA[{i - 1, j - 1}],
circA[{i, j - 1}]
]
],
DiscretizeRegion@RegionUnion[
circA[{i - 1, j + 1}],
circA[{i, j + 1}]
]
],
{i, 2, 6},
{j, 2, 6}
] // Flatten // Show 