前段時間看到生成凸包的Graham算法，查了一些資料，沒看到c#版的，於是想動手寫一個，該程序草草完成，其中有值得優化的地方，諸位可自行改正。
Graham算法的原理：
（1）給定一個點集P={P0,P1,.....Pn),找出點集中Y值最小的點，如果Y值最小的點有多個，可選擇其中X值最小的 www.2cto.com
（2）假設上一步找出的點位P0，現在對剩下的點進行極角排序，按逆時針（本程序是這樣，其他的方式是一樣的）極角從小到大排序，假定排序結果集為{p1,p2,p3,p4,....pn}。注意，這裡不一定有N個點，因為可能存在幾個點的極角相同，這時取據P0距最遠的點，也或者你設定了一個極角阈值，當兩個點的極角差小於該阈值時，取據P0距離遠的點，這時生成的凸包有一點點誤差，誤差的大小取決於你設置的阈值。本程序沒采用阈值，所以生成的凸包理論上不存在誤差。
在進行極角排序時，不需要真的算法每個點的極角（注意，這裡的極角是該點與P0相對於X軸的夾角），只需要使用向量叉積來判斷即可，這個過程我使用了鏈表來存儲排序結果，因為這個過程會進行頻繁的插入。
（3）將p0,p1,p3入棧，p0,p1這兩個點肯定在凸包上（原因很簡單，p0不用解釋了吧，p1因為它是極角最小的且為第一個點），p2則不一定在凸包上，然後進行循環(for(int i=2;i<n;i++），判定棧頂的下一個點，棧頂點，及p[i]點這三點組成的折線段是否向左轉，如果是的話，則p[i]入棧；否則，當前位於棧頂的點不在凸包上，彈棧（該過程進行回朔,確保之前的所有點轉向正確，這個步驟很重要，不然不能生成凸多邊形），最後返回棧即可。
下面是整個代碼：
算法部分：
class ConvexAogrithm
{
private List<PointF> nodes;
private Stack<PointF> sortedNodes;
public PointF[] sor_nodes;
public ConvexAogrithm(List<PointF> points)
{
nodes = points;
}
private double DistanceOfNodes(PointF p0, PointF p1)
{
if (p0.IsEmpty || p1.IsEmpty)
return 0.0;
return Math.Sqrt((p1.X - p0.X) * (p1.X - p0.X) + (p1.Y - p0.Y) * (p1.Y - p0.Y));
}
public void GetNodesByAngle( out PointF p0)
{
LinkedList<PointF> list_node = new LinkedList<PointF>();
p0 = GetMinYPoint();
LinkedListNode<PointF> node = new LinkedListNode<PointF>(nodes[0]);
list_node.AddFirst(node);
for (int i = 1; i < nodes.Count; i++)
{
int direct = IsClockDirection(p0, node.Value, nodes[i]);
if (direct == 1)
{
list_node.AddLast(nodes[i]);
node = list_node.Last;
//node.Value = nodes[i];
}
else if (direct == -10)
{
list_node.Last.Value = nodes[i];
//node = list_node.Last
//node.Value = nodes[i];
}
else if (direct == 10)
continue;
else if (direct == -1)
{
LinkedListNode<PointF> temp = node.Previous;
while (temp != null && IsClockDirection(p0, temp.Value, nodes[i]) == -1)
{
temp = temp.Previous;
}
if (temp == null)
{
list_node.AddFirst(nodes[i]);
continue;
}
if (IsClockDirection(p0, temp.Value, nodes[i]) == -10)
temp.Value = nodes[i];
else if (IsClockDirection(p0, temp.Value, nodes[i]) == 10)
continue;
else
list_node.AddAfter(temp, nodes[i]);
}
}
sor_nodes = list_node.ToArray();
sortedNodes = new Stack<PointF>();
sortedNodes.Push(p0);
sortedNodes.Push(sor_nodes[0]);
sortedNodes.Push(sor_nodes[1]);
for (int i = 2; i<sor_nodes.Length; i++)
{
PointF p2 = sor_nodes[i];
PointF p1 = sortedNodes.Pop();
PointF p0_sec = sortedNodes.Pop();
sortedNodes.Push(p0_sec);
sortedNodes.Push(p1);
if (IsClockDirection1(p0_sec, p1, p2) == 1)
{
sortedNodes.Push(p2);
continue;
}
while (IsClockDirection1(p0_sec, p1, p2) != 1)
{
sortedNodes.Pop();
p1 = sortedNodes.Pop();
p0_sec = sortedNodes.Pop();
sortedNodes.Push(p0_sec);
sortedNodes.Push(p1);
}
sortedNodes.Push(p2);
}
}
private int IsClockDirection1(PointF p0, PointF p1, PointF p2)
{
PointF p0_p1 = new PointF(p1.X - p0.X, p1.Y - p0.Y);
PointF p0_p2 = new PointF(p2.X - p0.X, p2.Y - p0.Y);
return (p0_p1.X * p0_p2.Y - p0_p2.X * p0_p1.Y) > 0 ? 1 : -1;
}
private PointF GetMinYPoint()
{
PointF succNode;
float miny=nodes.Min(r=>r.Y);
IEnumerable<PointF> pminYs = nodes.Where(r => r.Y == miny);
PointF[] ps = pminYs.ToArray();
if (pminYs.Count() > 1)
{
succNode = pminYs.Single(r => r.X == pminYs.Min(t => t.X));
nodes.Remove(succNode);
return succNode;
}
else
{
nodes.Remove(ps[0]);
return ps[0];
}
}
private int IsClockDirection(PointF p0, PointF p1, PointF p2)
{
PointF p0_p1 = new PointF(p1.X-p0.X,p1.Y-p0.Y) ;
PointF p0_p2 = new PointF(p2.X - p0.X, p2.Y - p0.Y);
if ((p0_p1.X * p0_p2.Y - p0_p2.X * p0_p1.Y) != 0)
return (p0_p1.X * p0_p2.Y - p0_p2.X * p0_p1.Y) > 0 ? 1 : -1;
else
return DistanceOfNodes(p0, p1) > DistanceOfNodes(p0, p2) ? 10 : -10;
}
public Stack<PointF> SortedNodes
{
get { return sortedNodes; }
}
}
界面部分，供測試使用：
public partial class Form1 : Form
{
private List<PointF> nodes;
private Graphics g;
private Pen pen;
public Form1()
{
InitializeComponent();
g=this.panel1.CreateGraphics();
g.TranslateTransform(0f,this.panel1.Height);
g.ScaleTransform(1f,-1f);
pen = new Pen(Color.Blue);
}
private void button2_Click(object sender, EventArgs e)
{
g.Clear(panel1.BackColor);
nodes = new List<PointF>();
nodes.Clear();
Random rand = new Random();
Point p = new Point(); ;
for (int i = 0; i < 1000; i++)
{
p.X = rand.Next(10, panel1.Width - 9);
p.Y = rand.Next(10, panel1.Height - 9);
nodes.Add(p);
DrawCircle(p);
}
}
private void DrawCircle(Point p)
{
g.DrawEllipse(pen, p.X - 4, p.Y - 4, 8, 8);
g.FillEllipse(Brushes.Blue, p.X - 4, p.Y - 4, 8, 8);
}
private void button1_Click(object sender, EventArgs e)
{
ConvexAogrithm ca = new ConvexAogrithm(nodes);
PointF p;
ca.GetNodesByAngle(out p);
//PointF[] ps = ca.sor_nodes;
//float[] psangle=new float[ps.Length];
//for (int i = 0; i < psangle.Length; i++)
// psangle[i] = CalcAngle(p, ps[i]);
g.DrawEllipse(pen, p.X - 8, p.Y - 8, 16,16);
g.FillEllipse(Brushes.Blue, p.X - 8, p.Y - 8, 16, 16);
Stack<PointF> p_nodes = ca.SortedNodes;
pen = new Pen(Color.Black, 2.0f);
g.SmoothingMode = SmoothingMode.HighQuality;
pen.LineJoin = LineJoin.Round;
g.DrawPolygon(pen, p_nodes.ToArray());
}
private float CalcAngle(PointF p1,PointF p2)
{
float angle = (float)(Math.Atan(Math.Abs(p2.Y - p1.Y + 0.0) / Math.Abs(p2.X - p1.X + 0.0)) * 180 / Math.PI);
if ((p2.Y - p1.Y + 0.0) / (p2.X - p1.X + 0.0) < 0)
angle = 180 - angle;
return angle;
}
}
上面注釋的為我當初測試排序極角的結果使用的,psangle是排序的結果，調試可看到角度從小到大。1000個點凸包如下（那個大圓點是P0，我為標記使用，無其他含義）：
200個點的凸包如下：
