将凹壳算法转换为 c#答案

【问题标题】：Translating concave hull algorithm to c#将凹壳算法转换为 c#
【发布时间】：2013-05-06 21:24:19
【问题描述】：

所以我试图翻译这里找到的凹壳算法：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;
    }

}

【问题讨论】：

使用下面的代码并将new Polygon(...).Contains(..) 替换为光线套管算法有效。我很快就会 os 代码。
前几天我用标记答案和光线投射算法测试了代码，但它实际上不起作用。你做了哪些改变让它发挥作用？
好吧，我认为它有效。我将进一步调查我的数据以检查发生了什么。也许光线投射已关闭？
没有。问题是离开递归的条件。
我去看看。我最近没有太多时间回到这个话题。

标签： c# geometry concave-hull

【解决方案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 的角度，以-pi 和+pi 之间的弧度表示。不要混合度数和弧度。我的建议是始终使用弧度，并且仅在必要时转换为度数以实现人类可读的输出。

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 将顶点从point 与自身以及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 度。

【讨论】：

这个答案理论上是正确的。我不明白的是为什么要限制小于180度的角度。由于论文规定角度是按右旋排序的，即大于180度的角度是可能的，除非该顶点是先前计算的候选顶点。
@KenKin 我可能是错的，但这会计算两个角度之间的锐差，这两个角度永远不会 > 180 度。我相信通过取出AngleDifference 中的前两行，它将计算右转差异，也许还有b - a。

【解决方案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();
}

根据代码，如果您在相同距离处有多个顶点，则函数仅返回一个。由于 Dictionary 使用唯一键。

顺便说一句，有人完成了吗？

【讨论】：

我认为没有人这样做。当它得到回答时，我已经开始做其他工作了。

【解决方案3】：

我现在没有时间阅读这篇论文，但根据我对 conVEX hull 算法的了解，我假设您正在围绕特定方向的点寻找下一个要链接的点。

如果是这种情况，“角度”将是船体最近线段的角度，您希望按与该线的角度对点进行排序。因此，您想要计算一条线（在船体上）和一组线（从当前点到正在考虑的其他点）之间的角度。计算出的角度是正还是负取决于你是顺时针还是逆时针。要计算角度，请查看以下内容：

Calculating the angle between two lines without having to calculate the slope? (Java)

然后按角度排序。

【讨论】：

在这种情况下，OP 真的意味着凹壳。
这很酷，我只是说我对凸包的了解可能（或可能不）适用于这种情况。
好的，来自论文“SortByAngle[kNearestPoints,currentPoint,prevAngle] ► 将候选（邻居）按右转降序排序”，所以我相信我的回答基本上是正确的。跨度>

【解决方案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;
    }

【讨论】：