# How to quickly find all grid cells covered by a geometric region?

I'm trying to compute a complete set of unit grid cells ("pixels") which a geometric region covers at least partially. I have constructed a method that produces correct results:

With[{reg = Polygon@RandomReal[{0, 10}, {3, 2}]},
Graphics[{reg, Blue, Opacity[1/2],
Rectangle[{x, y}] /.
Quiet@Solve[
Resolve[Exists[{xp, yp},
RegionMember[RegionIntersection[Rectangle[{x, y}], reg], {xp, yp}]],
Reals], {x, y}, Integers]}]]


The problem with this implementation is that it's absurdly slow - producing this result takes several seconds for just a small triangle.

How correct results (preferably for arbitrary regions) could be produced more efficiently? Before suggesting use of Rasterize, please note that non-antialised Rasterize does not always produce identical results with code above.

• Not sure if this is what you are looking for, but maybe it'll be useful in one way or another. Or maybe not. Sep 14, 2016 at 18:30
• @corey979 I fear it's a solution to a different problem... :o Sep 14, 2016 at 18:36
• Is your region always a polygon or do you need to work with more general region objects? Sep 14, 2016 at 18:39
• @george2079 I would fancy a solution which would work at least for all (finite) semialgebraic sets, but even non-self-intersecting polygons are a start. Sep 14, 2016 at 18:40
• Do I understand correctly that input is a list of points-vertices of a polygon on 2d plane plus their connectivity? Disconnected sets and non convex sets are also allowed? May 15, 2022 at 9:20

This may seem a bit crude but it should be fast and robust - finely discretize the region and look for the mesh vertices to be in cells:

