PDA

View Full Version : Polygon's



mcmichael13
04-21-2009, 06:21 AM
Hey guys, first post here, looking to learn a lot from the site :)

First question:

So I have several sets of points which I am going to use to create several polygon "bins". I was wondering how to load these point sets to define my polygons in vba? After I get the polygon's loaded, I will use "Point-In-Polygon" techniques to test which bin a specified point (in excel) is located in. Right now they are using the summation of angles to verify the existence of a point in a particular polygon, but the code has bugs and is not correct all the time. So I figured that using the "ray" technique would be more accurate. That being said, I just need some help in learning syntax, I've been reading a lot and using the help extensively. Thanks in advance.

Rob

mcmichael13
04-21-2009, 10:48 AM
anyone? please?

mdmackillop
04-21-2009, 12:39 PM
Welcome to VBAX
Never come across these things. Post an example with your code and I'll see if I can follow things. Add any comments etc that might help to your workbook. Use Manage Attachments in the GoAdvanced reply section.
What version of Excel? (in case it makes a difference)

mcmichael13
04-21-2009, 12:51 PM
This is the code as it stands now, looks like its a modified code:


Option Explicit
'****************************************************************
' Routines as posted on Planet Source Code
'
'
'This is a very useful class to store polygon's vertices and
'check if the polygon is convex or not and if a point is inside
'the polygon or not
'By Raul Fragoso on 14-March-2002
'****************************************************************
'Declares PI
Public Const PI As Double = 3.14159265358979


Public Function udfPointInPolygon(ByVal x As Range, ByVal y As Range, Vertices As Range) As Boolean

Dim intIndex As Integer
Dim sngTotalAngle As Single
Dim intVerticeCount As Integer

intVerticeCount = Vertices.Rows.Count

' Get the angle between the point and the
' first and last vertices.
sngTotalAngle = GetAngle( _
Vertices.Cells(intVerticeCount, 1).Value, _
Vertices.Cells(intVerticeCount, 2).Value, _
x, y, _
Vertices.Cells(1, 1).Value, _
Vertices.Cells(1, 2).Value)

' Add the angles from the point to each other
' pair of vertices.
For intIndex = 1 To intVerticeCount - 1
sngTotalAngle = sngTotalAngle + GetAngle(Vertices.Cells(intIndex, 1).Value, _
Vertices.Cells(intIndex, 2).Value, _
x, y, _
Vertices.Cells(intIndex + 1, 1).Value, _
Vertices.Cells(intIndex + 1, 2).Value)
Next

' The total angle should be 2 * PI or -2 * PI if
' the point is in the polygon and close to zero
' if the point is outside the polygon.
udfPointInPolygon = (Abs(sngTotalAngle) > 4)
End Function
' Return the angle ABC.
' Return a value between PI and -PI.
' Note that the value is the opposite of what you might
' expect because Y coordinates increase downward.
Private Function GetAngle(ByVal Ax As Single, ByVal Ay As Single, ByVal Bx As Single, ByVal By As Single, ByVal Cx As Single, ByVal Cy As Single) As Single
Dim dot_product As Single
Dim cross_product As Single

' Get the dot product and cross product.
dot_product = DotProduct(Ax, Ay, Bx, By, Cx, Cy)
cross_product = CrossProductLength(Ax, Ay, Bx, By, Cx, Cy)

' Calculate the angle.
GetAngle = ATan2(cross_product, dot_product)
End Function
' Return the dot product AB · BC.
' Note that AB · BC = |AB| * |BC| * Cos(theta).
Private Function DotProduct( _
ByVal Ax As Single, ByVal Ay As Single, _
ByVal Bx As Single, ByVal By As Single, _
ByVal Cx As Single, ByVal Cy As Single _
) As Single

Dim BAx As Single
Dim BAy As Single
Dim BCx As Single
Dim BCy As Single

' Get the vectors' coordinates.
BAx = Ax - Bx
BAy = Ay - By
BCx = Cx - Bx
BCy = Cy - By

' Calculate the dot product.
DotProduct = BAx * BCx + BAy * BCy
End Function


' Return the angle with tangent opp/hyp. The returned
' value is between PI and -PI.
Private Function ATan2(ByVal opp As Single, ByVal adj As Single) As Single

Dim sngAngle As Single

' Get the basic angle.
If Abs(adj) < 0.0001 Then
sngAngle = PI / 2
Else
sngAngle = Abs(Atn(opp / adj))
End If

' See if we are in quadrant 2 or 3.
If adj < 0 Then
' angle > PI/2 or angle < -PI/2.
sngAngle = PI - sngAngle
End If

' See if we are in quadrant 3 or 4.
If opp < 0 Then
sngAngle = -sngAngle
End If

' Return the result.
ATan2 = sngAngle

End Function
' Return the cross product AB x BC.
' The cross product is a vector perpendicular to AB
' and BC having length |AB| * |BC| * Sin(theta) and
' with direction given by the right-hand rule.
' For two vectors in the X-Y plane, the result is a
' vector with X and Y components 0 so the Z component
' gives the vector's length and direction.
Private Function CrossProductLength( _
ByVal Ax As Single, ByVal Ay As Single, _
ByVal Bx As Single, ByVal By As Single, _
ByVal Cx As Single, ByVal Cy As Single _
) As Single

Dim sngBAx As Single
Dim sngBAy As Single
Dim sngBCx As Single
Dim sngBCy As Single

' Get the vectors' coordinates.
sngBAx = Ax - Bx
sngBAy = Ay - By
sngBCx = Cx - Bx
sngBCy = Cy - By

