Le portail francophone de la géomatique

VERSION 5.00
Begin {C62A69F0-16DC-11CE-9E98-00AA00574A4F} frm_FindDiameter 
   Caption         =   "Find Polgon Diameters"
   ClientHeight    =   3120
   ClientLeft      =   45
   ClientTop       =   330
   ClientWidth     =   4710
   OleObjectBlob   =   "frm_FindDiameter.frx":0000
   StartUpPosition =   1  'CenterOwner
End
Attribute VB_Name = "frm_FindDiameter"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Public turncomplete As Boolean
Private Sub cbx_Sel_Layer_Change()
  Dim pMxDoc As IMxDocument
  Dim pLayer As ILayer
  Dim pFeatLayer As IFeatureLayer
  Dim pEnumLayers As IEnumLayer
  Dim pFeatClass As IFeatureClass
  Dim strLayer As String
  Dim pAV As IActiveView
  Dim pFields As IFields
  Dim pField As IField
  Dim pEnumField As IEnumFieldMap
  Dim myArr As IArray
  
  chkbx_CreateAz.Enabled = True
  chkbx_CreateCard.Enabled = True
   
  Set pMxDoc = Application.Document
  Set pAV = pMxDoc.FocusMap
  Set pEnumLayers = pMxDoc.FocusMap.Layers
  pEnumLayers.Reset
  Set pLayer = pEnumLayers.Next
  Do While Not pLayer Is Nothing
    If UCase(pLayer.Name) = UCase(cbx_Sel_Layer.Value) Then Exit Do
    Set pLayer = pEnumLayers.Next
  Loop
  If pLayer Is Nothing Then
      Exit Sub
  End If

  Set pFeatLayer = pLayer
  Set pFeatClass = pFeatLayer.FeatureClass
  IndexA = pFeatClass.FindField("Azimuth")
  If (IndexA < 0) Then
  chkbx_CreateAz.Enabled = True
  End If
  IndexA = pFeatClass.FindField("Cardinal")
    If (IndexA < 0) Then
  chkbx_CreateCard.Enabled = True
  End If
 End Sub

Private Sub chkbx_CreateAz_Click()
  If chkbx_CreateAz Then
    btn_OK.Enabled = True
  Else
    If chkbx_CreateCard = False Then
      btn_OK.Enabled = False
    End If
  End If
End Sub

Private Sub chkbx_CreateCard_Click()
  If chkbx_CreateCard Then
    btn_OK.Enabled = True
  Else
    If chkbx_CreateAz = False Then
      btn_OK.Enabled = False
    End If
  End If
End Sub
Private Sub UserForm_Initialize()
  Dim pDoc As IMxDocument
  Dim pELayer As IEnumLayer
  Dim pLayer As ILayer, pLayer2 As ILayer
  Dim pFLayer As IFeatureLayer
  Dim pFClass As IFeatureClass
  Dim pMap As IMap
  Dim i As Integer

  Set pDoc = ThisDocument
  Set pELayer = pDoc.ActiveView.FocusMap.Layers
  pELayer.Reset
  Set pLayer = pELayer.Next

  Do While Not pLayer Is Nothing
    If TypeOf pLayer Is IFeatureLayer Then
      Set pFLayer = pLayer
      If pFLayer.FeatureClass.ShapeType = esriGeometryPolygon Then
        cbx_Sel_Layer.AddItem pLayer.Name
      End If
    End If
    Set pLayer = pELayer.Next
  Loop

  If cbx_Sel_Layer.ListCount = 0 Then
    MsgBox "Tool needs polygon layer", vbCritical, "EXITING"
    Exit Sub
  End If

' Set controls on Form
  chkbx_CreateAz.Enabled = False
  chkbx_CreateCard.Enabled = False
  btn_OK.Enabled = False

 End Sub
Private Sub btn_OK_Click()
  Dim strLayer As String
  Dim strFile As String
  strLayer = cbx_Sel_Layer.SelText
  ' Set controls on Form
  Call ProcessSelections(strLayer)
  MsgBox "Done", vbOKOnly, "Successfully Completed"
  'frm_FindDiameter.Hide
