数学篇 cad.net 葛立恒凸包算法和面积最小包围盒
凸包
参考
安德鲁算法
分治法(其中nfox的项目实现的是分治法)
多边形快速凸包算法(Melkman‘s Algorithm)
还可以这看cpp的代码: https://www.cnblogs.com/VividBinGo/p/11637684.html
定义
凸包又叫凸多边形,本篇文章可能混用两种说法,形象的理解就是一些点(点集)用一根橡皮筋紧紧地包裹外边点.
如果知道了这个定义,那么还有:
用一个保鲜膜裹着三维点,求膜上点集.
用一个最小的球裹着三维点,求球球的中心点和直径.
这样就进入了一个叫拓扑学的学科上.......我的老天鹅.
我竟然搜了半天没看到可以直接拿来用的c#代码,还是小轩轩给我的....
葛立恒凸包
注意一下,如果点集形成的是正交矩形的情况时,
算出来的凸包会多一个点,可以进行后处理.
(你会发现代码实际上是右上角最大点开始的,其他的教程说明从最小点开始算,
这需要知道的是最小最大都是边界点,边界点必然是凸包的边点,用谁起算都可以)
包围盒
参考
因为没看懂的关系,所以我自己想了以下的粗糙方法实现
我还发现游戏行业还有快速平方根,模糊距离,快速包围盒等等的实现......而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.
得到轴向矩形
矩形其他角点只需要旋转点集得到,这样就是一个轴向矩形.
其余包围盒:
现在知道了ab段,接着循环bc,cd,de等等....即可得到每个包围矩形,然后就可以求得面积最大或者最小的包围盒.
代码
bug致谢: JieGou(杰狗)
子函数
db.Action委托无返回值写法
入门数学 PointV+凸度
曲线采样
曲线长度
RectAngle类
GetEntityPoint3ds
主函数
#if !HC2020
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.Runtime;
using Acap = Autodesk.AutoCAD.ApplicationServices.Application;
#else
using GrxCAD.DatabaseServices;
using GrxCAD.EditorInput;
using GrxCAD.Geometry;
using GrxCAD.Runtime;
using Acap = GrxCAD.ApplicationServices.Application;
#endif
using System;
using System.Collections.Generic;
using System.Linq;
using JoinBox.BasalMath;
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 class CmdTestGrahamClass
{
/// <summary>
/// 求凸包测试命令
/// </summary>
[CommandMethod("CmdTest_Graham")]
public void CmdTest_Graham()
{
var dm = Acap.DocumentManager;
var doc = dm.MdiActiveDocument;
var db = doc.Database;
var ed = doc.Editor;
ed.WriteMessage("\n****{惊惊连盒}求凸包,选择曲线:");
//定义选择集选项
var pso = new PromptSelectionOptions
{
RejectObjectsOnLockedLayers = true, //不选择锁定图层对象
AllowDuplicates = true, //不允许重复选择
};
var ssPsr = ed.GetSelection(pso);
if (ssPsr.Status != PromptStatus.OK)
return;
db.Action(tr => {
var getPts = new List<PointV>();
var ents = ssPsr.Value.GetObjectIds().ToEntity(tr);
ents.ForEach(ent => {
if (ent is Line line)
{
getPts.Add(line.StartPoint);
getPts.Add(line.EndPoint);
}
else if (ent is Polyline pl)
{
pl.GetEntityPoint3ds().ForEach(pt => { getPts.Add(pt); });
}
else if (ent is Curve curve)
{
var cs = new CurveSample(curve);
cs.GetSamplePoints.ForEach(pt => { getPts.Add(pt); });
}
else if (ent is DBPoint bPoint)
getPts.Add(bPoint.Position);
else
{
var entPosition = ent.GetType().GetProperty("Position");//反射获取属性
if (entPosition is not null)
{
var pt = (Point3d)entPosition.GetValue(null, null);
getPts.Add(pt);
}
}
});
//葛立恒方法
var ptList = GrahamConvexHull(getPts);
if (ptList == null)
return;
ed.WriteMessage("\n\r凸包对踵点最大距离:" + RotateCalipersMax(ptList));
ed.WriteMessage("\n\r凸包对踵点最小距离:" + RotateCalipersMin(ptList));
if (ptList.Length > 1)
{
var disList = new List<double>();
for (int i = 0; i < ptList.Length - 1; i++)
disList.Add(ptList[i].DistanceTo(ptList[i + 1]));
disList.Add(ptList[ptList.Length - 1].DistanceTo(ptList[0]));
ed.WriteMessage("\n\r凸包点集最大距离:" + disList.Max());
ed.WriteMessage("\n\r凸包点集最小距离:" + disList.Min());
}
//生成凸包边界
var bv = new List<BulgeVertexWidth>();
for (int i = 0; i < ptList.Length; i++)
bv.Add(new BulgeVertexWidth(ptList[i].ToPoint2d(), 0));
var pl = EntityAdd.AddPolyLineToEntity(bv);
tr.AddEntityToMsPs(db, pl);
var rects = Boundingbox(ptList);
//Test: 求凸包,修改数值以便测试
int creatAll = 0;
if (creatAll == 0)
{
//生成所有包围盒
int ColorIndex = 0;
foreach (var rec in rects)
{
var ent = rec.ToPolyLine();
ent.EntityRotate(rec.Angle, Vector3d.ZAxis, rec.MinPoint.ToPoint3d());
ent.ColorIndex = ++ColorIndex;
tr.AddEntityToMsPs(db, ent);
}
}
else
{
//生成计算面积最小的包围盒
var rec = rects.OrderBy(item => item.Area).First();
var ent = rec.ToPolyLine();
ent.EntityRotate(rec.Angle, Vector3d.ZAxis, rec.MinPoint.ToPoint3d());
tr.AddEntityToMsPs(db, ent);
}
});
}
/// <summary>
/// 求凸包_葛立恒算法,
/// 凸包多段线在正交的情况下会多一个点
/// </summary>
/// <param name="points"></param>
/// <returns>凸包点集(顺时针)</returns>
PointV[]? GrahamConvexHull(IEnumerable<PointV> points)
{
if (points == null)
throw new ArgumentNullException(nameof(points));
//最小最大都是边界点,边界点必然是凸包的边点
//因为Cosine排序方式是负角度,所以从右上角计算
PointV maxPt;//p0
//消重,排序,此处排序将导致点集从顺时针开始
var ptsOrder = points.Distinct(PointV.Distinct);
//按角度排序,同角度用距离排序
ptsOrder = ptsOrder.OrderBy(p => p.X).ThenBy(p => p.Y);
maxPt = ptsOrder.Last();
ptsOrder = ptsOrder.OrderByDescending(pn => Cosine(maxPt, pn)).ThenBy(p => maxPt.DistanceTo(p));
var pts = ptsOrder.ToArray();
bool whileBreak = false;
var stack = new Stack<PointV>();
stack.Push(maxPt);
stack.Push(pts[1]);
//从大到小遍历,第三个点开始
for (int i = 2; i < pts.Length; i++)
{
var pn = pts[i]; //第一次为 p2 ,相当于pn
var p1 = stack.Pop(); //第一次为 p1 ,相当于前一个点(返回并剔除顶部对象)
var p0 = stack.Peek();//第一次为 maxPt ,相当于后面一个点(查询顶部对象)
//为真表示要剔除
while (!whileBreak && p1.CrossAclockwise(p0, pn))
{
if (stack.Count > 1)
{
//删除顶部对象
stack.Pop();
//前后点交换,用于while循环,
//可参考 https://www.bilibili.com/video/BV1v741197YM 04:15
//栈顶就是回滚之后的,再次和qn进行向量叉乘,看看是不是叉乘方向,
//是就继续回滚,否则结束循环后加入栈中.
p1 = p0;
p0 = stack.Peek();
}
else
{
//栈少于1,就不再剔除顶部(保护栈中 maxPt 不剔除),结束循环
whileBreak = true;
}
}
stack.Push(p1);
stack.Push(pn);
whileBreak = false;
}
return RemoveBulgeSamll(stack);
}
/// <summary>
/// 剔除凸度太小的,以免凸包有多余的边点
/// </summary>
/// <param name="stack"></param>
/// <returns></returns>
PointV[]? RemoveBulgeSamll(IEnumerable<PointV> stack)
{
if (stack == null)
throw new ArgumentNullException(nameof(stack));
//旧方案没有考虑(尾-头0-头1)这样的情况,现在改过来
//可以用LinkedList代替
var link = new LoopList<PointV>(stack);
var node = link.First;
if (node == null)
return null;
for (int i = 0; i < link.Count; i++)
{
var bulge = MathHelper.GetArcBulge(
node!.Value.ToPoint2d(),
node!.Next!.Value.ToPoint2d(),
node!.Next!.Next!.Value.ToPoint2d());
if (Math.Abs(bulge) < 1e-6)
link.Remove(node.Next);//剔除中间
node = node.Next;
}
return link.ToArray();
}
/// <summary>
/// 角度p0和pn的角度
/// </summary>
/// <param name="p0">最靠近x轴的点(必然是凸包边界的点)</param>
/// <param name="pn"></param>
/// <returns></returns>
double Cosine(PointV p0, PointV pn)
{
double dis = p0.DistanceTo(pn);
//求角度(高/斜)==sin(角度)
//距离是0表示是自己和自己的距离,那么0不可以是除数,否则Nan
//如果1.0为sin周期函数最大,改成0则会令0度的点被忽略了
return dis == 0.0 ? 1.0 : (pn.Y - p0.Y) / dis;
}
/// <summary>
/// 有向包围盒
/// </summary>
/// <param name="ptsOrd">凸包点集(顺时针)</param>
/// <returns>返回每条边的包围盒</returns>
List<RectAngle> Boundingbox(IEnumerable<PointV> ptsOrd)
{
/*
* 最小包围盒(外接矩形)
* 重要的条件:凸多边形的最小周长(最小面积)外接矩形存在一条边和多边形的一条边重合;
* 角点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 ptList = ptsOrd.ToList();
var rects = new List<RectAngle>();
//因为利用三点获取第一个交点,
//所以为了最后的边界,需要加入集合前面,以使得成为一个环
ptList.Add(ptList[0]);
ptList.Add(ptList[1]);
ptList.Add(ptList[2]);
//边循环,ab..bc..cd..de..
int numB = 1;
for (int i = 0; i < ptList.Count - numB; i++)
{
double acLengthLast = -1;
var ptA = ptList[i];
var ptB = ptList[i + numB];
//此处循环是求长度求r1点,如果出现点乘距离收缩,就结束
for (int c_nmuber = i + 2; c_nmuber < ptList.Count; c_nmuber++)
{
//ac->ab点乘求矩形的一个角点(绕轴转点)
//角点距离如果出现缩小,就确定了这个点是r1
var acLength = ptA.DistanceTo(ptA.DotProduct(ptList[c_nmuber], ptB));
//第一次赋值
if (acLengthLast == -1)
{
acLengthLast = acLength;
}
else if (acLengthLast < acLength)
{
acLengthLast = acLength;
}
else
{
//此处点乘距离已经收缩,求得绕轴转点r1(也是矩形角点),结束循环.
PointV r1;
var ptCLast = ptList[c_nmuber - 1]; //边界点-1就是上次的
var ptC2 = ptA.DotProduct(ptCLast, ptB); //角点计算
if (ptA.DistanceTo(ptC2) < ptA.GetDistanceTo(ptB))
r1 = ptB;//叉乘点出现在ab线内,此时b才是矩形边点
else
r1 = ptC2;
rects.Add(CreatRect(ptList, ptA, r1));
break;
}
}
}
return rects;
}
/// <summary>
/// 得到一个轴向矩形和夹角
/// </summary>
/// <param name="pts"></param>
/// <param name="a_point"></param>
/// <param name="r1">绕轴转点</param>
/// <returns></returns>
RectAngle CreatRect(List<PointV> pts, PointV a_point, PointV r1)
{
var angle = r1.GetVectorTo(a_point).GetAngle2XAxis();
//旋转点后求一个轴向矩形
var ptsRo = new List<PointV>();
pts.ForEach(pt => {
ptsRo.Add(pt.RotateBy(-angle, PointV.ZAxis, r1));
});
var ox = ptsRo.OrderBy(pt => pt.X);
var oy = ptsRo.OrderBy(pt => pt.Y);
var minX = ox.First().X;
var minY = oy.First().Y;
var maxX = ox.Last().X;
var maxY = oy.Last().Y;
//RectAngle的旋转点是min点
return new RectAngle(new Rect(minX, minY, maxX, maxY), angle);
}
// 概念参考了这里,但是它代码好像有点问题
// http://www.cppblog.com/staryjy/archive/2009/11/19/101412.html
/// <summary>
/// 凸包对踵点最大的距离_旋转卡壳
/// </summary>
/// <returns></returns>
double RotateCalipersMax(IEnumerable<PointV> ptsOrd)
{
var pts = ptsOrd.ToList();
if (pts.Count < 2)
return 0;
else if (pts.Count == 2)
return pts[0].DistanceTo(pts[1]);
double result = 0;
int ps = 2; //2是从下下一个点开始,因为1开始会导致[0]和[0]判断距离
pts.Add(pts[0]);//但是下下一个点开始就表示没有判断到0->1线段,必须加入尾部判断
int first = pts.Count - ps;
for (int i = 0; i < first; i++) //点序是叉乘方向的.
{
//叉乘求面积:
//面积呈现为单峰函数(函数图像中间大两边小,函数从递增到递减),
//满足结束条件:
//面积最高峰的时候下一次判断即为>前者(函数递减) || 取余为0(即为遍历到最后了)
//叉乘求面积,A<B,表示求后者面积数大,直到后者面积数小,结束循环,求最大值长度即可
while (Math.Abs(pts[i].Cross(pts[i + 1], pts[ps])) <
Math.Abs(pts[i].Cross(pts[i + 1], pts[ps + 1])))
ps = (ps + 1) % first;//X%Y,X<Y返回X.取余为0(即为遍历到最后了)
//峰值时候求的三点距离
//第一次3->0和3->1
result = Math.Max(result, Math.Max(pts[ps].DistanceTo(pts[i]), pts[ps].DistanceTo(pts[i + 1])));
}
return result;
}
/// <summary>
/// 凸包对踵点最小的距离_旋转卡壳
/// </summary>
/// <returns></returns>
double RotateCalipersMin(IEnumerable<PointV> ptsOrd)
{
var pts = ptsOrd.ToList();
if (pts.Count < 2)
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(); //返回的长度
}
}
}
(完)
