多角形の内外判定

多角形の内外判定をする必要があって、pythonで作ってみた。

参考にしたのがmatlabのinpolygon関数（http://www.mathworks.co.jp/jp/help/matlab/ref/inpolygon.html)

だけど、諸事情でend pointが閉じていないものにした。

使ったのはCrossing Number Algorithmという一番単純なもので、多角形の辺が、ある点より上にある数が奇数であれば、その点は多角形の内側にあるというもの。

import matplotlib.pyplot as plt
import numpy as np

def main():
    n = 2000
    x = np.ndarray(n, dtype=float)
    y = np.ndarray(n, dtype=float)
    inside = np.ndarray(n, dtype=bool)
        
    x = (np.random.rand(n) - 0.5)*2.5
    y = (np.random.rand(n) - 0.5)*2.5
    inside[:] = False
        
    x1, y1 = [], []
    
    # define a polygon
    for i in range(5):
        x1.append(np.cos(i*2.*np.pi/5.))
        y1.append(np.sin(i*2.*np.pi/5.))
    
    # check the points inside the polygon
    for i, (sx, sy) in enumerate(zip(x, y)):
        inside[i] = inpolygon(sx, sy, x1, y1)
    
    # plot results
    x1.extend([x1[0]])
    y1.extend([y1[0]])
    plt.plot(x[inside], y[inside], 'ro')
    plt.plot(x[inside==False], y[inside==False], 'ko')
    plt.plot(x1, y1)
    plt.savefig("inpoly.png")
        
def inpolygon(sx, sy, x, y):
    ''' 
    x[:], y[:]: polygon
    sx, sy: point
    '''     
    np = len(x)
    inside = False
    for i1 in range(np): 
        i2 = (i1+1)%np
        if min(x[i1], x[i2]) < sx < max(x[i1], x[i2]):
            #a = (y[i2]-y[i1])/(x[i2]-x[i1])
            #b = y[i1] - a*x[i1]
            #dy = a*sx+b - sy
            #if dy >= 0:
            if (y[i1] + (y[i2]-y[i1])/(x[i2]-x[i1])*(sx-x[i1]) - sy) > 0:
                inside = not inside

    return inside

if __name__ == "__main__":
    main()

ためしてみると、うまくいった。こんな単純な計算でできるんだなあ。