Adds FloodFilling

Config/DefaultEngine.ini

@@ -122,8 +122,10 @@ r.DynamicGlobalIlluminationMethod = 1
 r.ReflectionMethod = 1
 
 r.Shadow.Virtual.Enable = 1
-r.Mobile.AntiAliasing=0
-r.AntiAliasingMethod=0
+r.Mobile.AntiAliasing=3
+r.AntiAliasingMethod=3
+vr.InstancedStereo=True
+r.MobileHDR=True
 
 [/Script/HardwareTargeting.HardwareTargetingSettings]
 TargetedHardwareClass = Mobile

Source/MetaCastBachelor/DensityField.cpp

@@ -7,12 +7,12 @@
 
 FDensityField::FDensityField(): XNum(0), YNum(0), ZNum(0), XStep(0), YStep(0), ZStep(0)
 {
-	GridNodes = TArray<FVoxel>();
+	VoxelList = TArray<FVoxel>();
 }
 
 void FDensityField::InitializeDensityField(const FVector& MinBounds, const FVector& MaxBounds, const int32 XAxisNum, const int32 YAxisNum, const int32 ZAxisNum)
 {
-	GridNodes.Empty();
+	VoxelList.Empty();
 
 	XStep = (MaxBounds.X - MinBounds.X) / (XAxisNum - 1);
 	YStep = (MaxBounds.Y - MinBounds.Y) / (YAxisNum - 1);
@@ -26,7 +26,7 @@ void FDensityField::InitializeDensityField(const FVector& MinBounds, const FVect
 			{
 				FVector Position(MinBounds.X + X * XStep, MinBounds.Y + Y * YStep, MinBounds.Z + Z * ZStep);
 				FVector GridPosition(X, Y, Z);
-				GridNodes.Add(FVoxel(Position, GridPosition));
+				VoxelList.Add(FVoxel(Position, GridPosition));
 			}
 		}
 	}
