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.

$\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$
    – corey979
    Sep 14, 2016 at 18:30
  • $\begingroup$ @corey979 I fear it's a solution to a different problem... :o $\endgroup$
    – kirma
    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$
    – george2079
    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$
    – kirma
    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$
    – yarchik
    May 15, 2022 at 9:20

5 Answers 5

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

$\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$
    – rcollyer
    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$ 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$
    – george2079
    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 &]}]]
$\endgroup$
4
  • 2
    $\begingroup$ MinMax /@ Transpose@reg[[1]] could be replaced with more general-purpose RegionBounds[reg], BTW. $\endgroup$
    – kirma
    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$
    – Jason B.
    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$
    – kirma
    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$ 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]] &]]}]]
$\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$
    – user484
    Sep 14, 2016 at 19:58
  • $\begingroup$ @Rahul Hmm, didn't know, or maybe just remember that one! $\endgroup$
    – kirma
    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$
    – kirma
    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

$\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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.