I honestly have no idea where to begin with this problem. In summary, I have a 2D coarse grid with an intersecting line. For an easy example, let's assume it's a 4x4 grid.

Coarse Grid with intersect.

I wish to pass through each grid block multiple times and subsequently divide each grid block in half that either has the line passing through it (or possibly through one of its vertices).

1st pass 3rd pass 4th pass


After a set number of passes, I want to assign a property to each grid block that is intersected. (example, set color to blue).

Assign properties

The only requirement is the algorithm provide the final set of refined grid blocks have at least one shared side until the termination points. The final output would be a list of vertices and their cartesian coordinates.

Thanks in advance.

UPDATE: The grid blocks (which will eventually total in the millions) will have CPU intensive calculations performed on them. So the ideal solution will only refine the grid block size close to the line while leaving the areas away from the line in their "coarse" state.

  • $\begingroup$ Is this a homework question by any chance? $\endgroup$ – Sjoerd C. de Vries Jul 25 '12 at 1:53
  • 3
    $\begingroup$ @SjoerdC.deVries, Nope. An attempt to solve a problem the industry software gurus have yet to provide a solution to. If you think it's a homework problem due to my awesome Visio graphics, I thank you. $\endgroup$ – kale Jul 25 '12 at 1:57
  • 1
    $\begingroup$ Useful algorithms: stackoverflow.com/questions/99353/… $\endgroup$ – Sjoerd C. de Vries Jul 25 '12 at 2:04
  • $\begingroup$ I'd love to to this one; it's not very difficult either, but I'm on vacation and my iPad doesn't run mma yet. The approach would be a recursive algorithm that finds all squares in your grid through which the segment goes using the above link. Each square is then divided in four sub squares on which the same function is recursively run. An additional depth parameter is decreased with every call and is used as a stopping criterion. $\endgroup$ – Sjoerd C. de Vries Jul 25 '12 at 2:09
  • 1
    $\begingroup$ Possibly useful: stackoverflow.com/questions/6691491/… $\endgroup$ – Oleksandr R. Jul 25 '12 at 4:05

This is my implementation using Graphics primitives and rules. Here's the final result; the implementation details and edge cases follow.

enter image description here

1. General approach

First, we start with a single square and build up a test grid:

square = Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}];
grid = Graphics[{EdgeForm[Black], FaceForm[None], 
   Table[Transpose@First@square + {i, j} // Transpose // Polygon, {i, 4}, {j, 4}]}];

enter image description here

For use in a later step, I'm translating the square explicitly here instead of using Translate or GeometricTransform, which remains in the FullForm (Normal doesn't get rid of it either, although it does for simple cases)

Next, I define a function intersectQ that tells you whether a line segment intersects a polygon. This uses an internal function Graphics`Mesh`IntersectQ, which is only available in v8 onwards.

intersectQ[line_, poly_] := Or @@ With[
    {l = Developer`PartitionMap[Line, First@poly /. {x_, y__} :> {x, y, x}, 2, 1]}, 
    Graphics`Mesh`IntersectQ[{line, #}] & /@ l

The logic is simple — check to see if the line intersects any of the sides of the polygon (Note: make sure to read section 3 below).

Next, a function/rule that divides a square if a line passes through it (using the above test):

divide[line_] := p : Polygon[v_] :> 
        {c = Mean@v, s = (Most@# + Rest@#)/2 &[v /. {x_, y__} :> {x, y, x}] /.
             {x__, y_} :> {y, x, y}}, 
        Polygon@{v[[#]], s[[# + 1]], c, s[[#]]} & /@ Range@4
    ] /; intersectQ[line, p]

Finally, we apply this to the grid and an arbitrary line to get the refined grid:

refGrid = With[{line = Line[{{3/2, 5/2}, {9/2, 9/2}}]},
    Show[Nest[# /. divide[line] &, grid, 3], Graphics[{Red, line}]]]

You can nest this how many ever times you want.

2. Obtaining the coordinates of the top-right corners of the grids

To get the top-right vertex from each grid block (left image), you can do:

Cases[refGrid, Polygon[v_] :> 
    With[{c = Mean@v}, Select[v, ArcTan @@ (# - c) == π/4 &]], Infinity]~Flatten~1

enter image description here

or get them just for the intersecting grids by using the above with intersectQ (right image)

enter image description here enter image description here

3. Edge cases

ruebenko pointed out some edge cases wherein lines that overlap with an edge of the grid are not considered as intersecting the grid. This is what it looks like at present:

enter image description here

Notice that the third and fourth large grids in the first row are not divided properly.

One simple way suggested by ruebenko was to find the area of the triangle formed by the end-points of the line and each vertex of the polygon — if the area of any triangle is 0, then it is overlaps with an edge. This is good, but fails when the end-points of the line and the end-points of the edge are collinear, but the line and the edge don't overlap.

A mathematical solution to this can be obtained by recognizing that a line segment $L$ and a polygon $P$ are both convex sets. So a line can be conclusively said to be touching a polygon or intersecting it if

$$\exists\ \{x,y\}: \{x,y\}\in L \bigwedge \{x,y\}\in P$$

We can tell Mathematica exactly that to get our condition:

divideQ[line_, poly_] := divideQ[line, poly] = 
        Exists[{x, y}, 
            With[{P = Array[p, 4], Q = Array[q, 2]}, 
                Exists[P, {x, y} == (P.First@poly) && Total@P == 1 && 
                    (And @@ Thread[P >= 0])
                ] && 
                Exists[Q, {x, y} == (Q.First@line) && Total@Q == 1 && 
                    (And @@ Thread[Q >= 0])
        ], Reals

Now simply use the above condition instead of intersectQ in section 1 above and you're all set. I've used memoization so that the result for the same polygon-line pair is cached. Using this, the above edge case is taken care of:

enter image description here

While this is a bit slower than the earlier one, it is cleaner and correct for all cases. If one knows well in advance that there are no such edge cases, you could go with intersectQ. All the figures in this answer barring the incorrect one for the edge case were done using divideQ.

  • $\begingroup$ How did you get to Graphics`Mesh`IntersectQ? Makes me wonder (again) what kind of badly desired CAGD goodies are already available... $\endgroup$ – Yves Klett Jul 25 '12 at 8:20
  • $\begingroup$ @R.M, It seems your method is quicker than wxffles, works with multiple line segments (although if they cross, it seems to blow up), but doesn't quite provide as "crisp" of graphics. One question: what is the easiest way to grab the top right vertices coordinates from each grid block? $\endgroup$ – kale Jul 25 '12 at 12:43
  • $\begingroup$ @YvesKlett I don't remember... I found that (and other related functions) when looking for something else. They're all undocumented though, so you'll have to spend some time figuring out what the arguments are. Another one I've used a few times is Graphics`Mesh`PolygonArea@poly to compute the area of a polygon. $\endgroup$ – rm -rf Jul 25 '12 at 14:14
  • $\begingroup$ @R.M, it shows in notifications that you responded, but I'm not seeing it here... I can only see the first few words. Also, is there a way to Parallelize this for increased speed? $\endgroup$ – kale Jul 25 '12 at 15:59
  • $\begingroup$ @kalewallace I'm not sure about parallelizing off the top of my head (will give it a try later when I have some time). $\endgroup$ – rm -rf Jul 25 '12 at 16:08

Here's my go at it.

This tells you if two line segments intersect (unless they lie on the same line, in which case it fails horribly):

segmentsIntersect[{a_, b_}, {p_, q_}] :=
  Module[{s, t, soln},
    soln = NSolve[a + t (b - a) == p + s (q - p), {s, t}];
    If[Length@soln == 0, False, (0 <= s <= 1 && 0 <= t <= 1) /. soln[[1]]]];

This does the box-segment intersection using the above:

  segment : {a : {ax_, ay_}, b : {bx_, by_}},
  r : Rectangle[{minx_, miny_}, {maxx_, maxy_}]
] :=
If[minx <= ax <= maxx && miny <= ay <= maxy || 
  minx <= bx <= maxx && miny <= by <= maxy, True,
  segmentsIntersect[segment, {{minx, miny}, {minx, maxy}}] ||
  segmentsIntersect[segment, {{minx, miny}, {maxx, miny}}] ||
  segmentsIntersect[segment, {{maxx, maxy}, {minx, maxy}}] ||
  segmentsIntersect[segment, {{maxx, maxy}, {minx, maxy}}]];

These functions can be replaced with an efficient and correct algorithm or your choice. Also, here will be where you decide whether corners and edges are in or out.

This recursively refines a given box:

refine[segment_, r : Rectangle[{minx_, miny_}, {maxx_, maxy_}], n_] :=
  If[n > 0 && segmentIntersectsRectangle[segment, r], 
    With[{midx = (minx + maxx)/2, midy = (miny + maxy)/2}, {
      refine[segment, Rectangle[{minx, miny}, {midx, midy}], n - 1],
      refine[segment, Rectangle[{midx, miny}, {maxx, midy}], n - 1],
      refine[segment, Rectangle[{minx, midy}, {midx, maxy}], n - 1],
      refine[segment, Rectangle[{midx, midy}, {maxx, maxy}], n - 1]
    }], r];

This creates some graphics:

refinedGrid[segment : {{ax_, ay_}, {bx_, by_}}, n_] :=
    minx = Floor[Min[ax, bx]],
    miny = Floor[Min[ay, by]],
    maxx = Ceiling[Max[ax, bx]],
    maxy = Ceiling[Max[ay, by]]},
    Graphics[{FaceForm[], EdgeForm[Black],
      Table[refine[segment, Rectangle[{x, y}, {x + 1, y + 1}], n],
        {x, minx, maxx - 1}, {y, miny, maxy - 1}],
      Red, Line[segment]}]]

As so:

refinedGrid[{{0.5, 0.5}, {4.5, 2.5}}, 4]


Exercise for the reader:


  • $\begingroup$ Nice approach. However, I'd start with a predefined coarse grid and then refine from there, not have a varying grid based on length/direction of line. The method overall seems slower than R.M.'s, but as you show (and I'll have to study to figure out) that it works for non-linear shapes. $\endgroup$ – kale Jul 25 '12 at 12:46
  • $\begingroup$ The more I study your algorithm, the more I like it. I've made it use a set coarse grid and parallelized it to make it quicker as well. How much performance penalty is due to the NSolve? $\endgroup$ – kale Jul 25 '12 at 21:11
  • 3
    $\begingroup$ @kalewallace, I am sure you could improve on the performance of NSolve. You could code (e.g. compile) the symbolic solution to Reduce[a + t (b - a) == p + s (q - p), {s, t}] and just insert coords. $\endgroup$ – user21 Jul 26 '12 at 0:22

Does this work for you:

linesquares[coords_, n_] := Module[{image},
  image = Rasterize[Graphics@Line[coords], RasterSize -> 2^n, ImageSize -> 2^n];
  Position[ImageData@image /. {1., 1., 1.} -> "", {_, _, _}]

linesquares[{{1, 1}, {4, 3}}, 3]


{{6, 7}, {6, 8}, {7, 6}, {7, 7}, {7, 8}, {8, 4}, {8, 5}, {8, 6}, {8, 7}, {9, 3}, {9,4}, {9, 5}, {10, 1}, {10, 2}, {10, 3}, {10, 4}, {11,1}, {11, 2}}


And to render the output:

Graphics[Rectangle /@ linesquares[{{1, 1}, {4, 3}}, 3]]

Mathematica graphic

  • $\begingroup$ Nice, but--and as I will now update my question--all grid blocks will eventually have very CPU intensive calculations performed on them. So where your method, it seems, involves a immediately refined grid, I want to start with a coarse grid and then only refine the areas close to the line (to save CPU cycles). $\endgroup$ – kale Jul 25 '12 at 1:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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