Delaunay
using System;using System.Collections.Generic;using System.ComponentModel;using System.Data;using System.Drawing;using System.Text;using System.Windows.Forms;using ESRI.ArcGIS.Carto;using ESRI.ArcGIS.Display;using ESRI.ArcGIS.Geometry;using ESRI.ArcGIS.Geodatabase;namespace drawtest{ public partial class Form1 : Form { public struct Vertex { public double x; public double y; public int ID; } public struct Triangle { public int V1Index; public int V2Index; public int V3Index; } public static int maxVertex = 300; public static int maxTiangle = 500; Vertex[] Vertices = new Vertex[maxVertex]; Triangle[] TinTriangles = new Triangle[maxTiangle]; public int[,] edges = new int[2, maxTiangle * 3]; public Form1() { InitializeComponent(); } private bool InCircle(float xp, float yp, float x1, float y1, float x2, float y2, float x3, float y3, double xc, double yc, double r) { double eps, m1, m2, mx1, mx2, my1, my2, dx, dy, rsqr, drsqr; eps = 0.00000001; if (Math.Abs(y1 - y2) < eps && Math.Abs(y2 - y3) < eps) { MessageBox.Show("INCIRCUM - F - Points are coincident !!"); return false; } if (Math.Abs(y2 - y1) < eps) { m2 = (-(Convert.ToDouble(x3) - Convert.ToDouble(x2)) / (Convert.ToDouble(y3) - Convert.ToDouble(y2))); mx2 = Convert.ToDouble((x2 x3) / 2.0); my2 = Convert.ToDouble((y2 y3) / 2.0); xc = Convert.ToDouble((x2 x1) / 2.0); yc = Convert.ToDouble(m2 * (xc - mx2) my2); } else if (Math.Abs(y3 - y2) < eps) { m1 = (-(Convert.ToDouble(x2) - Convert.ToDouble(x1)) / (Convert.ToDouble(y2) - Convert.ToDouble(y1))); mx1 = Convert.ToDouble((x1 x2) / 2.0); my1 = Convert.ToDouble((y1 y2) / 2.0); xc = Convert.ToDouble((x3 x2) / 2.0); yc = Convert.ToDouble(m1 * (xc - mx1) my1); } else { m1 = (-(Convert.ToDouble(x2) - Convert.ToDouble(x1)) / (Convert.ToDouble(y2) - Convert.ToDouble(y1))); m2 = (-(Convert.ToDouble(x3) - Convert.ToDouble(x2)) / (Convert.ToDouble(y3) - Convert.ToDouble(y2))); mx1 = Convert.ToDouble((x1 x2) / 2.0); mx2 = Convert.ToDouble((x2 x3) / 2.0); my1 = Convert.ToDouble((y1 y2) / 2.0); my2 = Convert.ToDouble((y2 y3) / 2.0); xc = Convert.ToDouble((m1 * mx1 - m2 * mx2 my2 - my1) / (m1 - m2)); yc = Convert.ToDouble(m1 * (xc - mx1) my1); } dx = (Convert.ToDouble(x2) - Convert.ToDouble(xc)); dy = (Convert.ToDouble(y2) - Convert.ToDouble(yc)); rsqr = Convert.ToDouble(dx * dx dy * dy); r = Convert.ToDouble(Math.Sqrt(rsqr)); dx = Convert.ToDouble(xp - xc); dy = Convert.ToDouble(yp - yc); drsqr = Convert.ToDouble(dx * dx dy * dy); if (drsqr <= rsqr) { return true; } return false; } private void button1_Click(object sender, EventArgs e) { IFeatureLayer pFeatureLayer; IFeatureClass pFeatureClass; IFeature pFeature; IQueryFilter pQueryFilter = new QueryFilterClass(); IPointCollection pPntCollection = new PolygonClass(); IPointCollection pPntCollection1 = new PolygonClass(); object missing = Type.Missing; int EdgeNum = 0; int TriangleNum = 0; int i; pFeatureLayer = axMapControl1.Map.get_Layer(0) as IFeatureLayer; pQueryFilter.WhereClause = ""; pFeatureClass = pFeatureLayer.FeatureClass; int PatchNumber = pFeatureClass.FeatureCount(pQueryFilter); int index = 0; pPntCollection = pFeatureLayer.FeatureClass.GetFeature(0).Shape as IPointCollection; int vertexnum = pPntCollection.PointCount; for ( i = 0; i < pPntCollection.PointCount; i ) { Vertices[i 1].x = pPntCollection.get_Point(i).X; Vertices[i 1].y = pPntCollection.get_Point(i).Y; Vertices[i 1].ID = pPntCollection.get_Point(i).ID; } //生成Delaunay bool[] Complete = new bool[maxTiangle ]; double xMin, xMax, yMin, yMax; double xMid, yMid, dMax; double dx, dy; //一般变量 int j, k; double xc = 0; double yc = 0; double r = 0; bool IsinCirlce = false; //找到离散点的最大值和最小值，构建supertriangle xMin = Vertices[1].x; yMin = Vertices[1].y; xMax = xMin; yMax = yMin; for ( i = 2; i <= vertexnum ; i ) { if (Vertices[i].x < xMin) xMin = Vertices[i].x; if (Vertices[i].x > xMax) xMax = Vertices[i].x; if (Vertices[i].y < yMin) yMin = Vertices[i].y; if (Vertices[i].y > yMax) yMax = Vertices[i].y; } dx = xMax - xMin; dy = yMax - yMin; if (dx > dy) dMax = dx; else dMax = dy; xMid = (xMax xMin) / 2; yMid = (yMax yMin) / 2; //将包含所有点的supertriangle作为第一个三角形 Vertices[vertexnum 1].x = xMid - 2 * dMax; Vertices[vertexnum 1].y = yMid - dMax; Vertices[vertexnum 2].x = xMid; Vertices[vertexnum 2].y = yMid 2 * dMax; Vertices[vertexnum 3].x = xMid 2 * dMax; Vertices[vertexnum 3].y = yMid - dMax; TinTriangles[1].V1Index =vertexnum 1; TinTriangles[1].V2Index = vertexnum 2; TinTriangles[1].V3Index = vertexnum 3; Complete[1] = false; TriangleNum = 1; for ( i = 1; i <= vertexnum; i ) { EdgeNum = 0; j = 0; do { j = j 1; if (Complete[j] != true) { IsinCirlce = InCircle((float ) Vertices[i].x,(float ) Vertices[i].y,(float ) Vertices[TinTriangles[j].V1Index].x,(float ) Vertices[TinTriangles[j].V1Index].y,(float ) Vertices[TinTriangles[j].V2Index].x,(float)Vertices[TinTriangles[j].V2Index].y,(float ) Vertices[TinTriangles[j].V3Index].x,(float ) Vertices[TinTriangles[j].V3Index].y, xc, yc, r); } if (IsinCirlce) { edges [0,EdgeNum 1 ] = TinTriangles[j].V1Index; edges[1, EdgeNum 1] = TinTriangles[j].V2Index; edges[0, EdgeNum 2] = TinTriangles[j].V2Index; edges[1, EdgeNum 2] = TinTriangles[j].V3Index; edges[0, EdgeNum 3] = TinTriangles[j].V3Index; edges[1, EdgeNum 3] = TinTriangles[j].V1Index; EdgeNum = EdgeNum 3; TinTriangles[j].V1Index = TinTriangles[TriangleNum].V1Index; TinTriangles[j].V2Index = TinTriangles[TriangleNum].V2Index; TinTriangles[j].V3Index = TinTriangles[TriangleNum].V3Index; Complete[j] = Complete[TriangleNum]; j = j - 1; TriangleNum = TriangleNum - 1; } } while (j < TriangleNum); for (j = 1; j < EdgeNum; j ) { if (edges[0,j ] != 0 && edges[1,j ] != 0) { for (k = j 1; k <= EdgeNum; k ) { if (edges[0,k ] != 0 && edges[1,k ] != 0) { if (edges[0,j ] == edges[1,k ]) { if (edges[1,j ] == edges[0,k ]) { edges[0,j ] = 0; edges[1,j ] = 0; edges[0,k ] = 0; edges[1,k ] = 0; } } } } } } for (j = 1; j <= EdgeNum; j ) { if (edges[0,j ] != 0 && edges[1,j ] != 0) { TriangleNum = TriangleNum 1; TinTriangles[TriangleNum].V1Index = edges[0,j ]; TinTriangles[TriangleNum].V2Index = edges[1,j ]; TinTriangles[TriangleNum].V3Index = i; Complete[TriangleNum] = false; } } } i = 0; do { i = i 1; if (TinTriangles[i].V1Index > vertexnum || TinTriangles[i].V2Index > vertexnum || TinTriangles[i].V3Index > vertexnum) { TinTriangles[i].V1Index = TinTriangles[TriangleNum].V1Index; TinTriangles[i].V2Index = TinTriangles[TriangleNum].V2Index; TinTriangles[i].V3Index = TinTriangles[TriangleNum].V3Index; i = i - 1; TriangleNum = TriangleNum - 1; } } while (i < TriangleNum); int many = TriangleNum; //画Delaunay ILine pLine = new LineClass(); ILine pLine1 = new LineClass(); ILine pLine2 = new LineClass(); //IPointCollection pPntCollection1 = new PolylineClass(); IPoint pPoint1 = new PointClass(); IPoint pPoint2 = new PointClass(); IPoint pPoint3 = new PointClass(); IGraphicsContainer pContainer = axMapControl1.Map as IGraphicsContainer; index = 0; for (i = 1; i <= TriangleNum; i ) { pPoint1.X = Vertices[TinTriangles[i].V1Index].x; pPoint1.Y = Vertices[TinTriangles[i].V1Index].y; pPoint2.X = Vertices[TinTriangles[i].V2Index].x; pPoint2.Y = Vertices[TinTriangles[i].V2Index].y; pPoint3.X = Vertices[TinTriangles[i].V3Index].x; pPoint3.Y = Vertices[TinTriangles[i].V3Index].y; pLine.PutCoords(pPoint1, pPoint2); pLine1.PutCoords(pPoint2, pPoint3); pLine2.PutCoords(pPoint1, pPoint3); ISegmentCollection pSeg = new PolylineClass() as ISegmentCollection; pSeg.AddSegment((ISegment)pLine, ref missing, ref missing); IPolyline pPolyline = pSeg as IPolyline; ISimpleLineSymbol pSimpleLineSymbol = new SimpleLineSymbolClass(); pSimpleLineSymbol.Style = esriSimpleLineStyle.esriSLSDot; ILineElement pLineElement = new LineElementClass(); IElement pElement = pLineElement as IElement; pElement.Geometry = pPolyline as IGeometry; pContainer.AddElement(pLineElement as IElement, index ); ISegmentCollection pSeg1 = new PolylineClass() as ISegmentCollection; pSeg1.AddSegment((ISegment)pLine1, ref missing, ref missing); IPolyline pPolyline1 = pSeg1 as IPolyline; ISimpleLineSymbol pSimpleLineSymbol1 = new SimpleLineSymbolClass(); pSimpleLineSymbol1.Style = esriSimpleLineStyle.esriSLSDot; ILineElement pLineElement1 = new LineElementClass(); IElement pElement1 = pLineElement1 as IElement; pElement1.Geometry = pPolyline1 as IGeometry; pContainer.AddElement(pLineElement1 as IElement, index ); ISegmentCollection pSeg2 = new PolylineClass() as ISegmentCollection; pSeg2.AddSegment((ISegment)pLine2, ref missing, ref missing); IPolyline pPolyline2 = pSeg2 as IPolyline; ISimpleLineSymbol pSimpleLineSymbol2 = new SimpleLineSymbolClass(); pSimpleLineSymbol2.Style = esriSimpleLineStyle.esriSLSDot; ILineElement pLineElement2 = new LineElementClass(); IElement pElement2 = pLineElement2 as IElement; pElement2.Geometry = pPolyline2 as IGeometry; pContainer.AddElement(pLineElement2 as IElement, index ); } axMapControl1.CenterAt(pPoint2); Application.DoEvents(); axMapControl1.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGeography, pContainer as object, axMapControl1.ActiveView.Extent); MessageBox.Show("OK"); } }}
评论