将凸包算法翻译为C#

Question

我正在尝试翻译这里找到的凹壳算法： http://repositorium.sdum.uminho.pt/bitstream/1822/6429/1/ConcaveHull_ACM_MYS.pdf

(第65页)

我已经读完了整篇文章，但是我无法弄清楚如何实现sortByAngle和angle，我不确定应该在其中哪个方法。这是我目前的进展：

//Main method
public static Vertex[] ConcaveHull(Vertex[] points, int k = 3)
{
    if (k < 3)
        throw new ArgumentException("K is required to be 3 or more", "k");
    List<Vertex> hull = new List<Vertex>();
    //Clean first, may have lots of duplicates
    Vertex[] clean = RemoveDuplicates(points);
    if (clean.Length < 3)
        throw new ArgumentException("At least 3 dissimilar points reqired", "points");
    if (clean.Length == 3)//This is the hull, its already as small as it can be.
        return clean;
    if (clean.Length < k)
        throw new ArgumentException("K must be equal to or smaller then the amount of dissimilar points", "points");
    Vertex firstPoint = clean[0]; //TODO find mid point
    hull.Add(firstPoint);
    Vertex currentPoint = firstPoint;
    Vertex[] dataset = RemoveIndex(clean, 0);
    double previousAngle = 0;
    int step = 2;
    int i;
    while (((currentPoint != firstPoint) || (step == 2)) && (dataset.Length > 0))
    {
        if (step == 5)
            dataset = Add(dataset, firstPoint);
        Vertex[] kNearestPoints = nearestPoints(dataset, currentPoint, k);
        Vertex[] cPoints = sortByAngle(kNearestPoints, currentPoint, previousAngle);
        bool its = true;
        i = 0;
        while ((its) && (i < cPoints.Length))
        {
            i++;
            int lastPoint = 0;
            if (cPoints[0] == firstPoint)
                lastPoint = 1;
            int j = 2;
            its = false;
            while ((!its) && (j < hull.Count - lastPoint))
            {
                its = intersectsQ(hull[step - 1 - 1], cPoints[0], hull[step - i - j - 1], hull[step - j - 1]);
                j++;
            }
        }
        if (its)
        {
            return ConcaveHull(points, k + 1);
        }
        currentPoint = cPoints[0];
        hull.Add(currentPoint);
        previousAngle = angle(hull[step - 1], hull[step - 2]);
        dataset = RemoveIndex(dataset, 0);
        step++;
    }
    bool allInside = true;
    i = dataset.Length;
    while (allInside && i > 0)
    {
        allInside = new Polygon(dataset).Contains(currentPoint); //TODO havent finished ray casting yet.
        i--;
    }
    if (!allInside)
        return ConcaveHull(points, k + 1);
    return hull.ToArray();
}

private static Vertex[] Add(Vertex[] vs, Vertex v)
{
    List<Vertex> n = new List<Vertex>(vs);
    n.Add(v);
    return n.ToArray();
}

private static Vertex[] RemoveIndex(Vertex[] vs, int index)
{
    List<Vertex> removed = new List<Vertex>();
    for (int i = 0; i < vs.Length; i++)
        if (i != index)
            removed.Add(vs[i]);
    return removed.ToArray();
}

private static Vertex[] RemoveDuplicates(Vertex[] vs)
{
    List<Vertex> clean = new List<Vertex>();
    VertexComparer vc = new VertexComparer();
    foreach (Vertex v in vs)
    {
        if (!clean.Contains(v, vc))
            clean.Add(v);
    }
    return clean.ToArray();
}

private static Vertex[] nearestPoints(Vertex[] vs, Vertex v, int k)
{
    Dictionary<double, Vertex> lengths = new Dictionary<double, Vertex>();
    List<Vertex> n = new List<Vertex>();
    double[] sorted = lengths.Keys.OrderBy(d => d).ToArray();
    for (int i = 0; i < k; i++)
    {
        n.Add(lengths[sorted[i]]);
    }
    return n.ToArray();
}

private static Vertex[] sortByAngle(Vertex[] vs, Vertex v, double angle)
{
    //TODO
    return new Vertex[]{};
}