' Calculate the Z coordinate of the cross product.
CrossProductLength = sngBAx * sngBCy - sngBAy * sngBCx

End Function


This is the code I would like to implement, I'm just very new and don't know how to do it quite yet:


Private Type POINTAPI
X As Long
Y As Long
End Type

Private Function PtInPoly(Poly() As POINTAPI, _
ByVal Xray As Long, _
ByVal YofRay As Long) As Boolean
Dim X As Long
Dim Yintersect As Long
Dim PolyCount As Long
Dim NumSidesCrossed As Long
PolyCount = 1 + UBound(Poly) - LBound(Poly)
For X = LBound(Poly) To UBound(Poly)
If Poly(X).X > Xray Xor Poly((X + 1) Mod PolyCount).X > Xray Then
Yintersect = Y_at_X_Ray(Xray, Poly(X), Poly((X + 1) Mod
PolyCount))
If Yintersect > YofRay Then
NumSidesCrossed = NumSidesCrossed + 1
End If
End If
Next
If NumSidesCrossed Mod 2 Then PtInPoly = True
End Function

Private Function Y_at_X_Ray(ByVal Xray As Single, _
p1 As POINTAPI, _
p2 As POINTAPI) As Single
Dim m As Single
Dim b As Single
m = (p2.Y - p1.Y) / (p2.X - p1.X)
b = (p1.Y * p2.X - p1.X * p2.Y) / (p2.X - p1.X)
Y_at_X_Ray = m * Xray + b
End Function


I will take this code and loop it over each of the bins to classify which bin it is located in. Let me know if you need any other data or information.

The problem I am running into is that I have never used API at all, so when I am running an If statement for:

If PtInPoly(u, v, [CB]) Then
Bin = "CB"

This is where u is the x coordinate, v is the y, and CB is the bin range... but the debugger is asking for an array for the [CB].


If I am going about this all wrong, please let me know of a better way. As I'm at a dead end here.

Thanks,
Rob

mdmackillop
04-21-2009, 12:55 PM
Can you post a workbook containing data to run your code?

mcmichael13
04-21-2009, 01:07 PM
The information in the workbook is kind of proprietary, don't want to send it out for all to see, but perhaps I could pm it to you? Its a rather large data file as well.

Andy Pope
04-22-2009, 02:39 AM
This maybe of help.
http://www.andypope.info/fun/pointinarea.htm

mcmichael13
04-22-2009, 05:39 AM
hey andy thanks for the reply, problem is though... i'm already using that code and for whatever reason, the calculations aren't right. thats why i wanted to implement the other code to use the "ray-casting" technique.

here is a screen shot of the graph i have, these are the "bins", to which when i plot i point... the excel spreadsheet spits out on another column what bin the point is in... however this code is just blatantly wrong sometimes and gives me the wrong bin.

Thanks,

Rob

Andy Pope
04-22-2009, 05:52 AM
Can you email your workbook or create non confidential example of the problem.

andy AT andypope DOT info

mcmichael13
04-22-2009, 06:13 AM
email sent, thanks a lot

Andy Pope
04-22-2009, 07:00 AM
The problem is to do with the very small values.

If you multiple the xy values by a factor of 1000 then the function returns a much more accurate result.

mcmichael13
04-22-2009, 07:03 AM
i also just was messing with the line:

udfPointInPolygon = (Abs(sngTotalAngle) > .0001)

and found if i replace it with:

udfPointInPolygon = (Abs(sngTotalAngle) > 10)

it is much more accurate.

Thanks so much, i'll keep testing to make sure i'm getting accurate results all over.

mcmichael13
04-22-2009, 07:08 AM
this may be a dumb question, but how would i multiply the x and y values by 1000 in the code? i'm drawing a blank haha

edit: it seems that with editing that line of code... the M and L bins no longer work at all, so i would like to see what you did, see if that works

Andy Pope
04-22-2009, 07:30 AM
Here is one way. I have added an optional argument to the function, Factor.


Public Function udf_InPolygon(ByVal X As Range, ByVal Y As Range, Vertices As Range, Optional Factor As Single = 1) As Boolean

Dim intIndex As Integer
Dim sngTotalAngle As Single
Dim intVerticeCount As Integer

intVerticeCount = Vertices.Rows.Count

' Get the angle between the point and the
' first and last vertices.
sngTotalAngle = GetAngle( _
Vertices.Cells(intVerticeCount, 1).Value * Factor, _
Vertices.Cells(intVerticeCount, 2).Value * Factor, _
X * Factor, Y * Factor, _
Vertices.Cells(1, 1).Value * Factor, _
Vertices.Cells(1, 2).Value * Factor)

' Add the angles from the point to each other
' pair of vertices.
For intIndex = 1 To intVerticeCount - 1
sngTotalAngle = sngTotalAngle + GetAngle(Vertices.Cells(intIndex, 1).Value * Factor, _
Vertices.Cells(intIndex, 2).Value * Factor, _
X * Factor, Y * Factor, _
Vertices.Cells(intIndex + 1, 1).Value * Factor, _
Vertices.Cells(intIndex + 1, 2).Value * Factor)
Next

' The total angle should be 2 * PI or -2 * PI if
' the point is in the polygon and close to zero
' if the point is outside the polygon.
udf_InPolygon = (Abs(sngTotalAngle) > 0.000001)
End Function


and the formula would look like, where the factor is 100

udf_inpolygon(B9,C9,E18:F37,100)

mcmichael13
04-22-2009, 09:02 AM
Hey andy, thanks so much... i got the code working great with awesome accuracy, thanks for all your help.

Rob