diff --git a/Content/MetaPointMap.umap b/Content/MetaPointMap.umap
index f44a6fa3bd8788f0813c8f4c4351e9c5bf22240f..7e4733c7f9d52777e794102217da9ff21ea628ba 100644
Binary files a/Content/MetaPointMap.umap and b/Content/MetaPointMap.umap differ
diff --git a/Content/SelectionVolume.uasset b/Content/SelectionVolume.uasset
index b4a7c48878aadac755046e2d70262c2c42fa0168..f5d2c569e5837f227551627058837d7703770703 100644
Binary files a/Content/SelectionVolume.uasset and b/Content/SelectionVolume.uasset differ
diff --git a/Source/MetaCastBachelor/DensityField.cpp b/Source/MetaCastBachelor/DensityField.cpp
index f95293968bea14ae7527061b0f7c92618d168060..ab49740d352689c4683859bce08b5445cde2753b 100644
--- a/Source/MetaCastBachelor/DensityField.cpp
+++ b/Source/MetaCastBachelor/DensityField.cpp
@@ -5,12 +5,26 @@
 
 // INITIALIZATION FUNCTIONS
 
-FDensityField::FDensityField(): XNum(0), YNum(0), ZNum(0), XStep(0), YStep(0), ZStep(0)
+FDensityField::FDensityField(): XNum(0), YNum(0), ZNum(0), XStep(0), YStep(0), ZStep(0), VoxelPointLookupTable(nullptr)
 {
 	VoxelList = TArray<FVoxel>();
 }
 
