Point In Polygon

Point In Polygon Problem

The problem of finding whether a point is inside or outside a polygon of n-sides.

Angle method:

Th idea is, the angle made by a point A with the vertices of a polygon would be 360 degrees if the point is inside the polygon. The below picture shows the concept

The code implementation below:

            '------ Angle Method 
            Dim InnerAng As Double
            For Each Pt In _Grid_Points '--- Go through all Grid Points
                InnerAng = 0
                For i = 0 To (SPoints.Count - 1) Step +1 '--- Go through all lines for intersection
                    Scalar_Angle_Summation(Pt.Sc_X, Pt.Sc_Y, _SPoints(i).Sc_X, _SPoints(i).Sc_Y, _EPoints(i).Sc_X, _EPoints(i).Sc_Y, InnerAng)
                Next
                '--- check whether the intersection count is odd
                If Math.Round(Math.Abs(InnerAng), 3) <> 0 Then
                    Inside_Points.Add(Pt)
                End If
            Next
            Private Sub Scalar_Angle_Summation(ByVal APcoord1 As Integer, ByVal APcoord2 As Integer, _
                                     ByVal SPcoord1 As Integer, ByVal SPcoord2 As Integer, _
                                     ByVal EPcoord1 As Integer, ByVal EPcoord2 As Integer, _
                                     ByRef InnerAngle As Double)
            Dim Local_IProduct As Double
            Dim v1_1, v1_2 As Double
            Dim v2_1, v2_2 As Double
            Dim Normalize1, Normalize2 As Double
            v1_1 = SPcoord1 - APcoord1
            v1_2 = SPcoord2 - APcoord2
            Normalize1 = Math.Sqrt((v1_1 ^ 2) + (v1_2 ^ 2))
            v1_1 = v1_1 / Normalize1
            v1_2 = v1_2 / Normalize1
            v2_1 = EPcoord1 - APcoord1
            v2_2 = EPcoord2 - APcoord2
            Normalize2 = Math.Sqrt((v2_1 ^ 2) + (v2_2 ^ 2))
            v2_1 = v2_1 / Normalize2
            v2_2 = v2_2 / Normalize2
            Dim t_sin, t_cos As Double
            t_sin = v1_1 * v2_2 - v2_1 * v1_2
            t_cos = v1_1 * v2_1 + v1_2 * v2_2
            Local_IProduct = Math.Atan2(t_sin, t_cos) * (180 / Math.PI)
            InnerAngle = InnerAngle + Local_IProduct
        End Sub

But this method is not correct. Assume a polygon with 5 sides which makes a star like shape now the question is whether the points in the center highlighted area of the polygon is inside or outside of the polygon. Actually it is on the outside.

Internet gifted me the algorithm below

Ray Casting method:

Ray casting is also a simple algorithm which doesn't seems to have any limitations. From each point cast a ray (any direction should suffice but lets cast both horizontal and vertical ray). Now the number of intersection made will be odd if the point is inside and even if the point is on the outside. The picture below shows the concept.

The core idea behind the implementation is to find whether a ray intersect a line. Consider a line which is defined by 2 points S(Sx,Sy) & E(Ex,Ey) and we are casting an horizontal ray from A(Ax,Ay) to infinity (let say (1000,Ay)

So A can be defined as

Horizontal Ray Parameterization from point A(Ax,Ay):

X = Ax(1-t) + 1000t

Y = Ay

The line SE Parameterization

X = Sx(1-s) + Ex(s)

Y = Sy(1-s) + Ey(s)

Now we know y= Ay which is a constant because our ray is horizontal

so, s = (Ay – Sy)/(Ey – Sy)

if s< 0 or s>1 then the intersection point is outside the bounds of the line segment (so no intersection on the segment)

if not (ie., 0<s<1)

find X = Sx(1-s) + Ex(s) using our s

Now find t = (X – Ax)/(1000 - Sx)

if t< 0 or t>1 then the intersection point is outside the bounds of the ray (so no intersection on the segment)other wise there is intersection

The code implementation below:

          '------ Ray Casting Method 
            Dim Horiz_Ray_Intersect As Integer
            Dim Verti_Ray_Intersect As Integer
            For Each Pt In _Grid_Points '--- Go through all Grid Points
                Horiz_Ray_Intersect = 0
                Verti_Ray_Intersect = 0
                For i = 0 To SPoints.Count - 1 Step +1 '--- Go through all lines for intersection
                    Ray_intersection_Algorithm(Pt.Sc_X, Pt.Sc_Y, _SPoints(i).Sc_X, _SPoints(i).Sc_Y, _EPoints(i).Sc_X, _EPoints(i).Sc_Y, Horiz_Ray_Intersect)
                    Ray_intersection_Algorithm(Pt.Sc_Y, Pt.Sc_X, _SPoints(i).Sc_Y, _SPoints(i).Sc_X, _EPoints(i).Sc_Y, _EPoints(i).Sc_X, Verti_Ray_Intersect)
                Next
                '--- check whether the intersection count is odd
                If Horiz_Ray_Intersect Mod 2 <> 0 And Verti_Ray_Intersect Mod 2 <> 0 Then
                    Inside_Points.Add(Pt)
                End If
            Next
          Private Sub Ray_intersection_Algorithm(ByVal APcoord1 As Integer, ByVal APcoord2 As Integer, _
                                               ByVal SPcoord1 As Integer, ByVal SPcoord2 As Integer, _
                                               ByVal EPcoord1 As Integer, ByVal EPcoord2 As Integer, _
                                               ByRef Intersection_Count As Integer)
            '------ Polygon sides
            '------ Pcoord1 = SPcoord1(1-s) + EPcoord1(s)
            '------ Pcoord2 = SPcoord2(1-s) + EPcoord2(s)
            '------ our Ray Coordinates
            '------ Rcoord1 = APcoord1(1-s) + 1000(s)
            '------ Rcoord2 = APcoord2
            If EPcoord2 - EPcoord1 = 0 Then 'To avoid 0/0 error
                Exit Sub
            End If
            Dim s As Double
            s = (APcoord2 - SPcoord2) / (EPcoord2 - SPcoord2)
            '--- continue if fall between the line bounds (ie.,0<s<1)
            If s > 0 And s < 1 Then
                '---- Second coordinate, remember the first coordinate intersect at AP
                Dim Icoord1 As Double
                Icoord1 = (SPcoord1 * (1 - s)) + (EPcoord1 * (s))
                '------ A Ray to infinity (let say 1000)
                '------ coord = AP(1-t) + 1000t
                Dim t As Double
                t = (Icoord1 - APcoord1) / (1000 - APcoord1)
                If t > 0 And t < 1 Then
                    Intersection_Count = Intersection_Count + 1
                End If
            End If
        End Sub

Suppose we have a cubic line, my solution is to segment it into multiple lines, which is better a better method than solving the cubic equation.