8
$\begingroup$

Suppose I have a set of "rectangles" and a set of "points":

SeedRandom[5]
rects = RandomReal[10, {5, 2, 2}]
pts = RandomReal[10, {5, 2}]

{{{0.00790584, 0.650192}, {9.89555, 9.68768}}, {{2.00866, 8.19521}, {0.897634, 9.70701}}, {{2.2991, 6.12503}, {0.96816, 5.48855}}, {{1.32548, 2.32332}, {7.76135, 5.50949}}, {{0.586896, 9.60602}, {0.982487, 0.343521}}}

{{8.06562, 4.39186}, {1.42284, 0.27687}, {0.794711, 8.59505}, {2.42136, 8.42835}, {5.54556, 7.21645}}

I want to find out which points are members of their associated rectangles. One slow possibility is to use RegionMember:

MapThread[RegionMember[Rectangle@@#1, #2]&, {rects, pts}]

{True, RegionMember[ Rectangle[{2.00866, 8.19521}, {0.897634, 9.70701}], {1.42284, 0.27687}], RegionMember[ Rectangle[{2.2991, 6.12503}, {0.96816, 5.48855}], {0.794711, 8.59505}], False, RegionMember[ Rectangle[{0.586896, 9.60602}, {0.982487, 0.343521}], {5.54556, 7.21645}]}

This approach doesn't work because RegionMember needs the first Rectangle coordinate to be strictly smaller than the second Rectangle coordinate. It also unpacks the rects variable. I would like a function inRange that returns 1 if the point is in the rectangle, and 0 otherwise, and I want to avoid unpacking. For the above example:

inRange[rects, pts]

should return:

{1, 0, 0, 0, 0}

A compiled solution is acceptable, but I would prefer a version that works for both packed arrays and mixed data types.

An alternate version using vectors only is fine, e.g.:

inRange[x1, y1, x2, y2, x, y]

A 1D version would also be interesting, with intervals instead of rectangles. The simplest approach in that case would be to use IntervalMemberQ, but that would cause unpacking.

Improve this question
$\endgroup$
1
  • $\begingroup$ Also I think we can reduce this to a 1D problem since your regions are rectangles. Take the 1D intervals for x coordinates and the 1D intervals for the y coordinates. The final answer would be the indices that satisfy both 1D cases. $\endgroup$
    – Greg Hurst
    Commented Jan 21, 2018 at 18:50

2 Answers 2

Reset to default
5
$\begingroup$

How about

inRect[{{xmin_, ymin_}, {xmax_, ymax_}}, {px_, py_}] := 
  Boole[xmin <= px <= xmax && ymin <= py <= ymax];
inRange[rects_, points_] := inRect @@@ Transpose[{rects, points}];

inRange[rects, pts]
(* {1, 0, 0, 0, 0} *)

A compiled solution has the advantage, that it can determine the correct level of listable:

inRangeC = Compile[{{r, _Real, 2}, {p, _Real, 1}},
  Boole[r[[1, 1]] <= p[[1]] <= r[[2, 1]] && 
    r[[1, 2]] <= p[[2]] <= r[[2, 2]]],
  Parallelization -> True,
  RuntimeAttributes -> {Listable}
  ]

inRangeC[rects, pts]

Edit

It was not clear from the original post that the assumption that xmin <= xmax and ymin<=ymax does not hold. This can be incorporated easily in the same style as Carl does with UnitStep.

inRangeC = Compile[{{r, _Real, 2}, {p, _Real, 1}},
  If[
    (-p[[1]] + r[[1, 1]]) (p[[1]] - r[[2, 1]]) < 0 || 
    (-p[[2]] + r[[1, 2]]) (p[[2]] - r[[2, 2]]) < 0,
    0,
    1
   ],
  Parallelization -> True,
  RuntimeAttributes -> {Listable}
]

While I don't agree with the way Carl measures the running time, I agree that a pure high-level implementation that is so short is preferable. But if you have data in a certain format, which seems to be a list of rectangles and a list of points, and your function needs to change the layout to operate on it, you need to include this change of layout in the running time (same is true for the compilation time, but this needs to be done only once).

Improve this answer
$\endgroup$
3
  • $\begingroup$ Your solution assumes that xmin is less than xmax, but my example for rects doesn't satisfy this constraint. You would need to include xmin>=px>=xmax as well to allow for arbitrarily specified rectangles. $\endgroup$
    – Carl Woll
    Commented Jan 21, 2018 at 20:48
  • $\begingroup$ Yes, this is true. I assumed you were using the standard definition of Rectangle in your real application and that I can count on this assumption. Although Rectangle works with wrongly defined regions, the Details section in its documentation (third line) clearly states the condition I have used. It is a simple fix though $\endgroup$
    – halirutan
    Commented Jan 21, 2018 at 20:54
  • $\begingroup$ I have an application in mind where the rectangles are constructed in such a way that the first point is not guaranteed to be the lower left corner. I would prefer to avoid pre-processing the rectangles to get them in the right form. $\endgroup$
    – Carl Woll
    Commented Jan 21, 2018 at 21:02
4
$\begingroup$

The fastest method I can think of is to use UnitStep in the vector approach:

inRange[x1_, y1_, x2_, y2_, x_, y_] := UnitStep[(x1-x)(x-x2)] UnitStep[(y1-y)(y-y2)]

A comparison with @halirutan's Compile approach (modified to allow rectangle specifications where the first point is not necessarily the lower left coordinate):

rects = RandomReal[10, {10^6, 2, 2}];
pts = RandomReal[10, {10^6, 2}];

x1 = rects[[All, 1, 1]];
y1 = rects[[All, 1, 2]];
x2 = rects[[All, 2, 1]];
y2 = rects[[All, 2, 2]];

x = pts[[All, 1]];
y = pts[[All, 2]];

r1 = inRange[x1, y1, x2, y2, x, y]; //RepeatedTiming
r2 = inRangeC[rects, pts]; //RepeatedTiming

r1===r2

{0.0192, Null}

{0.061, Null}

True

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ (+1) Using Subtract instead of Plus gives about a 20% speedup on my machine. $\endgroup$
    – Greg Hurst
    Commented Jan 21, 2018 at 23:35

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.