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.

Google Sites

Report abuse