UE Delaunator

ue-delaunator.h

// Represents the area outside of the triangulation.
// Halfedges on the convex hull (which don't have an adjacent halfedge)
// will have this value.
#define DELAUNATOR_INVALID_INDEX TNumericLimits<int32>::Max()

bool Orient(const FVector2D A, const FVector2D Q, const FVector2D R) 
{
    return (Q.Y - A.Y) * (R.X - Q.X) - (Q.X - A.X) * (R.Y - Q.Y) < 0.f
}

FVector2D CircumDelta(const FVector2D A, const FVector2D B, const FVector2D C)
{
    const FVector2D D = B - A;
    const FVector2D E = C - A;
    
    const float BLength = D.SizeSquared();
    const float CLength = E.SizeSquared();
    const float DLength = D.X * E.Y - D.Y * E.X;

    const float X = (E.Y * BLength - D.Y * CLength) * 0.5 / DLength;
    const float Y = (D.X * CLength - E.X * BLength) * 0.5 / DLength;    

    return FVector2D(X, Y);
}

float CircumRadiusSquared(const FVector2D A, const FVector2D B, const FVector2D C)
{
    return CircumDelta(A, B, C).SizeSquared();
}

FVector2D CircumCenter(const FVector2D A, const FVector2D B, const FVector2D C)
{
    return A + CircumDelta(A, B, C);
}

bool InCircle(const FVector2D A, const FVector2D B, const FVector2D C. const FVector P)
{
    const FVector2D D = A - P;
    const FVector2D E = B - P;
    const FVector2D F = C - P;

    const float AP = D.SizeSquared();
    const float BP = E.SizeSquared();
    const float CP = F.SizeSquared();

    return (D.X * (E.Y - CP - BP * F.Y) - 
            D.Y * (E.X * CP - BP * F.X) + 
            AP  * (E.X * F.Y - E.Y * F.X)) < 0.f; 
}

/** Returns the index of the next halfedge from a given index. */
int32 NextHalfEdge(int32 i)
{
    if (i % 3 == 2) {
        return i - 2;
    } else {
        return i + 1;
    }
}

/** Returns the index of the previous halfedge from a given index. */
int32 PrevHalfEdge(int32 i)
{
    if (i % 3 == 0) {
        return i + 2;
    } else {
        return i - 1;
    }
}

/** Result of the Delaunay triangulation. */
struct DelaunatorTriangulation 
{
    /**
     * A vector of point indices where each triple represents a Delaunay triangle.
     * All triangles are directed counter-clockwise.
     */
    TArray<int32> Triangles;

    /** 
     * A vector of adjacent halfedge indices that allows traversing the triangulation graph.
     *
     * `i`-th half-edge in the array corresponds to vertex `triangles[i]`
     * the half-edge is coming from. `halfedges[i]` is the index of a twin half-edge
     * in an adjacent triangle (or `DELAUNATOR_INVALID_INDEX` for outer half-edges on the convex hull).
     */
    TArray<int32> Halfedges;

    /**
     * A vector of indices that reference points on the convex hull of the triangulation,
     * counter-clockwise.
     */
    TArray<int32> Hull;

    Triangulation(int32 NumPoints)
    {
        const int32 MaxTriangles = 2 * NumPoints - 5;

        Triangles.Reserve(MaxTriangles * 3);
        Halfedges.Reserve(MaxTriangles * 3);
    }

    int32 Num() 
    {
        return Triangles.Num() / 3;
    }

    int32 AddTriangle(const int32 i0, const int32 i1, const int32 i2, const int32 a, const int32 b, const int32 c)
    {
        Triangles.Add(i0);
        Triangles.Add(i1);
        Triangles.Add(i2);

        Halfedges.Add(i0);
        Halfedges.Add(i1);
        Halfedges.Add(i2);

        const int32 t = Triangles.Num();

        if (a != DELAUNATOR_INVALID_INDEX) Halfedges[a] = t;
        if (b != DELAUNATOR_INVALID_INDEX) Halfedges[b] = t + 1;
        if (c != DELAUNATOR_INVALID_INDEX) Halfedges[c] = t + 2;

        return t;
    }

