Efektywność

Teoretycznie półprostą można przeprowadzić w dowolnym kierunku, jednak ze względu na uciążliwe testowanie przecięć w ogólnych warunkach warto uprościć nieco zadanie. Najlepiej gdy półprosta jest pozioma lub pionowa, wtedy problem praktycznie upraszcza się do jednego wymiaru.

Odcinki A oraz B są odrzucane od razu jako nieprzecinające półprostej, podobnie odcinek C jest natychmiast klasyfikowany jako przecinający półprostą. Dopiero odcinek D wymaga bardziej zaawansowanych testów.

Poniższy kod powinien rozwiać wszelkie wątpliwości.

# (x,y) - współrzędna punktu P
# (x0,y0), (x1,y1) - początek i koniec odcinka

def edge_cross_halfline(x,y, x0,y0, x1,y1):

        # odcinek poziomy
        if y1==y0:
                return 0

        # odcinek A
        if y0 > y xor y1 > 0: # (y0>y and y1>y) or (y0<y and y1<y)
                return 0

        # odcinek B
        if x0 < x and x1 < x:
                return 0

        # odcinek C
        if x0 >= x and x1 >= x:
                return 1

        # odcinek D
        # 1. liczymy punkt przecięcia (xp, yp)
        t  = (y-y0)/float(y1-y0)
        xp = x0 + t*(x1-x0)      # yp=y

        # 2. klasyfikujemy
        if xp >= x:
                return 1
        else:
                return 0

Warunki brzegowe

Testując przynależność punktu do wielokąta operacje dokonywane są na poziomie krawędzi (odcinków). Gdy półprosta przechodzi przez wierzchołek, to policzone zostaną dwa przecięcia, bowiem punkt ten należy do dwóch krawędzi z niego wychodzących.

Uwaga: zakłada się, że punkt nie należy do brzegu wielokąta!

Tak więc punkty P1 (2 mod 2 = 0) i 2 (3 mod 2 = 1) zostaną prawidłowo sklasyfikowane, okaże się natomiast, że punkt P3 (2 mod 2 = 0) leży poza wielokątem.

Rozwiązaniem tego problemu jest przejęcie, że punkt należy tylko do jednej z krawędzi; przy czym krawędź do której nie należy znajduje się powyżej półprostej.

I tak dla punktu P1 zostaną policzone dwa przecięcia, dla punktu P2 oraz P3 zliczone zostanie jedno przecięcie.

Poniżej wymagane zmiany w funkcji edge_cross_halfline.

# sortowanie: y0 < y1
if y0 > y1:
        x0,y0, x1,y1 = x1,y1, x0,y0

# odcinek A
if y0 < y and y1 <= y: # odcinek powyżej
        return 0

if y0 > y and y1 > y:  # odcinek poniżej półprostej
        return 0
# ...

Funkcja point_in_polygon

Funkcja przejmuje wielokąt opisany listą krawędzi i wierzchołków. Każda krawędź to 2-elementowa krotka zawierająca indeksy do tablicy wierzchołków. Każdy wierzchołek to również krotka dwuelementowa która zawiera współrzędne punktu.

def point_in_polygon(x,y, edges_list, vertex_list):
        count = 0

        for e in edges_list:
                a,b = e

                x0, y0 = vertex_list[a]
                x1, y1 = vertex_list[b]

                count += edge_cross_halfline(x,y, x0,y0, x1,y1)

        return count%2 == 1

Punkt wewnątrz trójkąta

Trójkąt jest oczywiście figurą wypukłą, dlatego można z powodzeniem zastosować przedstawioną wyżej funkcję. Można również wprost odwołać się do tych ledwie trzech prostych AB, BC, CA:

def point_in_triangle(x,y, xa,ya, xb,yb, xc,yc):

        def side(x,y, x0,y0, x1,y1):
                dx = x1-x0 # przesunięcie odcinka
                dy = y1-y0 # o wektor (-x0,-y0)
                 x =  x-x0
                 y =  y-y0

                # f(x,y) = ax+by = dy*x - dx*y
                return dy*x - dx*y > 0

        sideAB = side(x,y, xa,ya, xb,yb)
        sideBC = side(x,y, xb,yb, xc,yc)
        sideCA = side(x,y, xc,yc, xa,ya)

        return (sideAB == sideBC) and (sideAB == sideCA) and (sideBC == sideCA)
        # albo:
        # not ((sideAB xor sideBC) or (sideAB xor sideCA))

W przypadku używania zmiennych całkowitych w funkcji side rozsądnym zdaje się użycie MMX-owego rozkazu pmaddwd.

def point_in_triangle(x,y, x0,y0, x1,y1, x2,y2):

        x  -= x2 # przesunięcie wszystkich punktów
        x0 -= x2 # o wektor (-x2,-y2)
        x1 -= x2

        y  -= y2
        y0 -= y2
        y1 -= y2

        # rozwiązanie układu równań
        # x = t0 x0 + t1 x1
        # y = t0 y0 + t1 y1

        # wyznacznik główny
        det = x0*y1 - y0*x1

        # ewentualnie, jeśli współrzędne mogą być podane
        # jako całkowite:
        det = float(x0*y1 - y0*x1)

        if  abs(det) < epsilon: # brak rozwiązań
                return False

        t0 = ( x*y1 -  y*x1)/det
        t1 = (x0*y  - y0*x )/det

        if t0 < 0 or t1 < 0 or t0+t1 > 1:
                return False
        else:
                return True

# punkty to trzyelementowe krotki (x,y,z)
def point_in_triangle3D(p, p0,p1,p2)
        x,y,z = p
        x0,y0,z0 = p0
        x1,y1,z1 = p1
        x2,y2,z2 = p2

        dx = max(x0,x1,x2) - min(x0,x1,x2) # wysokości
        dy = max(y0,y1,y2) - min(y0,y1,y2) # prostopadłościanu
        dz = max(z0,z1,z2) - min(z0,z1,z2)

        min_height = min(dx,dy,dz)

        if dx == min_height: # płaszczyzna YZ: pomijamy wsp. x
                return point_in_triangle(y,z, y0,z0, y1,z1, y2,z2)

        if dy == min_height: # płaszczyzna ZX: pomijamy wsp. y
                return point_in_triangle(z,x, x0,z0, x1,z1, x2,z2)

        if dz == min_height: # płaszczyzna XY: pomijamy wsp. z
                return point_in_triangle(x,y, x0,y0, x1,y1, x2,y2)