数学篇 cad.net 葛立恒凸包算法和面积最小包围盒

凸包 

参考

安德鲁算法

分治法(其中nfox的项目实现的是分治法)

多边形快速凸包算法(Melkman‘s Algorithm)

还可以这看cpp的代码: https://www.cnblogs.com/VividBinGo/p/11637684.html

  

定义

凸包又叫凸多边形,本篇文章可能混用两种说法,

形象的理解就是一些点(点集)用一根橡皮筋紧紧地包裹外边点.

如果知道了这个定义,那么还有:

    用一个保鲜膜裹着三维点,求膜上点集.

    用一个最小的球裹着三维点,求球球的中心点和直径.

这样就进入了一个叫拓扑学的学科上.......我的老天鹅.

我竟然搜了半天没看到可以直接拿来用的c#代码,还是小轩轩给我的....

葛立恒凸包

注意一下,如果点集形成的是正交矩形的情况时,

算出来的凸包会多一个点,可以进行后处理.

(你会发现代码实际上是右上角最大点开始的,其他的教程说明从最小点开始算,

这需要知道的是最小最大都是边界点,边界点必然是凸包的边点,用谁起算都可以)

实现方式

db.Action(tr => 的委托见:https://www.cnblogs.com/JJBox/p/12196287.html

#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.DatabaseServices;
using GrxCAD.EditorInput;
using GrxCAD.Geometry;
using GrxCAD.ApplicationServices;
using GrxCAD.Runtime;
using GrxCAD.Colors;
using GrxCAD.GraphicsInterface;
using Viewport = GrxCAD.DatabaseServices.Viewport;
#endif
using System.Collections.Generic;
using System.Linq;

namespace JoinBox.src.数学操作
{
    public class Graham
    {
        /*
            视频参考: 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]的每个点,利用叉积判断,
              若栈顶的两个点与它不构成"向左转(逆时针)"的关系,则将栈顶的点出栈,直至没有点需要出栈以后,将当前点进栈
            ④ 所有点处理完之后栈中保存的点就是凸包了。
        */

        /// <summary>
        /// 最靠近x轴的点(必然是凸包边界的点)
        /// </summary>
        private Point2d _p0;
         
        /// <summary>
        /// 角度p0和pn的角度
        /// </summary>
        /// <param name="pn"></param>
        /// <returns></returns>
        private double Cosine(Point2d pn)
        {
            double d = _p0.GetDistanceTo(pn);
            //距离是相同返回1.0表示true:那么求角度(高/斜)==sin(角度)
            return d == 0.0 ? 1.0 : (pn.Y - _p0.Y) / d;
        }

        /// <summary>
        /// 叉乘,顺时针方向为真,表示要剔除
        /// </summary>
        /// <param name="n"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <returns></returns>
        private bool Clockwise(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;
        }

        
        /// <summary>
        /// 葛立恒求凸包
        /// </summary>
        /// <param name="pts"></param>
        /// <returns></returns> 
        private Point2d[] ConvexHull(List<Point2d> pts)
        {
              //消重
            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<Point2d>();
            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/BV1v741197YM 04:15
                        //栈顶就是回滚之后的,再次和qn进行向量叉乘,看看是不是顺时针,是就继续回滚..
                        //否则结束循环后加入栈中.
                        q1 = q0;
                        q0 = stack.Peek();
                    }
                    else
                    {
                        tf = false;
                    }
                }
                stack.Push(q1);
                stack.Push(qn);
                tf = true;
            }

            return stack.ToArray();
        }

        /// <summary>
        /// 葛立恒求法
        /// </summary>
        [CommandMethod("Test_Gra")]
        public void Test_Gra()
        {
            Document doc = Application.DocumentManager.MdiActiveDocument;
            Editor ed = doc.Editor;
            Database db = doc.Database;//当前的数据库
            ed.WriteMessage("
****{惊惊盒子}葛立恒求凸包,选择曲线:");

            //定义选择集选项
            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 pts = new List<Point2d>();
                foreach (ObjectId id in ssPsr.Value.GetObjectIds())
                {
                    var ent = id.ToEntity(tr);
                    if (ent is Curve curve)
                    {
                        var cs = new CurveSample<Point2d>(curve, 300);
                        pts.AddRange(cs.GetSamplePoints);
                    }
                    else if (ent is DBPoint dbPt)
                    {
                        pts.Add(new Point2d(dbPt.Position.X, dbPt.Position.Y));
                    }
                }
                Point2d[] npts = ConvexHull(pts);

                var bv = new List<BulgeVertex>();
                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);
            });
        } 
    }
}