private static bool intersectsQ(Vertex v1, Vertex v2, Vertex v3, Vertex v4)
{
    return intersectsQ(new Edge(v1, v2), new Edge(v3, v4));
}

private static bool intersectsQ(Edge e1, Edge e2)
{
    double x1 = e1.A.X;
    double x2 = e1.B.X;
    double x3 = e2.A.X;
    double x4 = e2.B.X;

    double y1 = e1.A.Y;
    double y2 = e1.B.Y;
    double y3 = e2.A.Y;
    double y4 = e2.B.Y;

    var x = ((x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4)) / ((x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4));
    var y = ((x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4)) / ((x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4));
    if (double.IsNaN(x) || double.IsNaN(y))
    {
        return false;
    }
    else
    {
        if (x1 >= x2)
        {
            if (!(x2 <= x && x <= x1)) { return false; }
        }
        else
        {
            if (!(x1 <= x && x <= x2)) { return false; }
        }
        if (y1 >= y2)
        {
            if (!(y2 <= y && y <= y1)) { return false; }
        }
        else
        {
            if (!(y1 <= y && y <= y2)) { return false; }
        }
        if (x3 >= x4)
        {
            if (!(x4 <= x && x <= x3)) { return false; }
        }
        else
        {
            if (!(x3 <= x && x <= x4)) { return false; }
        }
        if (y3 >= y4)
        {
            if (!(y4 <= y && y <= y3)) { return false; }
        }
        else
        {
            if (!(y3 <= y && y <= y4)) { return false; }
        }
    }
    return true;
}

private static double angle(Vertex v1, Vertex v2)
{
    // TODO fix
    Vertex v3 = new Vertex(v1.X, 0);
    if (Orientation(v3, v1, v2) == 0)
        return 180;

    double b = EuclideanDistance(v3, v1);
    double a = EuclideanDistance(v1, v2);
    double c = EuclideanDistance(v3, v2);
    double angle = Math.Acos((Math.Pow(a, 2) + Math.Pow(b, 2) - Math.Pow(c, 2)) / (2 * a * b));

    if (Orientation(v3, v1, v2) < 0)
        angle = 360 - angle;

    return angle;
}

private static double EuclideanDistance(Vertex v1, Vertex v2)
{
    return Math.Sqrt(Math.Pow((v1.X - v2.X), 2) + Math.Pow((v1.Y - v2.Y), 2));
}

public static double Orientation(Vertex p1, Vertex p2, Vertex p)
{
    double Orin = (p2.X - p1.X) * (p.Y - p1.Y) - (p.X - p1.X) * (p2.Y - p1.Y);
    if (Orin > 0)
        return -1;//Left
    if (Orin < 0)
        return 1;//Right
    return 0;//Colinier
}

我知道这里有很多代码。但是我不确定是否可以在没有它的情况下展示上下文和我的内容。

其他类：

public class Polygon
{

    private Vertex[] vs;

    public Polygon(Vertex[] Vertexes)
    {
        vs = Vertexes;
    }

    public Polygon(Bounds bounds)
    {
        vs = bounds.ToArray();
    }

    public Vertex[] ToArray()
    {
        return vs;
    }

    public IEnumerable<Edge> Edges()
    {
        if (vs.Length > 1)
        {
            Vertex P = vs[0];
            for (int i = 1; i < vs.Length; i++)
            {
                yield return new Edge(P, vs[i]);
                P = vs[i];
            }
            yield return new Edge(P, vs[0]);
        }
    }

    public bool Contains(Vertex v)
    {
        return RayCasting.RayCast(this, v);
    }
}

public class Edge
{
    public Vertex A = new Vertex(0, 0);
    public Vertex B = new Vertex(0, 0);
    public Edge() { }
    public Edge(Vertex a, Vertex b)
    {
        A = a;
        B = b;
    }
    public Edge(double ax, double ay, double bx, double by)
    {
        A = new Vertex(ax, ay);
        B = new Vertex(bx, by);
    }
}

public class Bounds
{
    public Vertex TopLeft;
    public Vertex TopRight;
    public Vertex BottomLeft;
    public Vertex BottomRight;
    public Bounds() { }

    public Bounds(Vertex TL, Vertex TR, Vertex BL, Vertex BR)
    {
        TopLeft = TL;
        TopRight = TR;
        BottomLeft = BL;
        BottomRight = BR;
    }