-void FDensityField::InitializeDensityField(const FVector& MinBounds, const FVector& MaxBounds, const int32 XAxisNum, const int32 YAxisNum, const int32 ZAxisNum)
+void FDensityField::InitializeDensityField(const APointCloud* PointCloud, const FVector& MinBounds, const FVector& MaxBounds, const FIntVector AxisNumbers)
+{
+	InitializeDataStructures(MinBounds, MaxBounds, AxisNumbers.X, AxisNumbers.Y, AxisNumbers.Z);
+
+	constexpr float InfluenceRadius = 1.0f; // Size of a voxel
+	constexpr float Sigma = 1.0f;  // Standard deviation for the Gaussian kernel
+	
+	CalculateVoxelDensities(PointCloud, InfluenceRadius, Sigma);
+	CalculateAndStoreGradients();
+
+	VoxelPointLookupTable = new FVoxelPointLookupTable();
+	VoxelPointLookupTable->Initialize(PointCloud->PositionVectors, this);
+}
+
+void FDensityField::InitializeDataStructures(const FVector& MinBounds, const FVector& MaxBounds, const int32 XAxisNum, const int32 YAxisNum, const int32 ZAxisNum)
 {
 	VoxelList.Empty();
 
@@ -57,12 +71,12 @@ void FDensityField::CalculateVoxelDensities(const APointCloud* PointCloud, const
 	// Iterate over each particle
 	for (const FVector& ParticlePosition : PointCloud->PositionVectors)
 	{
-		const int32 StartX = FMath::Max(0, static_cast<int32>((ParticlePosition.X - InfluenceRadius - PointCloud->MinBounds.X) / XStep));
-	    const int32 EndX = FMath::Min(XNum - 1, static_cast<int32>((ParticlePosition.X + InfluenceRadius - PointCloud->MinBounds.X) / XStep));
-	    const int32 StartY = FMath::Max(0, static_cast<int32>((ParticlePosition.Y - InfluenceRadius - PointCloud->MinBounds.Y) / YStep));
-	    const int32 EndY = FMath::Min(YNum - 1, static_cast<int32>((ParticlePosition.Y + InfluenceRadius - PointCloud->MinBounds.Y) / YStep));
-	    const int32 StartZ = FMath::Max(0, static_cast<int32>((ParticlePosition.Z - InfluenceRadius - PointCloud->MinBounds.Z) / ZStep));
-	    const int32 EndZ = FMath::Min(ZNum - 1, static_cast<int32>((ParticlePosition.Z + InfluenceRadius - PointCloud->MinBounds.Z) / ZStep));
+		const int32 StartX = FMath::Max(0, FMath::FloorToInt((ParticlePosition.X - InfluenceRadius - PointCloud->MinBounds.X) / XStep));
+	    const int32 EndX = FMath::Min(XNum - 1, FMath::FloorToInt((ParticlePosition.X + InfluenceRadius - PointCloud->MinBounds.X) / XStep));
+	    const int32 StartY = FMath::Max(0, FMath::FloorToInt((ParticlePosition.Y - InfluenceRadius - PointCloud->MinBounds.Y) / YStep));
+	    const int32 EndY = FMath::Min(YNum - 1, FMath::FloorToInt((ParticlePosition.Y + InfluenceRadius - PointCloud->MinBounds.Y) / YStep));
+	    const int32 StartZ = FMath::Max(0, FMath::FloorToInt((ParticlePosition.Z - InfluenceRadius - PointCloud->MinBounds.Z) / ZStep));
+	    const int32 EndZ = FMath::Min(ZNum - 1, FMath::FloorToInt((ParticlePosition.Z + InfluenceRadius - PointCloud->MinBounds.Z) / ZStep));
 
 	    for (int32 Z = StartZ; Z <= EndZ; ++Z)
 	    {
@@ -76,9 +90,11 @@ void FDensityField::CalculateVoxelDensities(const APointCloud* PointCloud, const
 	                    FVector NodePosition = VoxelList[Index].GetVoxelPosition();
 	                    const float Distance = FVector::Dist(NodePosition, ParticlePosition);
 
-	                    if (Distance < InfluenceRadius)
+	                    if (Distance <= InfluenceRadius)
 	                    {
-	                        const double Weight = FUtilities::GaussianKernel(Distance, Sigma);
+	                        //const double Weight = FUtilities::GaussianKernel(Distance, Sigma);
+	                        const double Weight = FUtilities::SmoothingKernel(Distance, InfluenceRadius);
+	                        //const double Weight = InfluenceRadius - Distance;
 	                        VoxelList[Index].AddToVoxelDensity(Weight);
 	                        VoxelList[Index].SetClosePointsNumber(VoxelList[Index].GetClosePointsNumber() + 1.0);
 	                    }
@@ -249,19 +265,19 @@ void FDensityField::DrawDebugVoxelDensity(const UWorld* World, const AActor* Act
 }
 
 
-
 //	WORLD POSITION CONVERSION FUNCTIONS
 
 int32 FDensityField::WorldPositionToIndex(const FVector &Position) const {
 	// Convert world position to grid position
 	const FIntVector3 GridPosition = WorldToGridPosition(Position);
 
-	// Check if the grid position is valid
-	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
+	// Check if the calculated grid coordinates are within the bounds
+	if (!IsValidGridPosition(GridPosition.X, GridPosition.Y, GridPosition.Z))
+	{
+		//UE_LOG(LogTemp, Warning, TEXT("Position out of grid bounds, returning invalid index."));
+		return -1; 
 	}
-
+	
 	return GridPositionToIndex(GridPosition.X, GridPosition.Y, GridPosition.Z);
 }
 
@@ -277,14 +293,7 @@ FIntVector3 FDensityField::WorldToGridPosition(const FVector& Position) const
 	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)
-	{
-		UE_LOG(LogTemp, Warning, TEXT("Position out of grid bounds, returning invalid grid coordinates."));
-		return FIntVector3(-1, -1, -1);  // Return an invalid grid position if out of bounds
-	}
-
+	
 	return FIntVector3(GridX, GridY, GridZ);
 }
 
diff --git a/Source/MetaCastBachelor/DensityField.h b/Source/MetaCastBachelor/DensityField.h
index 3f4ff1008b2ee812e351e5e8972346fda050215f..c854dcdecc9bbb3e0fbf913733f6f46cbbdee95e 100644
--- a/Source/MetaCastBachelor/DensityField.h
+++ b/Source/MetaCastBachelor/DensityField.h
@@ -1,6 +1,7 @@
 #pragma once
 
 #include "CoreMinimal.h"
+#include "VoxelPointLookupTable.h"
 #include "DensityField.generated.h"
 
 class APointCloud;
@@ -33,24 +34,6 @@ public:
     void AddToVoxelDensity(const double Dis) { VoxelDensity += Dis; }
 };
 
-USTRUCT(BlueprintType)
-struct FLUTUnit
-{
-    GENERATED_BODY()
-
-    TArray<int32> LUTUnit;
-
-    void AddToLUT(const int32 TargetInt)
-    {
-        LUTUnit.Add(TargetInt);
-    }
-
-    TArray<int32> GetLTUnit() const
-    {
-        return LUTUnit;
-    }
-};
-
 class FDensityField
 {
 	//	VARIABLES
@@ -58,17 +41,20 @@ class FDensityField
     int32 XNum, YNum, ZNum;
     float XStep, YStep, ZStep;
 	FVector GridOrigin;
+
+	void CalculateVoxelDensities(const APointCloud* PointCloud, float InfluenceRadius, float Sigma);
+	void CalculateVoxelDensitiesByNumber(const APointCloud* PointCloud, float InfluenceRadius);
+	void CalculateAndStoreGradients();
+	void InitializeDataStructures(const FVector& MinBounds, const FVector& MaxBounds, int32 XAxisNum, int32 YAxisNum, int32 ZAxisNum);
 	
 public:
+	FVoxelPointLookupTable* VoxelPointLookupTable;
 	mutable FCriticalSection DataGuard;
 	// CONSTRUCTOR
 	FDensityField();
 
 	// INITIALIZATION FUNCTIONS
-	void InitializeDensityField(const FVector& MinBounds, const FVector& MaxBounds, const int32 XAxisNum, const int32 YAxisNum, const int32 ZAxisNum);
-    void CalculateVoxelDensities(const APointCloud* PointCloud, float InfluenceRadius, float Sigma);
-	void CalculateVoxelDensitiesByNumber(const APointCloud* PointCloud, float InfluenceRadius);
-	void CalculateAndStoreGradients();
+	void InitializeDensityField(const APointCloud* PointCloud, const FVector& MinBounds, const FVector& MaxBounds, const FIntVector AxisNumbers);
 
 	//	DEBUG FUNCTIONS
     void DrawDebugVoxelDensity(const UWorld* World, const AActor* Actor, float Duration, float Scale, float DensityThreshold) const;
@@ -97,7 +83,6 @@ public:
 	bool IsValidWorldPosition(const FVector& Position) const;
 	TArray<int32> IndexToVoxelNeighbors(const int32 Index) const;
 	bool IsValidIndex(int32 Index) const;
-	bool IsValid();
 	void SetVoxelDensityByIndex(int32 Index, double Density);
     void SetVoxelGradientByIndex(int32 Index, const FVector& Gradient);
 };
diff --git a/Source/MetaCastBachelor/MetaPoint.cpp b/Source/MetaCastBachelor/MetaPoint.cpp
index 4a71e325d0e0f39d79a93d784a5eba1efab578e7..cde81c999f8b30ab9f8e5808b402b6e4983f317c 100644
--- a/Source/MetaCastBachelor/MetaPoint.cpp
+++ b/Source/MetaCastBachelor/MetaPoint.cpp
@@ -1,10 +1,11 @@
 #include "MetaPoint.h"
 #include "Utilities.h"
+#include "Generators/MarchingCubes.h"
 
 UMetaPoint::UMetaPoint() : LocalMaximumIndex(0), MyDensityField(nullptr), MetaPointThreshold(0), World(nullptr)
 {
 	// Initialize the Marching Cubes algorithm
-	MyMarchingCubes = new MarchingCubes();
+	//MyMarchingCubes = new FMarchingCubes();
 
 	// Create the procedural mesh component
 	ProceduralMesh = CreateDefaultSubobject<UProceduralMeshComponent>(TEXT("GeneratedMesh"));
@@ -24,6 +25,13 @@ void UMetaPoint::BeginPlay()
 	ProceduralMesh->SetMaterial(0, SelectionVolumeMat);
 	ProceduralMesh->SetMaterial(1, SelectionVolumeMat);
 	ProceduralMesh->SetMaterial(2, SelectionVolumeMat);
+	ProceduralMesh->bCastDynamicShadow = false;
+	ProceduralMesh->bRenderCustomDepth = true;
+	ProceduralMesh->SetCastShadow(false);
+	ProceduralMesh->SetCollisionEnabled(ECollisionEnabled::NoCollision);
+	ProceduralMesh->SetMobility(EComponentMobility::Movable);
+	ProceduralMesh->SetVisibility(true);
+
 }
 
 void UMetaPoint::TickComponent(float DeltaTime, ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction)
@@ -31,6 +39,8 @@ void UMetaPoint::TickComponent(float DeltaTime, ELevelTick TickType, FActorCompo
 	Super::TickComponent(DeltaTime, TickType, ThisTickFunction);
 
 	AccumulatedTime += DeltaTime;
+
+	ProceduralMesh->SetVisibility(Select);
 }
 
 void UMetaPoint::EraseParticles(const FVector& InputPosition)
@@ -42,8 +52,8 @@ void UMetaPoint::HandleMetaSelectReleased(const FInputActionInstance& Instance)
 {
 	Super::HandleMetaSelectReleased(Instance);
 
-	ProceduralMesh->ClearAllMeshSections();
-	ProceduralMesh->SetVisibility(false);
+	//ProceduralMesh->ClearAllMeshSections();
+	//ProceduralMesh->SetVisibility(false);
 
 	AsyncTask(ENamedThreads::Type::AnyBackgroundHiPriTask, [this]()
 	{
@@ -72,6 +82,8 @@ void UMetaPoint::InitMetaPointSelection()
 	// Convert the local maximum back to world coordinates
 	LocalMaximumIndex = MyDensityField->WorldPositionToIndex(LocalMaximum);
 
+	SelectStartPosition = LocalMaximum;
+
 	MetaPointThreshold = FUtilities::InterpolateDensityAtPosition(MyDensityField, SelectionLocalPosition);
 
 	// Initialize Visited array to false for all voxels
@@ -87,6 +99,8 @@ void UMetaPoint::InitMetaPointSelection()
 
 	ProceduralMesh->ClearAllMeshSections();
 	World = GetWorld();
+
+	LastHandInput = SelectionWorldPosition;
 }
 
 void UMetaPoint::SelectParticles(const FVector& InputPosition)
@@ -98,6 +112,15 @@ void UMetaPoint::SelectParticles(const FVector& InputPosition)
 		AsyncTask(ENamedThreads::AnyBackgroundThreadNormalTask, [this]()
 		{
 			const FVector SelectionWorldPosition = SelectionObject->GetComponentLocation();
+			
+			if(FVector::Distance(LastHandInput, SelectionWorldPosition) < MinimalHandMovement)
+			{
+				AccumulatedTime = 0.0f;
+				return;
+			}
+			
+			LastHandInput = SelectionWorldPosition;
+						
 			// Convert the world position of the selection object to the local position relative to the point cloud
 			const FVector SelectionLocalPosition = MyPointCloud->PointCloudVisualizer->GetComponentTransform().InverseTransformPosition(SelectionWorldPosition);
 
@@ -105,19 +128,19 @@ void UMetaPoint::SelectParticles(const FVector& InputPosition)
 			LocalMaximum = FUtilities::FollowGradientToMaximum(MyDensityField, SelectionLocalPosition);
 
 			// Convert the local maximum back to world coordinates
-			LocalMaximumIndex = MyDensityField->WorldPositionToIndex(LocalMaximum);
+			//LocalMaximumIndex = MyDensityField->WorldPositionToIndex(LocalMaximum);
 
 			MetaPointThreshold = FUtilities::InterpolateDensityAtPosition(MyDensityField, SelectionLocalPosition);
-			FloodedIndices = FUtilities::FloodFilling(MyDensityField, LocalMaximumIndex, MetaPointThreshold);
-			
-			GenerateVoxelMesh(FloodedIndices);
+
+			const float MaxDistance = FVector::Distance(SelectStartPosition, SelectionLocalPosition);
+			FloodedIndices = FUtilities::FloodFilling_2(MyDensityField, LocalMaximumIndex, MetaPointThreshold, MaxDistance);
 			
-			AccumulatedTime = 0.0f;
+			GenerateVoxelMeshSmooth(FloodedIndices);
 		});
 	}
 }
 
-void UMetaPoint::GenerateVoxelMesh(const TArray<int32> Voxels) const
+void UMetaPoint::GenerateVoxelMeshWithCubes(const TArray<int32> Voxels) const
 {
 	TArray<FVector> Vertices;
 	TArray<int32> Triangles;
@@ -164,35 +187,79 @@ void UMetaPoint::GenerateVoxelMesh(const TArray<int32> Voxels) const
 	});
 }
 