End Sub
Public Sub ProcessSelections(sLayerName As String)
Dim pMxDoc As IMxDocument
Dim pLayer As ILayer
Dim pFeatLayer As IFeatureLayer
Dim pEnumLayers As IEnumLayer
Dim pFeatClass As IFeatureClass
Dim pFeatCursor As IFeatureCursor
Dim pFeat As IFeature
Dim nFeat As Integer
Dim pTopoOp As ITopologicalOperator
Dim pHull As IPolygon
Dim pPolygon As IPolygon
Dim pPtColl As IPointCollection
Dim PolyDiPts As IPointCollection
Dim PolyDiMP As IMultipoint
Dim pAV As IActiveView
Dim pGC As IGraphicsContainer
Dim pline As ILine
Dim getRadian As Double
Dim getAngle As Integer, iAriAngle As Integer, iGeoAngle As Integer
Dim pi As Double
Dim IndexA As Long
Dim pField As IFieldEdit
    
pi = 4 * Atn(1) 'Get Calculated Value of pi
Set pMxDoc = Application.Document
Set pAV = pMxDoc.FocusMap
Set pGC = pAV
Set pEnumLayers = pMxDoc.FocusMap.Layers
Set PolyDiMP = New Multipoint
Set PolyDiPts = PolyDiMP
' find the requested layer:
  Set pLayer = pEnumLayers.Next
  Do While Not pLayer Is Nothing
    If UCase(pLayer.Name) = UCase(sLayerName) Then Exit Do
    Set pLayer = pEnumLayers.Next
  Loop
  If pLayer Is Nothing Then
      Exit Sub
  End If
' get the feature cursor:
  Set pFeatLayer = pLayer
  Set pFeatClass = pFeatLayer.FeatureClass
  If chkbx_CreateAz Then
    IndexA = pFeatClass.FindField("Azimuth")
    If IndexA < 0 Then
      Set pField = New Field
      With pField
        .Type = esriFieldTypeInteger
        .Name = "Azimuth"
      End With
      pFeatClass.AddField pField
    End If
    Index_Az = pFeatClass.FindField("Azimuth")
  End If
  If chkbx_CreateCard Then
    IndexA = pFeatClass.FindField("Cardinal")
    If IndexA < 0 Then
      Set pField = New Field
      With pField
        .Type = esriFieldTypeInteger
        .Name = "Cardinal"
      End With
      pFeatClass.AddField pField
    End If
    Index_Card = pFeatClass.FindField("Cardinal")
  End If
Set pFeatCursor = pFeatClass.Search(Nothing, False) 'Get a Cursor for all polygons in layer

' loop through the features:

Set pFeat = pFeatCursor.NextFeature 'Get first polygon
If (pFeat.Shape.GeometryType = esriGeometryPolygon) Then  'make sue it is a polgyon
  nFeat = 0  'Initialize polygon counter
  Do While Not pFeat Is Nothing   'Loop through all polygons
    nFeat = nFeat + 1   'Keep track of number of polygon processed
    Set pFeatPoints = pFeat.Shape   'Get the points that make up a polygon
    Set pPolygon = pFeat.ShapeCopy
    Set pTopoOp = pFeatPoints   'Assign points to a Topological Operator which allows for constructing new geometries
    Set pHull = pTopoOp.ConvexHull  'Get the convex hull of the point collection and assign to a polygon
    Set PolyDiPts = GetPolyDiPts(pHull) 'Pass the convex hull polgyon to a function to retrieve antipodal points of maximum distance apart

    Dim Apt As IPoint
    Dim Bpt As IPoint
    Set Apt = New Point
    Set Bpt = New Point

    Set pline = New Line    'Create empty line object
    PolyDiPts.QueryPoint 0, Apt 'Assign first antipodal point to a new point
    PolyDiPts.QueryPoint 1, Bpt 'Assgin second antipodal point to a new Point
    pline.PutCoords Apt, Bpt    'Create line from antipodal points
    getRadian = pline.Angle     'Retrieve angle of line object from a property of that object
    If getRadian < 0 Then   'Convert radians to Polar
      iAriAngle = 360 - (((getRadian * 360) * -1) / (2 * pi))
    Else
      iAriAngle = (getRadian * 360) / (2 * pi)
    End If
    
    getAngle = 45 - (iAriAngle - 45)
    If getAngle < 0 Then    'Convert polar to Azimuth
      iGeoAngle = getAngle + 360
    Else
      iGeoAngle = getAngle
    End If
    If chkbx_CreateAz Then
      pFeat.Value(Index_Az) = iGeoAngle 'Set Attribute field of polygon to Azimuth
      pFeat.Store
    End If
    If chkbx_CreateCard Then
      pFeat.Value(Index_Card) = GetCard(pPolygon) 'Set Attribute field of polygon to Azimuth
      pFeat.Store
    End If
    Set pFeat = pFeatCursor.NextFeature 'Select next polygon to process
  Loop