    public Vertex[] ToArray()
    {
        return new Vertex[] { TopLeft, TopRight, BottomRight, BottomLeft };
    }

}

public class Vertex
{
    public double X = 0;
    public double Y = 0;
    public Vertex() { }
    public Vertex(double x, double y)
    {
        X = x;
        Y = y;
    }

    public static Vertex[] Convert(string vs)
    {
        vs = vs.Replace("[", "");
        vs = vs.Replace("]", "");
        string[] spl = vs.Split(';');
        List<Vertex> nvs = new List<Vertex>();
        foreach (string s in spl)
        {
            try
            {
                nvs.Add(new Vertex(s));
            }
            catch
            {

            }
        }
        return nvs.ToArray();
    }

    public static string Stringify(Vertex[] vs)
    {
        string res = "[";
        foreach (Vertex v in vs)
        {
            res += v.ToString();
            res += ";";
        }
        res = res.RemoveLastCharacter();
        res += "]";
        return res;
    }

    public static string ToString(Vertex[] array)
    {
        string res = "[";
        foreach (Vertex v in array)
            res += v.ToString() + ",";
        return res.RemoveLastCharacter() + "]";
    }

    /*
    //When x < y return -1
    //When x == y return 0
    //When x > y return 1
    public static int Compare(Vertex x, Vertex y)
    {
        //To find lowest
        if (x.X < y.X)
        {
            return -1;
        }
        else if (x.X == y.X)
        {
            if (x.Y < y.Y)
            {
                return -1;
            }
            else if (x.Y == y.Y)
            {
                return 0;
            }
            else
            {
                return 1;
            }
        }
        else
        {
            return 1;
        }
    }
    */
    public static int CompareY(Vertex a, Vertex b)
    {
        if (a.Y < b.Y)
            return -1;
        if (a.Y == b.Y)
            return 0;
        return 1;
    }

    public static int CompareX(Vertex a, Vertex b)
    {
        if (a.X < b.X)
            return -1;
        if (a.X == b.X)
            return 0;
        return 1;
    }

    public double distance (Vertex b){
        double dX = b.X - this.X;
        double dY = b.Y - this.Y;
        return Math.Sqrt((dX*dX) + (dY*dY));
    }

    public double slope (Vertex b){
        double dX = b.X - this.X;
        double dY = b.Y - this.Y;
        return dY / dX;
    }

    public static int Compare(Vertex u, Vertex a, Vertex b)
    {
        if (a.X == b.X && a.Y == b.Y) return 0;

        Vertex upper = new Vertex();
        Vertex p1 = new Vertex();
        Vertex p2 = new Vertex();
        upper.X = (u.X + 180) * 360;
        upper.Y = (u.Y + 90) * 180;
        p1.X = (a.X + 180) * 360;
        p1.Y = (a.Y + 90) * 180;
        p2.X = (b.X + 180) * 360;
        p2.Y = (b.Y + 90) * 180;
        if(p1 == upper) return -1;
        if(p2 == upper) return 1;

        double m1 = upper.slope(p1);
        double m2 = upper.slope(p2);

        if (m1 == m2)
        {
            return p1.distance(upper) < p2.distance(upper) ? -1 : 1;
        }

        if (m1 <= 0 && m2 > 0) return -1;

        if (m1 > 0 && m2 <= 0) return -1;

        return m1 > m2 ? -1 : 1;
    }

    public static Vertex UpperLeft(Vertex[] vs)
    {
        Vertex top = vs[0];
        for (int i = 1; i < vs.Length; i++)
        {
            Vertex temp = vs[i];
            if (temp.Y > top.Y || (temp.Y == top.Y && temp.X < top.X))
            {
                top = temp;
            }
        }
        return top;
    }

}

- FabianCook

---

我几天前使用标记答案和光线投射算法测试了代码，但实际上它并没有起作用。你做了哪些更改让它能够工作？ - Ken Kin

不，问题在于递归结束的条件。 - Ken Kin

你好。你能分享一份包含顶点/边缘/多边形类的完整代码吗？ - Nickon

完成，添加了类。 - FabianCook