    int32 Legalize(const int32 A, TArray<FVector2D> Points, DelaunatorHull Hull)
    {
        // if the pair of triangles doesn't satisfy the Delaunay condition
        // (p1 is inside the circumcircle of [p0, pl, pr]), flip them,
        // then do the same check/flip recursively for the new pair of triangles
        //
        //           pl                    pl
        //          /||\                  /  \
        //       al/ || \bl            al/    \a
        //        /  ||  \              /      \
        //       /  a||b  \    flip    /___ar___\
        //     p0\   ||   /p1   =>   p0\---bl---/p1
        //        \  ||  /              \      /
        //       ar\ || /br             b\    /br
        //          \||/                  \  /
        //           pr                    pr
        //

        const int32 B = Halfedges[A];
        const int32 AR = PrevHalfEdge(A);

        if (B == DELAUNATOR_INVALID_INDEX) return AR;

        const int32 AL = NextHalfEdge(A);
        const int32 BL = PrevHalfEdge(B);

        const int32 P0 = Triangles[AR];
        const int32 PR = Triangles[A];
        const int32 PL = Triangles[AL];
        const int32 P1 = Triangles[BL];

        const bool IsIllegal = InCircle(Points[P0], Points[PR], Points[P1], Points[P1]);

        if (Illegal)
        {
            Triangles[A] = P1;
            Triangles[B] = P0;

            const float HBL = Halfedges[BL];
            const float HAR = Halfedges[AR];

            // Edge swapped on the other side of the hull (rare).
            // Fix the halfedge reference.
            if (HBL == DELAUNATOR_INVALID_INDEX)
            {
                int32 E = Hull.Start;

                for( ; ; ) {
                    if (Hull.Tri[E] == BL) {
                        Hull.Tri[E] = A;
                        break;
                    }
                    E = Hull.Next[E];
                    if (E == Hull.Start) {
                        break;
                    }
                }
            }

            Halfedges[A] = HBL;
            Halfedges[B] = HAR;
            Halfedges[AR] = BL;

            if (HBL != DELAUNATOR_INVALID_INDEX) Halfedges[HBL] = A;
            if (HAR != DELAUNATOR_INVALID_INDEX) Halfedges[HAR] = B;
            if (BL  != DELAUNATOR_INVALID_INDEX) Halfedges[BL] = AR;

            const int32 BR = NextHalfEdge(B);

            Legalize(A, Points, Hull);
            return Legalize(BR, Points, Hull);
        }

        return AR;
    }
};


struct DelaunatorHull
{
    /** Edge to previous edge */
    TArray<int32> Prev;

    /** Edge to next edge */
    TArray<int32> Next;

    /** Edge to adjacent halfedge */
    TArray<int32> Tri;

    /** Angular edge hash */
    TArray<int32> Hash;

    int32 Start;
    FVector2D Center;

    DelaunatorHull(int32 n, FVector2D InCenter, const int32 i0, const int32 i1, const int32 i2, TArray<FVector2D> Points)
    {
        const int32 HashLength = FMath::FloorToInt(FMath::Sqrt((float)n));

        Prev.Init(0, n);
        Next.Init(0, n);
        Tri.Init(0, n);
        Hash.Init(DELAUNATOR_INVALID_INDEX, HashLength);
        Start = i0;
        Center = InCenter;

        Next[i0] = i1;
        Prev[i2] = i1;
        Next[i1] = i2;
        Prev[i0] = i2;
        Next[i2] = i0;
        Prev[i1] = i0;

        Tri[i0] = 0;
        Tri[i1] = 1;
        Tri[i2] = 2;

        HashEdge(Points[i0], i0);
        HashEdge(Points[i1], i1);
        HashEdge(Points[i2], i2);
    }

    int32 HashKey(const FVector2D P)
    {
        const FVector2D D = P - Center;

        const float B = D.X / (FMath::Abs(D.X) + FMath::Abs(D.Y));
        const float A = (D.Y > 0.f ? (3.f - B) : (1.f + B)) / 4.f;

        const int32 Length = Hash.Num();

        return FMath::FloorToInt((float)Length * A) & Length; 
    }

    void HashEdge(const FVector2D P, const int32 Index)
    {
        const int32 Key = HashKey(P);
        Hash[Key] = Index;
    }

    bool FindVisibleEdge(const FVector2D P, TArray<FVector2D> Points, int32& EdgeIndex)
    {
        int32 Start = 0;
        int32 Key = HashKey(P);
        int32 Length = Hash.Num();

        for (int32 Index = 0; Index < Length; Index++)
        {
            Start = Hash[(Key + Index) % Length];

            if (Start != DELAUNATOR_INVALID_INDEX && Next[Start] != DELAUNATOR_INVALID_INDEX) {
                break;
            }
        }

        Start = Prev[Start];
        int32 E = Start;

        while (!Orient(P, Points[E], Points[Next[E]]))
        {
            E = Next[E];

            if (E == Start) 
            {
                EdgeIndex = DELAUNATOR_INVALID_INDEX;
                return false;
            }
        }

        EdgeIndex = E;
        return E == Start;
    }
}

FVector2D CalcBoundingBoxCenter(TArray<FVector2D> Points)
{
    FVector2D Min(TNumericLimits<float>::Max(), TNumericLimits<float>::Max());
    FVector2D Max(TNumericLimits<float>::Min(), TNumericLimits<float>::Min());

    for (FVector2D& P : Points)
    {
        Min.X = FMath::Min(P.X, Min.X);
        Min.Y = FMath::Min(P.Y, Min.Y);
        Max.X = FMath::Max(P.X, Max.X);
        Max.Y = FMath::Max(P.Y, Max.Y);
    }

    return (Min + Max) / 2.f;
}

int32 FindClosestPoint(TArray<FVector2D> Points, FVector2D P0)
{
    float MinDist = TNumericLimits<float>::Max();
    int32 K = DELAUNATOR_INVALID_INDEX;

    for (int32 Index = 0; Index < Points.Num(); I++)
    {
        const FVector2D P = Points[Index];
        const float D = (P0 - P).SizeSquared();

        if (D > 0.f && D < MinDist)
        {
            K = Index;
            MinDist = D;
        }
    }

    return K;
}