8

I am recently working on a project which has an algorithm to tell whether a point is inside an area.

The area looks like this:

{"state": "Washington", "border": [[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]]}

It is easy if the area is a rectangle. However, the area is irregular. One of the idea I have is to check whether a point is in the 'inner' side of every line of the area. However, the performance is not good. Any idea?

SmAc
  • 105
  • 2
  • 5
  • 1
    You could tesselate the semi-random shapes into triangles, which would make it so much easier to decide. – Agoston Horvath May 05 '17 at 20:14
  • 2
    You can just delegate to existing implementations, but if you want to do it yourself, there are [plenty of resources regarding known algorithms](https://www.google.com/search?q=point+inside+polygon). – user2357112 May 05 '17 at 20:15

3 Answers3

12

First of all, very interesting question!! Although it might be a duplicated question, but I am still going to post another workable answer different from that post to encourage this new guy.

The idea is to use the sum of angles to decide whether the target is inside or outside. If the target is inside an area, the sum of angle form by the target and every two border points will be 360. If the target is outside, the sum will not be 360. The angle has direction. If the angle going backward, the angle is a negative one. This is just like calculating the winding number.

The provided input data [[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]] is clockwise (you can check google map). Therefore, my code assume that the positive angle is clockwise one.

Refer this for the idea: enter image description here

The following is the python code that implements it.

def isInside(self, border, target):
    degree = 0
    for i in range(len(border) - 1):
        a = border[i]
        b = border[i + 1]

        # calculate distance of vector
        A = getDistance(a[0], a[1], b[0], b[1]);
        B = getDistance(target[0], target[1], a[0], a[1])
        C = getDistance(target[0], target[1], b[0], b[1])

        # calculate direction of vector
        ta_x = a[0] - target[0]
        ta_y = a[1] - target[1]
        tb_x = b[0] - target[0]
        tb_y = b[1] - target[1]

        cross = tb_y * ta_x - tb_x * ta_y
        clockwise = cross < 0

        # calculate sum of angles
        if(clockwise):
            degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
        else:
            degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))

    if(abs(round(degree) - 360) <= 3):
        return True
    return False
Junbang Huang
  • 1,927
  • 19
  • 26
  • Very great idea! So if inside, the sum of triangles is 360, and if it is not, the sum is not 360. One thing I don't understand is abs(round(degree) - 360) <= 3. why 3? – SmAc May 05 '17 at 20:31
  • This answer is wrong. It only applies to convex polygons; for concave polygons, the angles could double back on themselves. – user2357112 May 05 '17 at 20:33
  • @user2357112 hmm, can you make it clear? – Junbang Huang May 05 '17 at 20:35
  • @SmAc yes, 3 is a threshold about how accurate you want to be. If you are very confident about my idea, you can try to make it 0.1 instead of 3. The smaller the more accurate. – Junbang Huang May 05 '17 at 20:37
  • @user2357112 I don't see how it is wrong, if I understand you correcly. By double back you mean the angle move backward? I already take that into consideration, different backward angle has a negative value. By the way, I really want to see the problem. – Junbang Huang May 05 '17 at 20:40
  • It's probably worth pointing out that this is essentially computing the [winding number](https://en.wikipedia.org/wiki/Winding_number), which has some interesting properties, in particular it is zero outside of a closed curve, and a positive integer inside. @user2357112 it's not necessary that the interior be convex; see e.g. the illustrations in the Wikipedia article. – Robert Dodier May 05 '17 at 20:41
  • @JunbangHuang: I was thinking of a case like [this](http://imgur.com/a/4dYsK). I suppose it works if you treat the angles as signed, but you should be more explicit about that. Neither your text nor your diagram say anything about handling such cases. – user2357112 May 05 '17 at 20:42
  • @user2357112 make sense. I will update it – Junbang Huang May 05 '17 at 20:43
  • @RobertDodier Thank you for your link. It would be a lot easier for me to explain. – Junbang Huang May 05 '17 at 20:44
  • Note that your code currently makes an assumption about which direction the polygon's edges wind in. If the angles add up to -360 degrees instead of 360 degrees, your program will consider the point to be outside the polygon. – user2357112 May 05 '17 at 20:52
  • @user2357112 I make that assumption based on the input data provided. I just checked that out on my google map before I wrote code, the follow a certain direction, that is why I also made that assumption. – Junbang Huang May 05 '17 at 20:56
1

Sounds like a good use case for a modified convex hull algorithm.

  1. Start out with all points inside of the "area" (since no hull is yet created).
  2. Then, as you progress through your chosen algorithm (e.g. Graham scan is O(nlog(n)) performance), if the point is chosen as part of the convex hull, it is no longer inside of the "area" -- (i.e. a point that comprises the convex hull is not is not part of the final answer).
  3. Repeat until you have created the convex hull. The leftover points which are not part of the hull is therefore inside the area, which is the answer you are looking for.
ifma
  • 3,673
  • 4
  • 26
  • 38
  • Doesn't work, because [a point might lie inside the convex hull and yet outside the polygon](http://imgur.com/a/AjvWp). – user2357112 May 05 '17 at 20:48
  • @user2357112 ahh that is an interesting corner case I failed to account for =) – ifma May 06 '17 at 15:31
0

Here is a pseudo code:

var point = [x,y] // Point to evaluate

int i, j, c = 0
int size = border.count
for (i = 0, j = size-1; i < size; j = i++) {
    var ii = border[i]
    var jj = border[j]
    if ( ((ii.[1] > point.[1]) != (jj.[1] > point.[1])) &&
        (point.[0] < (jj.[0]-ii.[0]) * (point.[1]-ii.[1]) / (jj.[1]-ii.[1]) + ii.[0])) {
        c = !c
    }
}
return c
Camilo Ortegón
  • 3,414
  • 3
  • 26
  • 34
  • 1
    can you explain a bit? – SmAc May 05 '17 at 20:33
  • Sorry I got it from otherwhere, I'm using it in a Tangram app that I did for iOS, but I'm not sure if it works only on convex polygons. I only know that is a DP algorithm which means that it skips many steps. – Camilo Ortegón May 05 '17 at 21:09