r = Polygon@RandomReal[{0, 10}, {3, 2}]
gridorigin = Floor[{Min[r[[1, All, 1]]], Min[r[[1, All, 2]]]}];
gridspace = {1, 1};
cc = Union[(({Round[#[[1]], gridspace[[1]]],
Round[#[[2]], gridspace[[2]]]} + gridorigin +
gridspace/2) &@ (# - gridorigin - gridspace/2) ) & /@
MeshCoordinates[DiscretizeRegion[r, MaxCellMeasure -> .01]]];
RegionPlot[r, PlotRange -> All,
Prolog -> {{FaceForm[Blue],
Rectangle[# - gridspace/2, # + gridspace/2] & /@ cc},
Line[{gridorigin + {#, 0}, gridorigin + {#, 10}}] & /@
Range[0, 8, gridspace[[1]]],
Line[{gridorigin + {0, #}, gridorigin + {10, #}}] & /@
Range[0, 8, gridspace[[2]]]} , AspectRatio -> Automatic]


it should work with anything that DiscritizeRegion can handle:

r = ImplicitRegion[x^2 - y^2 <= 1 || x^2 + y^2 == 4, {x, y}];
gridorigin = {-5, -5};
gridspace = {.25, .5};
cc = Union[(({Round[#[[1]], gridspace[[1]]],
Round[#[[2]], gridspace[[2]]]} + gridorigin +
gridspace/2) &@ (# - gridorigin - gridspace/2) ) & /@
MeshCoordinates[
DiscretizeRegion[r, {{-4, 4}, {-4, 4}}, MaxCellMeasure -> .01]]];
RegionPlot[r, PlotRange -> {{-4, 4}, {-4, 4}},
Prolog -> {{FaceForm[Blue],
Rectangle[# - gridspace/2, # + gridspace/2] & /@ cc},
Line[{gridorigin + {#, 0}, gridorigin + {#, 10}}] & /@
Range[0, 10, gridspace[[1]]],
Line[{gridorigin + {0, #}, gridorigin + {10, #}}] & /@
Range[0, 10, gridspace[[2]]]} , AspectRatio -> Automatic]


note the region here includes an open circle, so the "holes" are correct. here are the mesh vertices:

• I was trying to use MinValue with ManhattanDistance, which amounts to the same thing, but this is far and away faster. +1 Sep 14, 2016 at 20:16
• Does your figure 2 show the limitation of a numerical method (potential inaccuracy) by the asymmetry of the blue highlighting in the centre? Numerical routes are likely the fastest though if one is not concerned with small chances of small errors. Sep 15, 2016 at 7:12
• @Quantum_Oli, yes in the second example the region precisely mathematically hits a grid line at a tangent point so its a round off issue whether that adjacent box gets counted or not. You might discard mesh points that fall within some tolerance of the grid lines to remedy that issue. Sep 15, 2016 at 18:35

A slightly quicker method, no time to explain atm. It might be sped up by reworking the Area[RegionIntersection[...]] bit.

With[{reg = Polygon@RandomReal[{0, 10}, {3, 2}]},
{{xMin, xMax}, {yMin, yMax}} = Floor[MinMax /@ Transpose@reg[[1]]];
Graphics[{reg, Blue, Opacity[1/2],
Select[Flatten[
Table[Rectangle[{x, y}, {x + 1, y + 1}], {x, xMin, xMax}, {y, yMin, yMax}]],
Area[RegionIntersection[reg, #]] != 0 &]}]]

• MinMax /@ Transpose@reg[[1]] could be replaced with more general-purpose RegionBounds[reg], BTW. Sep 14, 2016 at 18:55
• +1 Beat me to it. You can replace Rectangle[{x,y}, {x+1, y+1}] with Rectangle[{x,y}]. This is sometimes faster than the OP code, but sometimes not. Here's my version, and two examples, one where the Select version is faster, and one where the Resolve version is faster. Sep 14, 2016 at 19:06
• It would seem it can be faster to check for nonzero area (or more specifically, existence of points belonging to the region) by using Resolve-Exists style test than using Area (which computes an useless value). Still not very fast, I'd assume this problem to be resolved at least for polygons with millisecond-level timing. Sep 14, 2016 at 19:42
• Thank you for the feedback guys, always learning! It looks like @kirma's approach below combining this rough approach but cleverly pre-evaluating the Select predicate is much faster than this (and the fastest non-discretizing approach so far). Sep 15, 2016 at 7:14

This answer combines my own code from the question and the idea in @Quantum_Oli's answer. It is already quite much faster than either, at least for polygons:

With[{reg = Polygon@RandomReal[{0, 10}, {3, 2}]},
Graphics[{reg, Blue, Opacity[1/2],
Rectangle /@
Select[Flatten[CoordinateBoundsArray@Floor@RegionBounds@reg, 1],
Apply[
Evaluate[
Resolve[Exists[{xp, yp},
RegionMember[
RegionIntersection[Rectangle[{#1, #2}], reg], {xp, yp}]],
Reals]] &]]}]]

• One can also check whether RegionDimension[...] > -Infinity, but I don't know if that's faster or slower than the Resolve[Exists[...]] method.
– user484
Sep 14, 2016 at 19:58
• @Rahul Hmm, didn't know, or maybe just remember that one! Sep 14, 2016 at 19:59
• @Rahul Now I remember why I don't use it. RegionDimension doesn't evaluate for parametric regions, thus making pre-evaluation of Select predicate impossible. Probably it is better for more complicated regions, while my answer is better for simpler ones - at least triangles. Sep 14, 2016 at 20:06

This is pretty fast, but one probably needs to be careful with MaxCellMeasure. I use DensityHistogram of the MeshCoordinates to pull out the squares.

Block[{reg = Polygon@RandomReal[{0, 10}, {3, 2}], pix},
pix = Cases[
DensityHistogram[
MeshCoordinates[
TriangulateMesh[DiscretizeGraphics[reg],
MaxCellMeasure -> 0.01]], {{1}, {1}},
PerformanceGoal -> "Speed"] //
InputForm, _Rectangle, \[Infinity]];
Graphics[{Opacity[0.5, Blue], EdgeForm[Black], pix, reg}]]


A method based on RegionDilation introduced in v12.3:

With[{reg = Polygon@RandomReal[{0, 10}, {3, 2}]},
Graphics[{reg, Blue, Opacity[1/2], Rectangle /@
Select[Flatten[CoordinateBoundsArray@Floor@RegionBounds@reg, 1],
RegionMember[RegionDilation[reg, Rectangle[{-1, -1}]]]]}]]


This method is several hundred times faster for triangles than my original method on the question.

Another method based on computing a symbolic selection function based on RegionDisjoint (introduced in v11.1). This is about ten times slower for triangles than the method above, but still faster than the original:

With[{reg = Polygon@RandomReal[{0, 10}, {3, 2}]},
Graphics[{reg, Blue, Opacity[1/2], Rectangle /@
Select[Flatten[CoordinateBoundsArray@Floor@RegionBounds@reg, 1],
Apply@
Function[{x, y}, Evaluate[! RegionDisjoint[reg, Rectangle[{x, y}]]]]]}]]


These methods are usable on other regions too, say

reg = Annulus[RandomReal[{-1, 1}, 2], Accumulate@RandomReal[{2, 5}, 2]]