凸包
参考
分治法(其中nfox的项目实现的是分治法)
定义
凸包又叫凸多边形,本篇文章可能混用两种说法,
形象的理解就是一些点(点集)用一根橡皮筋紧紧地包裹外边点.
如果知道了这个定义,那么还有:
用一个保鲜膜裹着三维点,求膜上点集.
用一个最小的球裹着三维点,求球球的中心点和直径.
这样就进入了一个叫拓扑学的学科上.......我的老天鹅.
我竟然搜了半天没看到可以直接拿来用的c#代码,还是小轩轩给我的....
葛立恒凸包
注意一下,如果点集形成的是正交矩形的情况时,
算出来的凸包会多一个点,可以进行后处理.
(你会发现代码实际上是右上角最大点开始的,其他的教程说明从最小点开始算,
这需要知道的是最小最大都是边界点,边界点必然是凸包的边点,用谁起算都可以)
运行命令test_gra
#if !HC2020
usingAutodesk.AutoCAD.ApplicationServices;usingAutodesk.AutoCAD.DatabaseServices;usingAutodesk.AutoCAD.EditorInput;usingAutodesk.AutoCAD.Geometry;usingAutodesk.AutoCAD.Runtime;#else
usingGrxCAD.DatabaseServices;usingGrxCAD.EditorInput;usingGrxCAD.Geometry;usingGrxCAD.ApplicationServices;usingGrxCAD.Runtime;usingGrxCAD.Colors;usingGrxCAD.GraphicsInterface;using Viewport =GrxCAD.DatabaseServices.Viewport;#endif
usingSystem.Collections.Generic;usingSystem.Linq;
namespaceJoinBox.src.数学操作
{public classGraham
{
///
///最靠近x轴的点(必然是凸包边界的点)///
privatePoint2d _p0;///
///角度p0和pn的角度///
///
///
private doubleCosine(Point2d pn)
{double d =_p0.GetDistanceTo(pn);//距离是相同返回1.0表示true:那么求角度(高/斜)==sin(角度)
return d == 0.0 ? 1.0 : (pn.Y - _p0.Y) /d;
}///
///叉乘,顺时针方向为真,表示要剔除///
///
///
///
///
private boolClockwise(Point2d n, Point2d a, Point2d b)
{return ((a.X - b.X) * (a.Y - n.Y) - (a.X - n.X) * (a.Y - b.Y)) > -1e-6;
}///
///葛立恒求凸包///
///
///
private Point2d[] ConvexHull(Listpts)
{
//消重
pts = pts.Distinct().ToList();//靠近x轴的点
_p0 = pts.OrderBy(p => p.X).ThenBy(p =>p.Y).Last();//按角度及距离排序
pts = pts.OrderByDescending(p => Cosine(p)).ThenBy(p =>_p0.GetDistanceTo(p)).ToList();var stack = new Stack();
stack.Push(_p0);//顶部加入对象
stack.Push(pts[1]);bool tf = true;//遍历所有的点,因为已经角度顺序,所有是有序遍历.从第三个点开始
for (int i = 2; i < pts.Count; i++)
{
Point2d qn= pts[i]; //第一次为p2,相当于pn
Point2d q1 = stack.Pop(); //第一次为p1,相当于前一个点,删除顶部对象(相当于点回退)
Point2d q0 = stack.Peek();//第一次为_p0,相当于后面一个点,查询顶部对象
while (tf && Clockwise(qn, q1, q0))//顺时针方向为真,表示要剔除
{if (stack.Count > 1)//保护栈中_p0不剔除
{
stack.Pop();//删除顶部对象(相当于删除前一个点进行回退)//前后点交换,用于while循环,//可参考https://www.bilibili.com/video/BV1v741197YM04:15//栈顶就是回滚之后的,再次和qn进行向量叉乘,看看是不是顺时针,是就继续回滚..//否则结束循环后加入栈中.
q1 =q0;
q0=stack.Peek();
}else{
tf= false;
}
}
stack.Push(q1);
stack.Push(qn);
tf= true;
}returnstack.ToArray();
}///
///葛立恒求法///
[CommandMethod("Test_Gra")]public voidTest_Gra()
{
Document doc=Application.DocumentManager.MdiActiveDocument;
Editor ed=doc.Editor;
Database db= doc.Database;//当前的数据库
ed.WriteMessage("n****{惊惊盒子}葛立恒求凸包,选择曲线:");//定义选择集选项
var pso = newPromptSelectionOptions
{
RejectObjectsOnLockedLayers= true, //不选择锁定图层对象
AllowDuplicates = true, //不允许重复选择
};var ssPsr = ed.GetSelection(pso);//手选 这里输入al会变成all,无法删除ssget的all关键字
if (ssPsr.Status !=PromptStatus.OK)
{return;
}
db.Action(tr=>{var pts = new List();foreach (ObjectId id inssPsr.Value.GetObjectIds())
{var ent =id.ToEntity(tr);if (ent isCurve curve)
{var cs = new CurveSample(curve, 300);
pts.AddRange(cs.GetSamplePoints);
}else if (ent isDBPoint dbPt)
{
pts.Add(newPoint2d(dbPt.Position.X, dbPt.Position.Y));
}
}
Point2d[] npts=ConvexHull(pts);var bv = new List();for (int i = 0; i < npts.Length; i++)
{
bv.Add(new BulgeVertex(npts[i], 0));
}
Entity pl= EntityAdd.AddPolyLineToEntity(0, bv);
tr.AddEntityToMsPs(db, pl);
});
}
}
}
View Code
cad曲线采样
namespaceJoinBox.src.数学操作
{public class CurveSplit where TPoint : struct, IFormattable
{
Curve _curve {get; set; }double _fixedValue { get; set; }///
///求间隔点///
/// 曲线
/// 定值分割
public CurveSplit(Curve curve, doublefixedValue)
{
_curve=curve;
_fixedValue=fixedValue;
}///
///获取定值分割的曲线集合///
///
public ListGetSplitCurves
{get{var lstCurves = new List();//算曲线长度
double totalLength =_curve.GetLength();//若少于定值,则直接返回这长度
if (totalLength <_fixedvalue>
{
lstCurves.Add(_curve);returnlstCurves;
}double addLength = 0;var pt3dCol = newPoint3dCollection();//这段代码通过定值采集点
while (addLength
{//求起点到长度的点
pt3dCol.Add(_curve.GetPointAtDist(addLength));
addLength+=_fixedValue;
}if (addLength !=totalLength)
{
pt3dCol.Add(_curve.GetPointAtDist(totalLength));
}//通过点集,分割曲线,可能有精度问题......
var splits =_curve.GetSplitCurves(pt3dCol);foreach (var item insplits)
{
lstCurves.Add((Curve)item);
}
splits.Dispose();//释放
returnlstCurves;
}
}
}public class CurveSample where TPoint : struct, IFormattable
{
Curve _curve {get; set; }int _numSample { get; set; } = 1;///
///求采样///
/// 曲线
/// 采样份数
public CurveSample(Curve curve, intnumSample)
{
_curve=curve;
_numSample=numSample;
}///
///曲线采样///
///
public ListGetSamplePoints
{get{if (_numSample == 0)
{throw new System.Exception("NumSample参数不能为0");
}// https://www.cnblogs.com/luludongxu/p/5669729.html
//泛型构造传参//泛型是point2d时候参数是XY,我传入的obj数组有3个值,竟然能自动忽略到末尾的z
Type T = typeof(TPoint);var length =_curve.GetLength();var fixedValue = length /_numSample;var cs = new CurveSplit(_curve, fixedValue);var spls =cs.GetSplitCurves;var pts = new List();
pts.Add((TPoint)Activator.CreateInstance(T, spls[0].StartPoint.ToArray()));//起点
foreach (var item inspls)
{
pts.Add((TPoint)Activator.CreateInstance(T, item.EndPoint.ToArray()));//间隔点,尾点
}returnpts;
}
}
}
}
View Code
包围盒
参考
因为没看懂的关系,所以我自己想了以下的粗糙方法实现....毕竟我真的看不懂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 static List Boundingbox(Listpts)
{
var recs = new List();//因为利用三点获取第一个交点,//所以为了最后的边界,需要加入集合前面,以使得成为一个环
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.GetDistanceTo(DotProduct(a_point, pts[c_nmuber], b_point));//第一次赋值
if (ac_lengthLast == -1)
{
ac_lengthLast=ac_length;
}else if (ac_lengthLast
{
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();foreach (var pt inpts)
{
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 = newPoint2d(r1.X, maxY).RotateBy(angle, r1);
r3= newPoint2d(maxX, maxY).RotateBy(angle, r1);
r4= newPoint2d(maxX, r1.Y).RotateBy(angle, r1);
recs.Add(newRectangular(r1, r2, r3, r4));
}break;
}
}
}returnrecs;
}///
///面积最小包围盒///
/// 点集
///
public static Rectangular BoundingboxAreaMin(Listpts)
{var recs =Boundingbox(pts);return recs.OrderBy(rec => rec.Area).ToArray()[0];
}public classRectangular
{publicPoint2d R1;publicPoint2d R2;publicPoint2d R3;publicPoint2d R4;///
///矩形类///
/// 矩形的原点,依照它来旋转
///
///
///
publicRectangular(Point2d r1, Point2d r2, Point2d r3, Point2d r4)
{
R1=r1;
R2=r2;
R3=r3;
R4=r4;
R4=r4;
}///
///面积///
public doubleArea
{get{var x =R1.GetDistanceTo(R4);var y =R1.GetDistanceTo(R2);return x *y;
}
}
}
View Code
子函数
2: X轴到向量的弧度,cad的获取的弧度是1PI,所以转换为2PI(上小,下大)
///
///X轴到向量的弧度,cad的获取的弧度是1PI,所以转换为2PI(上小,下大)///
/// 向量
/// 弧度
public static double GetAngle2XAxis(thisVector3d ve)
{//世界重合到用户 Vector3d.XAxis->两点向量
double alz = Vector3d.XAxis.GetAngleTo(ve);//观察方向不要设置
alz = ve.Y > 0 ? alz : Pi2 - alz; //逆时针为正,如果-负值控制正反
alz = Math.Abs(Pi2 - alz) < 1e-10 ? 0: alz;returnalz;
}
View Code
(完)
