Thiessen Polygon

2018年12月02日

Thiessen polygon

In mathematics, a Thiessen polygon(Voronoi diagram) is a way to divide a space into a number of regions. A set of control points (called seeds, sites, or stations) is specified beforehand and for each seed points there will be a corresponding closed region consisting of all points closer to that seed point than to any other. The regions are called Voronoi cells. It is dual to the Delaunay triangulation.

It is named after Georgy Voronoy, and is also called a Voronoi tessellation, a Voronoi decomposition, a Voronoi partition, or a Dirichlet tessellation. Voronoi diagrams can be found in many fields in science, enginerring , technology, art, even in medicine ,and they have found numerous theoretical algorithms and practical applications.[

The simplest case

In the simplest and most familiar case (shown in Figure 1 and Figure 2), are given a finite set of points {p1,p2,…,pi,….pn} in the Euclidean plane. In this case, each seed point pi is simply a point and its corresponding Voronoi region R, consisting of every point whose distance to pi is less than or equal to its distance to any other site point. Each such cell is obtained from the intersection of half-spaces, and hence it is a convex polygon. The segments of the Voronoi diagram are all the points in the plane that are equidistant to the two nearest sites. The Voronoi vertices (nodes) are the points equidistant to three (or more) sites.

Figure 3

Figure 4

Figure 5

Build Thiessen Polygons

Thiessen (Voronoi) polygons define individual areas of influence around each of a set of points. Thiessen polygons are polygons whose boundaries define the area that is closest to each point relative to all other points. They are mathematically defined by the perpendicular bisectors of the lines between all points

Inputs:

2D Points data

Outputs:

New polygon

Notes :

The process of create Thiessen polygongs:

(1)Collects and input the points data (remove duplicate points)

(2)Create a Convex Hull (map boundary board)

(3)Creates a TIN structure(Trianglation)

Figure 1

(4)Find out external points and internal points.

Figure 2

(5)Create external Thiessen polygons and internal Thiessen polygons

Figure 3

(6)Computing the areas of each Thiessen polygon if need.

Figure 4

Public Function ThiessenPolygonInternal(ByVal canvas As PictureBox, ByVal Gr As Graphics, ByVal mPen As Pen, ByVal mBrush As Brush, ByVal IptGiven As Integer, Optional ByRef isFill As Boolean = True) As String

Dim j, i, K As Short

Dim indpts(2) As Short

Dim CollsTriInd As New Collection

Dim CollsTriIndNew As New Collection

Dim stArray() As String

Dim indTriangles() As Short

Dim nCount As Short

Dim iptA As String

Dim iptB As String

Dim stSideTest, stSideA, stSideBeg As String

Dim stTest(3) As String

Dim areaTpt As Single

CollsTriInd.Clear()

CollsTriIndNew.Clear()

For i = 1 To nTriangles

ReDim Preserve Triangles(i).ipts(0 To 2)

If Triangles(i).ipts(0) <> IptGiven And Triangles(i).ipts(1) <> IptGiven And Triangles(i).ipts(2) <> IptGiven Then GoTo 100 '

If Triangles(i).ipts(0) = IptGiven Then

CollsTriInd.Add(Trim(i & "," & Triangles(i).ipts(1) & "," & Triangles(i).ipts(2)))

GoTo 100

End If

If Triangles(i).ipts(1) = IptGiven Then

CollsTriInd.Add(Trim(i & "," & Triangles(i).ipts(2) & "," & Triangles(i).ipts(0)))

GoTo 100

End If

If Triangles(i).ipts(2) = IptGiven Then

CollsTriInd.Add(Trim(i & "," & Triangles(i).ipts(0) & "," & Triangles(i).ipts(1)))

GoTo 100

End If

100: Next i

nCount = 0

stArray = Split(CollsTriInd.Item(1), ",", -1, CompareMethod.Text)

ReDim Preserve indTriangles(nCount)

indTriangles(nCount) = Val(stArray(0))

iptA = stArray(1)

iptB = stArray(2)

stSideBeg = Trim(IptGiven & "," & iptA)

stSideTest = Trim(IptGiven & "," & iptB)

CollsTriIndNew.Add(CollsTriInd.Item(1))

CollsTriInd.Remove((1))

nCount = nCount + 1

On Error Resume Next

200: If CollsTriInd.Count >= 1 Then

For i = 1 To CollsTriInd.Count

If stSideTest = stSideBeg Then Exit For

stArray = Split(CollsTriInd.Item(i), ",", -1, CompareMethod.Text)

iptA = stArray(1)

iptB = stArray(2)

stTest(0) = Trim(IptGiven & "," & iptA)

stTest(1) = Trim(iptA & "," & IptGiven)

stTest(2) = Trim(IptGiven & "," & iptB)

stTest(3) = Trim(iptB & "," & IptGiven)

If stSideTest = stTest(0) Or stSideTest = stTest(1) Then

CollsTriIndNew.Add(CollsTriInd.Item(i))

CollsTriInd.Remove((i))

ReDim Preserve indTriangles(nCount)

indTriangles(nCount) = Val(stArray(0))

stSideTest = stTest(2)

nCount = nCount + 1

GoTo 200

End If

If stSideTest = stTest(2) Or stSideTest = stTest(3) Then

CollsTriIndNew.Add(CollsTriInd.Item(i))

CollsTriInd.Remove((i))

ReDim Preserve indTriangles(nCount)

indTriangles(nCount) = Val(stArray(0))

stSideTest = stTest(0)

nCount = nCount + 1

GoTo 200

End If

Next i

End If

'?????????????????

Dim PtsCirCen() As PointF ', Ptmid As ptXy

Dim nptPtCen As Integer

nptPtCen = 0

For i = 1 To CollsTriIndNew.Count()

stArray = Split(CollsTriIndNew.Item(i), ",", -1, CompareMethod.Text)

ReDim Preserve PtsCirCen(nptPtCen)

PtsCirCen(nptPtCen) = CircumCenter(Triangles(Val(stArray(0))).ipts(0), Triangles(Val(stArray(0))).ipts(1), Triangles(Val(stArray(0))).ipts(2))

nptPtCen = nptPtCen + 1

Next i

ReDim Preserve PtsCirCen(nptPtCen)

PtsCirCen(nptPtCen) = PtsCirCen(0)

'求Polygon是否與圖界有交點

Dim tPLinexy As PLineXy

Dim ptsi() As PointF

Dim iPoints() As PointF

Dim nAns As Short

Dim ptsIntF() As PointF

For i = 0 To UBound(PtsCirCen)

ReDim Preserve tPLinexy.pts(i)

tPLinexy.pts(i) = PtsCirCen(i)

Next i

nAns = PolylinePolylineIntReal(tPLinexy, MapBDLPlXy, iPoints)

'------------------------

'圖界交點

Dim ptCenTpt As PointF

Dim nPtAdd As Short, nptIntF As Integer

Dim ptsIntFB() As PointF

Dim LineTptA, LineTptB As LineXy

ReDim Preserve LineTptA.pts(1), LineTptB.pts(1)

If nAns <> -1 Then

Erase ptsIntF

ReDim Preserve iPoints(UBound(iPoints))

nptIntF = 0

For j = 0 To UBound(iPoints)

ReDim Preserve ptsIntF(nptIntF)

ptsIntF(nptIntF) = iPoints(j)

nptIntF = nptIntF + 1

Next j

'圖界交點是否有轉折點

LineTptA.pts(0) = ptsIntF(0)

LineTptA.pts(1) = ptsIntF(UBound(ptsIntF))

nPtAdd = 0

For i = 0 To UBound(MapBDLPlXy.pts) - 1

LineTptB.pts(0).X = ptsXyz(IptGiven).X

LineTptB.pts(0).Y = ptsXyz(IptGiven).Y

LineTptB.pts(1) = MapBDLPlXy.pts(i)

Erase iPoints

nAns = LineLineIntReal(LineTptA, LineTptB, iPoints)

If nAns <> -1 Then

ReDim Preserve iPoints(UBound(iPoints))

ReDim Preserve ptsIntFB(nPtAdd)

ptsIntFB(nPtAdd) = MapBDLPlXy.pts(i)

nPtAdd = nPtAdd + 1

End If

Next i

End If

nptIntF = nptIntF - 1

nPtAdd = nPtAdd - 1

'綜合所有控制點

'(1)圖界內所有外心

Dim angTpt As Single

Erase ptsi

nCount = 0

If nptPtCen >= 0 Then

For j = 0 To nptPtCen

If isPtInRegion(MapBDLPlXy.pts, PtsCirCen(j), angTpt) Then

ReDim Preserve ptsi(nCount)

ptsi(nCount) = PtsCirCen(j)

nCount = nCount + 1

End If

Next j

End If

'(2)圖界點

If nptIntF >= 0 Then

For j = 0 To UBound(ptsIntF)

ReDim Preserve ptsi(nCount)

ptsi(nCount) = ptsIntF(j)

nCount = nCount + 1

Next j

End If

'(3)轉折點

If nPtAdd >= 0 Then

For j = 0 To nPtAdd

ReDim Preserve ptsi(nCount)

ptsi(nCount) = ptsIntFB(j)

nCount = nCount + 1

Next j

End If

' 繞測站依角度大小排列

ptCenTpt.X = ptsXyz(IptGiven).X

ptCenTpt.Y = ptsXyz(IptGiven).Y

areaTpt = System.Math.Abs(Area_afterReorder(ptCenTpt, ptsi))

Erase PtsCirCen

For j = 0 To UBound(ptsi)

ReDim Preserve PtsCirCen(j)

PtsCirCen(j) = ptsi(j) '以PtsCirCen(j) 代替ptsI(j)

Next j

If isFill Then

mBrush = New SolidBrush(RandColors(IptGiven Mod 23))

Call FillRegionXY(Gr, mBrush, PtsCirCen)

Else

Gr.DrawPolygon(New Pen(Color.Blue, 2), PtsCirCen)

End If

ReDim Preserve ThiessenXyzs(nThiessen)

With ThiessenXyzs(nThiessen)

.indPt = IptGiven

.coverArea = areaTpt

.ptStaXyz.X = ptsXyz(IptGiven).X

.ptStaXyz.Y = ptsXyz(IptGiven).Y

.ptStaXyz.z = ptsXyz(IptGiven).z

For i = 0 To UBound(PtsCirCen)

ReDim Preserve ThiessenXyzs(nThiessen).PLgon.pts(i)

.PLgon.pts(i) = PtsCirCen(i)

Next i

Frmdocument.List1.Items.Add("sta_ " & IptGiven & " ;area= " & areaTpt)

sumAreaOfAll = sumAreaOfAll + areaTpt

nThiessen = nThiessen + 1

End With

canvas.Refresh()

End Function

Public Function ThiessenPolygonExternal(ByVal canvas As PictureBox, ByVal Gr As Graphics, ByVal mPen As Pen, ByVal mBrush As Brush, ByVal IptGiven As Integer, Optional ByRef isFill As Boolean = True) As String

Dim j, i As Short

Dim indpts(2) As Short

Dim CollsTriInd As New Collection

Dim CollsTriIndNew As New Collection

Dim stArray() As String

Dim nCount As Short

Dim stTest(3) As String

Dim LineOutA, LineOutB As LineXy

ReDim LineOutA.pts(0 To 1), LineOutB.pts(0 To 1)

Dim ptsCenKey(1) As PointF

Dim PtMid As PointF

Dim isPtOnBDL As Boolean

Dim IndKeyTriAng As Short

Dim LineTptA, LineTpt, LineTptB As LineXy

ReDim Preserve LineTptA.pts(0 To 1), LineTpt.pts(0 To 1), LineTptB.pts(0 To 1)

Dim PtsCirCen() As PointF

Dim nAns As Short

Dim iPoints() As PointF

Dim ptsIntF() As PointF

Dim areaTpt As Single

Dim Lmin, lTest As Single

Dim ptsVert(1) As PointF

Dim IndKeyIpt As Short

Dim indMin As Short

Dim sideStrGiven, sidestrEnd As String

Dim ptsKey() As PointF

'求任意一個外邊

For i = 1 To nTriangles

If Triangles(i).ipts(0) = IptGiven Or Triangles(i).ipts(1) = IptGiven Or Triangles(i).ipts(2) = IptGiven Then

For j = 0 To 2

If Triangles(i).ipts(j) <> IptGiven Then '不等於iptGiven 之點編碼

IndKeyIpt = Triangles(i).ipts(j)

'求最外邊

PtMid.X = (ptsXyz(IptGiven).X + ptsXyz(IndKeyIpt).X) / 2.0#

PtMid.Y = (ptsXyz(IptGiven).Y + ptsXyz(IndKeyIpt).Y) / 2.0#

isPtOnBDL = isPtOnOutBdl(PtMid)

'000000

If isPtOnBDL = True Then

sideStrGiven = Trim(IptGiven & "," & IndKeyIpt)

IndKeyTriAng = i

Exit For

End If

'00000

Next j

End If

Next i

Dim IndTriAngs() As Short

Call From1SideFindAllTriAng(IptGiven, sideStrGiven, IndKeyTriAng, sidestrEnd, IndTriAngs)

nCount = 0

For i = 0 To UBound(IndTriAngs)

ReDim Preserve PtsCirCen(nCount)

PtsCirCen(nCount) = CircumCenter(Triangles(IndTriAngs(i)).ipts(0), Triangles(IndTriAngs(i)).ipts(1), Triangles(IndTriAngs(i)).ipts(2))

nCount = nCount + 1

Next i

'求各外心連線兩外段是否需要延長

Dim ptsCenTmp() As PointF

'求各外心連線與圖邊界之交點

nCount = 0

Dim mPLxy As PLineXy

Dim nIntCen_Map_1 As Short

For i = 0 To UBound(PtsCirCen)

ReDim Preserve mPLxy.pts(nCount)

mPLxy.pts(nCount) = PtsCirCen(i)

nCount = nCount + 1

Next i

nAns = PolylinePolylineIntReal(mPLxy, MapBDLPlXy, iPoints)

If nAns = -1 Then

nIntCen_Map_1 = 0

ptsXyz(IptGiven).indIntMapFore = -99

ptsXyz(IptGiven).indIntMapBack = -99

GoTo 100

End If

'?????????????????????/有交點

Dim nptInt As Integer

If nAns <> -1 Then

ReDim Preserve iPoints(UBound(iPoints))

nptInt = 0

For j = 0 To UBound(iPoints)

ReDim Preserve ptsIntF(nCount)

ptsIntF(nptInt) = iPoints(j)

nptInt = nptInt + 1

Next j

End If

100:

Dim isInRgnLL, isInRgnUU As Boolean

isInRgnLL = isPtInRegion(MapBDLPlXy.pts, PtsCirCen(0))

isInRgnUU = isPtInRegion(MapBDLPlXy.pts, PtsCirCen(UBound(PtsCirCen)))

Dim PtsExtA, PtsExtB As PointF

Dim nptExtA As Integer, nptExtB As Integer

nptExtA = 0 : nptExtB = 0

If isInRgnLL = True Then

stArray = Split(sideStrGiven, ",", -1, CompareMethod.Text)

PtMid.X = (ptsXyz(Val(stArray(0))).X + ptsXyz(Val(stArray(1))).X) / 2

PtMid.Y = (ptsXyz(Val(stArray(0))).Y + ptsXyz(Val(stArray(1))).Y) / 2

LineTptA.pts(0) = ptOnLine_GivenDist(PtMid, PtsCirCen(0), 2000)

LineTptA.pts(1) = ptOnLine_GivenDist(PtsCirCen(0), PtMid, 2000)

nAns = PolyLineLineIntReal(MapBDLPlXy, LineTptA, iPoints)

If nAns <> -1 Then

Lmin = 99999

indMin = -99

ReDim Preserve iPoints(UBound(iPoints))

For j = 0 To UBound(iPoints)

lTest = LenLineXy(PtsCirCen(0), iPoints(j))

If lTest < Lmin Then

Lmin = lTest

indMin = j

End If

Next j

PtsExtA = iPoints(indMin)

nptExtA = nptExtA + 1

End If

'---------------------------

If isInRgnUU = True Then

stArray = Split(sidestrEnd, ",", -1, CompareMethod.Text)

PtMid.X = (ptsXyz(Val(stArray(0))).X + ptsXyz(Val(stArray(1))).X) / 2

PtMid.Y = (ptsXyz(Val(stArray(0))).Y + ptsXyz(Val(stArray(1))).Y) / 2

LineTptB.pts(0) = ptOnLine_GivenDist(PtMid, PtsCirCen(UBound(PtsCirCen)), 2000)

LineTptB.pts(1) = ptOnLine_GivenDist(PtsCirCen(UBound(PtsCirCen)), PtMid, 2000)

nAns = PolyLineLineIntReal(MapBDLPlXy, LineTptB, iPoints)

If nAns <> -1 Then

Lmin = 99999

indMin = -99

ReDim Preserve iPoints(UBound(iPoints))

For j = 0 To UBound(iPoints)

lTest = LenLineXy(PtsCirCen(0), iPoints(j))

If lTest < Lmin Then

Lmin = lTest

indMin = j

End If

Next j

PtsExtB = iPoints(indMin)

nptExtB = nptExtB + 1

End If

nptInt = nptInt - 1

nptExtA = nptExtA - 1

nptExtB = nptExtB - 1

'********************************

'求兩圖界交點間是否夾變化點

nCount = 0

Dim ptsTp() As PointF

Erase ptsTp

If nptExtA >= 0 Then

ReDim Preserve ptsTp(nCount)

ptsTp(nCount) = PtsExtA

nCount = nCount + 1

End If

If nptExtB >= 0 Then

ReDim Preserve ptsTp(nCount)

ptsTp(nCount) = PtsExtB

nCount = nCount + 1

End If

If nptInt >= 0 Then

For i = 0 To nptInt

ReDim Preserve ptsTp(nCount)

ptsTp(nCount) = ptsIntF(i)

nCount = nCount + 1

End If

On Error Resume Next

If nCount >= 1 Then

Dim collPts As New Collection

collPts.Clear()

For i = 0 To nCount

collPts.Add(Trim(Round(ptsTp(i).X, 0) & "," & Round(ptsTp(i).Y, 0)))

Call SortStr_Duplicate(collPts)

If collPts.Count <> 2 Then MsgBox("error in collpts computing")

If collPts.Count = 2 Then

For i = 1 To collPts.Count

stArray = Split(collPts(i), ",", -1, CompareMethod.Text)

LineTptA.pts(i - 1).X = Val(stArray(0))

LineTptA.pts(i - 1).Y = Val(stArray(1))

End If

'----------------------

Dim ptRefer As PointF

ptRefer.X = ptsXyz(IptGiven).X

ptRefer.Y = ptsXyz(IptGiven).Y

Dim collsPtxy As New Collection

Dim collsMile As New Collection

collsPtxy.Clear()

collsMile.Clear()

Call ptXY_MilesOfPLine端(MapBDLPlXy, collsPtxy, collsMile)

ptXY_MilesOfPLineAdd1Pt(MapBDLPlXy, LineTptA.pts(0), "交", collsPtxy, collsMile)

ptXY_MilesOfPLineAdd1Pt(MapBDLPlXy, LineTptA.pts(1), "交", collsPtxy, collsMile)

ptXY_MilesOfPLineAdd1Pt(MapBDLPlXy, ptRefer, "參", collsPtxy, collsMile)

Dim keyMiles交() As Single, keyMiles參As Single, iCase As Integer

Dim nptKeyMiles交As Integer

nptKeyMiles交= 0

For i = 1 To collsPtxy.Count

stArray = Split(collsPtxy(i), ",", -1, CompareMethod.Text)

If stArray(0) = "交" Then

ReDim Preserve keyMiles交(nptKeyMiles交)

keyMiles交(nptKeyMiles交) = collsMile(i)

nptKeyMiles交= nptKeyMiles交+ 1

End If

If stArray(0) = "參" Then keyMiles參= collsMile(i)

Next i

If keyMiles參> keyMiles交(0) And keyMiles參< keyMiles交(1) Then iCase = 0

If keyMiles參> keyMiles交(1) Or keyMiles參< keyMiles交(0) Then iCase = 1

Dim nPtAdd As Integer, ptsAddOnMapBdl() As PointF

nPtAdd = 0

Select Case iCase

Case 0

For i = 1 To collsPtxy.Count

stArray = Split(collsPtxy(i), ",", -1, CompareMethod.Text)

If Val(collsMile(i)) > keyMiles交(0) And Val(collsMile(i)) < keyMiles交(1) And Abs(Val(collsMile(i)) - keyMiles參) > 1.0 Then

ReDim Preserve ptsAddOnMapBdl(nPtAdd)

ptsAddOnMapBdl(nPtAdd).X = Val(stArray(1))

ptsAddOnMapBdl(nPtAdd).Y = Val(stArray(2))

nPtAdd = nPtAdd + 1

End If

Next i

Case 1

For i = 1 To collsPtxy.Count

stArray = Split(collsPtxy(i), ",", -1, CompareMethod.Text)

If (Val(collsMile(i)) < keyMiles交(0) Or Val(collsMile(i)) > keyMiles交(1)) And Abs(Val(collsMile(i)) - keyMiles參) > 1.0 Then

ReDim Preserve ptsAddOnMapBdl(nPtAdd)

ptsAddOnMapBdl(nPtAdd).X = Val(stArray(1))

ptsAddOnMapBdl(nPtAdd).Y = Val(stArray(2))

nPtAdd = nPtAdd + 1

End If

Next i

End Select

nPtAdd = nPtAdd - 1

'綜合所有控制點

'(1)圖界內所有外心

Dim ptsI() As PointF

Erase ptsI

nCount = 0

For i = 0 To UBound(PtsCirCen)

If isPtInRegion(MapBDLPlXy.pts, PtsCirCen(i)) Then

ReDim Preserve ptsI(nCount)

ptsI(nCount) = PtsCirCen(i)

nCount = nCount + 1

End If

Next i

'(2)圖界點

'---------------------------

If nptExtA >= 0 Then

ReDim Preserve ptsI(nCount)

ptsI(nCount) = PtsExtA

nCount = nCount + 1

End If

'--------------------------

If nptExtB >= 0 Then

ReDim Preserve ptsI(nCount)

ptsI(nCount) = PtsExtB

nCount = nCount + 1

End If

'(3)轉折點

If nPtAdd >= 0 Then

For j = 0 To nPtAdd

ReDim Preserve ptsI(nCount)

ptsI(nCount) = ptsAddOnMapBdl(j)

nCount = nCount + 1

Next j

End If

'(4) 圖界交點

If nptInt >= 0 Then

For j = 0 To nptInt

ReDim Preserve ptsI(nCount)

ptsI(nCount) = ptsIntF(j)

nCount = nCount + 1

Next j

End If

' 繞測站依角度大小排列

Dim ptCenTpt As PointF

ptCenTpt.X = ptsXyz(IptGiven).X

ptCenTpt.Y = ptsXyz(IptGiven).Y

areaTpt = System.Math.Abs(Area_afterReorder(ptCenTpt, ptsI))

Erase PtsCirCen

For j = 0 To UBound(ptsI)

ReDim Preserve PtsCirCen(j)

PtsCirCen(j) = ptsI(j) '以PtsCirCen(j) 代替ptsI(j)

Next j

If isFill Then

mBrush = New SolidBrush(RandColors(IptGiven Mod 23)) '(ZColA(IptGiven + 3 Mod 255))

Call FillRegionXY(Gr, mBrush, ptsI)

Else

Gr.DrawPolygon(New Pen(Color.Blue, 2), ptsI)

End If

ReDim Preserve ThiessenXyzs(nThiessen)

With ThiessenXyzs(nThiessen)

.indPt = IptGiven

.coverArea = areaTpt

.ptStaXyz.X = ptsXyz(IptGiven).X

.ptStaXyz.Y = ptsXyz(IptGiven).Y

.ptStaXyz.z = ptsXyz(IptGiven).z

For i = 0 To UBound(PtsCirCen)

ReDim Preserve ThiessenXyzs(nThiessen).PLgon.pts(i)

.PLgon.pts(i) = PtsCirCen(i)

Next i

Frmdocument.List1.Items.Add("sta_ " & IptGiven & " ;area= " & areaTpt)

sumAreaOfAll = sumAreaOfAll + areaTpt

nThiessen = nThiessen + 1

End With

200: canvas.Refresh()

End Function

上次修改此網站的日期： 2018年11月25日