18
$\begingroup$

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]}]]

enter image description here

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.

Improve this question
$\endgroup$
6
  • $\begingroup$ Not sure if this is what you are looking for, but maybe it'll be useful in one way or another. Or maybe not. $\endgroup$ Commented Sep 14, 2016 at 18:30
  • $\begingroup$ @corey979 I fear it's a solution to a different problem... :o $\endgroup$ Commented Sep 14, 2016 at 18:36
  • $\begingroup$ Is your region always a polygon or do you need to work with more general region objects? $\endgroup$ Commented Sep 14, 2016 at 18:39
  • $\begingroup$ @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. $\endgroup$ Commented Sep 14, 2016 at 18:40
  • $\begingroup$ 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? $\endgroup$ Commented May 15, 2022 at 9:20

5 Answers 5

Reset to default
11
$\begingroup$

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]

enter image description here

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]

enter image description here

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

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ I was trying to use MinValue with ManhattanDistance, which amounts to the same thing, but this is far and away faster. +1 $\endgroup$ Commented Sep 14, 2016 at 20:16
  • $\begingroup$ 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. $\endgroup$ Commented Sep 15, 2016 at 7:12
  • $\begingroup$ @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. $\endgroup$ Commented Sep 15, 2016 at 18:35
9
$\begingroup$

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 &]}]]
Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ MinMax /@ Transpose@reg[[1]] could be replaced with more general-purpose RegionBounds[reg], BTW. $\endgroup$ Commented Sep 14, 2016 at 18:55
  • 1
    $\begingroup$ +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. $\endgroup$ Commented Sep 14, 2016 at 19:06
  • $\begingroup$ 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. $\endgroup$ Commented Sep 14, 2016 at 19:42
  • $\begingroup$ 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). $\endgroup$ Commented Sep 15, 2016 at 7:14
6
$\begingroup$

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]] &]]}]]
Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ One can also check whether RegionDimension[...] > -Infinity, but I don't know if that's faster or slower than the Resolve[Exists[...]] method. $\endgroup$ Commented Sep 14, 2016 at 19:58
  • $\begingroup$ @Rahul Hmm, didn't know, or maybe just remember that one! $\endgroup$ Commented Sep 14, 2016 at 19:59
  • $\begingroup$ @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. $\endgroup$ Commented Sep 14, 2016 at 20:06
5
$\begingroup$

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}]]

Mathematica graphics

Improve this answer
$\endgroup$
2
$\begingroup$

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]]

enter image description here

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.