格雷厄姆扫描问题点数很高
Spa*_*n01 5 c# algorithm wpf computational-geometry
我的格雷厄姆扫描算法有一个问题,当我的列表有很多分数,但每次都很好,点数很少.我做了一些截图:
工作:(300分)
不工作(5000分)
角度计算:
public static double angle(MyVector3D vec1, MyVector3D vec2)
{
return Math.Atan2(vec2.Y - vec1.Y, vec2.X - vec1.X) * 180 / Math.PI;
}
极角分选点取决于最大Y点:
//bubblesort
private void sortList()
{
double temp = 0.0;
MyVector3D tempVector = new MyVector3D();
for (int i = 1; i < points.Count; i++)
{
for (int j = 1; j < points.Count - 1; j++)
{
if (angles[j] < angles[j + 1])
{
temp = angles[j + 1];
tempVector = points[j + 1];
angles[j + 1] = angles[j];
points[j + 1] = points[j];
angles[j] = temp;
points[j] = tempVector;
}
}
}
ccw方法:
private double ccw(MyVector3D vec1, MyVector3D vec2, MyVector3D vec3)
{
// ccwTest = ((vec2.X - vec1.X) * (vec3.Y - vec1.Y)) - ((vec2.Y - vec1.Y) * (vec3.X - vec1.X));
return ((vec2.X - vec1.X) * (vec3.Y - vec1.Y)) - ((vec2.Y - vec1.Y) * (vec3.X - vec1.X));
}
格雷厄姆扫描算法:
for (int i = 2; i < points.Count; i++)
{
while (ccw(points[M - 1], points[M], points[i]) > 0)
{
if (M > 1)
{
points.RemoveAt(M);
M -= 1;
i--;
}
else if (i == points.Count - 1)
{
break;
}
else
{
i += 1;
}
}
//goodPoints.Add(points[M]);
//listBoxInfos.Items.Add("g" + (int)points[M].X + "," + (int)points[M].Y + "," + 0);
//listBoxInfos.Items.Add("ccw" + ccwTest);
M += 1;
}
我真的不知道为什么我的程序会在800多点上爆炸......这很难调试,因为算法在300,400,500点时非常有效.
Ty为信息.
维基百科上的算法被打破了.它不处理多个点彼此共线且具有最小点的情况.例如,以下测试用例将失败:
Vector3[] points3 = new Vector3[]
{
new Vector3( 1, 1, 0),
new Vector3( 5, 5, 0),
new Vector3( 3, 3, 0),
new Vector3( 2, 2, 0),
new Vector3( 1, 1, 0),
new Vector3( 1, 10, 0),
};
问题在于,当沿着点行进时,如果该点位于船体的最后两个点之间,则可能需要丢弃当前点而不是延长船体或替换船体上的最后一个点.(只有当点也与最小点共线时才会发生这种情况,否则前面的角度排序会阻止这种双后退.)在所示的伪代码中没有逻辑.
我还认为维基百科算法可能会严重处理浮点舍入错误.特别是检查ccw <= 0看起来有问题.
这是尝试清理维基百科算法.我不得不摆脱(模糊地说是kludgy)"sentinal point",因为如果所有的点都是水平对齐的,它将基本上随机选择:
public static IList<Vector3> GrahamScanCompute(IList<Vector3> initialPoints)
{
if (initialPoints.Count < 2)
return initialPoints.ToList();
// Find point with minimum y; if more than one, minimize x also.
int iMin = Enumerable.Range(0, initialPoints.Count).Aggregate((jMin, jCur) =>
{
if (initialPoints[jCur].Y < initialPoints[jMin].Y)
return jCur;
if (initialPoints[jCur].Y > initialPoints[jMin].Y)
return jMin;
if (initialPoints[jCur].X < initialPoints[jMin].X)
return jCur;
return jMin;
});
// Sort them by polar angle from iMin,
var sortQuery = Enumerable.Range(0, initialPoints.Count)
.Where((i) => (i != iMin)) // Skip the min point
.Select((i) => new KeyValuePair<double, Vector3>(Math.Atan2(initialPoints[i].Y - initialPoints[iMin].Y, initialPoints[i].X - initialPoints[iMin].X), initialPoints[i]))
.OrderBy((pair) => pair.Key)
.Select((pair) => pair.Value);
List<Vector3> points = new List<Vector3>(initialPoints.Count);
points.Add(initialPoints[iMin]); // Add minimum point
points.AddRange(sortQuery); // Add the sorted points.
int M = 0;
for (int i = 1, N = points.Count; i < N; i++)
{
bool keepNewPoint = true;
if (M == 0)
{
// Find at least one point not coincident with points[0]
keepNewPoint = !NearlyEqual(points[0], points[i]);
}
else
{
while (true)
{
var flag = WhichToRemoveFromBoundary(points[M - 1], points[M], points[i]);
if (flag == RemovalFlag.None)
break;
else if (flag == RemovalFlag.MidPoint)
{
if (M > 0)
M--;
if (M == 0)
break;
}
else if (flag == RemovalFlag.EndPoint)
{
keepNewPoint = false;
break;
}
else
throw new Exception("Unknown RemovalFlag");
}
}
if (keepNewPoint)
{
M++;
Swap(points, M, i);
}
}
// points[M] is now the last point in the boundary. Remove the remainder.
points.RemoveRange(M + 1, points.Count - M - 1);
return points;
}
static void Swap<T>(IList<T> list, int i, int j)
{
if (i != j)
{
T temp = list[i];
list[i] = list[j];
list[j] = temp;
}
}
public static double RelativeTolerance { get { return 1e-10; } }
public static bool NearlyEqual(Vector3 a, Vector3 b)
{
return NearlyEqual(a.X, b.X) && NearlyEqual(a.Y, b.Y);
}
public static bool NearlyEqual(double a, double b)
{
return NearlyEqual(a, b, RelativeTolerance);
}
public static bool NearlyEqual(double a, double b, double epsilon)
{
// See here: http://floating-point-gui.de/errors/comparison/
if (a == b)
{ // shortcut, handles infinities
return true;
}
double absA = Math.Abs(a);
double absB = Math.Abs(b);
double diff = Math.Abs(a - b);
double sum = absA + absB;
if (diff < 4*double.Epsilon || sum < 4*double.Epsilon)
// a or b is zero or both are extremely close to it
// relative error is less meaningful here
return true;
// use relative error
return diff / (absA + absB) < epsilon;
}
static double CCW(Vector3 p1, Vector3 p2, Vector3 p3)
{
// Compute (p2 - p1) X (p3 - p1)
double cross1 = (p2.X - p1.X) * (p3.Y - p1.Y);
double cross2 = (p2.Y - p1.Y) * (p3.X - p1.X);
if (NearlyEqual(cross1, cross2))
return 0;
return cross1 - cross2;
}
enum RemovalFlag
{
None,
MidPoint,
EndPoint
};
static RemovalFlag WhichToRemoveFromBoundary(Vector3 p1, Vector3 p2, Vector3 p3)
{
var cross = CCW(p1, p2, p3);
if (cross < 0)
// Remove p2
return RemovalFlag.MidPoint;
if (cross > 0)
// Remove none.
return RemovalFlag.None;
// Check for being reversed using the dot product off the difference vectors.
var dotp = (p3.X - p2.X) * (p2.X - p1.X) + (p3.Y - p2.Y) * (p2.Y - p1.Y);
if (NearlyEqual(dotp, 0.0))
// Remove p2
return RemovalFlag.MidPoint;
if (dotp < 0)
// Remove p3
return RemovalFlag.EndPoint;
else
// Remove p2
return RemovalFlag.MidPoint;
}
顺便说一下,你的算法在几个地方是n阶平方:
- bubblesort是O(N ^ 2)的阶数
- 在找到船体时移除点而不是将它们交换到列表的前面可以是O(N ^ 2),因为所有后续点都必须向下移动.
让我知道如果它解决了你的问题,我已经测试了一下但不完全.
根据这个： http: //en.wikipedia.org/wiki/Graham_scan和其他人，您的 Graham 扫描算法实现至少存在 2 个问题，我认为您对较低的数字“很幸运”：
1）您在外部 for 和 else 中都增加了 i ，即通常您会跳过测试所有其他点。
2）我不相信你删除失败点的方法，是的，该点不在船体“此处”上，但可能是船体上更远的点，您需要向下交换这些点或使用堆栈基于的方法。