using System;
using System.Collections.Generic;
public class UTest
{
static void Main()
{
Vec2[] points = new Vec2[ 20000 ];
Random random = new Random(123456);
for (int i = 0; i < points.Length; i++)
{
points[i] = new Vec2(random.Next(200000) / 100.0f, random.Next(200000) / 100.0f);
}
int startTick, endTick;
startTick = System.Environment.TickCount;
float divideconquer = (float)Math.Sqrt((double)MinDistanceSquared(points)); // 140ms
endTick = System.Environment.TickCount;
Console.WriteLine("Divide and conquer finishes in {0,6}ms, result={1}", endTick - startTick, divideconquer);
startTick = System.Environment.TickCount;
float bruteforce = (float)Math.Sqrt((double)BruteForceMinDistanceSquared(points)); // 13200ms
endTick = System.Environment.TickCount;
Console.WriteLine("Brute force method finishes in {0,6}ms, result={1}", endTick - startTick, bruteforce);
}
static float MinDistanceSquared(Vec2[] points)
{
// 如果点的数目少于一定数量，直接穷举最短距离
if (points.Length < 6) return BruteForceMinDistanceSquared(points);
// 将点按照x轴排序。并分成左边部分，和右边部分。
Array.Sort<Vec2>(points, Vec2.CompareByX);
int middleIdx = points.Length / 2;
Vec2[] leftPoints = new Vec2[middleIdx];
Vec2[] rightPoints = new Vec2[points.Length - middleIdx];
for (int i = 0; i < leftPoints.Length; i++)
{
leftPoints[i] = new Vec2(points[i].y, points[i].x);
}
for (int i = 0; i < rightPoints.Length; i++)
{
rightPoints[i] = new Vec2(points[i + middleIdx].y, points[i + middleIdx].x);
}
// 分而治之: 算出左边部分的最短距离和右边部分的最短距离。
float l = MinDistanceSquared(leftPoints);
float r = MinDistanceSquared(rightPoints);
float minDistance = Math.Min(l, r);
// 靠近中间线的点，有可能跨越中间线出现更短的距离。
// 但该中间区域有很大局限。假定最短距离是s，中线的x为X，则该区域限定在X-S ~ X+S之间
float middle = points[middleIdx].x;
List<Vec2> middleStrip = new List<Vec2>();
foreach (Vec2 p in points)
{
if ((p.x - middle) * (p.x - middle) < minDistance) middleStrip.Add(new Vec2(p.y, p.x));
}
// 不同于wiki的算法描述，这里再次采用分而治之的原则（代码简洁很多），算出中间地带的最短距离
float m = MinDistanceSquared(middleStrip.ToArray());
return Math.Min(minDistance, m);
}
static float BruteForceMinDistanceSquared(Vec2[] points)
{
float minDistance = float.MaxValue;
// 穷举法，O(n*n)
for (int i = 0; i < points.Length; i++)
{
for (int j = i + 1; j < points.Length; j++)
{
float distance = (points[i] - points[j]).LengthSquare();
if (distance < minDistance) minDistance = distance;
}
}
return minDistance;
}
}
struct Vec2
{
public float x, y;
public Vec2(float x, float y)
{
this.x = x; this.y = y;
}
public float LengthSquare()
{
return x * x + y * y;
}
public static int CompareByX(Vec2 v1, Vec2 v2)
{
return v1.x - v2.x > 0 ? 1 : (v1.x - v2.x == 0 ? 0 : -1);
}
public static Vec2 operator -(Vec2 v1, Vec2 v2)
{
return new Vec2(v1.x - v2.x, v1.y - v2.y);
}
}
发表于 @ 2008年07月22日 17:06:00|评论(loading...)|收藏