Else
  MsgBox ("Wrong type feature geometry")
  Exit Sub
End If
End Sub
Private Function GetCard(pPolygon As IPolygon) As Integer
    'Get envelope
    Dim pArea As IArea
    Dim myCntrd As IPoint
    Dim ratio As Double
    Dim pEnv As IEnvelope
    Set pArea = pPolygon
    Set myCntrd = pArea.Centroid
    Set pEnv = pPolygon.Envelope
    ratio = pEnv.Height / pEnv.Width
    If ratio > 1 Then
      GetCard = 0
      Else
      GetCard = 270
    End If
End Function
Private Function MakeMarkerElement(pPoint As IPoint) As IMarkerElement
Dim pElement As IElement
Dim pMxDoc As IMxDocument
Dim pDocDefaultSymbols As IDocumentDefaultSymbols
Set pElement = New MarkerElement
pElement.Geometry = pPoint
Set MakeMarkerElement = pElement
Set pMxDoc = ThisDocument
Set pDocDefaultSymbols = pMxDoc
MakeMarkerElement.Symbol = pDocDefaultSymbols.MarkerSymbol
Exit Function
End Function
Public Function GetPolyDiPts(pHull As IPolygon) As IPointCollection 'Passed convex hull polygon.  returns max antipodal pair
Dim pPtColl As IPointCollection
Dim AntiPodalNdex As Integer
Dim maxDist As Double
Dim myPtColl As IPointCollection
Dim MaxDiColl As IPointCollection
Dim MaxMP As IMultipoint
Set MaxDiColl = MaxMP   'Set to empty collection of multipoint object
Set myPtColl = MaxMP    'Set to empty collection of multipoint object
Set GetPolyDiPts = MaxMP 'Set to empty collection of multipoint object
Set pPtColl = pHull 'Get points that make up convex hull polygon
pPtColl.RemovePoints pPtColl.PointCount - 1, 1  'Remove repeated polygon closure point
If (pPtColl.PointCount = 3) Then    'If poygon is a triangle
  Set GetPolyDiPts = FindPtsFor3(pPtColl, 0, 1, 2) 'Call function that solves easy case. Pass points and index starting points. Returns antipodal pair.
Else    'If polygon is not triangle
  Dim over As Boolean   'Need to define
  Dim currentIndex As Integer   'Algorithm needs to track 3 points in point collection
  Dim prevIndex As Integer  'previous index number to current point in point collection
  Dim Dist As Long  'Need to define
  over = False
  currentIndex = 0  'Initialize current point index in point collection
  prevIndex = FindPrvCCWpt(pPtColl, currentIndex)   'Find the previous index value of current index point. Passes point collection and current index point. returns previous point index value
  AntiPodalNdex = FindAntiPndx(pPtColl, currentIndex, 1) 'Find the index value to the antipodal point match for the current index point
  maxDist = FindMaxFor3(pPtColl, currentIndex, prevIndex, AntiPodalNdex)    'Finds maxDist of Antipodal pair (pair could be previous index or current index and Antpodal index
  Set MaxDiColl = FindPtsFor3(pPtColl, currentIndex, prevIndex, AntiPodalNdex)  'Determine antipodal pair from two posibilities
  Dist = 0  'Initialize max distance
  While (over = False)  'Loop through other points in collection traversing convex hull
    prevIndex = currentIndex    'increment previous index
    currentIndex = FindNxtCCWpt(pPtColl, currentIndex)  'increment current index
    Dist = FindMaxFor3(pPtColl, currentIndex, prevIndex, AntiPodalNdex) 'We can now start at previous antipodal index value when determining next antipodal point
    Set myPtColl = FindPtsFor3(pPtColl, currentIndex, prevIndex, AntiPodalNdex) 'Get antipodal pair
    If (Dist > maxDist) Then    'Compare distances of antipodal pairs
      maxDist = Dist    'Assign max distance if greater
      Set MaxDiColl = myPtColl  'Assign point collection with max Antipodal pair
    End If
    If turncomplete Then over = True    'Assign boolean to indicate traversing convex hull is complete
  Wend
 Set GetPolyDiPts = MaxDiColl   'return point collection of antipodal points representing polygon diameter
 End If
 End Function