-void UMetaPoint::CreateMeshSections(const TArray<FVector>& Vertices, const TArray<int32>& Triangles, const TArray<FColor>& VertexColors, const int32 MaxVerticesPerSection) const
+void UMetaPoint::GenerateVoxelMeshSmooth(TArray<int32> Voxels)
 {
-	const int32 NumVertices = Vertices.Num();
-	const int32 NumTriangles = Triangles.Num() / 3;
-
-	for (int32 SectionIndex = 0, VertexIndex = 0, TriangleIndex = 0; VertexIndex < NumVertices && TriangleIndex < NumTriangles;
-		 ++SectionIndex, VertexIndex += MaxVerticesPerSection, TriangleIndex += MaxVerticesPerSection / 3)
+	AsyncTask(ENamedThreads::AnyBackgroundThreadNormalTask, [this, Voxels]()
 	{
-		const int32 NumVerticesInSection = FMath::Min(MaxVerticesPerSection, NumVertices - VertexIndex);
-		const int32 NumTrianglesInSection = FMath::Min(MaxVerticesPerSection / 3, NumTriangles - TriangleIndex);
-
-		TArray<FVector> SectionVertices;
-		TArray<FColor> SectionVertexColors;
-		SectionVertices.Reserve(NumVerticesInSection);
-		SectionVertexColors.Reserve(NumVerticesInSection);
-		for (int32 i = 0; i < NumVerticesInSection; ++i)
+		FVector Min(FLT_MAX, FLT_MAX, FLT_MAX);
+		FVector Max(-FLT_MAX, -FLT_MAX, -FLT_MAX);
+		
+		for(const int FloodedIndex : Voxels)
 		{
-			SectionVertices.Add(Vertices[VertexIndex + i]);
-			SectionVertexColors.Add(VertexColors[VertexIndex + i]);
+			Min = FVector::Min(Min, MyDensityField->IndexToVoxelPosition(FloodedIndex));
+			Max = FVector::Max(Max, MyDensityField->IndexToVoxelPosition(FloodedIndex));
 		}
 
-		TArray<int32> SectionTriangles;
-		SectionTriangles.Reserve(NumTrianglesInSection * 3);
-		for (int32 i = 0; i < NumTrianglesInSection * 3; ++i)
+		UE::Geometry::FMarchingCubes MarchingCubes;
+		MarchingCubes.Bounds = UE::Geometry::TAxisAlignedBox3(Min, Max);
+		MarchingCubes.CubeSize = MyDensityField->GetStep().X;
+		MarchingCubes.IsoValue = 0;
+
+		// Define the implicit function
+		MarchingCubes.Implicit = [this, Voxels](const FVector3d& Pos) {
+			const FVector PosConverted = FVector(Pos.X, Pos.Y, Pos.Z);
+			const int Index = MyDensityField->WorldPositionToIndex(PosConverted);
+			return Voxels.Contains(Index) ? 1 : -1;
+		};
+
+		MarchingCubes.bParallelCompute = true;
+		const UE::Geometry::FMeshShapeGenerator& Generator = MarchingCubes.Generate();
+		TArray<FVector3d> Vertices3d = Generator.Vertices;
+		TArray<UE::Geometry::FIndex3i> Triangles = Generator.Triangles;
+		TArray<FVector3f> Normals3F = Generator.Normals;
+
+		// Convert FVector3d to FVector
+		TArray<FVector> Vertices;
+		Vertices.Reserve(Generator.Vertices.Num());
+		for (const FVector3d& Vertex : Vertices3d)
 		{
-			SectionTriangles.Add(Triangles[TriangleIndex * 3 + i] - VertexIndex);
+			Vertices.Add(FVector(Vertex.X, Vertex.Y, Vertex.Z));
 		}
 
-		// Create or update the mesh section
-		ProceduralMesh->CreateMeshSection(SectionIndex, SectionVertices, SectionTriangles, TArray<FVector>(), TArray<FVector2D>(), SectionVertexColors, TArray<FProcMeshTangent>(), false);
-	}
+		// Convert FIndex3i to int32
+		TArray<int32> Indices;
+		Indices.Reserve(Generator.Triangles.Num() * 3);
+		for (const UE::Geometry::FIndex3i& Triangle : Triangles)
+		{
+			Indices.Add(Triangle.A);
+			Indices.Add(Triangle.B);
+			Indices.Add(Triangle.C);
+		}
+
+		// Convert FVector3f to FVector
+		TArray<FVector> Normals;
+		Normals.Reserve(Generator.Normals.Num());
+		for (const FVector3f& Normal : Normals3F)
+		{
+			Normals.Add(FVector(Normal.X, Normal.Y, Normal.Z));
+		}
+
+		TArray<FColor> VertexColors;
+		VertexColors.Reserve(Vertices.Num());
+		for (int32 i = 0; i < Vertices.Num(); ++i)
+		{
+			VertexColors.Add(FColor::Green);
+		}
+
+
+		AsyncTask(ENamedThreads::Type::GameThread, [this, Vertices, Indices, Normals, VertexColors]()
+		{
+			ProceduralMesh->CreateMeshSection(0, Vertices, Indices, Normals, TArray<FVector2D>(), VertexColors, TArray<FProcMeshTangent>(), false);
+			ProceduralMesh->SetVisibility(true);
+			AccumulatedTime = 0.0f;
+		});
+	});
+	
+
+	
 }
