# Identifying lattice squares that are intersected by a closed curve

Given a lattice of unit squares I would like to find the number of squares that a curve (in the example below the curve is a circle) passes through and the number of squares that totally reside inside the curve.

In the example the number of squares that the curve passes through is 60 and the number that are fully within the curve is 156. These are the two numbers I would like to generate.

I have created a graphic of the situation with the following code.

grid = Table[{x, y}, {x, -8, 8, 1}, {y, -8, 8, 1}];
centres = Table[{x, y}, {x, -7.5, 7.5, 1}, {y, -7.5, 7.5, 1}];
g1 = Graphics@Line[#] & /@ grid;
g2 = Graphics@Line[#] & /@ Transpose@grid;
g3 = Graphics[{Red, PointSize[0.005], Point /@ centres}];
g4 = Graphics[Circle[{0, 0}, Sqrt]];
g5 = Graphics[{Black, PointSize[0.0075], Point@{0, 0}}];
Show[g1, g2, g3, g4, g5, Frame -> False,
PlotRange -> {{-9, 9}, {-9, 9}}, AspectRatio -> 1]


which produces: I am not sure how to set such a thing up - I have looked in the documentation but I am not seeing something similar at all.

Any help would be much appreciated.

• Just to make it clear: Intersection is valid only if the curve passes through the square (therefore it has two intersections with squares' boundary), not if it only touches the square in one point, right? Mar 17 at 12:34
• @Domen Yes that is correct - just touching doesn't count. Which is why I made the radius of the circle in the example the square root of a prime of the form 4k+3 so it wasn't the sum of two squares. Mar 17 at 12:42

rectangles = Rectangle /@ Tuples[Range[-8, 7, 1], {2}]

Count[_?(Not @ RegionDisjoint[Circle[{0, 0}, Sqrt], #] &)] @ rectangles

60

Count[_?(RegionWithin[Disk[{0, 0}, Sqrt], #] &)] @ rectangles

156


Update: To exclude rectangles "just touching" the circle, we can use

Count[_?(FreeQ[_EmptyRegion | _Point]@
RegionIntersection[Circle[{0, 0}, Sqrt], #] &)] @ rectangles

60

Count[_?(FreeQ[_EmptyRegion | _Point]@
RegionIntersection[Circle[{0, 0}, Sqrt], #] &)] @ rectangles

52

• Thank you for that - nice and brief. But I do have a question - if I choose a radius of root 53 there are eight points (-7,-2),(-7,2),(7,-2),(7,2), (-2,-7),(-2,7),(2,-7) and (2,7) where the circle touches those points - but then the first count is 68 and not 52. It is counting the two squares that touch each of those points and thus adding 16 to the total. How can I reject those situations and just count the squares where the curve actually crosses into a square? Mar 17 at 15:20
• @1729taxi, please see the update.
– kglr
Mar 18 at 1:48
• Fantastic. This is the sort of thing I have never done in Mathematica - your help is much appreciated. I shall accept this answer. Mar 18 at 7:44
• If I make my square a curve (of side length 7) then the second returned number should be zero - _Point does not work then - I tried replacing with _Line and it just gives the result of 4 - do you know how to detect this situation but return the correct number of zero? Mar 19 at 11:46

## Edit

For r=Sqrt, there are 148 squares in the curve, and there are 56 outside the curve, and 16*16-148-56=52 on the curve.

Clear["Global*"];
lines0 =
Flatten[{Table[
Line[{{i, j}, {i, j} + {1, 0}}], {i, -8, 7}, {j, -8, 8}],
Table[Line[{{i, j}, {i, j} + {0, 1}}], {i, -8, 8}, {j, -8, 7}]}];
p = {0, 0};
r = Sqrt;
range1 = Disk[p, r];
range2 = Annulus[p, {r, 4 r}];
intersec[range_] := Module[{lines, pts, edges, g, faces, squares},
lines =
Cases[RegionIntersection[range, #] & /@ lines0, _Line, Infinity];
pts = DeleteDuplicates[Flatten[lines[[;; , 1]], {2, 1}]];
edges =
Table[UndirectedEdge @@ Flatten[FirstPosition[pts, #] & /@ d], {d,
lines[[;; , 1]]}];
g = Graph[edges,
faces = PlanarFaceList[g];
squares = Select[faces, Length[#] == 4 &];
{MeshRegion[pts, Polygon /@ squares], Length[squares]}];
Labeled[Graphics[{lines0, {FaceForm[], EdgeForm[Cyan],
range1}, {FaceForm[Yellow], EdgeForm[Red],
intersec[range2] // First}, {FaceForm[Green], EdgeForm[Red],
intersec[range1] // First}}], {Style[intersec[range1] // Last,
Green, 20], Style[intersec[range2] // Last, Darker@Yellow],
Style[2*8*2*8 - intersec[range1][] - intersec[range2][],
20]}] Original

lines0 =
Flatten[{Table[
Line[{{i, j}, {i, j} + {1, 0}}], {i, -8, 7}, {j, -8, 8}],
Table[Line[{{i, j}, {i, j} + {0, 1}}], {i, -8, 8}, {j, -8, 7}]}];
Manipulate[Module[{range, lines, pts, edges, g, faces, squares},
range = Disk[p, Sqrt];
(* range=Annulus[p,{3,Sqrt}]; *)
lines =
Cases[RegionIntersection[range, #] & /@ lines0, _Line, Infinity];
pts = DeleteDuplicates[Flatten[lines[[;; , 1]], {2, 1}]];
edges =
Table[UndirectedEdge @@ Flatten[FirstPosition[pts, #] & /@ d], {d,
lines[[;; , 1]]}];
g = Graph[edges,
` 