凸包
参考
安德鲁算法
分治法(其中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;
}
}
}
}
(完)