\ No newline at end of file
diff --git a/Source/MetaCastBachelor/MetaPoint.h b/Source/MetaCastBachelor/MetaPoint.h
index e221d596fa68eed248fa363f616f08c7d2d4eb3f..ecb85c63abdc3ef3818c920993345c2302df7df8 100644
--- a/Source/MetaCastBachelor/MetaPoint.h
+++ b/Source/MetaCastBachelor/MetaPoint.h
@@ -4,7 +4,7 @@
 #include "CoreMinimal.h"
 #include "ProceduralMeshComponent.h"
 #include "MetaCastBaseline.h"
-#include "MetaCastBachelor/MarchingCubes.h"
+#include "Generators/MarchingCubes.h"
 #include "MetaPoint.generated.h"
 
 UCLASS( ClassGroup=(Custom), meta=(BlueprintSpawnableComponent) )
@@ -12,7 +12,6 @@ class UMetaPoint : public UMetaCastBaseline
 {
 	GENERATED_BODY()
 	
-	MarchingCubes* MyMarchingCubes;
 	FVector WorldMaximum;
 	int32 LocalMaximumIndex;
 	FVector LocalMaximum;
@@ -21,6 +20,10 @@ class UMetaPoint : public UMetaCastBaseline
 	UPROPERTY()
 	UWorld* World;
 	float AccumulatedTime = 0.0;
+	FVector SelectStartPosition;
+	FVector LastHandInput;
+	UPROPERTY(EditAnywhere)
+	float MinimalHandMovement = 0.01;
 
 	UPROPERTY()
 	UProceduralMeshComponent* ProceduralMesh;
@@ -34,7 +37,7 @@ class UMetaPoint : public UMetaCastBaseline
 	TArray<int32> RevisitIndices;
 
 	UPROPERTY(EditAnywhere)
-	int MetaCastPerSecond = 10;
+	int MetaCastPerSecond = 1;
 	
 public:
 	UMetaPoint();
@@ -48,8 +51,8 @@ public:
 	virtual void HandleMetaSelectPressed(const FInputActionInstance& Instance) override;
 	void InitMetaPointSelection();
 
-	void GenerateVoxelMesh(const TArray<int32> Voxels) const;
-	void CreateMeshSections(const TArray<FVector>& Vertices, const TArray<int32>& Triangles, const TArray<FColor>& VertexColors, int32 MaxVerticesPerSection) const;
+	void GenerateVoxelMeshWithCubes(const TArray<int32> Voxels) const;
+	void GenerateVoxelMeshSmooth(const TArray<int32> Voxels);
 
 	virtual void TickComponent(float DeltaTime, ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction) override;
 };
diff --git a/Source/MetaCastBachelor/MarchingCubes.cpp b/Source/MetaCastBachelor/MyMarchingCubes.cpp
similarity index 86%
rename from Source/MetaCastBachelor/MarchingCubes.cpp
rename to Source/MetaCastBachelor/MyMarchingCubes.cpp
index b44eb56b33d467a77ff415b9b16986d622d1f18f..57862015ed6a0d7493b047df806aeb2b0054d0dc 100644
--- a/Source/MetaCastBachelor/MarchingCubes.cpp
+++ b/Source/MetaCastBachelor/MyMarchingCubes.cpp
@@ -1,16 +1,16 @@
 // Copyright 2017 Patrick Chadbourne (PChadbourne). All Rights Reserved.
 
-#include "MarchingCubes.h"
+#include "MyMarchingCubes.h"
 
-MarchingCubes::MarchingCubes()
+FMyMarchingCubes::FMyMarchingCubes()
 {
 }
 
-MarchingCubes::~MarchingCubes()
+FMyMarchingCubes::~FMyMarchingCubes()
 {
 }
 
-void MarchingCubes::Polyganize(FCell Cell, FMeshData* Data)
+void FMyMarchingCubes::Polyganize(FCell Cell, FMeshData* Data)
 {
     int CubeIndex = 0;
     FVector VertList[12];
@@ -72,12 +72,10 @@ void MarchingCubes::Polyganize(FCell Cell, FMeshData* Data)
     }
 }
 