谢谢你的课程。我有一个问题。它对你有效吗？因为我得到了1.4k分，但只输入了1k（看起来没问题，但是...）- 当我想要检查返回的分数时，所有分数都相同：S也许你完整的算法和上面的代码之间存在一些差异。我使用了这篇文章下面的人写的方法。 - Nickon

Answer 1

关于规范的一点提示：您应该以大写字母开头命名函数名称，以小写字母开头命名变量。在函数sortByAngle中，您同时引用了参数angle和函数angle。

假设Angle(...)仅仅是计算两个点之间的角度：

private static double Angle(Vertex v1, Vertex v2)
{
    return Math.Atan2(v2.Y - v1.Y, v2.X - v1.X);
}

将为您提供从v1到v2的角度，以弧度表示，范围在-π和+π之间。不要混淆度和弧度。我的建议是始终使用弧度，仅在需要人类可读输出时转换为度数。

private static Vertex[] SortByAngle(Vertex[] vs, Vertex v, double angle)
{
    List<Vertex> vertList = new List<Vertex>(vs);
    vertList.Sort((v1, v2) => AngleDifference(angle, Angle(v, v1)).CompareTo(AngleDifference(angle, Angle(v, v2))));
    return vertList.ToArray();
}

使用List.Sort函数将顶点按照与自身的角度差和angle从大到小排序，输入元组中v1和v2的顺序被交换以进行降序排序，即首先计算最大角度差。角度差的计算方式如下：

private static double AngleDifference(double a, double b)
{
    while (a < b - Math.PI) a += Math.PI * 2;
    while (b < a - Math.PI) b += Math.PI * 2;

    return Math.Abs(a - b);
}

前两行确保角度之间不超过180度。

- Elle

---

这个答案在理论上是正确的。我不明白的是为什么要限制角度小于180度。因为文件指定角度按右手拐弯排序，也就是说，大于180度的角度是可能存在的，除非顶点在之前的计算中已经被排除了。 - Ken Kin

@KenKin 我可能错了，但这个计算两个角度之间的锐角差异，其值永远不会大于180度。我认为通过删除 AngleDifference 中的前两行，它将计算正确的右转差异，也许还可以使用 b - a。 - Elle

Answer 2

您的代码出现了错误

private static Vertex[] nearestPoints(Vertex[] vs, Vertex v, int k)
{
    Dictionary<double, Vertex> lengths = new Dictionary<double, Vertex>();
    List<Vertex> n = new List<Vertex>();
    double[] sorted = lengths.Keys.OrderBy(d => d).ToArray();
    for (int i = 0; i < k; i++)
    {
        n.Add(lengths[sorted[i]]);
    }
    return n.ToArray();
}

根据代码的逻辑，如果有多个顶点距离相同时，函数只返回一个。因为字典使用唯一键。

顺便问一下，有人完成了吗？

- nredko

---

我不认为有人回答了。在问题得到解答的时候，我已经转而处理其他工作了。 - FabianCook

Answer 3

我现在没有时间阅读这篇论文，但是根据我对凸包算法的了解，我认为你正在按照特定方向绕着点寻找下一个要连接的点。

如果是这样的话，“角度”就是凸包最近线段的角度，你需要按照每个点与该线之间的角度对点进行排序。因此，你需要计算凸包上一条线段和当前点到其他考虑点的连线之间的角度。计算出的角度是正数还是负数取决于你是顺时针还是逆时针。可以参考以下内容来计算角度：

如何在不计算斜率的情况下计算两条直线之间的夹角？（Java）

然后按照角度进行排序即可。

- David Cummins

---

在这种情况下，OP实际上是指凹壳体。 - madth3

很酷，我只是想说我的凸包知识可能（或者可能不）适用于这种情况。 - David Cummins

好的，根据论文中的“SortByAngle[kNearestPoints,currentPoint,prevAngle] ► sort the candidates (neighbours) in descending order of right-hand turn”这句话，我相信我的答案基本上是正确的。 - David Cummins

Answer 4

这个怎么样？

private List<Vector> sortClockwiseFromCentroid(List<Vector> points, Vector center)
    {
        points = points.OrderBy(x => Math.Atan2(x.X - center.X, x.Y - center.Y)).ToList();
        return points;
    }

- GeoGecco