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.