Sadly I am not that good with math and I am working on a Delaunay triangulation of a 2D point field. So far I have the field and “points” are being triangulated. In the example 79 triangles were calculated for 19 points.
The problem is that too many triangles are being made and points are connected to incorrect neighbours.
I hope someone better at math can take a quick look at the code I have currently and help me with this.
References used:
https://github.com/Yonaba/delaunay/blob/master/delaunay.lua#L16
https://github.com/jbegaint/delaunay-cpp/blob/master/delaunay.hpp#L176
https://www.youtube.com/watch?v=4TtMmTOuqBk
Result:
.cpp:
#include "Delaunay.h"
TArray<FDTriangle> FDelaunay::Triangulate(TArray<FDPoint>& InPoints, int32 InDelaunayConvexMultiplier) const {
TArray<FDTriangle> Triangles;
int32 NPoints = InPoints.Num();
if (NPoints < 3) {
UE_LOG(LogActor, Error, TEXT("Triangulate needs at least 3 points."));
return Triangles;
}
if (NPoints == 3) {
Triangles.Add(FDTriangle(InPoints[0], InPoints[1], InPoints[2]));
return Triangles;
}
// Start (Bowyer Watson) Delaunay triangulation.
// Get the max amount of expected triangles.
int32 TrMax = NPoints * 4;
// Get the min / max dimensions of the grid containing the points.
float MinX = InPoints[0].X;
float MinY = InPoints[0].Y;
float MaxX = MinX;
float MaxY = MinY;
for (int32 i = 0; i < NPoints; i++) {
FDPoint& Point = InPoints[i];
Point.Id = i;
MinX = FMath::Min(MinX, Point.X);
MaxX = FMath::Max(MaxX, Point.X);
MinY = FMath::Min(MinY, Point.Y);
MaxY = FMath::Max(MaxY, Point.Y);
}
float Dx = ((MaxX - MinX) * InDelaunayConvexMultiplier);
float Dy = ((MaxY - MinY) * InDelaunayConvexMultiplier);
float DeltaMax = FMath::Max(Dx, Dy);
float MidX = (MinX + MaxX) * 0.5f;
float MidY = (MinY + MaxY) * 0.5f;
// Add Super Triangle. For simplicity add the generated Super points on top of the point array.
FDPoint SuP1 = FDPoint(MidX - 2.f * DeltaMax, MidY - DeltaMax, NPoints);
FDPoint SuP2 = FDPoint(MidX, MidY + 2.f * DeltaMax, NPoints + 1);
FDPoint SuP3 = FDPoint(MidX + 2.f * DeltaMax, MidY - DeltaMax, NPoints + 2);
InPoints.EmplaceAt(SuP1.Id, SuP1);
InPoints.EmplaceAt(SuP2.Id, SuP2);
InPoints.EmplaceAt(SuP3.Id, SuP3);
Triangles.Add(FDTriangle(SuP1, SuP2, SuP3));
// Use NPoints instead of InPoints.Num() when looping over the original points, to not include the initial super points.
for (int32 i = 0; i < NPoints; i++) {
TArray<FDEdge> Edges;
// For each point, look which triangles their CircumCircle contains this point
// , in which case we don't want the triangle because it is not a Delaunay triangle.
for (int32 j = Triangles.Num() - 1; j >= 0; j--) {
FDTriangle CurTriangle = Triangles[j];
if (CurTriangle.IsInCircumCircle(InPoints[i])) {
// For each edge, add if not yet present by IsSimilar
auto AddEdgeIfNotSimilar = [&Edges](const FDEdge& InEdge) {
bool bAdd = true;
for (const FDEdge& EdgeX : Edges) {
if (EdgeX.IsSimilar(InEdge)) {
bAdd = false;
break;
}
}
if (bAdd) {
Edges.Add(InEdge);
}
};
AddEdgeIfNotSimilar(CurTriangle.E1);
AddEdgeIfNotSimilar(CurTriangle.E2);
AddEdgeIfNotSimilar(CurTriangle.E3);
// Remove the unwanted triangle
Triangles.RemoveAt(j);
}
}
for (int32 j = 0; j < Edges.Num(); j++) {
if (Triangles.Num() > TrMax) {
UE_LOG(LogActor, Error, TEXT("Made more triangles than required. "));
}
Triangles.Add(FDTriangle(Edges[j].P1, Edges[j].P2, InPoints[i]));
}
}
// Remove triangles using the initial super points.
for (int i = Triangles.Num() - 1; i >= 0; i--) {
FDTriangle Triangle = Triangles[i];
if (Triangle.P1.Id >= NPoints
|| Triangle.P2.Id >= NPoints
|| Triangle.P3.Id >= NPoints
) {
Triangles.RemoveAt(i);
}
}
// Remove the initial super points.
InPoints.Remove(SuP1);
InPoints.Remove(SuP2);
InPoints.Remove(SuP3);
return Triangles;
}
.h:
#pragma once
#include "CoreMinimal.h"
struct FDPoint {
float X;
float Y;
int32 Id;
// Initialize
FDPoint(float InX, float InY, int32 InId)
: X (InX)
, Y (InY)
, Id (InId)
{}
// Functions
float GetDistSqr(const FVector2D& InPoint) const {
return FVector2D::DistSquared(InPoint, FVector2D(X, Y));
}
float GetDist(const FVector2D& InPoint) const {
return FMath::Sqrt(GetDistSqr(InPoint));
}
//bool IsInCircle(const FVector2D& InCenter, float InRadius) const {
// return FMath::Square(InCenter.X - X) + FMath::Square(InCenter.Y - Y) <= FMath::Square(InRadius);
//}
bool IsNearlyEqual(const FDPoint& Other) const {
return (FMath::IsNearlyEqual(X, Other.X)
&& FMath::IsNearlyEqual(Y, Other.Y)
);
}
// Operators
bool operator==(const FDPoint& InOther) const {
return (
Id == InOther.Id
);
}
};
struct FDEdge {
FDPoint P1;
FDPoint P2;
// Initialize
FDEdge(const FDPoint& InP1, const FDPoint& InP2)
: P1 (InP1)
, P2 (InP2)
{}
// Functions
float GetLength() const {
return FVector2D::Distance(FVector2D(P1.X, P1.Y), FVector2D(P2.X, P2.Y));
}
//FVector2D GetMidPoint() const {
// return FVector2D(
// (P1.X + P2.X) / 2.f
// , (P1.Y + P2.Y) / 2.f
// );
//}
//float GetSlope() const {
// return ((P2.Y - P1.Y) / (P2.X / P1.X));
//}
//float GetSlopePerpendicular() const {
// return (-1 * (1 / GetSlope()));
//}
bool IsSimilar(const FDEdge& InEdge) const {
// Returns true if P1 and P2 match in any direction
return ((P1.IsNearlyEqual(InEdge.P1) && P2.IsNearlyEqual(InEdge.P2))
|| (P2.IsNearlyEqual(InEdge.P1) && P1.IsNearlyEqual(InEdge.P2))
);
}
bool IsNearlyEqual(const FDEdge& InEdge) const {
return (P1.IsNearlyEqual(InEdge.P1) && P2.IsNearlyEqual(InEdge.P2));
}
// Operators
bool operator==(const FDEdge& InOther) const {
return (
P1 == InOther.P1
&& P2 == InOther.P2
);
}
};
struct FDTriangle {
FDPoint P1;
FDPoint P2;
FDPoint P3;
FDEdge E1;
FDEdge E2;
FDEdge E3;
// Initialize
FDTriangle(const FDPoint& InP1, const FDPoint& InP2, const FDPoint& InP3)
: P1 (InP1)
, P2 (InP2)
, P3 (InP3)
, E1 (FDEdge(InP1, InP2))
, E2 (FDEdge(InP2, InP3))
, E3 (FDEdge(InP3, InP1))
{}
// Functions
float QuatCross(float InA, float InB, float InC) const {
return FMath::Sqrt(
(InA + InB + InC)
* (InA + InB - InC)
* (InA - InB + InC)
* (-InA + InB + InC)
);
}
//float CrossProduct(const FDPoint& InP1, const FDPoint& InP2, const FDPoint& InP3) const {
// float X1 = InP2.X - InP1.X;
// float X2 = InP3.X - InP2.X;
// float Y1 = InP2.Y - InP1.Y;
// float Y2 = InP3.Y - InP2.Y;
// return X1 * Y2 - Y1 * X2;
//}
//bool IsFlatAngle(const FDPoint& InP1, const FDPoint& InP2, const FDPoint& InP3) const {
// return CrossProduct(InP1, InP2, InP3) == 0;
//}
//bool IsCW() const {
// return (CrossProduct(P1, P2, P3) < 0);
//}
//bool IsCCW() const {
// return (CrossProduct(P1, P2, P3) > 0);
//}
FVector GetSidesLength() const {
return FVector(E1.GetLength(), E2.GetLength(), E3.GetLength());
}
FVector2D GetCenter() const {
return FVector2D(
(P1.X + P2.X + P3.X) / 3.f
, (P1.Y + P2.Y + P3.Y) / 3.f
);
}
float GetCircumRadius() const {
FVector SidesLength = GetSidesLength();
float Cross = QuatCross(SidesLength.X, SidesLength.Y, SidesLength.Z);
// Safe division
return (Cross != 0.f ? (SidesLength.X * SidesLength.Y * SidesLength.Z) / Cross : 0.f);
}
FVector2D GetCircumCenter() const {
//FVector2D ME1 = E1.GetMidPoint();
//FVector2D ME2 = E2.GetMidPoint();
//float SE1 = E1.GetSlope();
//float SE2 = E2.GetSlope();
//float PE1 = E1.GetSlopePerpendicular();
//float PE2 = E2.GetSlopePerpendicular();
float D = (P1.X * (P2.Y - P3.Y) + P2.X * (P3.Y - P1.Y) + P3.X * (P1.Y - P2.Y)) * 2;
float X = ((P1.X * P1.X + P1.Y * P1.Y) * (P2.Y - P3.Y) + (P2.X * P2.X + P2.Y * P2.Y) * (P3.Y - P1.Y) + (P3.X * P3.X + P3.Y * P3.Y) * (P1.Y - P2.Y));
float Y = ((P1.X * P1.X + P1.Y * P1.Y) * (P3.X - P2.X) + (P2.X * P2.X + P2.Y * P2.Y) * (P1.X - P3.X) + (P3.X * P3.X + P3.Y * P3.Y) * (P2.X - P1.X));
return (D != 0.f ? FVector2D(X / D, Y / D) : FVector2D(0.f, 0.f));
}
float GetArea() const {
FVector SidesLength = GetSidesLength();
return (QuatCross(SidesLength.X, SidesLength.Y, SidesLength.Z) / 4.f);
}
bool IsInCircumCircle(const FDPoint& InPoint) const {
FVector2D Center = GetCircumCenter();
return (FMath::Square(Center.X - InPoint.X)
+ FMath::Square(Center.Y - InPoint.Y)
<= FMath::Square(GetCircumRadius())
);
}
bool IsNearlyEqual(const FDTriangle& InTriangle) const {
return (P1.IsNearlyEqual(InTriangle.P1)
&& P2.IsNearlyEqual(InTriangle.P2)
&& P3.IsNearlyEqual(InTriangle.P3)
);
}
// Operators
bool operator==(const FDTriangle& Other) const {
return (
P1 == Other.P1
&& P2 == Other.P2
&& P3 == Other.P3
);
}
public:
};
class FDelaunay {
private:
protected:
public:
// When small (~1) produces convex-hull, when large, produces concave hulls. Defaults to 1000.
//int32 ConvexMultiplier = 1;
private:
protected:
public:
TArray<FDTriangle> Triangulate(TArray<FDPoint>& InPoints, int32 InDelaunayConvexMultiplier) const;
};
Usage snippet:
#include "Delaunay.h"
#include "Kismet/KismetSystemLibrary.h"
#include "CoreMinimal.h"
TArray<FDPoint> Points;
Points.Add(10.f, 20.f, -1);
Points.Add(15.f, 3.f, -1);
Points.Add(6.f, 12.f, -1);
Points.Add(30.f, 40.f, -1);
Points.Add(25.f, 17.f, -1);
const TArray<FDTriangle> Triangles = Delaunay.Triangulate(Points, 1);
for (const FDTriangle& TriangleX : Triangles) {
for (int32 i = 0; i < 3; i++) {
const FDEdge EdgeX = (
i == 0 ? TriangleX.E1
: i == 1 ? TriangleX.E2
: i == 2 ? TriangleX.E3
// Invalid
: FDEdge(FDPoint(0.f, 0.f, -1), FDPoint(0.f, 0.f, -1))
);
// For each edge per triangle per room:
UKismetSystemLibrary::DrawDebugLine(
this
, FVector(EdgeX.P1.X, EdgeX.P1.Y, 0.f)
, FVector(EdgeX.P2.X, EdgeX.P2.Y, 0.f)
, FColor::Silver
, 5.f
, 8.f
);
}
}
