AreEngine 求最小面積的外接矩形，非IEnvelope，表達不清楚了
阿新 • • 發佈：2018-12-16
1，總是會得到一些奇奇怪怪的要求，求一個面對象的外接最小面積的矩形，和ArcToolBox中的Mininum Bounding Geometry功能下的RECTANGLE_BY_AREA想似。具體看下圖：
區別如上圖所示：IEnvelope 得到的是下圖所示，需要的是第一種
（只是記錄一下，可以解決問題，但不是最優方法，程式碼冗餘量大，待解決）
基本思路：獲取所有邊，以其中的某一條邊入手，遍歷所有點到這條邊的距離，取最大距離的點（pointA），以此點為準畫一條平行於第一條線的線，作為第二條邊，再以此點（pointA）畫一條垂直於第一條邊的線作為輔助線，再次遍歷所有點到此輔助線的距離，得到最大距離的點（pointB），以pointB畫一條垂直第一條邊的線作為第三條邊，再次遍歷所有點到第三條邊的距離，取得最大距離的點pointC,以pointC畫一條垂直第一條邊的線，作為第四條邊。然後相鄰兩條線求交點，得到四個點，用四個點構造矩形。計算面積，以此，將所有的矩形都得到然後，比較面積求最小的就好了，好雞兒囉嗦我。
不喜就噴！接受反駁！
實現方法如下：
1 private void Mininum() 2 { 3 IFeatureClass pFClass = (this.axMapControl1.get_Layer(0) as IFeatureLayer).FeatureClass; 4 IFeatureCursor pFCursor = pFClass.Search(null, false); 5 IFeature pFeature = pFCursor.NextFeature();6 while (pFeature != null) 7 { 8 //求凸包 9 IGeometry pGeometry = pFeature.Shape; 10 ITopologicalOperator pTOperator = pGeometry as ITopologicalOperator; 11 IGeometry pNEWGeometry = pTOperator.ConvexHull();12 //每個凸包邊的集合 13 ISegmentCollection pSCollection = pNEWGeometry as ISegmentCollection; 14 //每個凸包點的集合 15 IPointCollection pPCollection = pNEWGeometry as IPointCollection; 16 //用來篩選最小的面積的一個界值，只要足夠大，多少都可以 17 double area = 99999999999999; 18 //儲存當前遍歷的要素（feature）的最小面積矩形物件。 19 IGeometry pGeoMetry = null; 20 for (int i = 0; i < pSCollection.SegmentCount; i++) 21 { 22 double DIS = 0; 23 double DIS3 = 0; 24 double DIS4 = 0; 25 26 IPoint RESUPoint=new PointClass(); 27 IPoint Resu3point = new PointClass(); 28 IPoint Resu4point = new PointClass(); 29 ISegment pPolyline_Segment = pSCollection.get_Segment(i); 30 IPoint FromPoint = pPolyline_Segment.FromPoint; 31 IPoint ToPoint = pPolyline_Segment.ToPoint; 32 //y=kx+b下邊兩行為第一條邊的表示式的引數值 33 double k = (FromPoint.Y - ToPoint.Y) / (FromPoint.X - ToPoint.X); 34 double b = FromPoint.Y - FromPoint.X * k; 35 for (int j = 0; j < pPCollection.PointCount; j++) 36 { 37 IPoint eachPoint = pPCollection.get_Point(j); 38 double distanse = dis(eachPoint, k, b); 39 if (distanse > DIS) 40 { 41 DIS = distanse; 42 RESUPoint = eachPoint; 43 } 44 } 45 //第二條邊的引數值 46 double k1 = k; 47 double b1 = RESUPoint.Y - k1 * RESUPoint.X; 48 //求第二個點到第一條邊的投影點及所在直線的引數（輔助線） 49 double k2 = -1 / k; 50 double b2 = RESUPoint.Y - k2 * RESUPoint.X; 51 //遍歷第三條邊所在的位置 52 for (int m = 0; m < pPCollection.PointCount; m++) 53 { 54 IPoint pPoint = pPCollection.get_Point(m); 55 double distance = dis(pPoint, k2, b2); 56 if (distance > DIS3) 57 { 58 DIS3 = distance; 59 Resu3point = pPoint; 60 } 61 } 62 //第三條的引數值 63 double k3 = k2; 64 double b3 = Resu3point.Y - k3 * Resu3point.X; 65 //------------開始求第四條邊 66 for (int n = 0; n < pPCollection.PointCount; n++) 67 { 68 IPoint pPoint=pPCollection.get_Point(n); 69 double distance = dis(pPoint,k3,b3); 70 if(distance>DIS4) 71 { 72 DIS4=distance; 73 Resu4point = pPoint; 74 } 75 } 76 //第四條邊的引數值； 77 double k4 = k3; 78 double b4 = Resu4point.Y-k4*Resu4point.X; 79 80 //------------四個交點的值 81 //pointA 82 double x_a=(b4-b)/(k-k4); 83 double y_a = k * x_a + b; 84 //pointB 85 double x_b = (b4 - b1) / (k1 - k4); 86 double y_b = k1 * x_b + b1; 87 //pointC 88 double x_c = (b3 - b1) / (k1 - k3); 89 double y_c = k1 * x_c + b1; 90 //pointD 91 double x_d = (b3 - b) / (k - k3); 92 double y_d = k * x_d + b; 93 //構造矩形 94 IPoint pointa = new PointClass(); 95 pointa.PutCoords(x_a, y_a); 96 IPoint pointb = new PointClass(); 97 pointb.PutCoords(x_b, y_b); 98 IPoint pointc = new PointClass(); 99 pointc.PutCoords(x_c, y_c); 100 IPoint pointd = new PointClass(); 101 pointd.PutCoords(x_d, y_d); 102 IPoint pointa1 = new PointClass(); 103 pointa1.PutCoords(x_a, y_a); 104 IPointCollection RECTANG = new PolygonClass(); 105 RECTANG.AddPoint(pointa); 106 RECTANG.AddPoint(pointb); 107 RECTANG.AddPoint(pointc); 108 RECTANG.AddPoint(pointd); 109 RECTANG.AddPoint(pointa1); 110 IPolygon pPolygon = RECTANG as IPolygon; 111 IGeometry PPPPPPPPGeo = pPolygon as IGeometry; 112 //************************* 113 //axMapControl1.FlashShape(PPPPPPPPGeo); 114 //注意此處！因構造矩形的時候，順時針逆時針不能確定，因此從ArcGIS來說這個面積有可能為負值，並不是取絕對值就可以，外環與內環的區別，因此用數學方法。 115 double a1=Math.Sqrt((pointb.Y-pointa.Y)*(pointb.Y-pointa.Y)+(pointb.X-pointa.X)*(pointb.X-pointa.X)); 116 double c1=Math.Sqrt((pointc.Y-pointb.Y)*(pointc.Y-pointb.Y)+(pointc.X-pointb.X)*(pointc.X-pointb.X)); 117 double pArea = a1 * c1; 118 if (Math.Abs(pArea) < area) 119 { 120 area = Math.Abs(pArea); 121 pGeoMetry = PPPPPPPPGeo; 122 } 123 } 124 axMapControl1.FlashShape(pGeoMetry); 125 pFeature = pFCursor.NextFeature(); 126 } 127 System.Runtime.InteropServices.Marshal.ReleaseComObject(pFCursor); 128 }
1 public double dis(IPoint pPoint, double k, double b) 2 { 3 double numerator = System.Math.Abs(k * pPoint.X - pPoint.Y + b); 4 double Denominator = System.Math.Sqrt(k * k + 1); 5 double dis=numerator/Denominator; 6 return Math.Abs(dis); 7 }
弱弱的再囉嗦一句：這個功能在arcpy中很好實現
arcpy.MinimumBoundingGeometry_management("parks.shp", "c:/output/output.gdb/parks_mbg", "RECTANGLE_BY_AREA", "NONE")
就一句就可以了，只是這個東西速度慢的要死，再次使用最小得到的矩形的又得選擇。
OK！BBBBBBBBBBBBBBBBBBBBBBBBBBB完了，