Public Function Euclidean_Dist(Firstpt As IPoint, Secondpt As IPoint) As Double 'Passed points calculates distance. Returns distance
Euclidean_Dist = (Sqr(((Firstpt.X - Secondpt.X) ^ 2) + (Firstpt.Y - Secondpt.Y) ^ 2))   'Pythagorean Theorem
End Function
Public Function FindPtsFor3(pPtColl As IPointCollection, ndex1 As Integer, ndex2 As Integer, ndex3 As Integer) As IPointCollection  'Passed polygon point collection and index values of points. Returns antipodal pair
Dim dist_12 As Double
Dim dist_23 As Double
Dim dist_31 As Double
Dim pt1 As IPoint
Dim pt2 As IPoint
Dim pt3 As IPoint
Dim Chk As Integer
Dim myMP As IMultipoint
Dim myPtColl As IPointCollection
Set myMP = New Multipoint
Set myPtColl = myMP
Set pt1 = New esriGeometry.Point
Set pt2 = New esriGeometry.Point
Set pt3 = New esriGeometry.Point
pPtColl.QueryPoint ndex1, pt1   'Get point in point collection from passed index value
pPtColl.QueryPoint ndex2, pt2   'Get point in point collection from passed index value
pPtColl.QueryPoint ndex3, pt3   'Get point in point collection from passed index value
dist_12 = Euclidean_Dist(pt1, pt2)  'Call function to get distance between points
dist_23 = Euclidean_Dist(pt2, pt3)  'Call function to get distance between points
dist_31 = Euclidean_Dist(pt3, pt1)  'Call function to get distance between points
Dim maxDist As Double
maxDist = dist_12   'Assign initial maximum distance
Chk = 1
If (dist_23 > maxDist) Then 'Check to see which distance is greatest
  Chk = 2
ElseIf (dist_31 > maxDist) Then
  Chk = 3
End If
Select Case Chk 'Depending on Case assign max antipodal pair to point collection
  Case 1
    myPtColl.AddPoint pt1
    myPtColl.AddPoint pt2
  Case 2
    myPtColl.AddPoint pt2
    myPtColl.AddPoint pt3
  Case 3
    myPtColl.AddPoint pt3
    myPtColl.AddPoint pt1
  End Select
  Set FindPtsFor3 = myPtColl    'Assign point collection to returning value
End Function
Public Function FindMaxFor3(pPtColl As IPointCollection, ndex1 As Integer, ndex2 As Integer, ndex3 As Integer) As Double    '
Dim dist_23 As Double
Dim dist_31 As Double
Dim pt1 As IPoint
Dim pt2 As IPoint
Dim pt3 As IPoint
Dim Chk As Integer
Set pt1 = New esriGeometry.Point
Set pt2 = New esriGeometry.Point
Set pt3 = New esriGeometry.Point
pPtColl.QueryPoint ndex1, pt1   'Get points from their index position
pPtColl.QueryPoint ndex2, pt2
pPtColl.QueryPoint ndex3, pt3
dist_12 = Euclidean_Dist(pt1, pt2)  'Calculate distance between points
dist_23 = Euclidean_Dist(pt2, pt3)
dist_31 = Euclidean_Dist(pt3, pt1)
Dim maxDist As Double
maxDist = dist_12   'Initialize maxDist value
If (dist_23 > maxDist) Then maxDist = dist_23   'Determine maximum distance
If (dist_31 > maxDist) Then maxDist = dist_31
FindMaxFor3 = maxDist   'Return maximum distance
End Function
Public Function FindNxtCCWpt(pPtColl As IPointCollection, ndex As Integer) As Integer   'passes point collection and point index value
' Find Next Counter ClockWise Point
If (ndex >= pPtColl.PointCount - 1) Then    'If finished traversing convex hull
  turncomplete = True   'Notify calling routines that traversing is complete
   FindNxtCCWpt = 0
Else
  FindNxtCCWpt = ndex + 1
End If
End Function
Public Function FindPrvCCWpt(pPtColl As IPointCollection, ndex As Integer) As Integer   'Passed point collection and current point index value.  Returns previous point index value
'Find Previous Count Clockwise Point FindPrvCCWpt
If (ndex <= 0) Then 'If at the 0 index point then
  FindPrvCCWpt = pPtColl.PointCount - 1 'Get last point in point collection index value
Else    'if index not 0
  FindPrvCCWpt = ndex - 1   'Get the previous point index value