@@ -46,10 +46,10 @@ void FDensityField::CalculateVoxelDensities(const APointCloud* PointCloud, const
 	UE_LOG(LogTemp, Log, TEXT("Starting density calculation."));
 
 	// Clear existing densities and counts
-	for (auto& Node : GridNodes)
+	for (auto& Node : VoxelList)
 	{
-	    Node.SetNodeDensity(0.0);
-	    Node.SetEnclosedParticleNum(0.0);  // Reset the counter
+	    Node.SetVoxelDensity(0.0);
+	    Node.SetClosePointsNumber(0.0);  // Reset the counter
 	}
 
 	UE_LOG(LogTemp, Log, TEXT("Cleared previous densities."));
@@ -70,17 +70,17 @@ void FDensityField::CalculateVoxelDensities(const APointCloud* PointCloud, const
 	        {
 	            for (int32 X = StartX; X <= EndX; ++X)
 	            {
-	                const int32 Index = GridPositionToIndex(Z, Y, X);
-	                if (Index >= 0 && Index < GridNodes.Num())
+	                const int32 Index = GridPositionToIndex(X, Y, Z);
+	                if (Index >= 0 && Index < VoxelList.Num())
 	                {
-	                    FVector NodePosition = GridNodes[Index].GetNodePosition();
+	                    FVector NodePosition = VoxelList[Index].GetVoxelPosition();
 	                    const float Distance = FVector::Dist(NodePosition, ParticlePosition);
 
 	                    if (Distance < InfluenceRadius)
 	                    {
 	                        const double Weight = FUtilities::GaussianKernel(Distance, Sigma);
-	                        GridNodes[Index].NodeDensityPlusDis(Weight);
-	                        GridNodes[Index].SetEnclosedParticleNum(GridNodes[Index].GetEnclosedParticleNum() + 1.0);
+	                        VoxelList[Index].AddToVoxelDensity(Weight);
+	                        VoxelList[Index].SetClosePointsNumber(VoxelList[Index].GetClosePointsNumber() + 1.0);
 	                    }
 	                }
 	            }
@@ -88,11 +88,11 @@ void FDensityField::CalculateVoxelDensities(const APointCloud* PointCloud, const
 	    }
 	}
 	double MaxDensity = 0.0;
-	for (const auto& Node : GridNodes)
+	for (const auto& Node : VoxelList)
 	{
-		if (Node.GetNodeDensity() > MaxDensity)
+		if (Node.GetVoxelDensity() > MaxDensity)
 		{
-			MaxDensity = Node.GetNodeDensity();
+			MaxDensity = Node.GetVoxelDensity();
 		}
 	}
 	UE_LOG(LogTemp, Log, TEXT("Maximum density found: %f"), MaxDensity);	
@@ -114,17 +114,17 @@ void FDensityField::CalculateAndStoreGradients()
 		{
 			for (int32 X = 0; X < XNum; X++)
 			{
-				const int32 Index = GridPositionToIndex(Z, Y, X);
-				if (!GridNodes.IsValidIndex(Index))
+				const int32 Index = GridPositionToIndex(X, Y, Z);
+				if (!VoxelList.IsValidIndex(Index))
 					continue;
 
 				// Use one-sided differences at the boundaries
-				const double DensityX1 = X < XNum - 1 ? GridNodes[GridPositionToIndex(Z, Y, X + 1)].GetNodeDensity() : GridNodes[Index].GetNodeDensity();
-				const double DensityX0 = X > 0 ? GridNodes[GridPositionToIndex(Z, Y, X - 1)].GetNodeDensity() : GridNodes[Index].GetNodeDensity();
-				const double DensityY1 = Y < YNum - 1 ? GridNodes[GridPositionToIndex(Z, Y + 1, X)].GetNodeDensity() : GridNodes[Index].GetNodeDensity();
-				const double DensityY0 = Y > 0 ? GridNodes[GridPositionToIndex(Z, Y - 1, X)].GetNodeDensity() : GridNodes[Index].GetNodeDensity();
-				const double DensityZ1 = Z < ZNum - 1 ? GridNodes[GridPositionToIndex(Z + 1, Y, X)].GetNodeDensity() : GridNodes[Index].GetNodeDensity();
-				const double DensityZ0 = Z > 0 ? GridNodes[GridPositionToIndex(Z - 1, Y, X)].GetNodeDensity() : GridNodes[Index].GetNodeDensity();
+				const double DensityX1 = X < XNum - 1 ? VoxelList[GridPositionToIndex(X + 1, Y, Z)].GetVoxelDensity() : VoxelList[Index].GetVoxelDensity();
+				const double DensityX0 = X > 0 ? VoxelList[GridPositionToIndex(X - 1, Y, Z)].GetVoxelDensity() : VoxelList[Index].GetVoxelDensity();
+				const double DensityY1 = Y < YNum - 1 ? VoxelList[GridPositionToIndex(X, Y + 1, Z)].GetVoxelDensity() : VoxelList[Index].GetVoxelDensity();
+				const double DensityY0 = Y > 0 ? VoxelList[GridPositionToIndex(X, Y - 1, Z)].GetVoxelDensity() : VoxelList[Index].GetVoxelDensity();
+				const double DensityZ1 = Z < ZNum - 1 ? VoxelList[GridPositionToIndex(X, Y, Z + 1)].GetVoxelDensity() : VoxelList[Index].GetVoxelDensity();
+				const double DensityZ0 = Z > 0 ? VoxelList[GridPositionToIndex(X, Y, Z - 1)].GetVoxelDensity() : VoxelList[Index].GetVoxelDensity();
 
 				FVector Gradient(
 					(DensityX1 - DensityX0) / (X == 0 || X == XNum - 1 ? XStep : 2 * XStep),
@@ -132,7 +132,7 @@ void FDensityField::CalculateAndStoreGradients()
 					(DensityZ1 - DensityZ0) / (Z == 0 || Z == ZNum - 1 ? ZStep : 2 * ZStep)
 				);
 
-				GridNodes[Index].SetNodeGradient(Gradient);
+				VoxelList[Index].SetVoxelGradient(Gradient);
 			}
 		}
 	}
@@ -151,18 +151,18 @@ void FDensityField::DrawDebugVoxelDensity(const UWorld* World, const AActor* Act
 	}
 
 	double MinDensity = DBL_MAX, MaxDensity = 0;
-	for (const FVoxel& Node : GridNodes)
+	for (const FVoxel& Node : VoxelList)
 	{
-		const double LogDensity = log10(Node.GetNodeDensity() + 1);  // Avoiding log(0)
+		const double LogDensity = log10(Node.GetVoxelDensity() + 1);  // Avoiding log(0)
 		if (LogDensity < MinDensity) MinDensity = LogDensity;
 		if (LogDensity > MaxDensity) MaxDensity = LogDensity;
 	}
 
 	const FTransform ActorTransform = Actor->GetActorTransform();
 
-	for (const FVoxel& Node : GridNodes)
+	for (const FVoxel& Node : VoxelList)
 	{
-		const double Density = Node.GetNodeDensity();
+		const double Density = Node.GetVoxelDensity();
 
 		const double LogDensity = log10(Density + 1);  // Normalize using logarithmic scale
 		const float NormalizedDensity = static_cast<float>((LogDensity - MinDensity) / (MaxDensity - MinDensity));
@@ -172,7 +172,7 @@ void FDensityField::DrawDebugVoxelDensity(const UWorld* World, const AActor* Act
 			FLinearColor VoxelColor = FLinearColor::LerpUsingHSV(FLinearColor::Green, FLinearColor::Red, NormalizedDensity);
 
 			FVector VoxelSize = FVector(Scale * XStep, Scale * YStep, Scale * ZStep) / 2.0f;
-			FVector WorldPosition = ActorTransform.TransformPosition(Node.GetNodePosition());
+			FVector WorldPosition = ActorTransform.TransformPosition(Node.GetVoxelPosition());
 			FVector ScaledVoxelSize = ActorTransform.TransformVector(VoxelSize);
 
 			DrawDebugBox(World, WorldPosition, 0.95 * ScaledVoxelSize, VoxelColor.ToFColor(true), false, Duration, 0, 0.1f);
@@ -185,30 +185,22 @@ void FDensityField::DrawDebugVoxelDensity(const UWorld* World, const AActor* Act
 
 //	CONVERSION FUNCTIONS
 
-int32 FDensityField::GridPositionToIndex(const int32 Z, const int32 Y, const int32 X) const
+int32 FDensityField::GridPositionToIndex(const int32 X, const int32 Y, const int32 Z) const
 {
     return Z * XNum * YNum + Y * XNum + X;
 }
 
-int32 FDensityField::WorldPositionToIndex(const FVector &Position, const APointCloud * PointCloud) const {
-	if (!PointCloud) {
-		UE_LOG(LogTemp, Warning, TEXT("PointCloud is null, cannot convert position to box index."));
-		return -1;
-	}
-
+int32 FDensityField::WorldPositionToIndex(const FVector &Position) const {
 	// Convert world position to grid position
-	FVector GridPosition = WorldToGridPosition(Position);
+	const FVector GridPosition = WorldToGridPosition(Position);
 
 	// Check if the grid position is valid
-	if (GridPosition.X < 0 || GridPosition.Y < 0 || GridPosition.Z < 0) {
+	if (!IsValidGridPosition(GridPosition.X, GridPosition.Y, GridPosition.Z)) {
 		UE_LOG(LogTemp, Warning, TEXT("Position out of grid bounds or invalid grid settings, cannot find box index."));
 		return -1;  // Grid position was invalid
 	}
 
-	// Calculate the linear index from grid position
-	const int32 Index = static_cast<int32>(GridPosition.Z) * XNum * YNum + static_cast<int32>(GridPosition.Y) * XNum + static_cast<int32>(GridPosition.X);
-
-	return Index;
+	return GridPositionToIndex(GridPosition.X, GridPosition.Y, GridPosition.Z);
 }
 
 FVector FDensityField::WorldToGridPosition(const FVector& Position) const
@@ -220,9 +212,9 @@ FVector FDensityField::WorldToGridPosition(const FVector& Position) const
 	}
 
 	// Calculate the grid coordinates by subtracting the grid origin and dividing by the grid step sizes
-	int32 GridX = static_cast<int32>((Position.X - GridOrigin.X) / XStep);
-	int32 GridY = static_cast<int32>((Position.Y - GridOrigin.Y) / YStep);
-	int32 GridZ = static_cast<int32>((Position.Z - GridOrigin.Z) / ZStep);
+	const int32 GridX = FMath::FloorToInt((Position.X - GridOrigin.X) / XStep);
+	const int32 GridY = FMath::FloorToInt((Position.Y - GridOrigin.Y) / YStep);
+	const int32 GridZ = FMath::FloorToInt((Position.Z - GridOrigin.Z) / ZStep);
 
 	// Check if the calculated grid coordinates are within the bounds
 	if (GridX < 0 || GridX >= XNum || GridY < 0 || GridY >= YNum || GridZ < 0 || GridZ >= ZNum)
@@ -234,7 +226,7 @@ FVector FDensityField::WorldToGridPosition(const FVector& Position) const
 	return FVector(GridX, GridY, GridZ);
 }
 
-bool FDensityField::IsValidIndex(const int32 XBin, const int32 YBin, const int32 ZBin) const
+bool FDensityField::IsValidGridPosition(const int32 XBin, const int32 YBin, const int32 ZBin) const
 {
 	// Check if the provided indices are within the grid bounds
 	return (XBin >= 0 && XBin < XNum) &&
@@ -242,103 +234,99 @@ bool FDensityField::IsValidIndex(const int32 XBin, const int32 YBin, const int32
 		   (ZBin >= 0 && ZBin < ZNum);
 }
 
-
-//	GETTER AND SETTER FUNCTIONS
-
-FVector FDensityField::GetGridOrigin() const
-{
-	return GridOrigin;
-}
-
-FVector FDensityField::GetNodeGradient(const int32 Index) const
+FVector FDensityField::IndexToVoxelGradient(const int32 Index) const
 {
-    if (GridNodes.IsValidIndex(Index))
+	if (VoxelList.IsValidIndex(Index))
 	{
-        return GridNodes[Index].GetNodeGradient();
+		return VoxelList[Index].GetVoxelGradient();
 	}
 	return FVector::ZeroVector;
 }
 
-FVector FDensityField::GetNodePosition(const int32 Index) const
+FVector FDensityField::IndexToVoxelPosition(const int32 Index) const
 {
-    if (GridNodes.IsValidIndex(Index))
+	if (VoxelList.IsValidIndex(Index))
 	{
-        return GridNodes[Index].GetNodePosition();
+		return VoxelList[Index].GetVoxelPosition();
 	}
-    return FVector::ZeroVector;
+	return -FVector::One();
 }
 
-FVector FDensityField::GetNodeGridPos(const int32 Index) const
+FVector FDensityField::IndexToGridPosition(const int32 Index) const
 {
-    if (GridNodes.IsValidIndex(Index))
+	if (IsValidIndex(Index))
 	{
-        return GridNodes[Index].GetNodeGridPos();
+		return VoxelList[Index].GetVoxelGridPos();
 	}
-    return FVector::ZeroVector;
+	return -FVector::One();
 }
 
-double FDensityField::GetNodeDensity(const int32 Index) const
+
+//	GETTER AND SETTER FUNCTIONS
+
+FVector FDensityField::GetGridOrigin() const
 {
-    if (GridNodes.IsValidIndex(Index))
+	return GridOrigin;
+}
+
+double FDensityField::GetVoxelDensityFromIndex(const int32 Index) const
+{
+    if (VoxelList.IsValidIndex(Index))
     {
-        return GridNodes[Index].GetNodeDensity();
+        return VoxelList[Index].GetVoxelDensity();
     }
     return 0.0;
 }
 
-int32 FDensityField::GetNodeNum() const
+int32 FDensityField::GetVoxelNumber() const
 {
-    return GridNodes.Num();
+    return VoxelList.Num();
 }
 
-void FDensityField::SetNodeDensity(const int32 Index, const double Density)
+void FDensityField::SetVoxelDensity(const int32 Index, const double Density)
 {
-    if (GridNodes.IsValidIndex(Index))
+    if (VoxelList.IsValidIndex(Index))
     {
-        GridNodes[Index].SetNodeDensity(Density);
+        VoxelList[Index].SetVoxelDensity(Density);
     }
 }
 
-void FDensityField::SetNodeGradient(const int32 Index, const FVector& Gradient)
+void FDensityField::SetVoxelGradient(const int32 Index, const FVector& Gradient)
 {
-    if (GridNodes.IsValidIndex(Index))
+    if (VoxelList.IsValidIndex(Index))
     {
-        GridNodes[Index].SetNodeGradient(Gradient);
+        VoxelList[Index].SetVoxelGradient(Gradient);
     }
 }
 
-double FDensityField::GetDensity(const int32 XBin, const int32 YBin, const int32 ZBin) const
+double FDensityField::GetVoxelDensityFromGridPosition(const int32 XBin, const int32 YBin, const int32 ZBin) const
 {
 	// First, check if the provided indices are valid
-	if (IsValidIndex(XBin, YBin, ZBin))
+	if (IsValidGridPosition(XBin, YBin, ZBin))
 	{
 		// Convert 3D grid coordinates to a linear index
-		const int32 Index = GridPositionToIndex(ZBin, YBin, XBin);
+		const int32 Index = GridPositionToIndex(XBin, YBin, ZBin);
 		
-		// Ensure the index is within the valid range
-		if (GridNodes.IsValidIndex(Index))
-		{
-			// Return the density at the calculated index
-			return GridNodes[Index].GetNodeDensity();
-		}
+		return GetVoxelDensityFromIndex(Index);
 	}
+	
 	// Return 0 if the indices are out of bounds or the index is invalid
 	return 0.0;
 }
 
-FVector FDensityField::GetGradient(const int32 XBin, const int32 YBin, const int32 ZBin) const
+FVector FDensityField::GetVoxelGradientFromGridPosition(const int32 XBin, const int32 YBin, const int32 ZBin) const
 {
 	// First, check if the provided indices are valid
-	if (IsValidIndex(XBin, YBin, ZBin))
+	if (IsValidGridPosition(XBin, YBin, ZBin))
 	{
 		// Convert 3D grid coordinates to a linear index
-		const int32 Index = GridPositionToIndex(ZBin, YBin, XBin);
+		const int32 Index = GridPositionToIndex(XBin, YBin, ZBin);
 		
 		// Ensure the index is within the valid range
-		if (GridNodes.IsValidIndex(Index))
+		if (VoxelList.IsValidIndex(Index))
 		{
 			// Return the gradient at the calculated index
-			return GridNodes[Index].GetNodeGradient();
+			return VoxelList[Index].GetVoxelGradient();
 		}
 	}
 	// Return 0 if the indices are out of bounds or the index is invalid
@@ -362,13 +350,45 @@ int32 FDensityField::GetZNum() const {
 	return ZNum;
 }
 
-bool FDensityField::IsValidPosition(const FVector& Position) const
+bool FDensityField::IsValidWorldPosition(const FVector& Position) const
 {
 	// Convert the world position to grid coordinates
 	const FVector GridPosition = WorldToGridPosition(Position);
 
 	// Check if the grid position is within the valid bounds
-	return (GridPosition.X >= 0 && GridPosition.X < XNum) &&
-		   (GridPosition.Y >= 0 && GridPosition.Y < YNum) &&
-		   (GridPosition.Z >= 0 && GridPosition.Z < ZNum);
+	return IsValidGridPosition(GridPosition.X, GridPosition.Y, GridPosition.Z);
+}
+
+TArray<int32> FDensityField::GetNeighborsFromIndex(const int32 Index) const
+{
+	TArray<int32> Neighbors;
+	if (VoxelList.IsValidIndex(Index))
+	{
+		const int32 X = VoxelList[Index].GetVoxelGridPos().X;
+		const int32 Y = VoxelList[Index].GetVoxelGridPos().Y;
+		const int32 Z = VoxelList[Index].GetVoxelGridPos().Z;
+
+		// List all potential neighbors
+		TArray<FIntVector> PotentialNeighbors = {
+			FIntVector(X+1, Y, Z), FIntVector(X-1, Y, Z),
+			FIntVector(X, Y+1, Z), FIntVector(X, Y-1, Z),
+			FIntVector(X, Y, Z+1), FIntVector(X, Y, Z-1)
+		};
+
+		for (const FIntVector& NeighborPos : PotentialNeighbors)
+		{
+			if (NeighborPos.X >= 0 && NeighborPos.X < XNum &&
+				NeighborPos.Y >= 0 && NeighborPos.Y < YNum &&
+				NeighborPos.Z >= 0 && NeighborPos.Z < ZNum)
+			{
+				int NeighborIndex = GridPositionToIndex(NeighborPos.Z, NeighborPos.Y, NeighborPos.Z);
+				Neighbors.Add(NeighborIndex);
+			}
+		}
+	}
+	return Neighbors;
+}
+
+bool FDensityField::IsValidIndex(const int32 Index) const {
+	return VoxelList.IsValidIndex(Index);
 }

Source/MetaCastBachelor/DensityField.h

@@ -15,22 +15,22 @@ private:
 	FVector GridPosition;
     double VoxelDensity;
     FVector VoxelGradient;
-	double EnclosedParticleNumber;
+	double ClosePointsNumber;
 
 public:
-    FVoxel() : WorldPosition(FVector::ZeroVector), GridPosition(FVector::ZeroVector), VoxelDensity(0.0), VoxelGradient(FVector::ZeroVector), EnclosedParticleNumber(0.0) {}
-    FVoxel(const FVector &InPosition, const FVector &InGridPos) : WorldPosition(InPosition), GridPosition(InGridPos), VoxelDensity(0.0), VoxelGradient(FVector::ZeroVector), EnclosedParticleNumber(0.0) {}
-	
-    double GetNodeDensity() const { return VoxelDensity; }
-	FVector GetNodeGridPos() const { return GridPosition; }
-    FVector GetNodePosition() const { return WorldPosition; }
-    FVector GetNodeGradient() const { return VoxelGradient; }
-	double GetEnclosedParticleNum() const { return EnclosedParticleNumber; }
-	
-	void SetNodeDensity(const double Density) { this->VoxelDensity = Density; }
-    void SetNodeGradient(const FVector &Gradient) { this->VoxelGradient = Gradient; }
-	void SetEnclosedParticleNum(const double Dis) { EnclosedParticleNumber = Dis; }
-    void NodeDensityPlusDis(const double Dis) { VoxelDensity += Dis; }
+    FVoxel() : WorldPosition(FVector::ZeroVector), GridPosition(FVector::ZeroVector), VoxelDensity(0.0), VoxelGradient(FVector::ZeroVector), ClosePointsNumber(0.0) {}
+    FVoxel(const FVector &InPosition, const FVector &InGridPos) : WorldPosition(InPosition), GridPosition(InGridPos), VoxelDensity(0.0), VoxelGradient(FVector::ZeroVector), ClosePointsNumber(0.0) {}
+	
+    double GetVoxelDensity() const { return VoxelDensity; }
+	FVector GetVoxelGridPos() const { return GridPosition; }
+    FVector GetVoxelPosition() const { return WorldPosition; }
+    FVector GetVoxelGradient() const { return VoxelGradient; }
+	double GetClosePointsNumber() const { return ClosePointsNumber; }
+	
+	void SetVoxelDensity(const double Density) { this->VoxelDensity = Density; }
+    void SetVoxelGradient(const FVector &Gradient) { this->VoxelGradient = Gradient; }
+	void SetClosePointsNumber(const double Dis) { ClosePointsNumber = Dis; }
+    void AddToVoxelDensity(const double Dis) { VoxelDensity += Dis; }
 };
 
 USTRUCT(BlueprintType)
@@ -54,7 +54,7 @@ struct FLUTUnit
 class FDensityField
 {
 	//	VARIABLES
-    TArray<FVoxel> GridNodes;
+    TArray<FVoxel> VoxelList;
     int32 XNum, YNum, ZNum;
     float XStep, YStep, ZStep;
 	FVector GridOrigin;
@@ -72,27 +72,27 @@ public:
     void DrawDebugVoxelDensity(const UWorld* World, const AActor* Actor, float Duration, float Scale, float DensityThreshold) const;
 
 	// CONVERSION FUNCTIONS
-	int32 WorldPositionToIndex(const FVector &Position, const APointCloud * PointCloud) const;
+	int32 WorldPositionToIndex(const FVector& Position) const;
 	FVector WorldToGridPosition(const FVector& Position) const;
-	bool IsValidIndex(int32 XBin, int32 YBin, int32 ZBin) const;
-	int32 GridPositionToIndex(int32 Z, int32 Y, int32 X) const;
+	bool IsValidGridPosition(int32 XBin, int32 YBin, int32 ZBin) const;
+	int32 GridPositionToIndex(int32 X, int32 Y, int32 Z) const;
 	
 	// GETTER AND SETTER FUNCTIONS
 	FVector GetGridOrigin() const;
-    FVector GetNodeGradient(int32 Index) const;
-    FVector GetNodePosition(int32 Index) const;
-    FVector GetNodeGridPos(int32 Index) const;
-	int32 GetNodeNum() const;
-    double GetNodeDensity(int32 Index) const;
-	double GetDensity(int32 XBin, int32 YBin, int32 ZBin) const;
-	FVector GetGradient(int32 XBin, int32 YBin, int32 ZBin) const;
+    FVector IndexToVoxelGradient(int32 Index) const;
+    FVector IndexToVoxelPosition(int32 Index) const;
+    FVector IndexToGridPosition(int32 Index) const;
+	int32 GetVoxelNumber() const;
+    double GetVoxelDensityFromIndex(int32 Index) const;
+	double GetVoxelDensityFromGridPosition(int32 XBin, int32 YBin, int32 ZBin) const;
+	FVector GetVoxelGradientFromGridPosition(int32 XBin, int32 YBin, int32 ZBin) const;
 	FVector GetStep() const;
 	int32 GetXNum() const;
 	int32 GetYNum() const;
 	int32 GetZNum() const;
-	bool IsValidPosition(const FVector& Position) const;
-
-
-	void SetNodeDensity(int32 Index, double Density);
-    void SetNodeGradient(int32 Index, const FVector& Gradient);
+	bool IsValidWorldPosition(const FVector& Position) const;
+	TArray<int32> GetNeighborsFromIndex(const int32 Index) const;
+	bool IsValidIndex(int32 Index) const;
+	void SetVoxelDensity(int32 Index, double Density);
+    void SetVoxelGradient(int32 Index, const FVector& Gradient);
 };

Source/MetaCastBachelor/MetaCastBaseline.cpp

@@ -130,14 +130,12 @@ void UMetaCastBaseline::AttachToHand(const float DeltaTime) const
 	if (MyPointCloud && LeftHandComponent)
 	{
 		// Get the target position and rotation from the hand component, applying the offset
-		FVector TargetPosition = LeftHandComponent->GetComponentLocation() + PositionOffset;
-		FRotator TargetRotation = LeftHandComponent->GetComponentRotation() + RotationOffset;
+		const FVector TargetPosition = LeftHandComponent->GetComponentLocation() + PositionOffset;
+		const FRotator TargetRotation = LeftHandComponent->GetComponentRotation() + RotationOffset;
 
 		// Smoothly lerp the point cloud's position and rotation to the target values
-		FVector NewPosition =
-			FMath::VInterpTo(MyPointCloud->GetActorLocation(), TargetPosition, DeltaTime, PositionLerpSpeed);
-		FRotator NewRotation =
-			FMath::RInterpTo(MyPointCloud->GetActorRotation(), TargetRotation, DeltaTime, RotationLerpSpeed);
+		const FVector NewPosition = FMath::VInterpTo(MyPointCloud->GetActorLocation(), TargetPosition, DeltaTime, PositionLerpSpeed);
+		const FRotator NewRotation = FMath::RInterpTo(MyPointCloud->GetActorRotation(), TargetRotation, DeltaTime, RotationLerpSpeed);
 
 		// Apply the new position and rotation to the point cloud
 		MyPointCloud->SetActorLocation(NewPosition);
@@ -146,8 +144,7 @@ void UMetaCastBaseline::AttachToHand(const float DeltaTime) const
 }
 
 // Called every frame
-void UMetaCastBaseline::TickComponent(const float DeltaTime, const ELevelTick TickType,
-                                      FActorComponentTickFunction* ThisTickFunction)
+void UMetaCastBaseline::TickComponent(const float DeltaTime, const ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction)
 {
 	Super::TickComponent(DeltaTime, TickType, ThisTickFunction);
 

Source/MetaCastBachelor/MetaPoint.cpp

 #include "MetaPoint.h"
-#include "MaterialDomain.h"
 #include "Utilities.h"
 
-UMetaPoint::UMetaPoint() : MyProceduralMesh(nullptr) {
+UMetaPoint::UMetaPoint() : Index(0), MyDensityField(nullptr), Threshold(0), MyProceduralMesh(nullptr)
+{
 	// Initialize the Marching Cubes algorithm
 	MyMarchingCubes = new MarchingCubes();
 }
@@ -12,10 +12,7 @@ void UMetaPoint::BeginPlay()
 	Super::BeginPlay();
 }
 
-void UMetaPoint::SelectParticles(const FVector& InputPosition)
-{
 
-}
 
 void UMetaPoint::EraseParticles(const FVector& InputPosition)
 {
@@ -36,7 +33,7 @@ void UMetaPoint::HandleMetaSelectPressed(const FInputActionInstance& Instance)
 
 
 // Method to perform gradient ascent to find the local maximum density starting from a given position.
-void UMetaPoint::FindAndMarkLocalMaximum() const
+void UMetaPoint::FindAndMarkLocalMaximum()
 {
 	const UWorld* World = GetWorld();
 	
@@ -48,19 +45,46 @@ void UMetaPoint::FindAndMarkLocalMaximum() const
 	const FVector WorldPosition = SelectionObject->GetComponentLocation();
 	UE_LOG(LogTemp, Log, TEXT("World Position: %s"), *WorldPosition.ToString());
 
+	MyDensityField = MyPointCloud->MyDensityField;
+
 	// Convert the world position of the selection object to the local position relative to the point cloud
-	const FVector StartPosition = MyPointCloud->GetTransform().InverseTransformPosition(WorldPosition);
+	const FVector StartPosition = MyPointCloud->PointCloudVisualizer->GetComponentTransform().InverseTransformPosition(WorldPosition);
 
 	// Perform gradient ascent to find the local maximum starting from the converted local position
-	const FVector LocalMaximum = FUtilities::FollowGradientToMaximum(StartPosition, MyPointCloud->MyDensityField);
+	LocalMaximum = FUtilities::FollowGradientToMaximum(MyDensityField, StartPosition);
+	UE_LOG(LogTemp, Log, TEXT("Local maximum found at world position: %s"), *WorldMaximum.ToString());
 
 	// Convert the local maximum back to world coordinates
-	const FVector WorldMaximum = MyPointCloud->GetActorTransform().TransformPosition(LocalMaximum);
+	Index = MyDensityField->WorldPositionToIndex(LocalMaximum);
+
+	Threshold = FUtilities::InterpolateDensityAtPosition(MyDensityField, LocalMaximum);
+	TestingThresholdFactor = 1;
+}
+
+void UMetaPoint::TickComponent(float DeltaTime, ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction)
+{
+	Super::TickComponent(DeltaTime, TickType, ThisTickFunction);
+
+	TestingThresholdFactor -= DeltaTime * 0.1f;
+}
+
+void UMetaPoint::SelectParticles(const FVector& InputPosition)
+{
+	const UWorld* World = GetWorld();
+	WorldMaximum = MyPointCloud->PointCloudVisualizer->GetComponentTransform().TransformPosition(LocalMaximum);
 	
 	// Draw a debug sphere at the location of the local maximum in the world
-	DrawDebugSphere(World, WorldMaximum, SelectionRadius / 2, 32, FColor::Red, false, 10);
-	UE_LOG(LogTemp, Log, TEXT("Local maximum found at world position: %s"), *WorldMaximum.ToString());
+	DrawDebugSphere(World, WorldMaximum, SelectionRadius / 2, 32, FColor::Red, false, 0);
+	auto Voxels = FUtilities::FloodFilling(MyDensityField, Index, Threshold * TestingThresholdFactor);
 
+	for (const int32 Index2 : Voxels)
+	{
+		MyPointCloud->DrawVoxel(Index2, 0);
 	}
+}
+
+
+
+
 
 

Source/MetaCastBachelor/MetaPoint.h

@@ -13,13 +13,20 @@ class UMetaPoint : public UMetaCastBaseline
 	GENERATED_BODY()
 	
 	MarchingCubes* MyMarchingCubes;
-	
+	FVector WorldMaximum;
+	int32 Index;
+	FVector LocalMaximum;
+	FDensityField* MyDensityField;
+	float Threshold;
 	
 public:
 
 	UPROPERTY(EditAnywhere, BlueprintReadWrite)
 	UProceduralMeshComponent* MyProceduralMesh;
 
+	UPROPERTY(EditAnywhere)
+	float TestingThresholdFactor = 2.0;
+	
 	UMetaPoint();
 	
 	virtual void BeginPlay() override;
@@ -29,6 +36,8 @@ public:
 
 	virtual void HandleMetaSelectReleased(const FInputActionInstance& Instance) override;
 	virtual void HandleMetaSelectPressed(const FInputActionInstance& Instance) override;
-	void FindAndMarkLocalMaximum() const;
+	void FindAndMarkLocalMaximum();
+
+	virtual void TickComponent(float DeltaTime, ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction) override;
 };
 

Source/MetaCastBachelor/PointCloud.cpp

@@ -91,7 +91,7 @@ void APointCloud::SetupDensityFieldFromPointCloud() const
 	constexpr int32 ZAxisNum = 100; // Number of divisions in the grid along the Z-axis
 	MyDensityField->InitializeDensityField(MinBounds, MaxBounds, XAxisNum, YAxisNum, ZAxisNum);
 
-	constexpr float InfluenceRadius = 5.0f; // Half the size of a voxel
+	constexpr float InfluenceRadius = 0.5f; // Half the size of a voxel
 	constexpr float Sigma = 3.0f;  // Standard deviation for the Gaussian kernel
 	MyDensityField->CalculateVoxelDensities(this, InfluenceRadius, Sigma);
 
@@ -99,10 +99,9 @@ void APointCloud::SetupDensityFieldFromPointCloud() const
 
 	if (const UWorld* World = GetWorld())
 	{
-		MyDensityField->DrawDebugVoxelDensity(World, this, 100.0f, 1.0f, 0.5f);
+		//MyDensityField->DrawDebugVoxelDensity(World, this, 100.0f, 1.0f, 0.8f);
 	}
 }
-
 void APointCloud::UpdateSelection()
 {
 	for (int32 i = 0; i < PositionVectors.Num(); i++)
@@ -119,6 +118,25 @@ void APointCloud::UpdateSelection()
 	PointCloudVisualizer->SetInputAndConvert2(PositionVectors, PointColors);
 }
 
+void APointCloud::DrawVoxel(const int Index, float Time) const
+{
+	if(!MyDensityField->IsValidIndex(Index))
+	{
+		return;
+	}
+
+	const FVector VoxelHalfExtents = MyDensityField->GetStep() / 2.0f;
+	const FVector VoxelWorldPosition = PointCloudVisualizer->GetComponentTransform().TransformPosition(MyDensityField->IndexToVoxelPosition(Index));
+	const FVector WorldScale = PointCloudVisualizer->GetComponentScale();
+
+	// Scale the voxel size by the component's world scale to ensure it reflects the actual dimensions
+	FVector ScaledVoxelSize = GetActorRightVector();
+
+	// Draw the debug box with the correct world space dimensions
+	DrawDebugBox(GetWorld(), VoxelWorldPosition, VoxelHalfExtents, this->GetActorQuat(), FColor::Green, false, Time, 0, 0.1f);
+
+}
+
 // Called every frame
 void APointCloud::Tick(float DeltaTime)
 {

Source/MetaCastBachelor/PointCloud.h

@@ -13,7 +13,6 @@ class APointCloud : public AActor
 	
 public:
 	FDensityField* MyDensityField;
-	
 	TArray<FVector> PositionVectors;
 	TArray<bool> DefaultFlags;
 	TArray<bool> SelectionFlags;
@@ -59,6 +58,7 @@ public:
 	virtual void Tick(float DeltaTime) override;
 
 	void UpdateSelection();
+	void DrawVoxel(const int Index, float Time) const;
 
 	UFUNCTION(BlueprintCallable)
 	void ReadPointCloudFromFile(FFilePath FileNamePoints, FFilePath FileNameFlags);

Source/MetaCastBachelor/Utilities.cpp

@@ -26,128 +26,94 @@ void FUtilities::QuickSort(TArray<float>& Array, const int32 Begin, const int32
 	}
 }
 
-double FUtilities::InterpolateDensityAtPosition(const FVector& Position, const FDensityField* DensityField)
+double FUtilities::InterpolateDensityAtPosition(const FDensityField* DensityField, const FVector& Position)
 {
     if (!DensityField)
         return 0.0;
 
-    // Convert position into grid coordinate space
     const FVector GridPos = (Position - DensityField->GetGridOrigin()) / DensityField->GetStep();
 
     const int32 XBin = FMath::FloorToInt(GridPos.X);
     const int32 YBin = FMath::FloorToInt(GridPos.Y);
     const int32 ZBin = FMath::FloorToInt(GridPos.Z);
 
-    // Calculate fractional components for interpolation
-    const double XFraction = GridPos.X - XBin;
-    const double YFraction = GridPos.Y - YBin;
-    const double ZFraction = GridPos.Z - ZBin;
-
-    // Ensure indices are within the valid range
-    if (!DensityField->IsValidIndex(XBin, YBin, ZBin) ||
-        !DensityField->IsValidIndex(XBin + 1, YBin, ZBin) ||
-        !DensityField->IsValidIndex(XBin, YBin + 1, ZBin) ||
-        !DensityField->IsValidIndex(XBin, YBin, ZBin + 1) ||
-        !DensityField->IsValidIndex(XBin + 1, YBin + 1, ZBin) ||
-        !DensityField->IsValidIndex(XBin + 1, YBin, ZBin + 1) ||
-        !DensityField->IsValidIndex(XBin, YBin + 1, ZBin + 1) ||
-        !DensityField->IsValidIndex(XBin + 1, YBin + 1, ZBin + 1))
-    {
-        return 0.0;
-    }
+    const float XFraction = GridPos.X - XBin;
+    const float YFraction = GridPos.Y - YBin;
+    const float ZFraction = GridPos.Z - ZBin;
 
-    // Fetch the densities from the surrounding voxel corners
-    const double D000 = DensityField->GetDensity(XBin, YBin, ZBin);
-    const double D100 = DensityField->GetDensity(XBin + 1, YBin, ZBin);
-    const double D010 = DensityField->GetDensity(XBin, YBin + 1, ZBin);
-    const double D001 = DensityField->GetDensity(XBin, YBin, ZBin + 1);
-    const double D110 = DensityField->GetDensity(XBin + 1, YBin + 1, ZBin);
-    const double D101 = DensityField->GetDensity(XBin + 1, YBin, ZBin + 1);
-    const double D011 = DensityField->GetDensity(XBin, YBin + 1, ZBin + 1);
-    const double D111 = DensityField->GetDensity(XBin + 1, YBin + 1, ZBin + 1);
+    const double D000 = DensityField->GetVoxelDensityFromGridPosition(XBin, YBin, ZBin);
+    const double D100 = DensityField->GetVoxelDensityFromGridPosition(XBin + 1, YBin, ZBin);
+    const double D010 = DensityField->GetVoxelDensityFromGridPosition(XBin, YBin + 1, ZBin);
+    const double D001 = DensityField->GetVoxelDensityFromGridPosition(XBin, YBin, ZBin + 1);
+    const double D101 = DensityField->GetVoxelDensityFromGridPosition(XBin + 1, YBin, ZBin + 1);
+    const double D011 = DensityField->GetVoxelDensityFromGridPosition(XBin, YBin + 1, ZBin + 1);
+    const double D110 = DensityField->GetVoxelDensityFromGridPosition(XBin + 1, YBin + 1, ZBin);
+    const double D111 = DensityField->GetVoxelDensityFromGridPosition(XBin + 1, YBin + 1, ZBin + 1);
 
-    // Interpolate along x for each of the pairs
-    const double D00 = D000 * (1 - XFraction) + D100 * XFraction;
-    const double D01 = D001 * (1 - XFraction) + D101 * XFraction;
-    const double D10 = D010 * (1 - XFraction) + D110 * XFraction;
-    const double D11 = D011 * (1 - XFraction) + D111 * XFraction;
+    // Trilinear interpolation
+    const double DX00 = FMath::Lerp(D000, D100, XFraction);
+    const double DX01 = FMath::Lerp(D001, D101, XFraction);
+    const double DX10 = FMath::Lerp(D010, D110, XFraction);
+    const double DX11 = FMath::Lerp(D011, D111, XFraction);
 
-    // Interpolate along y
-    const double D0 = D00 * (1 - YFraction) + D10 * YFraction;
-    const double D1 = D01 * (1 - YFraction) + D11 * YFraction;
+    const double DY0 = FMath::Lerp(DX00, DX10, YFraction);
+    const double DY1 = FMath::Lerp(DX01, DX11, YFraction);
 
-    // Final interpolation along z
-    const double Density = D0 * (1 - ZFraction) + D1 * ZFraction;
+    const double DZ = FMath::Lerp(DY0, DY1, ZFraction);
 
-    return Density;
+    return DZ;
 }
 
-
-FVector FUtilities::InterpolateGradientAtPosition(const FVector& Position, const FDensityField* DensityField)
+FVector FUtilities::InterpolateGradientAtPosition(const FDensityField* DensityField, const FVector& Position)
 {
 	if (!DensityField)
 		return FVector::ZeroVector;
 
-    // Convert position into grid coordinate space
-    const FVector GridPos = (Position - DensityField->GetGridOrigin()) / DensityField->GetStep();
-
-    const int32 XBin = FMath::FloorToInt(GridPos.X);
-    const int32 YBin = FMath::FloorToInt(GridPos.Y);
-    const int32 ZBin = FMath::FloorToInt(GridPos.Z);
-
-    // Calculate fractional components for interpolation
-    const double XFraction = GridPos.X - XBin;
-    const double YFraction = GridPos.Y - YBin;
-    const double ZFraction = GridPos.Z - ZBin;
-
-    // Ensure indices are within the valid range
-    if (!DensityField->IsValidIndex(XBin, YBin, ZBin) ||
-        !DensityField->IsValidIndex(XBin + 1, YBin, ZBin) ||
-        !DensityField->IsValidIndex(XBin, YBin + 1, ZBin) ||
-        !DensityField->IsValidIndex(XBin, YBin, ZBin + 1) ||
-        !DensityField->IsValidIndex(XBin + 1, YBin + 1, ZBin) ||
-        !DensityField->IsValidIndex(XBin + 1, YBin, ZBin + 1) ||
-        !DensityField->IsValidIndex(XBin, YBin + 1, ZBin + 1) ||
-        !DensityField->IsValidIndex(XBin + 1, YBin + 1, ZBin + 1))
-    {
-        return FVector::ZeroVector;
-    }
-
-    // Fetch the gradients from the surrounding voxel corners
-    const FVector G000 = DensityField->GetGradient(XBin, YBin, ZBin);
-    const FVector G100 = DensityField->GetGradient(XBin + 1, YBin, ZBin);
-    const FVector G010 = DensityField->GetGradient(XBin, YBin + 1, ZBin);
-    const FVector G001 = DensityField->GetGradient(XBin, YBin, ZBin + 1);
-    const FVector G110 = DensityField->GetGradient(XBin + 1, YBin + 1, ZBin);
-    const FVector G101 = DensityField->GetGradient(XBin + 1, YBin, ZBin + 1);
-    const FVector G011 = DensityField->GetGradient(XBin, YBin + 1, ZBin + 1);
-    const FVector G111 = DensityField->GetGradient(XBin + 1, YBin + 1, ZBin + 1);
-
-    // Interpolate along x for each of the pairs
-    const FVector G00 = G000 * (1 - XFraction) + G100 * XFraction;
-    const FVector G01 = G001 * (1 - XFraction) + G101 * XFraction;
-    const FVector G10 = G010 * (1 - XFraction) + G110 * XFraction;
-    const FVector G11 = G011 * (1 - XFraction) + G111 * XFraction;
+	FVector GridPos = (Position - DensityField->GetGridOrigin()) / DensityField->GetStep();
+
+	int32 XBin = FMath::FloorToInt(GridPos.X);
+	int32 YBin = FMath::FloorToInt(GridPos.Y);
+	int32 ZBin = FMath::FloorToInt(GridPos.Z);
+
+	// Fractional parts to interpolate
+	float XFraction = GridPos.X - XBin;
+	float YFraction = GridPos.Y - YBin;
+	float ZFraction = GridPos.Z - ZBin;
+
+	// Fetch gradients from corners of the cell
+	FVector G000 = DensityField->GetVoxelGradientFromGridPosition(XBin, YBin, ZBin);
+	FVector G100 = DensityField->GetVoxelGradientFromGridPosition(XBin + 1, YBin, ZBin);
+	FVector G010 = DensityField->GetVoxelGradientFromGridPosition(XBin, YBin + 1, ZBin);
+	FVector G001 = DensityField->GetVoxelGradientFromGridPosition(XBin, YBin, ZBin + 1);
+	FVector G110 = DensityField->GetVoxelGradientFromGridPosition(XBin + 1, YBin + 1, ZBin);
+	FVector G101 = DensityField->GetVoxelGradientFromGridPosition(XBin + 1, YBin, ZBin + 1);
+	FVector G011 = DensityField->GetVoxelGradientFromGridPosition(XBin, YBin + 1, ZBin + 1);
+	FVector G111 = DensityField->GetVoxelGradientFromGridPosition(XBin + 1, YBin + 1, ZBin + 1);
+
+	// Interpolate along x for each of the yz pairs
+	FVector G00 = FMath::Lerp(G000, G100, XFraction);
+	FVector G01 = FMath::Lerp(G001, G101, XFraction);
+	FVector G10 = FMath::Lerp(G010, G110, XFraction);
+	FVector G11 = FMath::Lerp(G011, G111, XFraction);
 
 	// Interpolate along y
-    const FVector G0 = G00 * (1 - YFraction) + G10 * YFraction;
-    const FVector G1 = G01 * (1 - YFraction) + G11 * YFraction;
+	FVector G0 = FMath::Lerp(G00, G10, YFraction);
+	FVector G1 = FMath::Lerp(G01, G11, YFraction);
 
 	// Final interpolation along z
-    const FVector Gradient = G0 * (1 - ZFraction) + G1 * ZFraction;
+	FVector Gradient = FMath::Lerp(G0, G1, ZFraction);
 
 	return Gradient;
 }
 
-
 float FUtilities::GaussianKernel(const float Distance, const float Sigma) {
 	return FMath::Exp(-0.5 * FMath::Square(Distance / Sigma));
 }
 
-FVector FUtilities::FollowGradientToMaximum(const FVector& StartPosition, const FDensityField* DensityField)
+FVector FUtilities::FollowGradientToMaximum(const FDensityField* DensityField, const FVector& StartPosition)
 {
 	FVector CurrentPosition = StartPosition;
-	double CurrentDensity = InterpolateDensityAtPosition(CurrentPosition, DensityField);
+	double CurrentDensity = InterpolateDensityAtPosition(DensityField, CurrentPosition);
 
 	int Iterations = 0;
 	constexpr int MaxIterations = 100;  // Prevent infinite loops
@@ -155,22 +121,21 @@ FVector FUtilities::FollowGradientToMaximum(const FVector& StartPosition, const
 
 	while (Iterations < MaxIterations)
 	{
-		FVector Gradient = InterpolateGradientAtPosition(CurrentPosition, DensityField);
+		FVector Gradient = InterpolateGradientAtPosition(DensityField, CurrentPosition);
 		if (Gradient.IsNearlyZero()) {
 			UE_LOG(LogTemp, Log, TEXT("Gradient is zero at position %s"), *CurrentPosition.ToString());
 			break;  // Gradient is zero, likely at a local maximum
 		}
 
 		// Move towards the direction of the gradient
-		FVector NextPosition = CurrentPosition + StepSize * Gradient.GetSafeNormal();
-		if (!DensityField->IsValidPosition(NextPosition)) {
+		FVector NextPosition = CurrentPosition + (StepSize * Gradient.GetSafeNormal());
+		if (!DensityField->IsValidWorldPosition(NextPosition)) {
 			UE_LOG(LogTemp, Warning, TEXT("Next position out of bounds: %s"), *NextPosition.ToString());
 			break;  // Next position is out of valid bounds
 		}
 
-		const double NewDensity = InterpolateDensityAtPosition(NextPosition, DensityField);
-		UE_LOG(LogTemp, Log, TEXT("Iteration %d: CurrentPos = %s, Gradient = %s, NextPos = %s, NewDensity = %f, CurrentDensity = %f"),
-			Iterations, *CurrentPosition.ToString(), *Gradient.ToString(), *NextPosition.ToString(), NewDensity, CurrentDensity);
+		const double NewDensity = InterpolateDensityAtPosition(DensityField, NextPosition);
+		//UE_LOG(LogTemp, Log, TEXT("Iteration %d: CurrentPos = %s, Gradient = %s, NextPos = %s, NewDensity = %f, CurrentDensity = %f"), Iterations, *CurrentPosition.ToString(), *Gradient.ToString(), *NextPosition.ToString(), NewDensity, CurrentDensity);
 
 		if (NewDensity <= CurrentDensity) {  // Check if density has increased
 			UE_LOG(LogTemp, Log, TEXT("Density did not increase; stopping at position %s with density %f"), *CurrentPosition.ToString(), CurrentDensity);
@@ -185,3 +150,51 @@ FVector FUtilities::FollowGradientToMaximum(const FVector& StartPosition, const
 	return CurrentPosition;
 }
 
+TArray<int32> FUtilities::FloodFilling(const FDensityField* DensityField, const int32 StartIndex, const float Threshold)
+{
+	UE_LOG(LogTemp, Warning, TEXT("T: %f"), Threshold);
+	
+	TArray<int32> VisitedIndices;
+	if (!DensityField) {
+		UE_LOG(LogTemp, Warning, TEXT("DensityField is null."));
+		return VisitedIndices;
+	}
+
+	if (!DensityField->IsValidIndex(StartIndex)) {
+		UE_LOG(LogTemp, Warning, TEXT("Start index %d is invalid."), StartIndex);
+		return VisitedIndices;
+	}
+
+	TArray<int32> Stack;
+	Stack.Push(StartIndex);
+	TArray<bool> Visited;
+	Visited.Init(false, DensityField->GetVoxelNumber());
+
+	while (Stack.Num() > 0)
+	{
+		int32 CurrentIndex = Stack.Pop();
+		if (!Visited[CurrentIndex])
+		{
+			Visited[CurrentIndex] = true;
+			VisitedIndices.Add(CurrentIndex);
+
+			TArray<int32> Neighbors = DensityField->GetNeighborsFromIndex(CurrentIndex);
+			UE_LOG(LogTemp, Log, TEXT("Visiting Node %d at Position %s with %d neighbors."), CurrentIndex, *DensityField->IndexToGridPosition(CurrentIndex).ToString(), Neighbors.Num());
+            
+			for (int32 NeighborIndex : Neighbors)
+			{
+				if (DensityField->GetVoxelDensityFromIndex(NeighborIndex) > Threshold && !Visited[NeighborIndex]) {
+					Stack.Push(NeighborIndex);
+				}
+			}
+		}
+	}
+
+	// Logging final details
+	UE_LOG(LogTemp, Log, TEXT("Flood filling completed from start index %d. Total voxels selected: %d"), StartIndex, VisitedIndices.Num());
+	for (const int32 idx : VisitedIndices) {
+		UE_LOG(LogTemp, Log, TEXT("Selected Voxel Index: %d at Grid Position: %s"), idx, *DensityField->IndexToGridPosition(idx).ToString());
+	}
+
+	return VisitedIndices;
+}

Source/MetaCastBachelor/Utilities.h

@@ -6,8 +6,9 @@ class FUtilities
 {
 public:
 	static void QuickSort(TArray<float>& Array, int32 Begin, int32 End);
-	static double InterpolateDensityAtPosition(const FVector& Position, const FDensityField* DensityField);
-	static FVector InterpolateGradientAtPosition(const FVector& Position, const FDensityField* DensityField);
+	static double InterpolateDensityAtPosition(const FDensityField* DensityField, const FVector& Position);
+	static FVector InterpolateGradientAtPosition(const FDensityField* DensityField, const FVector& Position);
 	static float GaussianKernel(float Distance, float Sigma);
-	static FVector FollowGradientToMaximum(const FVector& StartPosition, const FDensityField* DensityField);
+	static FVector FollowGradientToMaximum(const FDensityField* DensityField, const FVector& StartPosition);
+	static TArray<int32> FloodFilling(const FDensityField* DensityField, const int32 StartIndex, const float Threshold);
 };