Code:
Imports Autodesk.AutoCAD.ApplicationServices
Imports Autodesk.AutoCAD.DatabaseServices
Imports Autodesk.AutoCAD.Runtime
Imports Autodesk.AutoCAD.EditorInput
Imports Autodesk.AutoCAD.GeometryPublic Class Triangulate2
Public Shared Function circum(ByVal x1 As Double, ByVal y1 As Double, ByVal x2 As Double, ByVal y2 As Double, ByVal x3 As Double, ByVal y3 As Double, _
ByRef xc As Double, ByRef yc As Double, ByRef r As Double) As Boolean
' Calculation of circumscribed circle coordinates and
' squared radius
Const eps As Double = 0.000001
Const big As Double = 1000000000000.0
Dim result As Boolean = True
Dim m1 As Double, m2 As Double, mx1 As Double, mx2 As Double, my1 As Double, my2 As Double, _
dx As Double, dy As Double
If (Math.Abs(y1 - y2) < eps) AndAlso (Math.Abs(y2 - y3) < eps) Then
result = False
xc = x1
yc = y1
r = big
Else
If Math.Abs(y2 - y1) < eps Then
m2 = -(x3 - x2) / (y3 - y2)
mx2 = (x2 + x3) / 2
my2 = (y2 + y3) / 2
xc = (x2 + x1) / 2
yc = m2 * (xc - mx2) + my2
ElseIf Math.Abs(y3 - y2) < eps Then
m1 = -(x2 - x1) / (y2 - y1)
mx1 = (x1 + x2) / 2
my1 = (y1 + y2) / 2
xc = (x3 + x2) / 2
yc = m1 * (xc - mx1) + my1
Else
m1 = -(x2 - x1) / (y2 - y1)
m2 = -(x3 - x2) / (y3 - y2)
If Math.Abs(m1 - m2) < eps Then
result = False
xc = x1
yc = y1
r = big
Else
mx1 = (x1 + x2) / 2
mx2 = (x2 + x3) / 2
my1 = (y1 + y2) / 2
my2 = (y2 + y3) / 2
xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2)
yc = m1 * (xc - mx1) + my1
End If
End If
End If
dx = x2 - xc
dy = y2 - yc
r = dx * dx + dy * dy
Return result
End Function
Friend Enum OutputObjectType
PolyFaceMesh = 1
Solid3d = 2
SubDMesh = 4
All = 7
End Enum
Friend Shared Sub TriangulatePoints(ByVal objType As OutputObjectType, ByVal maxpoints As Integer)
Dim doc As Document = Application.DocumentManager.MdiActiveDocument
Dim db As Database = doc.Database
Dim ed As Editor = doc.Editor
Dim createSubDMesh As Boolean = (objType And OutputObjectType.SubDMesh) > 0
Dim createPolyFaceMesh As Boolean = (objType And OutputObjectType.PolyFaceMesh) > 0
Dim createSolid3d As Boolean = (objType And OutputObjectType.Solid3d) > 0
Dim tvs As TypedValue() = {New TypedValue(0, "LWPOLYLINE")}
Dim sf As New SelectionFilter(tvs)
Dim pso As New PromptSelectionOptions()
pso.MessageForAdding = vbLf & "Select polylines:"
pso.AllowDuplicates = False
Dim psr As PromptSelectionResult = ed.GetSelection(pso, sf)
If psr.Status = PromptStatus.Error Then
Return
End If
If psr.Status = PromptStatus.Cancel Then
Return
End If
Dim ss As SelectionSet = psr.Value
Dim npts As Integer = ss.Count
Dim zref As Double = 0.0
If createSolid3d Then
Dim ps As PromptDoubleResult = ed.GetDouble(vbLf & "Enter Z coordinate of reference plane:")
If ps.Status <> PromptStatus.OK Then
Return
End If
zref = ps.Value
End If
Dim i As Integer, j As Integer, k As Integer, ntri As Integer, ned As Integer, nouted As Integer, _
status1 As Integer = 0, status2 As Integer = 0
Dim status As Boolean
' Point coordinates
Dim ptx As Double() = New Double(maxpoints + 2) {}
Dim pty As Double() = New Double(maxpoints + 2) {}
Dim ptz As Double() = New Double(maxpoints + 2) {}
' Triangle definitions
Dim pt1 As Integer() = New Integer(maxpoints * 2) {}
Dim pt2 As Integer() = New Integer(maxpoints * 2) {}
Dim pt3 As Integer() = New Integer(maxpoints * 2) {}
' Circumscribed circle
Dim cex As Double() = New Double(maxpoints * 2) {}
Dim cey As Double() = New Double(maxpoints * 2) {}
Dim rad As Double() = New Double(maxpoints * 2) {}
Dim xmin As Double, ymin As Double, xmax As Double, ymax As Double, dx As Double, dy As Double, _
xmid As Double, ymid As Double
Dim ed1 As Integer() = New Integer(maxpoints * 2) {}
Dim ed2 As Integer() = New Integer(maxpoints * 2) {}
Dim outed1 As Integer() = Nothing
If createSolid3d Then
outed1 = New Integer(maxpoints) {}
End If
'+++
Dim myExt As New Extents3d
'/+++
Dim idarray As ObjectId() = ss.GetObjectIds()
Dim tr As Transaction = db.TransactionManager.StartTransaction()
Dim myOutlines As New List(Of Curve)
Using tr
Dim Id As ObjectId
Dim a1 As Integer = 0
For Each Id In idarray
Dim v3d As Polyline = tr.GetObject(Id, OpenMode.ForRead)
myOutlines.Add(v3d)
Dim vertex As Point3d
Dim a2 As Integer = 0
myExt.AddExtents(v3d.GeometricExtents)
For a2 = 0 To v3d.NumberOfVertices - 1
vertex = v3d.GetPoint3dAt(a2)
ptx(a1 + a2) = vertex.X
pty(a1 + a2) = vertex.Y
ptz(a1 + a2) = vertex.Z
Next
npts = a1 + a2
a1 = a2
Next
tr.Commit()
End Using
If status2 > 0 Then
ed.WriteMessage(vbLf & "Ignored {0} point(s) with same coordinates.", status2)
End If
' Supertriangle
xmin = ptx(0)
xmax = xmin
ymin = pty(0)
ymax = ymin
For i = 0 To npts - 1
If ptx(i) < xmin Then
xmin = ptx(i)
End If
If ptx(i) > xmax Then
xmax = ptx(i)
End If
If pty(i) < xmin Then
ymin = pty(i)
End If
If pty(i) > xmin Then
ymax = pty(i)
End If
Next
dx = xmax - xmin
dy = ymax - ymin
xmid = (xmin + xmax) / 2
ymid = (ymin + ymax) / 2
i = npts
ptx(i) = xmid - (90 * (dx + dy)) - 100
pty(i) = ymid - (50 * (dx + dy)) - 100
ptz(i) = 0
pt1(0) = i
i += 1
ptx(i) = xmid + (90 * (dx + dy)) + 100
pty(i) = ymid - (50 * (dx + dy)) - 100
ptz(i) = 0
pt2(0) = i
i += 1
ptx(i) = xmid
pty(i) = ymid + 100 * (dx + dy + 1)
ptz(i) = 0
pt3(0) = i
ntri = 1
circum(ptx(pt1(0)), pty(pt1(0)), ptx(pt2(0)), pty(pt2(0)), ptx(pt3(0)), pty(pt3(0)), _
cex(0), cey(0), rad(0))
' Main loop
For i = 0 To npts - 1
ned = 0
xmin = ptx(i)
ymin = pty(i)
j = 0
While j < ntri
dx = cex(j) - xmin
dy = cey(j) - ymin
If ((dx * dx) + (dy * dy)) < rad(j) Then
ed1(ned) = pt1(j)
ed2(ned) = pt2(j)
ned += 1
ed1(ned) = pt2(j)
ed2(ned) = pt3(j)
ned += 1
ed1(ned) = pt3(j)
ed2(ned) = pt1(j)
ned += 1
ntri -= 1
pt1(j) = pt1(ntri)
pt2(j) = pt2(ntri)
pt3(j) = pt3(ntri)
cex(j) = cex(ntri)
cey(j) = cey(ntri)
rad(j) = rad(ntri)
j -= 1
End If
j += 1
End While
For j = 0 To ned - 2
For k = j + 1 To ned - 1
If (ed1(j) = ed2(k)) AndAlso (ed2(j) = ed1(k)) Then
ed1(j) = -1
ed2(j) = -1
ed1(k) = -1
ed2(k) = -1
End If
Next
Next
For j = 0 To ned - 1
If (ed1(j) >= 0) AndAlso (ed2(j) >= 0) Then
pt1(ntri) = ed1(j)
pt2(ntri) = ed2(j)
pt3(ntri) = i
status = circum(ptx(pt1(ntri)), pty(pt1(ntri)), ptx(pt2(ntri)), pty(pt2(ntri)), ptx(pt3(ntri)), pty(pt3(ntri)), _
cex(ntri), cey(ntri), rad(ntri))
If Not status Then
status1 += 1
End If
ntri += 1
End If
Next
Next
' Removal of outer triangles
i = 0
nouted = 0
While i < ntri
If (pt1(i) >= npts) OrElse (pt2(i) >= npts) OrElse (pt3(i) >= npts) Then
If createSolid3d Then
If (pt1(i) >= npts) AndAlso (pt2(i) < npts) AndAlso (pt3(i) < npts) Then
ed1(nouted) = pt2(i)
ed2(nouted) = pt3(i)
nouted += 1
End If
If (pt2(i) >= npts) AndAlso (pt1(i) < npts) AndAlso (pt3(i) < npts) Then
ed1(nouted) = pt3(i)
ed2(nouted) = pt1(i)
nouted += 1
End If
If (pt3(i) >= npts) AndAlso (pt1(i) < npts) AndAlso (pt2(i) < npts) Then
ed1(nouted) = pt1(i)
ed2(nouted) = pt2(i)
nouted += 1
End If
End If
ntri -= 1
pt1(i) = pt1(ntri)
pt2(i) = pt2(ntri)
pt3(i) = pt3(ntri)
cex(i) = cex(ntri)
cey(i) = cey(ntri)
rad(i) = rad(ntri)
i -= 1
End If
i += 1
End While
If createSolid3d Then
outed1(0) = 0
For i = 1 To nouted - 1
For j = 1 To nouted - 1
If ed2(outed1(i - 1)) = ed1(j) Then
outed1(i) = j
j = nouted
End If
Next
Next
outed1(nouted) = 0
End If
'+++
'punkt der garantiert auserhalb ist
Dim outerPoint As Point3d = myExt.MaxPoint
'+++
tr = db.TransactionManager.StartTransaction()
Using tr
Dim bt As BlockTable = DirectCast(tr.GetObject(db.BlockTableId, OpenMode.ForRead, False), BlockTable)
Dim btr As BlockTableRecord = DirectCast(tr.GetObject(bt(BlockTableRecord.ModelSpace), OpenMode.ForWrite, False), BlockTableRecord)
If createPolyFaceMesh Then
Dim pfm As New PolyFaceMesh()
btr.AppendEntity(pfm)
tr.AddNewlyCreatedDBObject(pfm, True)
For i = 0 To npts - 1
Dim vert As New PolyFaceMeshVertex(New Point3d(ptx(i), pty(i), ptz(i)))
pfm.AppendVertex(vert)
tr.AddNewlyCreatedDBObject(vert, True)
Next
For i = 0 To ntri - 1
'hier avarage Point
Dim pts As New List(Of Point3d)()
pts.Add(New Point3d(ptx(pt1(i)), pty(pt1(i)), ptz(pt1(i))))
pts.Add(New Point3d(ptx(pt2(i)), pty(pt2(i)), ptz(pt2(i))))
pts.Add(New Point3d(ptx(pt3(i)), pty(pt3(i)), ptz(pt3(i))))
Dim myavPoint As New Point3d(pts.Average(Function(a) a.X), pts.Average(Function(a) a.Y), pts.Average(Function(a) a.Z))
'Anzahl der Schittpunkt zählen
Dim inters As Integer = Interections(myOutlines, myavPoint)
Dim tt = inters Mod 2
'wenn ungerade sind wir innerhalb
If tt <> 0 Then
Dim face As New FaceRecord(CShort(pt1(i) + 1), CShort(pt2(i) + 1), CShort(pt3(i) + 1), 0)
pfm.AppendFaceRecord(face)
tr.AddNewlyCreatedDBObject(face, True)
End If
Next
End If
If createSubDMesh OrElse createSolid3d Then
Dim vertarray As New Point3dCollection()
Dim facearray As New Int32Collection()
For i = 0 To npts - 1
vertarray.Add(New Point3d(ptx(i), pty(i), ptz(i)))
Next
If createSolid3d Then
For i = 0 To nouted - 1
vertarray.Add(New Point3d(ptx(ed1(outed1(i))), pty(ed1(outed1(i))), zref))
Next
End If
j = 0
For i = 0 To ntri - 1
facearray.Add(3)
facearray.Add(pt1(i))
facearray.Add(pt2(i))
facearray.Add(pt3(i))
Next
If createSolid3d Then
For i = 0 To nouted - 1
facearray.Add(4)
k = outed1(i)
facearray.Add(ed1(k))
facearray.Add(ed2(k))
If i = nouted - 1 Then
facearray.Add(npts)
Else
facearray.Add(npts + i + 1)
End If
facearray.Add(npts + i)
Next
facearray.Add(nouted)
For i = 0 To nouted - 1
facearray.Add(npts + i)
Next
End If
Dim sdm As New SubDMesh()
sdm.SetDatabaseDefaults()
sdm.SetSubDMesh(vertarray, facearray, 0)
btr.AppendEntity(sdm)
tr.AddNewlyCreatedDBObject(sdm, True)
If createSolid3d Then
Dim sol As Solid3d = Nothing
Try
sol = sdm.ConvertToSolid(False, False)
btr.AppendEntity(sol)
tr.AddNewlyCreatedDBObject(sol, True)
Catch
ed.WriteMessage(vbLf & "Mesh was too complex to turn into a solid.")
End Try
If Not createSubDMesh Then
sdm.Erase()
End If
End If
End If
tr.Commit()
End Using
If status1 > 0 Then
ed.WriteMessage(vbLf & "Warning! {0} thin triangle(s) found!" & " Wrong result possible!", status1)
End If
Application.UpdateScreen()
End Sub
Public Shared Function Interections(ByVal crvs As List(Of Curve), ByVal pt As Point3d) As Integer
Dim myintersections As Integer
For Each crv As Curve In crvs
myintersections += InterctionsCount(crv, pt)
Next
Return myintersections
End Function
Enum IncidenceType
ToLeft = 0
ToRight = 1
ToFront = 2
Unknown
End Enum
Private Shared Function CurveIncidence(ByVal cur As Curve, ByVal param As Double, ByVal dir As Vector3d, ByVal normal As Vector3d) As IncidenceType
Dim deriv1 As Vector3d = cur.GetFirstDerivative(param)
If deriv1.IsParallelTo(dir) Then
' Need second degree analysis
Dim deriv2 As Vector3d = cur.GetSecondDerivative(param)
If (deriv2.IsZeroLength OrElse deriv2.IsParallelTo(dir)) Then
Return IncidenceType.ToFront
ElseIf (deriv2.CrossProduct(dir).DotProduct(normal) < 0) Then
Return IncidenceType.ToRight
Else
Return IncidenceType.ToLeft
End If
End If
If (deriv1.CrossProduct(dir).DotProduct(normal) < 0) Then
Return IncidenceType.ToLeft
Else
Return IncidenceType.ToRight
End If
End Function
Public Shared Function InterctionsCount(ByVal cur As Curve, ByVal testPt As Point3d) As Integer
If Not cur.Closed Then
Return False
End If
Dim ptOnCurve As Point3d = cur.GetClosestPointTo(testPt, False)
If Tolerance.Equals(testPt, ptOnCurve) Then
Return True
End If
' Check it's planar
Dim plane As Plane = cur.GetPlane
If Not cur.IsPlanar Then
Return False
End If
' Make the test ray from the plane
Dim normal As Vector3d = plane.Normal
Dim testVector As Vector3d = normal.GetPerpendicularVector
Dim ray As Ray = New Ray
ray.BasePoint = testPt
ray.UnitDir = testVector
Dim intersectionPoints As Point3dCollection = New Point3dCollection
' Fire the ray at the curve
cur.IntersectWith(ray, Intersect.OnBothOperands, intersectionPoints, 0, 0)
ray.Dispose()
Dim numberOfInters As Integer = intersectionPoints.Count
If (numberOfInters = 0) Then
Return False
End If
Dim nGlancingHits As Integer = 0
Dim epsilon As Double = 0.000002
' (trust me on this)
Dim i As Integer = 0
Do While (i < numberOfInters)
' Get the first point, and get its parameter
Dim hitPt As Point3d = intersectionPoints(i)
Dim hitParam As Double = cur.GetParameterAtPoint(hitPt)
Dim inParam As Double = (hitParam - epsilon)
Dim outParam As Double = (hitParam + epsilon)
Dim inIncidence As IncidenceType = CurveIncidence(cur, inParam, testVector, normal)
Dim outIncidence As IncidenceType = CurveIncidence(cur, outParam, testVector, normal)
If (((inIncidence = IncidenceType.ToRight) _
AndAlso (outIncidence = IncidenceType.ToLeft)) _
OrElse ((inIncidence = IncidenceType.ToLeft) _
AndAlso (outIncidence = IncidenceType.ToRight))) Then
nGlancingHits = (nGlancingHits + 1)
End If
i = (i + 1)
Loop
Return numberOfInters + nGlancingHits
End Function
End Class