End If
End Function
Public Function FindAntiPndx(pPtColl As IPointCollection, ndex As Integer, startndex As Integer) As Integer 'Passed the point collection, current point index value, index value to start process to find antipodal point
Dim a As Integer
Dim b As Integer
Dim c As Integer
Dim current As Integer
Dim area_abc As Double
Dim area_abcurrent As Double
Dim Apt As IPoint   'Establish point objects to work with
Dim Bpt As IPoint
Dim Cpt As IPoint
Set Apt = New esriGeometry.Point
Set Bpt = New esriGeometry.Point
Set Cpt = New esriGeometry.Point
b = ndex    'Value of current point index
c = startndex   'Value of starting index point for processing
a = FindPrvCCWpt(pPtColl, b)    'Find the previous point index value of current index
current = FindPrvCCWpt(pPtColl, c)  'find the previous point index of the starting index for processing
pPtColl.QueryPoint a, Apt   'Get the points from the point collection for processing
pPtColl.QueryPoint b, Bpt
pPtColl.QueryPoint c, Cpt
area_abc = AreaOfTriangle(Apt, Bpt, Cpt)    'Calculate the area the intial 3 points
pPtColl.QueryPoint current, Cpt 'Get index starting point
area_abcurrent = AreaOfTriangle(Apt, Bpt, Cpt)  'Calculate area
While (area_abc >= area_abcurrent)  'While traversing the points of a convex hull polygon
    current = c 'calculating area from current and previous point and next traversed point
    c = FindNxtCCWpt(pPtColl, c)    'When the area of the triangle starts decreasing you have found your
    pPtColl.QueryPoint c, Cpt   'antipodal pair of points
    area_abc = AreaOfTriangle(Apt, Bpt, Cpt)
    pPtColl.QueryPoint current, Cpt
    area_abcurrent = AreaOfTriangle(Apt, Bpt, Cpt)
 Wend
FindAntiPndx = current  'Return the index value of the antipodal point of the point collection
End Function
Public Function AreaOfTriangle(Apt As IPoint, Bpt As IPoint, Cpt As IPoint) As Double   'Pass three points return area
'Calculate area of triangle using the Heron's formula
Dim a As Double
Dim b As Double
Dim c As Double
Dim s As Double
a = Abs(Euclidean_Dist(Apt, Bpt))   'Get distance between points
b = Abs(Euclidean_Dist(Bpt, Cpt))
c = Abs(Euclidean_Dist(Cpt, Apt))
s = 0.5 * (a + b + c) ' Semi-perimeter
AreaOfTriangle = Sqr(s * ((s - a) * (s - b) * (s - c)))
End Function

Public Function Features2ConvexHull(sLayerName As String, sTextFileName As String, Optional nMaxFeatures As Double)
Dim pFeat As IFeature
Dim pLayer As ILayer
Dim pMxDoc As IMxDocument
Dim pEnumLayers As IEnumLayer
Dim i As Integer
Dim pFeatClass As IFeatureClass
Dim pFeatLayer As IFeatureLayer
Dim pFeatCursor As IFeatureCursor
Dim sY As String, sX As String, sFeatN As String, sFeatVerts As String
Dim strCenX As String, strCenY As String, strAngle As String
Dim lFileID As Long
Dim pFeatPoints As IPointCollection
Dim pArea As IArea
Dim pPoints As Point
Dim nFeat As Double
Dim pAV As IActiveView
Dim pGC As IGraphicsContainer

' get the document and map objects
Set pMxDoc = Application.Document
Set pAV = pMxDoc.FocusMap
Set pGC = pAV
Set pEnumLayers = pMxDoc.FocusMap.Layers

' find the requested layer:
  Set pLayer = pEnumLayers.Next
  Do While Not pLayer Is Nothing
    If UCase(pLayer.Name) = UCase(sLayerName) Then Exit Do
    Set pLayer = pEnumLayers.Next
  Loop
  
  If pLayer Is Nothing Then
      Exit Function
  End If
' get the feature cursor:
  Set pFeatLayer = pLayer
  Set pFeatClass = pFeatLayer.FeatureClass
  Set pFeatCursor = pFeatClass.Search(Nothing, False)
' loop through the features:
Set pFeat = pFeatCursor.NextFeature
If (pFeat.Shape.GeometryType = esriGeometryPolygon) Then
  Do While Not pFeat Is Nothing
    nFeat = nFeat + 1
    If (nFeat > nMaxFeatures) And (nMaxFeatures > 0) Then
      Exit Do
    End If
      Set pFeatPoints = pFeat.Shape
      Dim pTopoOp As ITopologicalOperator
      Set pTopoOp = pFeatPoints
      Dim pHull As IPolygon
      Set pHull = pTopoOp.ConvexHull
  Loop
Else
  MsgBox ("Wrong type feature geometry")
  Return
End If
End Function