# 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}, {x + 1, y + 1}] /.
Quiet@Solve[
Resolve[Exists[{xp, yp},
RegionMember[
RegionIntersection[Rectangle[{x, y}, {x + 1, y + 1}],
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. – corey979 Sep 14 '16 at 18:30
• @corey979 I fear it's a solution to a different problem... :o – kirma Sep 14 '16 at 18:36
• Is your region always a polygon or do you need to work with more general region objects? – george2079 Sep 14 '16 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. – kirma Sep 14 '16 at 18:40

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[#[], gridspace[]],
Round[#[], gridspace[]]} + 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[]],
Line[{gridorigin + {0, #}, gridorigin + {10, #}}] & /@
Range[0, 8, gridspace[]]} , 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[#[], gridspace[]],
Round[#[], gridspace[]]} + 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[]],
Line[{gridorigin + {0, #}, gridorigin + {10, #}}] & /@
Range[0, 10, gridspace[]]} , 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 – rcollyer Sep 14 '16 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. – Quantum_Oli Sep 15 '16 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. – george2079 Sep 15 '16 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[]];
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[] could be replaced with more general-purpose RegionBounds[reg], BTW. – kirma Sep 14 '16 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. – Jason B. Sep 14 '16 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. – kirma Sep 14 '16 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). – Quantum_Oli Sep 15 '16 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 '16 at 19:58
• @Rahul Hmm, didn't know, or maybe just remember that one! – kirma Sep 14 '16 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. – kirma Sep 14 '16 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}]] 