[推荐]几何算法相关
本帖最后由 作者 于 2010-7-5 7:52:27 编辑 <br /><br /> <p><font face="Verdana">一、从凸包开始（概念）</font></p><p>在ObjectArx.net论坛看到highflybird的帖子，关于凸包求解的Arx版本</p>
<p>原帖见<font face="Verdana"><a href="http://www.objectarx.net/forum.php?mod=viewthread&tid=1697">http://www.objectarx.net/forum.php?mod=viewthread&tid=1697</a></font></p>
<p>和highflybird的观点相同，几何算法对于Cad二次开发的初学者应该是一个门槛跨过去就是海阔天空：）</p>
<p>相关的链接在以后贴上</p>
<p> </p>
<p>二维凸包的概念是几何算法的开篇</p>
<p><font face="Verdana"><a href="http://baike.baidu.com/view/707209.htm">http://baike.baidu.com/view/707209.htm</a></font></p>
<p>概念：</p>
<p>1.1 点集Q的凸包(convex hull)是指一个最小凸多边形，满足Q中的点或者在多边形边上或者在其内。右图中由红色线段表示的多边形就是点集Q={p0,p1,...p12}的凸包。 </p>
<div class="spctrl"></div>
<p>1.2 一组平面上的点，求一个包含所有点的最小的凸多边形，这就是凸包问题了。这可以形象地想成这样：在地上放置一些不可移动的木桩，用一根绳子把他们尽量紧地圈起来，这就是凸包了。</p>
<p>
<p><font face="Verdana"></font></p>
<p></p>
<p>相关的算法网上很多，不过，起码NetApi里很少看到有人做</p> 好东西 值得学习
强！顶~！向大佬学习 过来膜拜大佬…… 本帖最后由 作者 于 2010-7-4 21:16:07 编辑
二、维护一个循环链表
理解算法后，先不要慌着看代码，最好自己先想想，脑袋里有点映像
这里首先贴上的是循环链表，我的感觉：这种方式解决点集的求解问题要更好些
当然，如果有更好的方法，欢迎一起来讨论
循环链表的概念可以看看这里：
http://student.zjzk.cn/course_ware/data_structure/web/xianxingbiao/xianxingbiao2.3.2.htm
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace TlsCad.Collections
{
/// <summary>
/// 环链表节点
/// </summary>
/// <typeparam name="T"></typeparam>
public class LoopListNode<T>
{
public T Value { get; set; }
/// <summary>
/// 上一个节点
/// </summary>
public LoopListNode<T> Previous {internal set; get; }
/// <summary>
/// 下一个节点
/// </summary>
public LoopListNode<T> Next {internal set; get; }
/// <summary>
///
/// </summary>
public LoopList<T> List {internal set; get; }
public LoopListNode(T value)
{
Value = value;
}
}
/// <summary>
/// 环链表
/// </summary>
/// <typeparam name="T"></typeparam>
public class LoopList<T> : IEnumerable<T>, IFormattable
{
/// <summary>
/// 节点数
/// </summary>
public int Count { get; private set; }
/// <summary>
/// 首节点
/// </summary>
public LoopListNode<T> First { get; private set; }
/// <summary>
/// 尾节点
/// </summary>
public LoopListNode<T> Last
{
get
{
if (First == null)
return null;
else
return First.Previous;
}
}
public bool SetFirst(LoopListNode<T> node)
{
if (Contains(node))
{
First = node;
return true;
}
return false;
}
/// <summary>
/// 交换两个节点的值
/// </summary>
/// <param name="node1"></param>
/// <param name="node2"></param>
public void Swap(LoopListNode<T> node1, LoopListNode<T> node2)
{
T value = node1.Value;
node1.Value = node2.Value;
node2.Value = value;
}
/// <summary>
/// 在首节点之前插入节点,并设置新节点为首节点
/// </summary>
/// <param name="value"></param>
/// <returns></returns>
public LoopListNode<T> AddFirst(T value)
{
LoopListNode<T> node = new LoopListNode<T>(value);
node.List = this;
if (Count == 0)
{
First = node;
First.Previous = First.Next = node;
}
else
{
LoopListNode<T> last = Last;
First.Previous = last.Next = node;
node.Next = First;
node.Previous = last;
First = node;
}
Count++;
return First;
}
/// <summary>
///在尾节点之后插入节点,并设置新节点为尾节点
/// </summary>
/// <param name="value"></param>
/// <returns></returns>
public LoopListNode<T> Add(T value)
{
LoopListNode<T> node = new LoopListNode<T>(value);
node.List = this;
if (Count == 0)
{
First = node;
First.Previous = First.Next = node;
}
else
{
LoopListNode<T> last = Last;
First.Previous = last.Next = node;
node.Next = First;
node.Previous = last;
}
Count++;
return Last;
}
/// <summary>
/// 删除首节点
/// </summary>
/// <returns></returns>
public bool RemoveFirst()
{
switch (Count)
{
case 0:
return false;
case 1:
First = null;
break;
default:
LoopListNode<T> last = Last;
First = First.Next;
First.Previous = last;
last.Next = First;
break;
}
Count--;
return true;
}
/// <summary>
/// 删除尾节点
/// </summary>
/// <returns></returns>
public bool RemoveLast()
{
switch (Count)
{
case 0:
return false;
case 1:
First = null;
break;
default:
LoopListNode<T> last = Last.Previous;
last.Next = First;
First.Previous = last;
break;
}
Count--;
return true;
}
/// <summary>
/// 删除节点
/// </summary>
/// <param name="node"></param>
/// <returns></returns>
public bool Remove(LoopListNode<T> node)
{
if (Contains(node))
{
if (Count == 1)
{
First = null;
}
else
{
if (node == First)
{
RemoveFirst();
}
else
{
node.Next.Previous = node.Previous;
node.Previous.Next = node.Next;
}
}
Count--;
return true;
}
return false;
}
public bool Contains(LoopListNode<T> node)
{
return node != null && node.List == this;
}
public bool Contains(T value)
{
LoopListNode<T> node = First;
if (node == null)
return false;
for (int i = 0; i < Count;i++ )
{
if (node.Value.Equals(value))
return true;
}
return false;
}
public LoopListNode<T> AddBefore(LoopListNode<T> node, T value)
{
if (node == First)
{
return AddFirst(value);
}
else
{
LoopListNode<T> tnode = new LoopListNode<T>(value);
node.Previous.Next = tnode;
tnode.Previous = node.Previous;
node.Previous = tnode;
tnode.Next = node;
Count++;
return tnode;
}
}
public LoopListNode<T> AddAfter(LoopListNode<T> node, T value)
{
LoopListNode<T> tnode = new LoopListNode<T>(value);
node.Next.Previous = tnode;
tnode.Next = node.Next;
node.Next = tnode;
tnode.Previous = node;
Count++;
return tnode;
}
/// <summary>
/// 链接两节点,并去除这两个节点间的所有节点
/// </summary>
/// <param name="from"></param>
/// <param name="to"></param>
public void LinkTo(LoopListNode<T> from, LoopListNode<T> to)
{
if (from != to && Contains(from) && Contains(to))
{
LoopListNode<T> node = from.Next;
bool isFirstChanged = false;
int number = 0;
while (node != to)
{
if (node == First)
isFirstChanged = true;
node = node.Next;
number++;
}
from.Next = to;
to.Previous = from;
if (number > 0 && isFirstChanged)
First = to;
Count -= number;
}
}
/// <summary>
/// 链接两节点,并去除这两个节点间的所有节点
/// </summary>
/// <param name="from"></param>
/// <param name="to"></param>
/// <param name="number"></param>
public void LinkTo(LoopListNode<T> from, LoopListNode<T> to, int number)
{
if (from != to && Contains(from) && Contains(to))
{
from.Next = to;
to.Previous = from;
First = to;
Count -= number;
}
}
/// <summary>
/// 链接两节点,并去除这两个节点间的所有节点
/// </summary>
/// <param name="from"></param>
/// <param name="to"></param>
/// <param name="number"></param>
public void LinkTo(LoopListNode<T> from, LoopListNode<T> to, int number, bool isFirstChanged)
{
if (from != to && Contains(from) && Contains(to))
{
from.Next = to;
to.Previous = from;
if (isFirstChanged)
First = to;
Count -= number;
}
}
#region IEnumerable<T> 成员
/// <summary>
/// 获取节点的查询器
/// </summary>
/// <param name="from"></param>
/// <returns></returns>
public IEnumerable<LoopListNode<T>> GetNodes(LoopListNode<T> from)
{
LoopListNode<T> node = from;
for (int i = 0; i < Count; i++)
{
yield return node;
node = node.Next;
}
}
/// <summary>
/// 获取节点的查询器
/// </summary>
/// <param name="from"></param>
/// <returns></returns>
public IEnumerable<LoopListNode<T>> GetNodes()
{
LoopListNode<T> node = First;
for (int i = 0; i < Count; i++)
{
yield return node;
node = node.Next;
}
}
/// <summary>
/// 获取节点值的查询器
/// </summary>
/// <param name="from"></param>
/// <returns></returns>
public IEnumerator<T> GetEnumerator()
{
LoopListNode<T> node = First;
for (int i = 0; i < Count; i++)
{
yield return node.Value;
node = node.Next;
}
}
IEnumerator<T> IEnumerable<T>.GetEnumerator()
{
return GetEnumerator();
}
#region IEnumerable 成员
System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator()
{
return GetEnumerator();
}
#endregion
#endregion
#region IFormattable 成员
public override string ToString()
{
string s = "( ";
foreach (T value in this)
{
s += value.ToString() + " ";
}
return s + ")";
}
string IFormattable.ToString(string format, IFormatProvider formatProvider)
{
return ToString();
}
#endregion
}
}
三、凸包类
这里放上的是三种计算凸包的算法实现
算法的关键是：
1、组成凸包的点应为左手系或右手系，判断方式是按行列式计算三点构成的面积
引：使用行列式(Determinant)来定义三角形面积，可以展开并化简为
Area(P1, P2, P3) = (x1y2 - x1y3 - x2y1 + x3y1 + x2y3 - x3y2) / 2
但请注意，上式并非单纯计算三角形的面积，而是「有向面积」(Signed Area)，它不仅告诉我们给定3点所构成三角形的大小，而且告诉我们这3点的相对方向。当Area(P1, P2, P3) 的值是负数时，这代表P1->P2->P3的走向是逆时针方向。反之，若 Area(P1, P2, P3)的值是正数，则代表这3点的走向是顺时针方向。若 Area(P1, P2, P3)为零，则代表这3点共线。因此Area这个函数在求凸包的运算中有极大的用途，可用来判断哪些点是凸包上的点，以及应用直线连结哪些点。事实上，本程序亦广泛运用这个函数。
概念相关http://home.pacific.net.hk/~kfzhou/Hulls.html
2、代码中实际实现的是左手系的凸包，核心是维护该循环链表始终为左手系
代码：
using System;
using System.Linq;
using System.Collections.Generic;
using Autodesk.AutoCAD.Geometry;
namespace TlsCad.Collections
{
public class ConvexHull2d : LoopList<Point2d>
{
#region Eval
internal void GrahamEval(List<Point2d> pnts)
{
//先找到一个基准点(XY均最小)
int i, n = 0;
for (i = 1; i < pnts.Count; i++)
{
if (pnts.Y < pnts.Y || (pnts.Y == pnts.Y && pnts.X < pnts.X))
n = i;
}
Point2d ptBase = pnts;
//按各点与基准点的极角和极长排序
var q2 =
from p in pnts
let v = p - ptBase
orderby v.Angle, v.Length
select p;
List<Point2d> ptlst = q2.ToList();
ptlst.Remove(ptBase);
//首先放入前三点
Add(ptBase);
if (ptlst.Count < 3)
{
ptlst.ForEach(p => Add(p));
return;
}
Add(ptlst);
Add(ptlst);
//如果共线
for (i = 2; i < ptlst.Count; i++)
{
if (GetArea(First.Next.Value, Last.Value, ptlst) == 0)
Last.Value = ptlst;
else
break;
};
//依次与链表末端比较
for (; i < ptlst.Count; i++)
{
//如果当前点在链表末端的顺时针方向
int num = 0;
LoopListNode<Point2d> node = Last;
while (IsClockWise(node.Previous.Value, node.Value, ptlst))
{
node = node.Previous;
num++;
}
LinkTo(node, First, num);
Add(ptlst);
}
}
internal void MelkmanEval(List<Point2d> pnts)
{
//按坐标排序,保证方向性
var q1 =
from p in pnts
orderby p.X
select p;
List<Point2d> ptlst = q1.ToList();
switch (ptlst.Count)
{
case 0:
return;
case 1:
Add(ptlst);
return;
default:
Add(ptlst);
Add(ptlst);
break;
}
//如果共线
int i = 2;
if (First.Value.X == Last.Value.X)
{
for (; i < ptlst.Count; i++)
{
if (ptlst.X == Last.Value.X)
{
double y = ptlst.Y;
if (y > Last.Value.Y)
Last.Value = ptlst;
else if (y < First.Value.Y)
First.Value = ptlst;
i++;
}
else
{
break;
}
}
}
for (; i < ptlst.Count; i++)
{
if (GetArea(First.Value, Last.Value, ptlst) == 0)
Last.Value = ptlst;
else
break;
}
if (i == ptlst.Count)
return;
//保证逆时针方向
if (IsClockWise(First.Value, Last.Value, ptlst))
Swap(First, Last);
AddFirst(ptlst);
//依次比较
for (i++; i < ptlst.Count; i++)
{
Point2d pnt = ptlst;
int num = 0;
LoopListNode<Point2d> from= First, to = First;
//做左链
while (IsClockWise(to.Next.Value, pnt, to.Value))
{
to = to.Next;
num++;
}
//做右链
while (IsClockWise(from.Previous.Value, from.Value, pnt))
{
from = from.Previous;
num++;
}
LinkTo(from, to, num - 1);
AddFirst(pnt);
}
}
internal void JarvisEval(List<Point2d> pnts)
{
switch (pnts.Count)
{
case 0:
return;
case 1:
Add(pnts);
return;
case 2:
Add(pnts);
Add(pnts);
return;
default:
Add(pnts);
Add(pnts);
break;
}
int i = 2;
List<Point2d> tpnts = new List<Point2d> { First.Value, Last.Value };
for (; i < pnts.Count; i++)
{
if (GetArea(First.Value, Last.Value, pnts) == 0)
tpnts.Add(pnts);
else
break;
}
var q1 =
from p in tpnts
orderby p.X, p.Y
select p;
First.Value = q1.First();
Last.Value = q1.Last();
if (i == pnts.Count)
return;
//保证逆时针方向
if (IsClockWise(First.Value, Last.Value, pnts))
Swap(First, Last);
AddFirst(pnts);
for (i++; i < pnts.Count; i++)
{
Point2d pnt = pnts;
Vector2d vec1 = Last.Value - pnt;
Vector2d vec2 = First.Value - pnt;
bool iscw1 = false;
bool iscw2 = IsClockWise(vec1, vec2);
LoopListNode<Point2d> from = null, to = null;
int num = 0;
foreach (var node in GetNodes())
{
vec1 = vec2;
vec2 = node.Next.Value - pnt;
iscw1 = iscw2;
iscw2 = IsClockWise(vec1, vec2);
if (iscw1)
{
if (iscw2)
num++;
else
to = node;
}
else if (iscw2)
{
from = node;
}
}
if (from != null)
{
LinkTo(from, to, num);
Add(pnt);
}
}
}
#endregion
public bool IsOutside(Point2d pnt)
{
foreach (var node in GetNodes())
{
if(IsClockWise(node.Value, node.Next.Value, pnt))
return true;
}
return false;
}
public bool IsInside(Point2d pnt)
{
foreach (var node in GetNodes())
{
if (IsClockWise(node.Value, pnt, node.Next.Value))
return true;
}
return false;
}
public bool IsOn(Point2d pnt)
{
foreach (var node in GetNodes())
{
using (var ls2d = new LineSegment2d(node.Value, node.Next.Value))
{
if (ls2d.IsOn(pnt))
return true;
}
}
return false;
}
//public double GetMaxDistance(out LoopList<Point2d> ptlst)
//{
//}
private static double GetArea(Point2d ptBase, Point2d pt1, Point2d pt2)
{
return (pt2 - ptBase).DotProduct((pt1 - ptBase).GetPerpendicularVector());
}
private static bool IsClockWise(Point2d ptBase, Point2d pt1, Point2d pt2)
{
return GetArea(ptBase, pt1, pt2) <= 0;
}
private static double GetArea(Vector2d vecBase, Vector2d vec)
{
return vec.DotProduct(vecBase.GetPerpendicularVector());
}
private static bool IsClockWise(Vector2d vecBase, Vector2d vec)
{
return GetArea(vecBase, vec) <= 0;
}
}
}
相关的测试代码和扩展函数
public void test21()
{
Document doc = Application.DocumentManager.MdiActiveDocument;
Editor ed = doc.Editor;
Random rand = new Random();
int num = Convert.ToInt32(ed.GetString("\ninput number of points:").StringResult);
using (DBTransaction tr = new DBTransaction())
{
tr.Database.Pdmode = 35;
tr.Database.Pdsize = -2;
var pnts =
Enumerable.Range(0, num).Select
(
i =>
{
Point3d pnt = new Point3d(rand.NextDouble() * 100, rand.NextDouble() * 100, 0);
DBPoint dpnt = new DBPoint(pnt);
return dpnt;
}
);
tr.OpenCurrentSpace(OpenMode.ForWrite);
tr.AddEntity(pnts);
}
}
/// <summary>
/// 这段测试演示凸包和最小包围圆的函数用法,
/// 测试前请打开新的图形文件(dwg)并画上点,
/// 这段代码不会产生任何实体
/// </summary>
public static void test22()
{
Document doc = Application.DocumentManager.MdiActiveDocument;
Editor ed = doc.Editor;
Database db = doc.Database;
int num = 0;
double total = 0;
List<Point2d> pnts;
Stopwatch watch = new Stopwatch();
watch.Reset();
watch.Start();
var resSel = ed.SelectAll(new ResultList { { 0, "point" } });
num = resSel.Value.Count;
watch.Stop();
ed.WriteMessage("\n总共{0}个点", num);
ed.WriteMessage("\n选择集耗用的时间:{0}(ms)", watch.ElapsedMilliseconds);
total += watch.ElapsedMilliseconds;
watch.Reset();
watch.Start();
using (Transaction tr = db.TransactionManager.StartTransaction())
{
pnts =
resSel.Value.GetObjectIds()
.Select(id => ((DBPoint)tr.GetObject(id, OpenMode.ForRead)).Position)
.Select(pt => new Point2d(pt.X, pt.Y)).ToList();
watch.Stop();
ed.WriteMessage("\n打开并转换为二维点耗时:{0}毫秒", watch.ElapsedMilliseconds);
total += watch.ElapsedMilliseconds;
watch.Reset();
watch.Start();
LoopList<Point2d> ptlst;
CircularArc2d ca2d = pnts.GetMinCircle(out ptlst);
watch.Stop();
ed.DrawVectors(ca2d.GetSamplePoints(1000), 5);
ed.WriteMessage("\n最小包围圆计算耗时:{0}毫秒\n通过了{1}点:\n{2}", watch.ElapsedMilliseconds, ptlst.Count, ptlst);
total += watch.ElapsedMilliseconds;
}
watch.Reset();
watch.Start();
ConvexHull2d pchMelkman = pnts.GetConvexHull();
watch.Stop();
ed.WriteMessage("\nMelkman计算凸包耗时:{0}毫秒 凸包点数:{1}", watch.ElapsedMilliseconds, pchMelkman.Count);
total += watch.ElapsedMilliseconds;
watch.Reset();
watch.Start();
ConvexHull2d pchJarvis = pnts.GetConvexHull2();
watch.Stop();
ed.WriteMessage("\nJarvis计算凸包耗时:{0}毫秒 凸包点数:{1}", watch.ElapsedMilliseconds, pchJarvis.Count);
total += watch.ElapsedMilliseconds;
watch.Reset();
watch.Start();
ConvexHull2d pchGraham = pnts.GetConvexHull3();
watch.Stop();
ed.WriteMessage("\nGraham计算凸包耗时:{0}毫秒 凸包点数:{1}", watch.ElapsedMilliseconds, pchGraham.Count);
total += watch.ElapsedMilliseconds;
ed.WriteMessage("\n总耗时:{0}毫秒", total);
ed.DrawVectors(pchMelkman, 4);
ed.DrawPoints(pchMelkman, 1, 1, 32);
ed.DrawPoints(pchJarvis, 2, 2, 32);
ed.DrawPoints(pchGraham, 3, 3, 32);
}
public void test23()
{
Stopwatch watch = new Stopwatch();
using (DBTransaction tr = new DBTransaction())
{
var resSel = tr.Editor.SelectAll();
if (resSel.Status == PromptStatus.OK)
{
var ss = resSel.Value;
int num = ss.Count;
watch.Start();
resSel.Value.ForEach(OpenMode.ForWrite, ent => ent.Erase());
watch.Stop();
tr.Editor.WriteMessage("\n{0} Entitys Erase Time:{1}(ms)", num, watch.ElapsedMilliseconds);
}
tr.ZoomWindow(new Point3d(-10, -10, 0), new Point3d(110, 110, 0));
}
}
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Runtime.InteropServices;
using System.Reflection;
using Autodesk.AutoCAD.Runtime;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.ApplicationServices;
namespace TlsCad.ExtendMethods
{
public static class EditorEx
{
/// <summary>
/// 2009版本提供了acedCmd函数的封装
/// </summary>
private static MethodInfo _runCommand =
typeof(Editor).GetMethod(
"RunCommand",
BindingFlags.NonPublic | BindingFlags.Instance);
/// <summary>
/// 反射调用AutoCad命令(2009版本以上)
/// </summary>
/// <param name="editor">Editor对象</param>
/// <param name="args">参数</param>
/// <returns></returns>
public static PromptStatus Command(this Editor editor, params object[] args)
{
return (PromptStatus)_runCommand.Invoke(editor, new object[] { args });
}
public static void DrawVectors(this Editor editor, IEnumerable<Point2d> pnts, short colorIndex)
{
var itor = pnts.GetEnumerator();
if (!itor.MoveNext())
return;
TypedValue tvFirst = new TypedValue((int)LispDataType.Point2d, itor.Current);
ResultBuffer rb =
new ResultBuffer
{
new TypedValue((int)LispDataType.Int16, colorIndex),
};
TypedValue tv1;
TypedValue tv2 = tvFirst;
while (itor.MoveNext())
{
tv1 = tv2;
tv2 = new TypedValue((int)LispDataType.Point2d, itor.Current);
rb.Add(tv1);
rb.Add(tv2);
}
rb.Add(tv2);
rb.Add(tvFirst);
editor.DrawVectors(rb, Matrix3d.Identity);
}
public static void DrawPoints(this Editor editor, IEnumerable<Point2d> pnts, short colorIndex, double radius, int numEdges)
{
foreach (Point2d pnt in pnts)
{
editor.DrawPoint(pnt, colorIndex, radius, numEdges);
}
}
public static void DrawPoint(this Editor editor, Point2d pnt, short colorIndex, double radius, int numEdges)
{
Vector2d vec = Vector2d.XAxis * radius;
double angle = Math.PI * 2 / numEdges;
List<Point2d> pnts = new List<Point2d>();
pnts.Add(pnt + vec);
for (int i = 1; i < numEdges; i++)
{
pnts.Add(pnt + vec.RotateBy(angle * i));
}
editor.DrawVectors(pnts, colorIndex);
}
}
}
四、最小包围圆
相关的讨论链接：
http://bbs.mjtd.com/forum.php?mod=viewthread&tid=55997
http://www.objectarx.net/forum.php?mod=viewthread&tid=4870
相关的概念和算法描述就直接贴过来了，：）
以下部分摘录自:汪卫，王文平，汪嘉业，求一个包含点集所有点的最小圆的算法，软件学报，2000，11（9）：1237-1240
求一个最小圆包含给定点集所有点的问题是人们在实践和理论上都十分感兴趣的一个问题.由于这个圆的圆心是到点集中最远点最近的一个点,因而在规划某些设施时很有实用价值.
这个圆心也可看成是点集的中心.此外,在图形学中,圆也常可取作边界盒,使用它可减少很多不必要的计算.
在空间数据库中可将该问题用于建立空间数据的索引以提高查询速度.这个问题看起来十分简单,但用直观的算法去解此问题,其复杂性可达O(n4),其中n为点集中点的数目.
有关此问题的讨论在计算几何的专著及论文中未见报道.
本文提出了一种新的算法并证明了这种算法的时间复杂性为O(|lg(d/R)|*n),其中R是所求的最小圆的半径,d为点集中不在圆周上但距圆周最近的点到圆周的距离.
1、算 法
第1步.在点集中任取3点A,B,C.
第2步.作一个包含A,B,C三点的最小圆.圆周可能通过这3点(如图1所示),也可能只通过其中两点,但包含第3点.后一种情况圆周上的两点一定是位于圆的一条直径的两端.
第3步.在点集中找出距离第2步所建圆圆心最远的点D.若D点已在圆内或圆周上,则该圆即为所求的圆,算法结束.否则,执行第4步.
第4步.在A,B,C,D中选3个点,使由它们生成的一个包含这4点的圆为最小.这3点成为新的A,B和C,返回执行第2步. 若在第4步生成的圆的圆周只通过A,B,C,D中的两点,则圆周上的两点取成新的A和B,从另两点中任取一点作为新的C.
2、算法正确性
本节要证明上述算法一定能终止,且最后一次求得的圆即是所要求的包含点集所有点的最小圆.
引理1.算法第4步所生成的圆的半径随着迭代过程递增. 证明:因为第4步每一次生成的圆是包含原来的A,B,C三点,又要加上圆外的一点,而上一次生成的圆是包含A,B,C的最小圆,因而新圆的半径一定比原来的圆半径要大.
定理1.上述算法是收敛的,且最后得到包含点集所有点的最小圆. 证明:因为在点集中任取3点或两点生成的包含这3点或两点的最小圆的个数是有限的.由引理1可知,算法进行过程中所生成的圆的半径是递增的,因而经过有限次迭代后,可求得这些最小圆中半径最大的一个.从算法第3步可知,只有当点集中所有点都在由3个点或两个点生成的最小圆内时,算法才结束,因而最后得到的圆一定是包含点集中所有点的最小圆.
代码
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using TlsCad.Collections;
using Autodesk.AutoCAD.Geometry;
namespace TlsCad.ExtendMethods
{
public static class PointEx
{
#region PointList
/// <summary>
/// 判断点是否属于圆
/// </summary>
/// <param name="ca2d"></param>
/// <param name="pnt"></param>
/// <returns></returns>
public static bool IsPartOf(this CircularArc2d ca2d, Point2d pnt)
{
return ca2d.IsOn(pnt) || ca2d.IsInside(pnt);
}
/// <summary>
/// 按两点返回最小包围圆
/// </summary>
/// <param name="pt1"></param>
/// <param name="pt2"></param>
/// <param name="ptlst">输出圆上的点</param>
/// <returns></returns>
public static CircularArc2d GetMinCircle(Point2d pt1, Point2d pt2, out LoopList<Point2d> ptlst)
{
ptlst = new LoopList<Point2d> { pt1, pt2 };
return
new CircularArc2d
(
(pt1 + pt2.GetAsVector()) / 2,
pt1.GetDistanceTo(pt2) / 2
);
}
/// <summary>
/// 按三点返回最小包围圆
/// </summary>
/// <param name="pt1"></param>
/// <param name="pt2"></param>
/// <param name="pt3"></param>
/// <param name="ptlst">输出圆上的点</param>
/// <returns></returns>
public static CircularArc2d GetMinCircle(Point2d pt1, Point2d pt2, Point2d pt3, out LoopList<Point2d> ptlst)
{
ptlst =
new LoopList<Point2d> { pt1, pt2, pt3 };
//遍历各点与下一点的向量长度,找到距离最大的两个点
double maxLength;
LoopListNode<Point2d> maxNode =
ptlst.GetNodes().FindByMaxKey
(
out maxLength,
node => node.Value.GetDistanceTo(node.Next.Value)
);
//以两点做最小包围圆
LoopList<Point2d> tptlst;
CircularArc2d ca2d =
GetMinCircle(maxNode.Value, maxNode.Next.Value, out tptlst);
//如果另一点属于该圆
if (ca2d.IsPartOf(maxNode.Previous.Value))
{
//返回
ptlst = tptlst;
return ca2d;
}
//否则按三点做圆
ptlst.SetFirst(maxNode);
ca2d = new CircularArc2d(pt1, pt2, pt3);
ca2d.SetAngles(0, Math.PI * 2);
return ca2d;
}
/// <summary>
/// 按四点返回最小包围圆
/// </summary>
/// <param name="pt1"></param>
/// <param name="pt2"></param>
/// <param name="pt3"></param>
/// <param name="pt4"></param>
/// <param name="ptlst">输出圆上的点</param>
/// <returns></returns>
public static CircularArc2d GetMinCircle(Point2d pt1, Point2d pt2, Point2d pt3, Point2d pt4, out LoopList<Point2d> ptlst)
{
LoopList<Point2d> iniptlst =
new LoopList<Point2d> { pt1, pt2, pt3, pt4 };
ptlst = null;
CircularArc2d ca2d = null;
//遍历C43的组合,环链表的优势在这里
foreach (LoopListNode<Point2d> firstNode in iniptlst.GetNodes())
{
//获取各组合下三点的最小包围圆
LoopListNode<Point2d> secondNode = firstNode.Next;
LoopListNode<Point2d> thirdNode = secondNode.Next;
LoopList<Point2d> tptlst;
CircularArc2d tca2d = GetMinCircle(firstNode.Value, secondNode.Value, thirdNode.Value, out tptlst);
//如果另一点属于该圆,并且半径小于当前值就把它做为候选解
if (tca2d.IsPartOf(firstNode.Previous.Value))
{
if (ca2d == null || tca2d.Radius < ca2d.Radius)
{
ca2d = tca2d;
ptlst = tptlst;
}
}
}
//返回直径最小的圆
return ca2d;
}
/// <summary>
/// 按点集返回最小包围圆
/// </summary>
/// <param name="pnts"></param>
/// <param name="ptlst">输出圆上的点</param>
/// <returns></returns>
public static CircularArc2d GetMinCircle(this List<Point2d> pnts, out LoopList<Point2d> ptlst)
{
//点数较小时直接返回
switch (pnts.Count)
{
case 0:
ptlst = new LoopList<Point2d>();
return null;
case 1:
ptlst = new LoopList<Point2d> { pnts };
return new CircularArc2d(pnts, 0);
case 2:
return GetMinCircle(pnts, pnts, out ptlst);
case 3:
return GetMinCircle(pnts, pnts, pnts, out ptlst);
case 4:
return GetMinCircle(pnts, pnts, pnts, pnts, out ptlst);
}
//按前三点计算最小包围圆
Point2d[] tpnts = new Point2d;
pnts.CopyTo(0, tpnts, 0, 3);
CircularArc2d ca2d = GetMinCircle(tpnts, tpnts, tpnts, out ptlst);
//找到点集中距离圆心的最远点为第四点
tpnts = pnts.FindByMaxKey(pnt => pnt.GetDistanceTo(ca2d.Center));
//如果最远点属于圆结束
while (!ca2d.IsPartOf(tpnts))
{
//如果最远点不属于圆,按此四点计算最小包围圆
ca2d = GetMinCircle(tpnts, tpnts, tpnts, tpnts, out ptlst);
//将结果作为新的前三点
if (ptlst.Count == 3)
{
tpnts = ptlst.Last.Value;
}
else
{
//如果计算的结果只有两个点
//选择上一步计算的圆心为基准,找出点集中圆心的最远点作为第三点
//似乎是算法的问题?这里第三点不能任意选择,否则无法收敛
//tpnts = pnts.GetMaxBy(pnt => pnt.GetDistanceTo(ca2d.Center));
//if (ca2d.IsPartOf(tpnts))
// return ca2d;
//看来是理解有误,第三点应该取另两点中距离圆心较远的点
//但按算法中描述的任选其中一点的话,还是无法收敛......
tpnts =
tpnts.Except(ptlst)
.FindByMaxKey(pnt => ca2d.Center.GetDistanceTo(pnt));
}
tpnts = ptlst.First.Value;
tpnts = ptlst.First.Next.Value;
//按此三点计算最小包围圆
ca2d = GetMinCircle(tpnts, tpnts, tpnts, out ptlst);
//找到点集中圆心的最远点为第四点
tpnts = pnts.FindByMaxKey(pnt => pnt.GetDistanceTo(ca2d.Center));
}
return ca2d;
}
public static ConvexHull2d GetConvexHull(this List<Point2d> pnts)
{
ConvexHull2d ch2d = new ConvexHull2d();
ch2d.MelkmanEval(pnts);
return ch2d;
}
public static ConvexHull2d GetConvexHull2(this List<Point2d> pnts)
{
ConvexHull2d ch2d = new ConvexHull2d();
ch2d.JarvisEval(pnts);
return ch2d;
}
public static ConvexHull2d GetConvexHull3(this List<Point2d> pnts)
{
ConvexHull2d ch2d = new ConvexHull2d();
ch2d.GrahamEval(pnts);
return ch2d;
}
#endregion
}
}
相关的扩展函数和效果(测试代码见4楼)
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace TlsCad.ExtendMethods
{
public static class LinqEx
{
/// <summary>
/// 按转换函数找出序列中最大键值的对应值
/// </summary>
/// <typeparam name="TValue"></typeparam>
/// <typeparam name="TKey"></typeparam>
/// <param name="enumerable"></param>
/// <param name="func"></param>
/// <returns></returns>
public static TValue FindByMaxKey<TValue, TKey>(this IEnumerable<TValue> enumerable, Func<TValue, TKey> func)
where TKey : IComparable<TKey>
{
var itor = enumerable.GetEnumerator();
if (!itor.MoveNext())
throw new ArgumentNullException();
TValue value = itor.Current;
TKey key = func(value);
while (itor.MoveNext())
{
TKey tkey = func(itor.Current);
if (tkey.CompareTo(key) > 0)
{
key = tkey;
value = itor.Current;
}
}
return value;
}
/// <summary>
/// 按转换函数找出序列中最大键值的对应值
/// </summary>
/// <typeparam name="TValue"></typeparam>
/// <typeparam name="TKey"></typeparam>
/// <param name="enumerable"></param>
/// <param name="maxResult">对应的最大键值</param>
/// <param name="func"></param>
/// <returns></returns>
public static TValue FindByMaxKey<TValue, TKey>(this IEnumerable<TValue> enumerable, out TKey maxResult, Func<TValue, TKey> func)
where TKey : IComparable<TKey>
{
var itor = enumerable.GetEnumerator();
if (!itor.MoveNext())
throw new ArgumentNullException();
TValue value = itor.Current;
TKey key = func(value);
while (itor.MoveNext())
{
TKey tkey = func(itor.Current);
if (tkey.CompareTo(key) > 0)
{
key = tkey;
value = itor.Current;
}
}
maxResult = key;
return value;
}
/// <summary>
/// 按转换函数找出序列中最小键值的对应值
/// </summary>
/// <typeparam name="TValue"></typeparam>
/// <typeparam name="TKey"></typeparam>
/// <param name="enumerable"></param>
/// <param name="maxResult">对应的最小键值</param>
/// <param name="func"></param>
/// <returns></returns>
public static TValue FindByMinKey<TValue, TKey>(this IEnumerable<TValue> enumerable, out TKey minKey, Func<TValue, TKey> func)
where TKey : IComparable<TKey>
{
var itor = enumerable.GetEnumerator();
if (!itor.MoveNext())
throw new ArgumentNullException();
TValue value = itor.Current;
TKey key = func(value);
while (itor.MoveNext())
{
TKey tkey = func(itor.Current);
if (tkey.CompareTo(key) < 0)
{
key = tkey;
value = itor.Current;
}
}
minKey = key;
return value;
}
/// <summary>
/// 按转换函数找出序列中最小键值的对应值
/// </summary>
/// <typeparam name="TValue"></typeparam>
/// <typeparam name="TKey"></typeparam>
/// <param name="enumerable"></param>
/// <param name="func"></param>
/// <returns></returns>
public static TValue FindByMinKey<TValue, TKey>(this IEnumerable<TValue> enumerable, Func<TValue, TKey> func)
where TKey : IComparable<TKey>
{
var itor = enumerable.GetEnumerator();
if (!itor.MoveNext())
throw new ArgumentNullException();
TValue value = itor.Current;
TKey key = func(value);
while (itor.MoveNext())
{
TKey tkey = func(itor.Current);
if (tkey.CompareTo(key) < 0)
{
key = tkey;
value = itor.Current;
}
}
return value;
}
/// <summary>
/// 按转换函数找出序列中最(小/大)键值的对应值
/// </summary>
/// <typeparam name="TValue"></typeparam>
/// <typeparam name="TKey"></typeparam>
/// <param name="enumerable"></param>
/// <param name="func"></param>
/// <returns></returns>
public static TValue[] FindByMKeys<TValue, TKey>(this IEnumerable<TValue> enumerable, Func<TValue, TKey> func)
where TKey : IComparable<TKey>
{
var itor = enumerable.GetEnumerator();
if (!itor.MoveNext())
throw new ArgumentNullException();
TValue[] values = new TValue;
values = values = itor.Current;
TKey minKey = func(values);
TKey maxKey = minKey;
while (itor.MoveNext())
{
TKey tkey = func(itor.Current);
if (tkey.CompareTo(minKey) < 0)
{
minKey = tkey;
values = itor.Current;
}
else if (tkey.CompareTo(maxKey) > 0)
{
maxKey = tkey;
values = itor.Current;
}
}
return values;
}
/// <summary>
/// 按比较器找出序列中最(小/大)键值的对应值
/// </summary>
/// <typeparam name="TValue"></typeparam>
/// <param name="enumerable"></param>
/// <param name="comparison"></param>
/// <returns></returns>
public static TValue[] FindByMKeys<TValue>(this IEnumerable<TValue> enumerable, Comparison<TValue> comparison)
{
var itor = enumerable.GetEnumerator();
if (!itor.MoveNext())
throw new ArgumentNullException();
TValue[] values = new TValue;
values = values = itor.Current;
while (itor.MoveNext())
{
if (comparison(itor.Current, values) < 0)
{
values = itor.Current;
}
else if (comparison(itor.Current, values) > 0)
{
values = itor.Current;
}
}
return values;
}
/// <summary>
/// 按转换函数找出序列中最(小/大)键值的对应键值
/// </summary>
/// <typeparam name="TValue"></typeparam>
/// <typeparam name="TKey"></typeparam>
/// <param name="enumerable"></param>
/// <param name="func"></param>
/// <returns></returns>
public static TKey[] FindMKeys<TValue, TKey>(this IEnumerable<TValue> enumerable, Func<TValue, TKey> func)
where TKey : IComparable<TKey>
{
var itor = enumerable.GetEnumerator();
if (!itor.MoveNext())
throw new ArgumentNullException();
TKey[] keys = new TKey;
keys = keys = func(itor.Current);
while (itor.MoveNext())
{
TKey tkey = func(itor.Current);
if (tkey.CompareTo(keys) < 0)
{
keys = tkey;
}
else if (tkey.CompareTo(keys) > 0)
{
keys = tkey;
}
}
return keys;
}
}
}
用于CAD的什么版本啊？我用CAD2008，无法编译通过，比如DBTransaction类就无法识别，是版本错误还是少引用了什么文件？ <p><font face="Verdana">DBTransaction:</font></p>
<p><font face="Verdana"><a href="http://bbs.mjtd.com/forum.php?mod=viewthread&tid=76123">http://bbs.mjtd.com/forum.php?mod=viewthread&tid=76123</a></font></p>
<p>不过这是个较老的版本，</p>
<p>tt3测试命令的</p>
<p>tr.ZoomWindow(new Point3d(-10, -10, 0), new Point3d(110, 110, 0));</p>
<p>还是不能用</p>
<p>不过，这是个纯测试代码，可以去掉</p>
<p>另外版本为：VS2008+ACad2008,上面的代码用到了很多C#3.0的语法，如果是VS2005，就要改很多了</p> <p>前来学习</p>
<p> </p> 顶！！！！！