RE (long): ... RE: rectangle intersection
- To: mathgroup at smc.vnet.net
- Subject: [mg36238] RE (long): [mg36199]... RE: [mg36093] rectangle intersection
- From: "Wolf, Hartmut" <Hartmut.Wolf at t-systems.com>
- Date: Wed, 28 Aug 2002 04:15:58 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Bobby,
here my solution(s) as promised. As you stated, my published solution was
only for rectangles oriented parallel to the X and Y axis, as this is
perhaps the predominant application.
The idea was, to decide whether all vertices of a rectangle lie to one side
of the other. This idea expands to the general case. And effectively it is
the same idea as of your suggestion, as well as of the improved version
after Garry Helzer, i.e. the "constructive" solutions, as opposed to those
involving more abstract mathematical reasoning as of David Park and Andrzej
Kozlowski. And these, esp. Andrzej's, are very interesting, as they show how
to derive a solution just from this kind of reasoning.
First some preliminaries: random rectangles
I had chosen to describe the rectangles as a list of counterclockwise
vertices. To generate them, I had used a similar description as David Park:
toRectangle[corner_, \[Theta]_, base_, side_] :=
Module[{vb = {Cos[\[Theta]], Sin[\[Theta]]}* base,
vs = {-Sin[\[Theta]], Cos[\[Theta]]}* side},
{corner, corner + vb, corner + vb + vs, corner + vs}]
makeRect[] :=
toRectangle[{Random[], Random[]}, \[Pi]/2 Random[], Random[],
Random[]]
rectLine[{pFirst_, pRest__}] := Line[{pFirst, pRest, pFirst}]
ShowRects[r1_, r2_] :=
Show[Graphics[{PointSize[0.05], Hue[.5], rectLine[r1],
MapIndexed[{Point[#1], GrayLevel[0], Text[First[#2], #1]} &,
r1],
Hue[0], rectLine[r2],
MapIndexed[{Point[#1], GrayLevel[0], Text[First[#2], #1]} &,
r2]}],
Background -> GrayLevel[.8], AspectRatio -> Automatic, PlotRange ->
All]
test cases:
rect[1] = toRectangle[{0, 0}, 0, 1, 1]; (* from David *)
rect[2] = toRectangle[{1, 1}*0.9, \[Pi]/4, 2, 1];
ShowRects[rect[1], rect[2]]
rx = rect[2] /. {x_?NumericQ, y_} :> {x, y} + {-0.1, 0.15};
ShowRects[rect[1], rx]
rxx = {{-.5, .2}, {1.5, .2}, {1.5, .8}, {-.5, .8}};
ShowRects[rect[1], rxx]
test data:
recs1000 = Table[makeRect[],{1000}];
Now here my solution. It is written such as to communicate the idea,
call it
elegant:
offSide[r2_][{p1_, p2_}] := And @@ ((p2 - p1).(# - p1) <= 0 &) /@
r2
rectOverlap[r1_, r2_] :=
=AC (Or @@ offSide[r2] /@ Partition[r1, 2, 1, {1, 1}] \[Or]
Or @@ offSide[r1] /@ Partition[r2, 2, 1, {1, 1}])
I came to this, when I tried to expand my solution for evenly oriented
rectangles to the general case. I first tried oblique coordinates, sighting
along the sides of the rectangles, then saw that giving them in the dual
base is computationally simpler, to recognize i.e. just the distance of a
point from the side's (straight line).
Correctness:
rectOverlap[rect[1], rect[2]]
True
rectOverlap[rect[1], rx]
False
rectOverlap[rect[1], rxx]
True
Performance:
(testr1 =
MapThread[
rectOverlap, {recs1000, RotateLeft[recs1000, 1]}]); // Timing
{3.635 Second, Null}
The costs for RotateLeft are negligable in Timing.
Garry Helzer had proposed a solution which wasn't correct (e.g. for rect[1},
rxx) as he noted himself. Bobby Treat however noticed that it can be fixed.
Here my version thereof. (You can see how similar it is, starting from
different reasoning.)
rightSide[a_, {b_, c_}] := Det[Prepend[#, 1] & /@ {a, b, c}] < 0
vertexExcluded[r1_, r2_] :=
=AC (Or @@ And @@@
Outer[rightSide[#2, #1] &, Partition[r2, 2, 1, {1, 1}], r1, 1]
\[Or]
Or @@ And @@@
Outer[rightSide[#2, #1] &, Partition[r1, 2, 1, {1, 1}], r2, 1])
(testr4 =
MapThread[
vertexExcluded, {recs1000, RotateLeft[recs1000, 1]}]); //
Timing
{3.646 Second, Null}
testr1 === testr4
True
of quite similar performance.
Bobby Treat gave a different implementation of this (as seen below) which is
marginally faster.
Now here Bobby's solution, rewritten in my style (which makes it slightly
faster):
extent[r1_, r2_] := {Min[#], Max[#]} & /@
((r1[[{1, 2}]] - r1[[{2, 3}]]).Transpose[r2])
cannotIntersect[{{min1_, max1_}, {min2_, max2_}}] :=
max2 < min1 || max1 < min2
intersects[r1_, r2_] :=
Catch[If[cannotIntersect[#], Throw[False]] & /@
Flatten[Transpose[
Join[Outer[extent, {r1}, {r1, r2}, 1],
Outer[extent, {r2}, {r2, r1}, 1]], {1, 3, 2}], 1]; True]
(testr5 =
MapThread[intersects, {recs1000, RotateLeft[recs1000, 1]}]); //
Timing
{1.502 Second, Null}
testr5 === testr1
True
Now, this here is a rewrite of my rectOverlap. It explicitely writes down
the tests for all sides. Such we come to use of the non-strict evaluation of
And[ ] and Or[ ]. This effectively corresponds to a Catch and Throw.
rectOverlap2[r1 : {p1_, p2_, p3_, p4_}, r2 : {q1_, q2_, q3_, q4_}] :=
! Or[
((p2 - p1).(# - p1) <= 0 &) /@ And @@ r2,
((p3 - p2).(# - p2) <= 0 &) /@ And @@ r2,
((p4 - p3).(# - p3) <= 0 &) /@ And @@ r2,
((p1 - p4).(# - p4) <= 0 &) /@ And @@ r2,
((q2 - q1).(# - q1) <= 0 &) /@ And @@ r1,
((q3 - q2).(# - q2) <= 0 &) /@ And @@ r1,
((q4 - q3).(# - q3) <= 0 &) /@ And @@ r1,
((q1 - q4).(# - q4) <= 0 &) /@ And @@ r1]
(testr2 =
MapThread[
rectOverlap2, {recs1000, RotateLeft[recs1000, 1]}]); // Timing
{1.502 Second, Null}
testr2 === testr1
True
Equal in performance to Bobby's solution. The use of And with ((p2 - p1).(#
- p1) <= 0 &) /@ And @@ r2, etc. is a bit tricky: it effectively prevents
the evaluation of (p2 - p1).(q1 - p1) <= 0 etc., such that this is evaluated
within And (non-standard evaluation) after mapping. And @@ ((p2 - p1).(# -
p1) <= 0 &) /@ r2 to the contrary is less performant. (You can see here how
clever the language is designed, to allow use of And as a container, until
the time comes to execute!)
This solution can be improved a little bit; we replace mapping by threading
and delay execution of the pure functions' body to the application of And
now by a different method (since we don't Map):
rectOverlap3[r1 : {p1_, p2_, p3_, p4_}, r2 : {q1_, q2_, q3_, q4_}] :=
! Or[
Block[{v = p2 - p1},
Or[And @@ Thread[Unevaluated[(Unevaluated[v.(# - p1) <= 0]
&)[r2]]],
And @@ Thread[Unevaluated[(Unevaluated[v.(# - p2) >= 0]
&)[r2]]]]],
Block[{v = p4 - p1},
Or[And @@ Thread[Unevaluated[(Unevaluated[v.(# - p1) <= 0]
&)[r2]]],
And @@ Thread[Unevaluated[(Unevaluated[v.(# - p4) >= 0]
&)[r2]]]]],
Block[{v = q2 - q1},
Or[And @@ Thread[Unevaluated[(Unevaluated[v.(# - q1) <= 0]
&)[r1]]],
And @@ Thread[Unevaluated[(Unevaluated[v.(# - q2) >= 0]
&)[r1]]]]],
Block[{v = q4 - q1},
Or[And @@ Thread[Unevaluated[(Unevaluated[v.(# - q1) <= 0]
&)[r1]]],
And @@ Thread[Unevaluated[(Unevaluated[v.(# - q4) >= 0]
&)[r1]]]]]]
(testr3 =
MapThread[
rectOverlap3, {recs1000, RotateLeft[recs1000, 1]}]); // Timing
{1.442 Second, Null}
testr3 === testr1
True
So this is fastest by a small margin, but has lost almost all elegance.
The
outer Unevaluated is necessary to prevent evaluation within Thread, the
inner does the essential trick noted above.
Hope this was of some interest,
Hartmut
>-----Original Message-----
>From: DrBob [mailto:majort at cox-internet.com]
To: mathgroup at smc.vnet.net
>Sent: Monday, August 26, 2002 10:16 AM
>Subject: [mg36238] [mg36199] RE: [mg36162] Fwd: [mg36140] Fwd: [mg36124] RE:
>[mg36093] rectangle intersection
>
>
>Garry,
>
>Also note your solution requires rectangle points to be in clockwise
>order (mine doesn't), but yours works for arbitrary convex polygons as
>written.
>
>Bobby
>
>-----Original Message-----
>From: DrBob [mailto:majort at cox-internet.com]
To: mathgroup at smc.vnet.net
>Subject: [mg36238] [mg36199] RE: [mg36162] Fwd: [mg36140] Fwd: [mg36124]
>RE: [mg36093]
>rectangle intersection
>
>Garry,
>
>Here's a solution using your "LeftSide" concept; it works perfectly but
>takes twice as much time as my solution. Both solutions look at every
>vertex of both rectangles, but mine uses two sides from each and yours
>requires looking at all four sides of each rectangle. I'd think yours
>should be a trifle faster than this, though. There may be efficiencies
>I'm missing (in both solutions).
>
>ClearAll[cis, rect, pickRect, extent, cannotIntersect, intersects,
>daveRect]
>cis[t_] := {Cos@t, Sin@t}
>rect[{pt : {_, _}, angle_, {len1_, len2_}}] := Module[{pt2},
> {pt, pt2 =
> pt + len1 cis[angle],
> pt2 - len2 cis[angle - Pi/2], pt - len2 cis[angle - Pi/2]}]
>daveRect := {{Random[], Random[]}, Random[] + Pi/2, {Random[],
>Random[]}}
>pickRect := rect@daveRect
>extent[r1_,
> r2_] := {Min@#, Max@#} & /@ ((Take[r1, 2] - r1[[{2,
>3}]]).Transpose@r2)
>cannotIntersect[{{min1_, max1_}, {min2_,
> max2_}}] := max2 < min1 || min2 > max1
>intersects[r1_, r2_] := Catch[
> If[cannotIntersect[#], Throw[False]] & /@
>Flatten[Transpose[Outer[extent, \
>{r1}, {r1, r2}, 1]~Join~Outer[extent, {r2}, {r2, r1}, 1], {1, 3, 2}],
>1];
> Throw[True]]
>
>ClearAll[leftSide,leftIntersects,sides]
>sides[a_List]:=Partition[Join[a,{First@a}],2,1]
>leftSide[{a_,b_},{{c_,d_},{e_,f_}}]:=-b c+a d+b e-d e-a f+c f>0
>leftSide[a:{{_,_}..},b:{{_,_},{_,_}}]:=leftSide[#,b]&/@a
>leftSide[a_List,b:{{{_,_},{_,_}}..}]:=leftSide[a,#]&/@b
>leftIntersects[a_,b_]:=!Or@@(And@@#&/@leftSide[a,sides@b])&&!
> Or@@(And@@#&/@leftSide[b,sides@a])
>
>davePairs={daveRect,daveRect}&/@Range[10000];
>rectanglePairs=Map[Reverse@rect[#]&,davePairs,{2}];
>Timing[right=intersects[Sequence@@#]&/@rectanglePairs;]
>Timing[test=leftIntersects[Sequence@@#]&/@rectanglePairs;]
>right\[Equal]test
>
>{3.187999999999999*Second, Null}
>{6.765000000000001*Second, Null}
>True
>
>Bobby Treat
>
>-----Original Message-----
>From: Garry Helzer [mailto:gah at math.umd.edu]
To: mathgroup at smc.vnet.net
>Subject: [mg36238] [mg36199] [mg36162] Fwd: [mg36140] Fwd: [mg36124] RE:
>[mg36093] rectangle
>intersection
>
>As Daniel Lichtblau pointed out, the statement below about vertices is
>nonsense. Consider two overlapping rectangles arranged as a cross. You
>need to compute intersections and test them instead of vertices.
>
>Begin forwarded message:
>
>> From: Garry Helzer <gah at math.umd.edu>
To: mathgroup at smc.vnet.net
>> Date: Fri Aug 23, 2002 12:25:13 AM US/Eastern
>> Subject: [mg36238] [mg36199] [mg36162] [mg36140] Fwd: [mg36124] RE:
>[mg36093] rectangle
>intersection
>>
>> Begin forwarded message:
>>
>> Dear colleagues,
>>
>> any hints on how to implement a very fast routine in Mathematica for
>> testing if two rectangles have an intersection area?
>> Thanks in advance
>> Frank Brand
>>
>>
>> Here is one approach.
>>
>> Given three points {x1,y1},{x2,y2},{x3,y3}, switch to homogenous
>> coordinates a={1,x1,y1}, b={1,x2,y2}, c={1,x3,y3}. Then
>> Sign[Det[{a,b,c}]] is +1 if and only if the point a lies on your left
>as
>> you walk along the line though b and c in the direction from b to c.
>> ( If the result is zero, then a lies on the line.)
>>
>> The value of the determinant is
>x2y3-x3y2-x1y3+x3y1+x1y2-x2y1, and the
>> speed of the algorithm depends essentially on how fast this quantity
>can
>> be computed. Suppose we write a function LeftSide[a,{b,c}] that
>computes
>> the sign of the determinant.
>>
>> Now let {p1,p2, . . ., pn} be a list of vertices (pi={xi,yi}) of a
>> convex polygon traced counterclockwise. Then a lies within or on the
>> boundary of the polygon if and only if none of the numbers
>> LeftSide[a,{pi,p(i+1)}] are -1. That is, if -1 does not appear in the
>> list LeftSide[a,#]&/@Partition[{p1,p2,. . .,pn,p1},2,1].
>>
>> Now use the fact that if two convex polynomials overlap, then some
>> vertex of one of them must lie inside or on the boundary of
>the other.
>>
>> If an overlap of positive area is required, then the check is that
>only
>> +1 appears--not that -1 does not appear.
>>
>> For two rectangles ( or parallelograms) this approach requires the
>> evaluation of 16 determinants, so it may be a bit expensive. If the
>> points have rational coordinates, then (positive) denominators may
be
>> cleared in the homogeneous coordinates and the computations can be
>done
>> in integer arithmetic, at the cost of at least three more
>> multiplications per determinant.
>>
>>
>Garry Helzer
>Department of Mathematics
>University of Maryland
>College Park, MD 20742
>301-405-5176
>gah at math.umd.edu
>>
>>
>
>
>
>
>
>
>