子函数:cad曲线采样

namespace JoinBox.src.数学操作
{
    public class CurveSplit<TPoint> where TPoint : struct, IFormattable
    {
        Curve _curve { get; set; }
        double _fixedValue { get; set; }

        /// <summary>
        /// 求间隔点
        /// </summary>
        /// <param name="curve">曲线</param>
        /// <param name="fixedValue">定值分割</param>
        public CurveSplit(Curve curve, double fixedValue)
        {
            _curve = curve;
            _fixedValue = fixedValue; 
        }

        /// <summary>
        /// 获取定值分割的曲线集合
        /// </summary>
        /// <returns></returns>
        public List<Curve> GetSplitCurves
        {
            get
            {
                var lstCurves = new List<Curve>();

                //算曲线长度
                double totalLength = _curve.GetLength();

                //若少于定值,则直接返回这长度
                if (totalLength < _fixedValue)
                {
                    lstCurves.Add(_curve);
                    return lstCurves;
                }

                double addLength = 0;
                var pt3dCol = new Point3dCollection();
                //这段代码通过定值采集点
                while (addLength < totalLength)
                {
                    //求起点到长度的点
                    pt3dCol.Add(_curve.GetPointAtDist(addLength));
                    addLength += _fixedValue;
                }
                if (addLength != totalLength)
                {
                    pt3dCol.Add(_curve.GetPointAtDist(totalLength));
                }

                //通过点集,分割曲线,可能有精度问题......
                var splits = _curve.GetSplitCurves(pt3dCol);
                foreach (var item in splits)
                {
                    lstCurves.Add((Curve)item);
                }
                splits.Dispose();//释放
                return lstCurves;
            }
        }
    }

    public class CurveSample<TPoint> where TPoint : struct, IFormattable
    {
        Curve _curve { get; set; }
        int _numSample { get; set; } = 1;

        /// <summary>
        /// 求采样
        /// </summary>
        /// <param name="curve">曲线</param>
        /// <param name="numSample">采样份数</param>
        public CurveSample(Curve curve, int numSample)
        {
            _curve = curve;
            _numSample = numSample;
        }
         
        /// <summary>
        /// 曲线采样
        /// </summary>
        /// <returns></returns>
        public List<TPoint> GetSamplePoints
        {
            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<TPoint>(_curve, fixedValue);
                var spls = cs.GetSplitCurves;

                var pts = new List<TPoint>();
                pts.Add((TPoint)Activator.CreateInstance(T, spls[0].StartPoint.ToArray()));//起点
                foreach (var item in spls)
                {
                    pts.Add((TPoint)Activator.CreateInstance(T, item.EndPoint.ToArray()));//间隔点,尾点
                }
                return pts;
            }
        }
    }
}

包围盒

参考

三维包围盒实现

凸包最小周长外接矩形

凸包最小面积外接矩形

因为没看懂的关系,所以我自己想了以下的粗糙方法实现....毕竟我真的看不懂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等等....即可得到每个包围矩形,然后就可以求得最大或者最小的包围盒.

嘻嘻.说明到此结束.

代码

        /// <summary>
        /// 有向包围盒
        /// </summary>
        /// <param name="pts">点集</param>
        /// <returns>返回每条边的包围盒</returns>
        public static List<Rectangular> Boundingbox(List<Point2d> pts)
        {
            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.GetDistanceTo(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>
        public static Rectangular BoundingboxAreaMin(List<Point2d> pts)
        {
            var recs = Boundingbox(pts);
            return recs.OrderBy(rec => rec.Area).ToArray()[0];
        }

        public 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.GetDistanceTo(R4);
                    var y = R1.GetDistanceTo(R2);
                    return x * y;
                }
            }
        }

子函数:旋转一个圆

1: 点乘叉乘的子函数参考 https://www.cnblogs.com/JJBox/p/14300098.html

2: X轴到向量的弧度(角度),cad的获取的弧度是1PI,所以转换为2PI(上小,下大) 

        /// <summary>
        /// X轴到向量的弧度,cad的获取的弧度是1PI,所以转换为2PI(上小,下大)
        /// </summary>
        /// <param name="ve">向量</param>
        /// <returns>弧度</returns>
        public static double GetAngle2XAxis(this Vector3d 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;
            return alz;
        }

(完)