-
-
-void MarchingCubes::ProcessGrid(TArray<float> Grid, int Width, FMeshData* Data, FVector A, FVector B, FVector C, FVector D)
+void FMyMarchingCubes::ProcessGrid(TArray<float> Grid, int Width, FMeshData* Data, FVector A, FVector B, FVector C, FVector D)
 {
-    int WidthSquared = Width * Width;
-    int GridSize = Width * Width * Width;
+	const int WidthSquared = Width * Width;
+	const int GridSize = Width * Width * Width;
     
     // Ensure the Grid array has enough elements
     if (Grid.Num() < GridSize)
@@ -94,14 +92,14 @@ void MarchingCubes::ProcessGrid(TArray<float> Grid, int Width, FMeshData* Data,
         {
             for (int k = 0; k < Width - 1; k++)
             {
-                int idx7 = i + Width * j + WidthSquared * k;
-                int idx4 = idx7 + 1;
-                int idx6 = idx7 + Width;
-                int idx5 = idx6 + 1;
-                int idx3 = idx7 + WidthSquared;
-                int idx0 = idx3 + 1;
-                int idx2 = idx3 + Width;
-                int idx1 = idx2 + 1;
+	            const int idx7 = i + Width * j + WidthSquared * k;
+                const int idx4 = idx7 + 1;
+                const int idx6 = idx7 + Width;
+                const int idx5 = idx6 + 1;
+                const int idx3 = idx7 + WidthSquared;
+                const int idx0 = idx3 + 1;
+                const int idx2 = idx3 + Width;
+                const int idx1 = idx2 + 1;
 
                 // Ensure all indices are within bounds
                 if (idx7 >= Grid.Num() || idx4 >= Grid.Num() || idx6 >= Grid.Num() || idx5 >= Grid.Num() ||
@@ -156,9 +154,7 @@ void MarchingCubes::ProcessGrid(TArray<float> Grid, int Width, FMeshData* Data,
     }
 }
 
-
-
-FVector MarchingCubes::VertexInterp(FVector P1, FVector P2, float P1Val, float P2Val, float Value)
+FVector FMyMarchingCubes::VertexInterp(FVector P1, FVector P2, float P1Val, float P2Val, float Value)
 {
 
 
@@ -179,7 +175,7 @@ FVector MarchingCubes::VertexInterp(FVector P1, FVector P2, float P1Val, float P
 	return(P);
 }
 
-FVector MarchingCubes::CalculateTriangularInterpolation(int i, int j, int k, FVector A, FVector B, FVector C, FVector D, int Width)
+FVector FMyMarchingCubes::CalculateTriangularInterpolation(int i, int j, int k, FVector A, FVector B, FVector C, FVector D, int Width)
 {
 	if (i + j < Width) {
 		FVector alpha = FMath::Lerp(A, B, (i*2.0f) / (float)Width);
diff --git a/Source/MetaCastBachelor/MarchingCubes.h b/Source/MetaCastBachelor/MyMarchingCubes.h
similarity index 99%
rename from Source/MetaCastBachelor/MarchingCubes.h
rename to Source/MetaCastBachelor/MyMarchingCubes.h
index ea31fc91287d0d4483bc195f8f50a31fbe30187c..d48a2ed843e217fb157133fdce6299ccf28d14a9 100644
--- a/Source/MetaCastBachelor/MarchingCubes.h
+++ b/Source/MetaCastBachelor/MyMarchingCubes.h
@@ -19,11 +19,11 @@ public:
 	TArray<FVector> N;
 };
 
-class MarchingCubes
+class FMyMarchingCubes
 {
 public:
-	MarchingCubes();
-	~MarchingCubes();
+	FMyMarchingCubes();
+	~FMyMarchingCubes();
 
 	void ProcessGrid(TArray<float> Grid, int Width, FMeshData* Data, FVector A, FVector B, FVector C, FVector D);
 
diff --git a/Source/MetaCastBachelor/PointCloud.cpp b/Source/MetaCastBachelor/PointCloud.cpp
index 747d186f28458c275465ed59f67bcbad0c023d37..d319ce91f39e29ecdb9443e723229bac0bd796cc 100644
--- a/Source/MetaCastBachelor/PointCloud.cpp
+++ b/Source/MetaCastBachelor/PointCloud.cpp
@@ -85,22 +85,9 @@ void APointCloud::UpdateBounds()
 void APointCloud::SetupDensityFieldFromPointCloud() const
 {
 	UE_LOG(LogTemp, Log, TEXT("Initializing DensityField!"));
-	
-	constexpr int32 XAxisNum = 100; // Number of divisions in the grid along the X-axis
-	constexpr int32 YAxisNum = 100; // Number of divisions in the grid along the Y-axis
-	constexpr int32 ZAxisNum = 100; // Number of divisions in the grid along the Z-axis
-	MyDensityField->InitializeDensityField(MinBounds, MaxBounds, XAxisNum, YAxisNum, ZAxisNum);
-
-	constexpr float InfluenceRadius = 2.0f; // Half the size of a voxel
-	constexpr float Sigma = 0.5f;  // Standard deviation for the Gaussian kernel
-	
-	MyDensityField->CalculateVoxelDensitiesByNumber(this, InfluenceRadius);
-	MyDensityField->CalculateAndStoreGradients();
 
-	if (const UWorld* World = GetWorld())
-	{
-		//MyDensityField->DrawDebugVoxelDensity(World, this, 100.0f, 1.0f, 0.8f);
-	}
+	const FIntVector AxisNumbers(100, 100, 100);
+	MyDensityField->InitializeDensityField(this, MinBounds, MaxBounds, AxisNumbers);
 }
 
 void APointCloud::UpdateSelection()
@@ -111,7 +98,7 @@ void APointCloud::UpdateSelection()
 		{
 			if (SelectionFlags[i])
 			{
-				PointColors[i] = FColor::Red;
+				PointColors[i] = FColor(0, 200, 0);
 			}else
 			{
 				PointColors[i] = FColor::Orange;
@@ -149,57 +136,31 @@ void APointCloud::DrawVoxel(const int Index, const float Time) const
 
 void APointCloud::ColorPointsInVoxels(const TArray<int32> VoxelIndices)
 {
+	FScopeLock Lock(&DataGuard);
 	if (!MyDensityField)
 	{
 		UE_LOG(LogTemp, Warning, TEXT("DensityField is not initialized."));
 		return;
 	}
 
-	TArray<FVector> MyPositionVectors = PositionVectors;
-
-	FScopeLock Lock(&DataGuard);
-	FScopeLock Lock2(&MyDensityField->DataGuard);
-	// Iterate through each index in VoxelIndices
 	for (const int32 VoxelIndex : VoxelIndices)
 	{
-		// Ensure the voxel index is valid
-		if (!MyDensityField || !MyDensityField->IsValidIndex(VoxelIndex))
-		{
-			UE_LOG(LogTemp, Warning, TEXT("Invalid voxel index: %d"), VoxelIndex);
-			continue;
-		}
-
-		// Get the bounds of the voxel
-		FVector VoxelMin, VoxelMax;
-		MyDensityField->IndexToVoxelBounds(VoxelIndex, VoxelMin, VoxelMax);
-
-		// Iterate over all points to check if they are within this voxel
-		for (int32 i = 0; i < MyPositionVectors.Num(); ++i)
+		const TArray<int32>& PointIndices = MyDensityField->VoxelPointLookupTable->GetPointsInVoxel(VoxelIndex);
+		for (const int32 PointIndex : PointIndices)
 		{
-			if(!MyPositionVectors.IsValidIndex(i))
+			if(SelectionFlags.IsValidIndex(PointIndex))
 			{
-				continue;
-			}
-			const FVector Point = MyPositionVectors[i];
-
-			// Check if the point is within the voxel bounds
-			if (Point.X >= VoxelMin.X && Point.X <= VoxelMax.X &&
-				Point.Y >= VoxelMin.Y && Point.Y <= VoxelMax.Y &&
-				Point.Z >= VoxelMin.Z && Point.Z <= VoxelMax.Z)
-			{
-				SelectionFlags[i] = true;
+				SelectionFlags[PointIndex] = true;
 			}
 		}
 	}
-
-	//UpdateSelection();
 }
 
+
 // Called every frame
 void APointCloud::Tick(float DeltaTime)
 {
 	Super::Tick(DeltaTime);
-
 	
 	UpdateSelection();
 }
diff --git a/Source/MetaCastBachelor/PointCloud.h b/Source/MetaCastBachelor/PointCloud.h
index 4770f1b4e6560b2f146abcf5a251f69d33c80ec4..6aae1269053c7110a6e69381a2e2bea7bb3600cf 100644
--- a/Source/MetaCastBachelor/PointCloud.h
+++ b/Source/MetaCastBachelor/PointCloud.h
@@ -13,9 +13,13 @@ class APointCloud : public AActor
 	
 public:
 	FDensityField* MyDensityField;
+	UPROPERTY()
 	TArray<FVector> PositionVectors;
+	UPROPERTY()
 	TArray<bool> DefaultFlags;
+	UPROPERTY()
 	TArray<bool> SelectionFlags;
+	UPROPERTY()
 	TArray<FColor> PointColors;
 	mutable FCriticalSection DataGuard;
 
diff --git a/Source/MetaCastBachelor/Utilities.cpp b/Source/MetaCastBachelor/Utilities.cpp
index 236264834ff8020d4e791220713fb65bcd11b198..91ab449b6ad2f9f62b3a0a4295633ab4c3af2e30 100644
--- a/Source/MetaCastBachelor/Utilities.cpp
+++ b/Source/MetaCastBachelor/Utilities.cpp
@@ -31,24 +31,21 @@ double FUtilities::InterpolateDensityAtPosition(const FDensityField* DensityFiel
     if (!DensityField)
         return 0.0;
 
-    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);
-
-    const float XFraction = GridPos.X - XBin;
-    const float YFraction = GridPos.Y - YBin;
-    const float ZFraction = GridPos.Z - ZBin;
-
-    const double D000 = DensityField->GridPositionToVoxelDensity(XBin, YBin, ZBin);
-    const double D100 = DensityField->GridPositionToVoxelDensity(XBin + 1, YBin, ZBin);
-    const double D010 = DensityField->GridPositionToVoxelDensity(XBin, YBin + 1, ZBin);
-    const double D001 = DensityField->GridPositionToVoxelDensity(XBin, YBin, ZBin + 1);
-    const double D101 = DensityField->GridPositionToVoxelDensity(XBin + 1, YBin, ZBin + 1);
-    const double D011 = DensityField->GridPositionToVoxelDensity(XBin, YBin + 1, ZBin + 1);
-    const double D110 = DensityField->GridPositionToVoxelDensity(XBin + 1, YBin + 1, ZBin);
-    const double D111 = DensityField->GridPositionToVoxelDensity(XBin + 1, YBin + 1, ZBin + 1);
+    const FVector GridPositionFraction = (Position - DensityField->GetGridOrigin()) / DensityField->GetStep();
+    const FIntVector3 GridPosition = DensityField->WorldToGridPosition(Position - DensityField->GetGridOrigin());
+	
+    const float XFraction = GridPositionFraction.X - GridPosition.X;
+    const float YFraction = GridPositionFraction.Y - GridPosition.Y;
+    const float ZFraction = GridPositionFraction.Z - GridPosition.Z;
+
+    const double D000 = DensityField->GridPositionToVoxelDensity(GridPosition.X, GridPosition.Y, GridPosition.Z);
+    const double D100 = DensityField->GridPositionToVoxelDensity(GridPosition.X + 1, GridPosition.Y, GridPosition.Z);
+    const double D010 = DensityField->GridPositionToVoxelDensity(GridPosition.X, GridPosition.Y + 1, GridPosition.Z);
+    const double D001 = DensityField->GridPositionToVoxelDensity(GridPosition.X, GridPosition.Y, GridPosition.Z + 1);
+    const double D101 = DensityField->GridPositionToVoxelDensity(GridPosition.X + 1, GridPosition.Y, GridPosition.Z + 1);
+    const double D011 = DensityField->GridPositionToVoxelDensity(GridPosition.X, GridPosition.Y + 1, GridPosition.Z + 1);
+    const double D110 = DensityField->GridPositionToVoxelDensity(GridPosition.X + 1, GridPosition.Y + 1, GridPosition.Z);
+    const double D111 = DensityField->GridPositionToVoxelDensity(GridPosition.X + 1, GridPosition.Y + 1, GridPosition.Z + 1);
 
     // Trilinear interpolation
     const double DX00 = FMath::Lerp(D000, D100, XFraction);
@@ -70,25 +67,22 @@ FVector FUtilities::InterpolateGradientAtPosition(const FDensityField* DensityFi
 		return FVector::ZeroVector;
 
 	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);
+	const FIntVector3 GridPosition = DensityField->WorldToGridPosition(Position - DensityField->GetGridOrigin());
 
 	// Fractional parts to interpolate
-	float XFraction = GridPos.X - XBin;
-	float YFraction = GridPos.Y - YBin;
-	float ZFraction = GridPos.Z - ZBin;
+	float XFraction = GridPos.X - GridPosition.X;
+	float YFraction = GridPos.Y - GridPosition.Y;
+	float ZFraction = GridPos.Z - GridPosition.Z;
 
 	// Fetch gradients from corners of the cell
-	FVector G000 = DensityField->GridPositionToVoxelGradient(XBin, YBin, ZBin);
-	FVector G100 = DensityField->GridPositionToVoxelGradient(XBin + 1, YBin, ZBin);
-	FVector G010 = DensityField->GridPositionToVoxelGradient(XBin, YBin + 1, ZBin);
-	FVector G001 = DensityField->GridPositionToVoxelGradient(XBin, YBin, ZBin + 1);
-	FVector G110 = DensityField->GridPositionToVoxelGradient(XBin + 1, YBin + 1, ZBin);
-	FVector G101 = DensityField->GridPositionToVoxelGradient(XBin + 1, YBin, ZBin + 1);
-	FVector G011 = DensityField->GridPositionToVoxelGradient(XBin, YBin + 1, ZBin + 1);
-	FVector G111 = DensityField->GridPositionToVoxelGradient(XBin + 1, YBin + 1, ZBin + 1);
+	FVector G000 = DensityField->GridPositionToVoxelGradient(GridPosition.X, GridPosition.Y, GridPosition.Z);
+	FVector G100 = DensityField->GridPositionToVoxelGradient(GridPosition.X + 1, GridPosition.Y, GridPosition.Z);
+	FVector G010 = DensityField->GridPositionToVoxelGradient(GridPosition.X, GridPosition.Y + 1, GridPosition.Z);
+	FVector G001 = DensityField->GridPositionToVoxelGradient(GridPosition.X, GridPosition.Y, GridPosition.Z + 1);
+	FVector G110 = DensityField->GridPositionToVoxelGradient(GridPosition.X + 1, GridPosition.Y + 1, GridPosition.Z);
+	FVector G101 = DensityField->GridPositionToVoxelGradient(GridPosition.X + 1, GridPosition.Y, GridPosition.Z + 1);
+	FVector G011 = DensityField->GridPositionToVoxelGradient(GridPosition.X, GridPosition.Y + 1, GridPosition.Z + 1);
+	FVector G111 = DensityField->GridPositionToVoxelGradient(GridPosition.X + 1, GridPosition.Y + 1, GridPosition.Z + 1);
 
 	// Interpolate along x for each of the yz pairs
 	FVector G00 = FMath::Lerp(G000, G100, XFraction);
@@ -110,6 +104,11 @@ float FUtilities::GaussianKernel(const float Distance, const float Sigma) {
 	return FMath::Exp(-0.5 * FMath::Square(Distance / Sigma));
 }
 
+float FUtilities::SmoothingKernel(const float Distance, const float Radius) {
+	const float Value = FMath::Max(0.0f, Radius * Radius - Distance * Distance);
+	return Value * Value * Value;
+}
+
 FVector FUtilities::FollowGradientToMaximum(const FDensityField* DensityField, const FVector& StartPosition)
 {
 	FVector CurrentPosition = StartPosition;
@@ -117,20 +116,20 @@ FVector FUtilities::FollowGradientToMaximum(const FDensityField* DensityField, c
 
 	int Iterations = 0;
 	constexpr int MaxIterations = 100;  // Prevent infinite loops
-	constexpr float StepSize = 0.1f;  // Define a reasonable step size
+	constexpr float StepSize = 0.005f;  // Define a reasonable step size
 
 	while (Iterations < MaxIterations)
 	{
 		FVector Gradient = InterpolateGradientAtPosition(DensityField, CurrentPosition);
 		if (Gradient.IsNearlyZero()) {
-			UE_LOG(LogTemp, Log, TEXT("Gradient is zero at position %s"), *CurrentPosition.ToString());
+			//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->IsValidWorldPosition(NextPosition)) {
-			UE_LOG(LogTemp, Warning, TEXT("Next position out of bounds: %s"), *NextPosition.ToString());
+			//UE_LOG(LogTemp, Warning, TEXT("Next position out of bounds: %s"), *NextPosition.ToString());
 			break;  // Next position is out of valid bounds
 		}
 
@@ -184,9 +183,9 @@ TArray<int32> FUtilities::FloodFilling(const FDensityField* DensityField, const
 			for (int32 NeighborIndex : Neighbors)
 			{
 				const double NeighborDensity = DensityField->IndexToVoxelDensity(NeighborIndex);
-				//UE_LOG(LogTemp, Log, TEXT("Lookign at Neighbor %d at Position %s with density: %f."), NeighborIndex, *DensityField->IndexToGridPosition(NeighborIndex).ToString(), NeighborDensity);
+				//UE_LOG(LogTemp, Log, TEXT("Looking at Neighbor %d at Position %s with density: %f."), NeighborIndex, *DensityField->IndexToGridPosition(NeighborIndex).ToString(), NeighborDensity);
 				
-				if (NeighborDensity > Threshold && !Visited[NeighborIndex]) {
+				if (NeighborDensity > Threshold * 1.1f && !Visited[NeighborIndex]) {
 					Stack.Push(NeighborIndex);
 					//UE_LOG(LogTemp, Log, TEXT("Pushing Neighbor %d to stack."), NeighborIndex);
 				}
@@ -203,55 +202,71 @@ TArray<int32> FUtilities::FloodFilling(const FDensityField* DensityField, const
 	return VisitedIndices;
 }
 
-void FUtilities::FloodFilling_2(const FDensityField* DensityField, const float Threshold, TArray<bool>& VisitedIndices, TArray<int32>& IndicesToVisit, TArray<int32>& FloodedIndices, TArray<int32>& PotentialRevisit)
+TArray<int32> FUtilities::FloodFilling_2(const FDensityField* DensityField, const int32 StartIndex, const float Threshold, const float MaxDistance)
 {
+	//UE_LOG(LogTemp, Warning, TEXT("Threshold for FloodFilling: %f"), Threshold);
+	
+	TArray<int32> VisitedIndices;
 	if (!DensityField) {
 		UE_LOG(LogTemp, Warning, TEXT("DensityField is null."));
-		return;
+		return VisitedIndices;
 	}
 
-	TArray<int32> CurrentPotentialRevisit;
+	if (!DensityField->IsValidIndex(StartIndex)) {
+		UE_LOG(LogTemp, Warning, TEXT("Start index %d is invalid."), StartIndex);
+		return VisitedIndices;
+	}
 
-	
-	constexpr int MaxIterations = 2000;
-	int IterationCount = 0;
-	
-	while (IndicesToVisit.Num() > 0)
+	TArray<int32> Stack;
+	Stack.Push(StartIndex);
+	TArray<bool> Visited;
+	Visited.Init(false, DensityField->GetVoxelNumber());
+
+	while (Stack.Num() > 0)
 	{
-		
-		IterationCount++;
-		if(IterationCount > MaxIterations)
+		int32 CurrentIndex = Stack.Pop();
+		if (!Visited[CurrentIndex])
 		{
-			while (IndicesToVisit.Num() > 0) {
-				CurrentPotentialRevisit.Add(IndicesToVisit.Pop());
-			}
-			break;
-		}
-		
-		int32 CurrentIndex = IndicesToVisit.Pop();
+			Visited[CurrentIndex] = true;
+			VisitedIndices.Add(CurrentIndex);
 
-		if(DensityField->IsValidIndex(CurrentIndex))
-		{
-			VisitedIndices[CurrentIndex] = true;
-			const double CurrentDensity = DensityField->IndexToVoxelDensity(CurrentIndex);
+			TArray<int32> Neighbors = DensityField->IndexToVoxelNeighbors(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)
+			{
+				const double NeighborDensity = DensityField->IndexToVoxelDensity(NeighborIndex);
+				//UE_LOG(LogTemp, Log, TEXT("Looking at Neighbor %d at Position %s with density: %f."), NeighborIndex, *DensityField->IndexToGridPosition(NeighborIndex).ToString(), NeighborDensity);
 
-			if (CurrentDensity >= Threshold) {
-				FloodedIndices.Add(CurrentIndex);
+				const float DistanceToStart = FVector::Dist(DensityField->IndexToVoxelPosition(StartIndex), DensityField->IndexToVoxelPosition(NeighborIndex));
+				
+				if (!Visited[NeighborIndex]){
+					if(DistanceToStart <= MaxDistance)
+					{
+						if(NeighborDensity > Threshold * 1.0f)
+						{
+							Stack.Push(NeighborIndex);
+						}
 
-				TArray<int32> Neighbors = DensityField->IndexToVoxelNeighbors(CurrentIndex);
-				for (int32 NeighborIndex : Neighbors)
-				{
-					if (!VisitedIndices[NeighborIndex])
+					}else
 					{
-						IndicesToVisit.Push(NeighborIndex);
+						const float Multiplier = 1 + FMath::Pow(DistanceToStart - MaxDistance, 2);
+						if(NeighborDensity > Threshold * Multiplier)
+						{
+							Stack.Push(NeighborIndex);
+						}
 					}
+					//UE_LOG(LogTemp, Log, TEXT("Pushing Neighbor %d to stack."), NeighborIndex);
 				}
-			} else {
-				CurrentPotentialRevisit.Add(CurrentIndex);
 			}
 		}
 	}
 
-	// Update PotentialRevisit with the new list from this run
-	PotentialRevisit = CurrentPotentialRevisit;
+	// 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;
 }
diff --git a/Source/MetaCastBachelor/Utilities.h b/Source/MetaCastBachelor/Utilities.h
index 6d406d74976619aeb28b4b9419e1dd45daf112f1..b453f8263106d69d8614f128c6a7a66b78f58143 100644
--- a/Source/MetaCastBachelor/Utilities.h
+++ b/Source/MetaCastBachelor/Utilities.h
@@ -9,7 +9,8 @@ public:
 	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 float SmoothingKernel(float Distance, float Radius);
 	static FVector FollowGradientToMaximum(const FDensityField* DensityField, const FVector& StartPosition);
 	static TArray<int32> FloodFilling(const FDensityField* DensityField, const int32 StartIndex, const float Threshold);
-	static void FloodFilling_2(const FDensityField* DensityField, const float Threshold, TArray<bool>& VisitedIndices, TArray<int32>& IndicesToVisit, TArray<int32>& FloodedIndices, TArray<int32>& PotentialRevisit);
+	static TArray<int32> FloodFilling_2(const FDensityField* DensityField, const int32 StartIndex, const float Threshold, const float MaxDistance);
 };
diff --git a/Source/MetaCastBachelor/VoxelPointLookupTable.cpp b/Source/MetaCastBachelor/VoxelPointLookupTable.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a9499c20304161e48985ac61d9cf723d90a450a4
--- /dev/null
+++ b/Source/MetaCastBachelor/VoxelPointLookupTable.cpp
@@ -0,0 +1,29 @@
+#include "VoxelPointLookupTable.h"
+#include "DensityField.h"
+
+
+FVoxelPointLookupTable::FVoxelPointLookupTable() : LinkedDensityField(nullptr) {}
+
+FVoxelPointLookupTable::~FVoxelPointLookupTable(){}
+
+void FVoxelPointLookupTable::Initialize(const TArray<FVector>& PositionVectors, const FDensityField* DensityField)
+{
+	LinkedDensityField = DensityField;
+	if (!LinkedDensityField) return;
+
+	for (int32 i = 0; i < PositionVectors.Num(); ++i)
+	{
+		const FVector& Point = PositionVectors[i];
+		int32 VoxelIndex = LinkedDensityField->WorldPositionToIndex(Point);
+		if (VoxelIndex != -1)
+		{
+			VoxelToPointIndicesMap.FindOrAdd(VoxelIndex).Add(i);
+		}
+	}
+}
+
+TArray<int32> FVoxelPointLookupTable::GetPointsInVoxel(const int32 VoxelIndex) const
+{
+	return VoxelToPointIndicesMap.Contains(VoxelIndex) ? VoxelToPointIndicesMap[VoxelIndex] : TArray<int32>();
+}
+
diff --git a/Source/MetaCastBachelor/VoxelPointLookupTable.h b/Source/MetaCastBachelor/VoxelPointLookupTable.h
new file mode 100644
index 0000000000000000000000000000000000000000..52076a383f10db918ba8b2d2209b593999648435
--- /dev/null
+++ b/Source/MetaCastBachelor/VoxelPointLookupTable.h
@@ -0,0 +1,26 @@
+#pragma once
+#include "CoreMinimal.h"
+
+class FDensityField;
+
+class FVoxelPointLookupTable
+{
+
+public:
+	FVoxelPointLookupTable();
+
+	virtual ~FVoxelPointLookupTable();
+
+	// Initializes the lookup table with position vectors and a reference to the density field
+	void Initialize(const TArray<FVector>& PositionVectors, const FDensityField* DensityField);
+
+	// Retrieves the array of point indices for a specific voxel index
+	TArray<int32> GetPointsInVoxel(int32 VoxelIndex) const;
+
+private:
+	// Maps voxel indices to arrays of point indices
+	TMap<int32, TArray<int32>> VoxelToPointIndicesMap;
+
+	// Reference to the density field for voxel bounds queries
+	const FDensityField* LinkedDensityField;
+};