凸包
參考
安德魯算法
分治法(其中nfox的項目實現的是分治法)
多邊形快速凸包算法(Melkman‘s Algorithm)
還可以這看cpp的代碼: https://www.cnblogs.com/VividBinGo/p/11637684.html
定義
凸包又叫凸多邊形,本篇文章可能混用兩種說法,形象的理解就是一些點(點集)用一根橡皮筋緊緊地包裹外邊點.
如果知道了這個定義,那么還有:
用一個保鮮膜裹着三維點,求膜上點集.
用一個最小的球裹着三維點,求球球的中心點和直徑.
這樣就進入了一個叫拓撲學的學科上.......我的老天鵝.
我竟然搜了半天沒看到可以直接拿來用的c#代碼,還是小軒軒給我的....
葛立恆凸包
注意一下,如果點集形成的是正交矩形的情況時,
算出來的凸包會多一個點,可以進行后處理.
(你會發現代碼實際上是右上角最大點開始的,其他的教程說明從最小點開始算,
這需要知道的是最小最大都是邊界點,邊界點必然是凸包的邊點,用誰起算都可以)
主函數:
#if !HC2020
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.Runtime;
#else
using GrxCAD.ApplicationServices;
using GrxCAD.DatabaseServices;
using GrxCAD.EditorInput;
using GrxCAD.Geometry;
using GrxCAD.Runtime;
#endif
using System.Collections.Generic;
using System.Linq;
using static System.Math;
using static JoinBox.MathHelper;
namespace JoinBox
{
/*
視頻參考: https://www.bilibili.com/video/BV1v741197YM
相關學習: https://www.cnblogs.com/VividBinGo/p/11637684.html
① 找到所有點中最左下角的點_p0(按 x 升序排列，如果 x 相同，則按 y 升序排列)
② 以_p0為基准求所有點的極角,並按照極角排序(按極角升序排列,若極角相同,則按距離升序排列),設這些點為p1,p2,……,pn-1
③ 建立一個棧,_p0,p1進棧，對於P[2..n-1]的每個點,利用叉積判斷,
若棧頂的兩個點與它不構成"向左轉(逆時針)"的關系,則將棧頂的點出棧,直至沒有點需要出棧以后,將當前點進棧
④ 所有點處理完之后棧中保存的點就是凸包了。
*/
public partial class Graham
{
/// <summary>
/// 最靠近x軸的點(必然是凸包邊界的點)
/// </summary>
Point3d _p0;
/// <summary>
/// 求凸包測試命令
/// </summary>
[CommandMethod("test_gra")]
public void test_gra()
{
var doc = Application.DocumentManager.MdiActiveDocument;
var ed = doc.Editor;
var db = doc.Database;//當前的數據庫
ed.WriteMessage("\n****{驚驚連盒}求凸包,選擇曲線:");
//定義選擇集選項
var pso = new PromptSelectionOptions
{
RejectObjectsOnLockedLayers = true, //不選擇鎖定圖層對象
AllowDuplicates = true, //不允許重復選擇
};
var ssPsr = ed.GetSelection(pso);//手選 這里輸入al會變成all,無法刪除ssget的all關鍵字
if (ssPsr.Status != PromptStatus.OK)
return;
db.Action(tr => {
var getPts = new List<Point3d>();
foreach (ObjectId id in ssPsr.Value.GetObjectIds())
{
var ent = id.ToEntity(tr);
if (ent is Curve curve)
{
var cs = new CurveSample(curve);
getPts.AddRange(cs.GetSamplePoints);
}
else if (ent is DBPoint bPoint)
getPts.Add(bPoint.Position);
else
{
var entPosition = ent.GetType().GetProperty("Position");//反射獲取屬性
if (entPosition != null)
{
var pt = (Point3d)entPosition.GetValue(null, null);
getPts.Add(pt);
}
}
}
//葛立恆方法
var pts = GrahamConvexHull(getPts).ToPoint2ds();
ed.WriteMessage("\n\r凸包對踵點最大距離:" + RotateCalipersMax(pts));
ed.WriteMessage("\n\r凸包對踵點最小距離:" + RotateCalipersMin(pts));
var psd = new PointsDistance<Point2d>();
ed.WriteMessage("\n\r凸包點集最大距離:" + psd.Min(pts));
ed.WriteMessage("\n\r凸包點集最小距離:" + psd.Max(pts));
var bv = new List<BulgeVertex>();
for (int i = 0; i < pts.Count(); i++)
{
bv.Add(new BulgeVertex(pts[i], 0));
}
Entity pl = EntityAdd.AddPolyLineToEntity(bv);
tr.AddEntityToMsPs(db, pl);
#if true3
var recs = Boundingbox(pts);
//生成所有的包圍盒,每條邊生成一個
int ColorIndex = 0;
foreach (var rec in recs)
{
bv = new List<BulgeVertex>
{
new BulgeVertex(rec.R1, 0),
new BulgeVertex(rec.R2, 0),
new BulgeVertex(rec.R3, 0),
new BulgeVertex(rec.R4, 0)
};
pl = EntityAdd.AddPolyLineToEntity(0, bv);
pl.ColorIndex = ++ColorIndex;
tr.AddEntityToMsPs(db, pl);
}
#endif
//生成計算面積最小的包圍盒
var recAreaMin = BoundingboxAreaMin(pts);
bv = new List<BulgeVertex>
{
new BulgeVertex(recAreaMin.R1, 0),
new BulgeVertex(recAreaMin.R2, 0),
new BulgeVertex(recAreaMin.R3, 0),
new BulgeVertex(recAreaMin.R4, 0)
};
pl = EntityAdd.AddPolyLineToEntity(bv);
tr.AddEntityToMsPs(db, pl);
});
}
/// <summary>
/// 角度p0和pn的角度
/// </summary>
/// <param name="pn"></param>
/// <returns></returns>
double Cosine(Point3d pn)
{
double d = _p0.DistanceTo(pn);
//距離是0表示是自己和自己的距離,那么0不可以是除數,否則Nan:求角度(高/斜)==sin(角度)
var angle = d == 0.0 ? 1.0 : (pn.Y - _p0.Y) / d;
//var angle = d == 0 ? 0 : (pn.Y - _p0.Y) / d; //0度會讓點被忽略了
return angle;
}
/// <summary>
/// 求凸包_葛立恆算法,出來的凸包做的多段線在正交的情況下會多點或者少點
/// </summary>
/// <param name="pts"></param>
/// <returns></returns>
Point3d[] GrahamConvexHull(IEnumerable<Point3d> pt2ds)
{
//消重,點排序
var pts = pt2ds.Distinct(new ToleranceDistinct()).OrderBy(p => p.X).ThenBy(p => p.Y).ToList();
//max右上角,因為負角度的問題,所以需要從右上角計算
_p0 = pts.Last();
//按角度及距離排序
pts = pts.OrderByDescending(p => Cosine(p)).ThenBy(p => _p0.DistanceTo(p)).ToList();
var stack = new Stack<Point3d>();
stack.Push(_p0);//頂部加入對象
stack.Push(pts[1]);
bool tf = true;
//遍歷所有的點,因為已經角度順序,所有是有序遍歷.從第三個點開始
for (int i = 2; i < pts.Count; i++)
{
Point3d qn = pts[i]; //第一次為p2,相當於pn
Point3d q1 = stack.Pop(); //第一次為p1,相當於前一個點,刪除頂部對象(相當於點回退)
Point3d q0 = stack.Peek();//第一次為_p0,相當於后面一個點,查詢頂部對象
//為真表示要剔除
while (tf && CrossAclockwise(q1, q0, qn))
{
if (stack.Count > 1)
{
stack.Pop();//刪除頂部對象(相當於刪除前一個點進行回退)
//前后點交換,用於while循環,
//可參考 https://www.bilibili.com/video/BV1v741197YM 04:15
//棧頂就是回滾之后的,再次和qn進行向量叉乘,看看是不是叉乘方向,是就繼續回滾
//否則結束循環后加入棧中.
q1 = q0;
q0 = stack.Peek();
}
else
{
//棧少於1,就不在剔除頂部.結束循環...
//保護棧中_p0不剔除
tf = false;
}
}
stack.Push(q1);
stack.Push(qn);
tf = true;
}
var npts = stack.ToList();
//過濾凸度過少的話,將點移除,以免凸包有多余的邊點.
for (int i = 0; i < npts.Count() - 2; i++)
{
var bu = GetArcBulge(npts[i], npts[i + 1], npts[i + 2]);
if (Abs(bu) < 1e-6)
{
npts.RemoveAt(i + 1);//移除中間
i--;
}
}
return npts.ToArray();
}
}
}
子函數:
cad曲線采樣
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.Geometry;
using System;
using System.Collections.Generic;
namespace JoinBox
{
public class CurveSplit
{
public List<Curve> Curves { get; set; }
/// <summary>
/// 獲取定值分割的曲線集合
/// </summary>
/// <param name="curve">曲線</param>
/// <param name="fixedValue">定值分割</param>
public CurveSplit(Curve curve, double fixedValue)
{
Curves = new();
//算曲線長度
double curveLength = curve.GetLength();
//若少於定值,則直接返回這條曲線,表示返回這段長度
if (curveLength < fixedValue)
{
Curves.Add(curve);
return;
}
var pts = new Point3dCollection();
//用來疊加長度
double overlyingLength = 0;
//定值采集點
while (overlyingLength < curveLength)
{
//求起點到長度的點
pts.Add(curve.GetPointAtDist(overlyingLength));
overlyingLength += fixedValue;
}
//最后沒有完全重合,加入尾巴點
if (overlyingLength - curveLength < 1e-6)
pts.Add(curve.GetPointAtDist(curveLength));
//通過點集,分割曲線
var splits = curve.GetSplitCurves(pts.Collection);
foreach (var item in splits)
{
var cuItem = (Curve)item;
Curves.Add(cuItem);
}
pts.Dispose();
splits.Dispose();//釋放
}
/// <summary>
/// 手動釋放生成出來的曲線,
/// 因為cad的Point3d沒有繼承,所以不能用 <see cref="IDisposable">進行釋放</see>
/// 否則提示:Forgot to call Dispose? (Autodesk.AutoCAD.DatabaseServices.Arc): DisposableWrapper
/// </summary>
public void Dispose()
{
Curves?.ForEach(cu => {
cu.Dispose();
});
}
}
public class CurveSample
{
Curve _curve { get; set; }
int _numSample { get; set; }
/// <summary>
/// 求采樣
/// </summary>
/// <param name="curve">曲線</param>
/// <param name="sampleNum">采樣份數</param>
public CurveSample(Curve curve, int sampleNum = 256)
{
_curve = curve;
_numSample = sampleNum;
}
/// <summary>
/// 曲線采樣(注意尾點是否缺少哦)
/// </summary>
/// <returns></returns>
public IEnumerable<Point3d> GetSamplePoints
{
get
{
if (_numSample == 0)
throw new System.Exception("NumSample參數不能為0");
var length = _curve.GetLength();
var fixedValue = length / _numSample;
var cs = new CurveSplit(_curve, fixedValue);
var curves = cs.Curves;
var pts = new List<Point3d>();
pts.Add(curves[0].StartPoint);//起點
foreach (var item in curves)
pts.Add(item.EndPoint);//間隔點,尾點
//末尾兩個點可能一樣,需要判斷去除
if (pts[pts.Count - 1] == pts[pts.Count - 2])
pts.RemoveAt(pts.Count - 1);
cs.Dispose();
return pts;
}
}
}
}
曲線長度
namespace JoinBox
{
public static partial class MathHelper
{
/// <summary>
/// 曲線長度
/// </summary>
/// <param name="curve">曲線</param>
/// <returns></returns>
public static double GetLength(this Curve curve)
{
return curve.GetDistanceAtParameter(curve.EndParam);
}
}
}
包圍盒
參考
因為沒看懂的關系,所以我自己想了以下的粗糙方法實現....畢竟我真的不想看cpp....
我還發現游戲行業還有快速平方根,模糊距離,快速包圍盒等等的實現......而cad只能用精確的包圍盒 😃
定義
包圍盒又叫外接矩形,如下圖形式:
先了解AABB叫軸向包圍盒,這是cad自帶函數可以求得的.(求最小點和最大點).
再是今次的主角OBB有向包圍盒.
包圍盒分二維和三維,由於三維在我的領域內沒啥用途,所以我只實現了二維包圍盒,
當然了,如果你已經掌握點乘和叉乘的要領,那么你可以根據二維來寫出三維來.
條件
條件1: 包圍盒不是出現在點集形成的圖形上,而是出現在點集的凸包上,看下圖可知:
條件2:包圍盒的矩形邊必然和凸包其中一條邊重合,所以利用凸包,就可以實現面積最小的包圍盒了.
實現方式
如圖所示:
先求矩形的第一個角點r1,沿着a->b向量(無限長的線),角點r1會是c點落在這a->b向量上,得到a到r1的距離長度.
再來d點可能比c點更遠,e點比d點更遠......於是乎需要循環......
當循環出現了本次距離比上次小,說明找到了最后的r1角點了......
用一些數學語言描述一下:
矩形角點r1==a->c點乘a->b得到r1角點,再a->d點乘a->b...輪詢,如果r1與a點距離開始進行縮小(點乘距離回歸),
那么表示r1確定,以及得到末尾點ptLast=c點,
最后明確得到:a->ptLast點乘a->b就是最近角點r1.
求剩余的點,如下圖所示:
矩形角點r2,r3,r4僅僅需要旋轉點集得到,
通過a->b與X軸弧度(角度),繞r1旋轉,然后進行坐標是數值交換即可得到一個ab段包圍盒矩形.
其余包圍盒:
現在知道了ab段,接着循環bc,cd,de等等....即可得到每個包圍矩形,然后就可以求得最大或者最小的包圍盒.
嘻嘻.說明到此結束.
代碼
public partial class Graham
{
/// <summary>
/// 有向包圍盒
/// </summary>
/// <param name="pts">點集</param>
/// <returns>返回每條邊的包圍盒</returns>
List<Rectangular> Boundingbox(List<Point2d> pts)
{
/*
最小包圍盒(外接矩形)
重要的條件:凸多邊形的最小周長(最小面積)外接矩形存在一條邊和多邊形的一條邊重合
角點r1==a->c叉乘a->b得到r1角點,再a->d叉乘a->b...輪詢,如果r1與a點距離開始進行縮小(叉乘距離回歸),
那么表示r1確定,以及得到末尾點ptLast圖上位c點,
最后明確得到:a->ptLast叉乘a->b就是最近角點r1
角點r2,r3,r4需要的僅僅是通過a->b向量與X軸角度繞r1旋轉得到
*/
var recs = new List<Rectangular>();
//因為利用三點獲取第一個交點,
//所以為了最后的邊界,需要加入集合前面,以使得成為一個環
pts.Add(pts[0]);
pts.Add(pts[1]);
pts.Add(pts[2]);
pts.Add(pts[3]);
//此處循環是作為邊,下次則是ab..bc..cd..de..循環下去
for (int i = 0; i < pts.Count() - 1; i++)
{
//矩形的點
Point2d r1;
Point2d r2;
Point2d r3;
Point2d r4;
//最后一次的點
Point2d c_point;
//上一次長度
double ac_lengthLast = -1;
Point2d a_point = pts[i];
Point2d b_point = pts[i + 1];
//此處循環是求長度求r1點,如果出現點乘距離收縮,就結束
for (int c_nmuber = i + 2; c_nmuber < pts.Count(); c_nmuber++)
{
//ac->ab點乘求矩形的一個角點
//角點距離如果出現縮小,就確定了這個點是j
var ac_length = a_point.DistanceTo(DotProduct(a_point, pts[c_nmuber], b_point));
//第一次賦值
if (ac_lengthLast == -1)
{
ac_lengthLast = ac_length;
}
else if (ac_lengthLast < ac_length)
{
ac_lengthLast = ac_length;
}
else
{
//此處條件是點乘距離已經收縮,求得r1點,最后會break結束循環.
c_point = pts[c_nmuber - 1]; //邊界點-1就是上次的
r1 = DotProduct(a_point, c_point, b_point);//角點計算
{
//根據角度旋轉所有的點
var v1 = r1.GetVectorTo(a_point);
var angle1 = v1.GetAngleTo(Vector2d.XAxis);
var angle = v1.GetAngle2XAxis();
//此處循環是求r2,r3,r4的點
//后面的點旋轉之后加入集合,再利用linq判斷Y軸X軸最大的
var houmiandesuoyoudian = new List<Point2d>();
foreach (var pt in pts)
{
houmiandesuoyoudian.Add(pt.RotateBy(-angle, r1));
}
var maxY = houmiandesuoyoudian.OrderByDescending(pt => pt.Y).ToList()[0].Y;
var maxX = houmiandesuoyoudian.OrderByDescending(pt => pt.X).ToList()[0].X;
//通過最大,計算角點,然后逆旋轉,回到原始圖形
r2 = new Point2d(r1.X, maxY).RotateBy(angle, r1);
r3 = new Point2d(maxX, maxY).RotateBy(angle, r1);
r4 = new Point2d(maxX, r1.Y).RotateBy(angle, r1);
recs.Add(new Rectangular(r1, r2, r3, r4));
}
break;
}
}
}
return recs;
}
/// <summary>
/// 面積最小包圍盒
/// </summary>
/// <param name="pts">點集</param>
/// <returns></returns>
Rectangular BoundingboxAreaMin(List<Point2d> pts)
{
var recs = Boundingbox(pts);
return recs.OrderBy(rec => rec.Area).ToArray()[0];
}
// 概念參考了這里,但是它代碼好像有點問題 http://www.cppblog.com/staryjy/archive/2009/11/19/101412.html
/// <summary>
/// 凸包對踵點最大的距離_旋轉卡殼
/// </summary>
/// <returns></returns>
double RotateCalipersMax(IEnumerable<Point2d> pt2ds)
{
var pts = pt2ds.ToList();
if (pts.Count == 0 || pts.Count == 1)
{
return 0;
}
else if (pts.Count == 2)
{
return pts[0].DistanceTo(pts[1]);
}
//返回的長度
double ans = 0;
int ps = 2; //2是從下下一個點開始,因為1開始會導致[0]和[0]判斷距離
pts.Add(pts[0]);//但是下下一個點開始就表示沒有判斷到0->1線段,必須加入尾部判斷
int p = pts.Count - ps;
for (int i = 0; i < p; i++) //點序是叉乘方向的.
{
//叉乘求面積,面積呈現為單峰函數(函數圖像中間大兩邊小,函數從遞增到遞減),
//滿足結束條件:面積最高峰的時候下一次判斷即為>前者(函數遞減) || 取余為0(即為遍歷到最后了)
//叉乘求面積,A<B,表示求后者面積數大,直到后者面積數小,結束循環,求最大值長度即可
while (Abs(Cross(pts[i], pts[i + 1], pts[ps])) <
Abs(Cross(pts[i], pts[i + 1], pts[ps + 1])))
{
ps = (ps + 1) % p;//X%Y,X<Y返回X.取余為0(即為遍歷到最后了)
}
//峰值時候求的三點距離
//第一次3->0和3->1
ans = Max(ans, Max(pts[ps].DistanceTo(pts[i]),
pts[ps].DistanceTo(pts[i + 1])
));
}
return ans;
}
/// <summary>
/// 凸包對踵點最小的距離_旋轉卡殼
/// </summary>
/// <returns></returns>
double RotateCalipersMin(IEnumerable<Point2d> pt2ds)
{
var pts = pt2ds.ToList();
if (pts.Count == 0 || pts.Count == 1)
{
return 0;
}
else if (pts.Count == 2)
{
return pts[0].DistanceTo(pts[1]);
}
var lstmin = new List<double>();
//點集順序是叉乘方向的.
//凸包對踵點最小的距離==非鄰點的最小點距,鄰邊點不計算
//計算方式是通過獲取循環數字,確定執行次數.
//從i從i+2點開始遞增,但是循環數字是遞減的,只因不需要重復計算0-2,2-0的情況.
//所以時間復雜度是常數項
var last = pts.Count - 1;
for (int i = 0; i < last; i++)
{
//循環次數 = 總數 - 1(前一個點鄰邊) - i(循環遞減)
int fornum = last - 1 - i;
//如果i==0,減多1次
if (i == 0)
{
fornum--;
}
int inumq = i + 1;//前一個點(鄰邊)
for (int j = 0; j < fornum; j++)
{
inumq++;//前前一個點
var dis = pts[i].DistanceTo(pts[inumq]);
lstmin.Add(dis);
}
}
return lstmin.Min(); //返回的長度
}
class Rectangular
{
public Point2d R1;
public Point2d R2;
public Point2d R3;
public Point2d R4;
/// <summary>
/// 矩形類
/// </summary>
/// <param name="r1">矩形的原點,依照它來旋轉</param>
/// <param name="r2"></param>
/// <param name="r3"></param>
/// <param name="r4"></param>
public Rectangular(Point2d r1, Point2d r2, Point2d r3, Point2d r4)
{
R1 = r1;
R2 = r2;
R3 = r3;
R4 = r4;
R4 = r4;
}
/// <summary>
/// 面積
/// </summary>
public double Area
{
get
{
var x = R1.DistanceTo(R4);
var y = R1.DistanceTo(R2);
return x * y;
}
}
}
}
(完)