diff --git a/Content/MagicWandMap.umap b/Content/MagicWandMap.umap
new file mode 100644
index 0000000000000000000000000000000000000000..432e39ca2127b23880bf35d2220201a73ca86de3
Binary files /dev/null and b/Content/MagicWandMap.umap differ
diff --git a/Content/MagicWandPawn.uasset b/Content/MagicWandPawn.uasset
new file mode 100644
index 0000000000000000000000000000000000000000..9472186a5940e9e6147f87a741f1b214548baeff
Binary files /dev/null and b/Content/MagicWandPawn.uasset differ
diff --git a/Content/MetaPointMap.umap b/Content/MetaPointMap.umap
index a90f205a3061d3582938745c1b4c3550588e57fa..e65a8b10eadfdfb9712a93fbbe490b559b6eb5ed 100644
Binary files a/Content/MetaPointMap.umap and b/Content/MetaPointMap.umap differ
diff --git a/Source/MetaCastBachelor/DensityField.cpp b/Source/MetaCastBachelor/DensityField.cpp
index 6c6af6b8a896b8c791f2ee8170ca6e2ceccfb4ba..0af1af8019c8539660da23f9d575f761fb0c33de 100644
--- a/Source/MetaCastBachelor/DensityField.cpp
+++ b/Source/MetaCastBachelor/DensityField.cpp
@@ -5,7 +5,7 @@
 
 // INITIALIZATION FUNCTIONS
 
-FDensityField::FDensityField(): XNum(0), YNum(0), ZNum(0), XStep(0), YStep(0), ZStep(0), VoxelPointLookupTable(nullptr)
+FDensityField::FDensityField(): XNum(0), YNum(0), ZNum(0), XStep(0), YStep(0), ZStep(0), MaxDensity(0), MinDensity(0), AvgDensity(0), VoxelPointLookupTable(nullptr)
 {
 	VoxelList = TArray<FVoxel>();
 }
@@ -17,8 +17,8 @@ void FDensityField::InitializeDensityField(const APointCloud* PointCloud, const
 	VoxelPointLookupTable = new FVoxelPointLookupTable();
 	VoxelPointLookupTable->Initialize(PointCloud->PositionVectors, this);
 	
-	CalculateVoxelDensities_2(PointCloud, PointCloud->InfluenceRadius);
-	//CalculateVoxelDensitiesByNumber(PointCloud, 4);
+	//CalculateVoxelDensities_2(PointCloud, PointCloud->InfluenceRadius);
+	CalculateVoxelDensitiesSimple(PointCloud, PointCloud->InfluenceRadius);
 	CalculateAndStoreGradients();
 }
 
@@ -101,7 +101,7 @@ void FDensityField::CalculateVoxelDensities(const APointCloud* PointCloud, const
 	        }
 	    }
 	}
-	double MaxDensity = 0.0;
+	MaxDensity = 0.0;
 	for (const auto& Node : VoxelList)
 	{
 		if (Node.GetVoxelDensity() > MaxDensity)
@@ -118,7 +118,7 @@ void FDensityField::CalculateVoxelDensities(const APointCloud* PointCloud, const
 
 }
 
-void FDensityField::CalculateVoxelDensitiesByNumber(const APointCloud* PointCloud, const float InfluenceRadius)
+void FDensityField::CalculateVoxelDensitiesSimple(const APointCloud* PointCloud, const float InfluenceRadius)
 {
     if (!PointCloud)
         return;
@@ -160,7 +160,8 @@ void FDensityField::CalculateVoxelDensitiesByNumber(const APointCloud* PointClou
 
                         if (Distance < InfluenceRadius)
                         {
-                            VoxelList[Index].AddToVoxelDensity(1.0);  // Simply increment the density count
+	                        const float Weight = FUtilities::SmoothingKernel(Distance, InfluenceRadius);
+                            VoxelList[Index].AddToVoxelDensity(Weight);  // Simply increment the density count
                         }
                     }
                 }
@@ -168,15 +169,21 @@ void FDensityField::CalculateVoxelDensitiesByNumber(const APointCloud* PointClou
         }
     }
 
-    double MaxDensity = 0.0;
+    MaxDensity = MIN_flt;
+    MinDensity = MAX_FLT;
+	AvgDensity = 0;
     for (const auto& Node : VoxelList)
     {
-        if (Node.GetVoxelDensity() > MaxDensity)
-        {
-            MaxDensity = Node.GetVoxelDensity();
-        }
+	    const float Dens = Node.GetVoxelDensity();
+    	MaxDensity = FMath::Max(Dens, MaxDensity);
+    	MinDensity = FMath::Min(Dens, MinDensity);
+    	AvgDensity += Dens;
     }
-    UE_LOG(LogTemp, Log, TEXT("Maximum density found: %f"), MaxDensity);
+	AvgDensity /= VoxelList.Num();
+	
+	UE_LOG(LogTemp, Log, TEXT("Maximum density found: %f"), MaxDensity);
+	UE_LOG(LogTemp, Log, TEXT("Minimum density found: %f"), MinDensity);
+	UE_LOG(LogTemp, Log, TEXT("Average density found: %f"), AvgDensity);
 
 	const double EndTime = FPlatformTime::Seconds();  // End timing
 	const double ElapsedTime = EndTime - StartTime;  // Calculate elapsed time
@@ -412,14 +419,15 @@ void FDensityField::CalculateAndStoreGradients()
 
 //	DEBUG FUNCTIONS
 
-void FDensityField::DrawDebugVoxelDensity(const UWorld* World, const AActor* Actor, const float Duration, const float Scale, const float DensityThreshold) const
+void FDensityField::DrawDebugVoxelDensity(const UWorld* World, const AActor* Actor, const float Duration, const float Scale, const float DensityThreshold)
 {
 	if (!World || !Actor)
 	{
 		return;
 	}
 
-	double MinDensity = DBL_MAX, MaxDensity = 0;
+	MinDensity = DBL_MAX;
+	MaxDensity = 0;
 	for (const FVoxel& Node : VoxelList)
 	{
 		const double LogDensity = log10(Node.GetVoxelDensity() + 1);  // Avoiding log(0)
diff --git a/Source/MetaCastBachelor/DensityField.h b/Source/MetaCastBachelor/DensityField.h
index 5a3cd269d9704f36933c56edfb8e44e5c564eaf4..1e86417bf5ac8d7b2b8bda3f21f85d75eebd9855 100644
--- a/Source/MetaCastBachelor/DensityField.h
+++ b/Source/MetaCastBachelor/DensityField.h
@@ -42,12 +42,18 @@ class FDensityField
     float XStep, YStep, ZStep;
 	FVector GridOrigin;
 
+	
+
 	void CalculateVoxelDensities(const APointCloud* PointCloud, float InfluenceRadius, float Sigma);
-	void CalculateVoxelDensitiesByNumber(const APointCloud* PointCloud, float InfluenceRadius);
+	void CalculateVoxelDensitiesSimple(const APointCloud* PointCloud, float InfluenceRadius);
 	void CalculateAndStoreGradients();
+	void DrawDebugVoxelDensity(const UWorld* World, const AActor* Actor, float Duration, float Scale, float DensityThreshold);
 	void InitializeDataStructures(const FVector& MinBounds, const FVector& MaxBounds, int32 XAxisNum, int32 YAxisNum, int32 ZAxisNum);
 	
 public:
+	double MaxDensity;
+	double MinDensity;
+	double AvgDensity;
 	FVoxelPointLookupTable* VoxelPointLookupTable;
 	mutable FCriticalSection DataGuard;
 	// CONSTRUCTOR
diff --git a/Source/MetaCastBachelor/MagicWand.cpp b/Source/MetaCastBachelor/MagicWand.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..360602aa66886594847bc8a45cf18126066a2a36
--- /dev/null
+++ b/Source/MetaCastBachelor/MagicWand.cpp
@@ -0,0 +1,315 @@
+#include "MagicWand.h"
+#include "Utilities.h"
+#include "Generators/MarchingCubes.h"
+
+UMagicWand::UMagicWand() : World(nullptr)
+{
+	// Create the procedural mesh component
+	ProceduralMesh = CreateDefaultSubobject<UProceduralMeshComponent>(TEXT("GeneratedMesh"));
+	SelectedClusterIndices.Reserve(100000);
+}
+
+void UMagicWand::BeginPlay()
+{
+	Super::BeginPlay();
+
+	World = GetWorld();
+	
+	if (!World) {
+		UE_LOG(LogTemp, Warning, TEXT("Invalid world provided."));
+	}
+
+	ProceduralMesh->AttachToComponent(MyPointCloud->PointCloudVisualizer, FAttachmentTransformRules::SnapToTargetIncludingScale);
+	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 UMagicWand::TickComponent(float DeltaTime, ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction)
+{
+	Super::TickComponent(DeltaTime, TickType, ThisTickFunction);
+
+	AccumulatedTime += DeltaTime;
+
+	ProceduralMesh->SetVisibility(Select);
+
+	if (Select)
+	{
+		PerformMagicWandSelection(SelectionObject->GetComponentLocation());
+	}
+}
+
+void UMagicWand::PerformMagicWandSelection(const FVector& InputPosition)
+{
+    if (MyPointCloud->SelectionFlags.Num() != MyPointCloud->PositionVectors.Num())
+    {
+        UE_LOG(LogTemp, Warning, TEXT("PerformMagicWandSelection: Positions and SelectionFlags array sizes do not match."));
+        return;
+    }
+
+    // Find the closest point to the input position as the seed
+    int32 SeedIndex = INDEX_NONE;
+    float MinDistance = FLT_MAX;
+
+    for (int32 i = 0; i < MyPointCloud->PositionVectors.Num(); i++)
+    {
+        FVector CurrentPoint = MyPointCloud->PositionVectors[i];
+        CurrentPoint = MyPointCloud->PointCloudVisualizer->GetComponentTransform().TransformPosition(CurrentPoint);
+
+        const float Distance = FVector::Dist(InputPosition, CurrentPoint);
+        if (Distance < MinDistance)
+        {
+            MinDistance = Distance;
+            SeedIndex = i;
+        }
+    }
+
+    if (SeedIndex != INDEX_NONE)
+    {
+        // Clear previous selection
+        SelectedClusterIndices.Empty();
+        for (bool& Flag : MyPointCloud->SelectionFlags)
+        {
+            Flag = false;
+        }
+
+        // Expand selection from the seed
+        ExpandSelection(SeedIndex);
+
+    	for (int32 i = 0; i < SelectedClusterIndices.Num(); i++)
+        {
+	        const int Index = SelectedClusterIndices[i];
+            MyPointCloud->SelectionFlags[Index] = true;
+        }
+    }
+}
+
+void UMagicWand::ExpandSelection(const int32 SeedIndex)
+{
+    TQueue<int32> ToProcess;
+    ToProcess.Enqueue(SeedIndex);
+    SelectedClusterIndices.Add(SeedIndex);
+    MyPointCloud->SelectionFlags[SeedIndex] = true;
+
+    while (!ToProcess.IsEmpty())
+    {
+        int32 CurrentIndex;
+        ToProcess.Dequeue(CurrentIndex);
+        FVector CurrentPoint = MyPointCloud->PositionVectors[CurrentIndex];
+        CurrentPoint = MyPointCloud->PointCloudVisualizer->GetComponentTransform().TransformPosition(CurrentPoint);
+
+        for (int32 i = 0; i < MyPointCloud->PositionVectors.Num(); i++)
+        {
+            if (!MyPointCloud->SelectionFlags[i])
+            {
+                FVector Point = MyPointCloud->PositionVectors[i];
+                Point = MyPointCloud->PointCloudVisualizer->GetComponentTransform().TransformPosition(Point);
+
+                if (FVector::Dist(CurrentPoint, Point) <= ProximityThreshold)
+                {
+                    ToProcess.Enqueue(i);
+                    SelectedClusterIndices.Add(i);
+                    MyPointCloud->SelectionFlags[i] = true;
+                }
+            }
+        }
+    }
+}
+
+
+void UMagicWand::EraseParticles(const FVector& InputPosition)
+{
+	
+}
+
+void UMagicWand::HandleMetaSelectReleased(const FInputActionInstance& Instance)
+{
+	Super::HandleMetaSelectReleased(Instance);
+
+	//ProceduralMesh->ClearAllMeshSections();
+	//ProceduralMesh->SetVisibility(false);
+
+
+	AsyncTask(ENamedThreads::Type::AnyBackgroundHiPriTask, [this]()
+	{
+		MyPointCloud->ColorPointsInVoxels(FloodedIndices);
+	});
+}
+
+void UMagicWand::HandleMetaSelectPressed(const FInputActionInstance& Instance)
+{
+	Super::HandleMetaSelectPressed(Instance);
+
+	UE_LOG(LogTemp, Warning, TEXT("PerformMagicWandSelection"));
+	InitMagicWandSelection();
+}
+
+void UMagicWand::InitMagicWandSelection()
+{	
+	const FVector SelectionWorldPosition = SelectionObject->GetComponentLocation();
+	MyDensityField = MyPointCloud->MyDensityField;
+
+	// 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);
+	
+	ProceduralMesh->ClearAllMeshSections();
+	World = GetWorld();
+	
+	AbortMarchingCubes = true;
+	AccumulatedTime = 0.0f;
+}
+
+void UMagicWand::SelectParticles(const FVector& InputPosition)
+{
+	if (AccumulatedTime >=  1 / EvaluationsPerSecond)
+	{
+		AccumulatedTime = -10000;
+
+		AsyncTask(ENamedThreads::AnyBackgroundThreadNormalTask, [this]()
+		{
+			const FVector SelectionWorldPosition = SelectionObject->GetComponentLocation();
+			
+			
+		});
+	}
+}
+
+void UMagicWand::GenerateVoxelMeshWithCubes(const TArray<int32> Voxels) const
+{
+	TArray<FVector> Vertices;
+	TArray<int32> Triangles;
+	TArray<FColor> VertexColors;
+
+	// Generate cube vertices and triangles for each voxel
+	for (const int32 VoxelIndex : Voxels)
+	{
+		if(!MyDensityField->IsValidIndex(VoxelIndex))
+		{
+			continue;
+		}
+
+		TArray<FVector> VoxelCorners = MyDensityField->IndexToVoxelCornersWorld(VoxelIndex);
+
+		// Add vertices for the voxel
+		const int32 BaseIndex = Vertices.Num();
+		Vertices.Append(VoxelCorners);
+
+		// Define the triangles (12 triangles for 6 faces of the cube)
+		// Each face of the cube has 2 triangles (6 faces * 2 triangles per face = 12 triangles)
+		Triangles.Append({
+			BaseIndex + 0, BaseIndex + 1, BaseIndex + 2, BaseIndex + 0, BaseIndex + 2, BaseIndex + 3, // Bottom face
+			BaseIndex + 4, BaseIndex + 6, BaseIndex + 5, BaseIndex + 4, BaseIndex + 7, BaseIndex + 6, // Top face
+			BaseIndex + 0, BaseIndex + 4, BaseIndex + 1, BaseIndex + 1, BaseIndex + 4, BaseIndex + 5, // Front face
+			BaseIndex + 1, BaseIndex + 5, BaseIndex + 2, BaseIndex + 2, BaseIndex + 5, BaseIndex + 6, // Right face
+			BaseIndex + 2, BaseIndex + 6, BaseIndex + 3, BaseIndex + 3, BaseIndex + 6, BaseIndex + 7, // Back face
+			BaseIndex + 3, BaseIndex + 7, BaseIndex + 0, BaseIndex + 0, BaseIndex + 7, BaseIndex + 4  // Left face
+		});
+		
+		// Add red color for each corner vertex
+		for (int32 i = 0; i < 8; ++i)
+		{
+			VertexColors.Add(FColor::Green);
+		}
+	}
+
+	// Create the mesh section
+	AsyncTask(ENamedThreads::Type::GameThread, [this, Vertices, Triangles, VertexColors]()
+	{
+		//CreateMeshSections(Vertices, Triangles, VertexColors, 1000);
+		ProceduralMesh->CreateMeshSection(0, Vertices, Triangles, TArray<FVector>(), TArray<FVector2D>(), VertexColors, TArray<FProcMeshTangent>(), false);
+		ProceduralMesh->SetVisibility(true);
+	});
+}
+
+void UMagicWand::GenerateVoxelMeshSmooth(const TArray<int32> Voxels)
+{
+	if(IsMarchingCubesRunning)
+		return;
+
+	IsMarchingCubesRunning = true;
+	
+	FScopeLock Lock(&ProceduralMeshGuard);
+	AsyncTask(ENamedThreads::AnyBackgroundThreadNormalTask, [this, Voxels]()
+	{
+		FVector Min(FLT_MAX, FLT_MAX, FLT_MAX);
+		FVector Max(-FLT_MAX, -FLT_MAX, -FLT_MAX);
+		
+		for(const int FloodedIndex : Voxels)
+		{
+			Min = FVector::Min(Min, MyDensityField->IndexToVoxelPosition(FloodedIndex));
+			Max = FVector::Max(Max, MyDensityField->IndexToVoxelPosition(FloodedIndex));
+		}
+
+		UE::Geometry::FMarchingCubes MarchingCubes;
+		MarchingCubes.Bounds = UE::Geometry::TAxisAlignedBox3(Min, Max);
+		MarchingCubes.CubeSize = MarchingCubeSize == 0 ? MyDensityField->GetStep().X : MarchingCubeSize;
+		MarchingCubes.CancelF = [this]() {
+			return AbortMarchingCubes;
+		};
+		MarchingCubes.IsoValue = 0;
+
+		AbortMarchingCubes = false;
+		// 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)
+		{
+			Vertices.Add(FVector(Vertex.X, Vertex.Y, Vertex.Z));
+		}
+
+		// 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]()
+		{
+			FScopeLock MeshLock(&ProceduralMeshGuard);
+			ProceduralMesh->CreateMeshSection(0, Vertices, Indices, Normals, TArray<FVector2D>(), VertexColors, TArray<FProcMeshTangent>(), false);
+			ProceduralMesh->SetVisibility(true);
+			AccumulatedTime = 0.0f;
+			IsMarchingCubesRunning = false;
+		});
+	});
+}
\ No newline at end of file
diff --git a/Source/MetaCastBachelor/MagicWand.h b/Source/MetaCastBachelor/MagicWand.h
new file mode 100644
index 0000000000000000000000000000000000000000..f20239e223f57bb4b4fa83fa02266b11195266b4
--- /dev/null
+++ b/Source/MetaCastBachelor/MagicWand.h
@@ -0,0 +1,57 @@
+
+#pragma once
+
+#include "CoreMinimal.h"
+#include "ProceduralMeshComponent.h"
+#include "MetaCastBaseline.h"
+#include "Generators/MarchingCubes.h"
+#include "MagicWand.generated.h"
+
+UCLASS( ClassGroup=(Custom), meta=(BlueprintSpawnableComponent) )
+class UMagicWand : public UMetaCastBaseline
+{
+	GENERATED_BODY()
+	
+	FDensityField* MyDensityField;
+	UPROPERTY()
+	UWorld* World;
+	bool IsMarchingCubesRunning = false;
+	bool AbortMarchingCubes = false;
+
+	// Critical section to protect shared variables
+	mutable FCriticalSection ProceduralMeshGuard;
+
+	UPROPERTY()
+	UProceduralMeshComponent* ProceduralMesh;
+	UPROPERTY(EditAnywhere)
+	UMaterialInterface* SelectionVolumeMat;
+	TArray<int32> FloodedIndices;
+	UPROPERTY(EditAnywhere)
+	float MarchingCubeSize = 0;
+	UPROPERTY(EditAnywhere)
+	int EvaluationsPerSecond = 10;
+	float AccumulatedTime = 0.0;
+	UPROPERTY(EditAnywhere)
+	float ProximityThreshold = 0.1f;
+	TArray<int32> SelectedClusterIndices;
+
+	
+public:
+	UMagicWand();
+	
+	virtual void BeginPlay() override;
+		
+	virtual void SelectParticles(const FVector& InputPosition) override;
+	virtual void EraseParticles(const FVector& InputPosition) override;
+
+	virtual void HandleMetaSelectReleased(const FInputActionInstance& Instance) override;
+	virtual void HandleMetaSelectPressed(const FInputActionInstance& Instance) override;
+	void InitMagicWandSelection();
+
+	void GenerateVoxelMeshWithCubes(const TArray<int32> Voxels) const;
+	void GenerateVoxelMeshSmooth(const TArray<int32> Voxels);
+
+	virtual void TickComponent(float DeltaTime, ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction) override;
+	void PerformMagicWandSelection(const FVector& InputPosition);
+	void ExpandSelection(int32 SeedIndex);
+};
\ No newline at end of file
diff --git a/Source/MetaCastBachelor/MetaPoint.cpp b/Source/MetaCastBachelor/MetaPoint.cpp
index 6b45ec3eb10b7998cb43a397de22dac85c24867f..2e908f081531254cac5f4c038fb546e1cd7c5f26 100644
--- a/Source/MetaCastBachelor/MetaPoint.cpp
+++ b/Source/MetaCastBachelor/MetaPoint.cpp
@@ -1,6 +1,7 @@
 #include "MetaPoint.h"
 #include "Utilities.h"
 #include "Generators/MarchingCubes.h"
+#include "UObject/GCObjectScopeGuard.h"
 
 UMetaPoint::UMetaPoint() : LocalMaximumIndex(0), MyDensityField(nullptr), MetaPointThreshold(0), World(nullptr)
 {
@@ -52,14 +53,13 @@ void UMetaPoint::HandleMetaSelectReleased(const FInputActionInstance& Instance)
 {
 	Super::HandleMetaSelectReleased(Instance);
 
-	//ProceduralMesh->ClearAllMeshSections();
-	//ProceduralMesh->SetVisibility(false);
+	ProceduralMesh->ClearAllMeshSections();
+	ProceduralMesh->SetVisibility(false);
 
-	MyPointCloud->ColorPointsInVoxels(FloodedIndices);
 
 	AsyncTask(ENamedThreads::Type::AnyBackgroundHiPriTask, [this]()
 	{
-		
+		MyPointCloud->ColorPointsInVoxels(FloodedIndices);
 	});
 }
 
@@ -70,6 +70,17 @@ void UMetaPoint::HandleMetaSelectPressed(const FInputActionInstance& Instance)
 	InitMetaPointSelection();
 }
 
+void UMetaPoint::HandleMetaEraseReleased(const FInputActionInstance& Instance)
+{
+	Super::HandleMetaEraseReleased(Instance);
+
+	//deselect all particles
+	for (int32 i = 0; i < MyPointCloud->SelectionFlags.Num(); ++i)
+	{
+		MyPointCloud->SelectionFlags[i] = false;
+	}
+}
+
 void UMetaPoint::InitMetaPointSelection()
 {	
 	const FVector SelectionWorldPosition = SelectionObject->GetComponentLocation();
@@ -114,7 +125,7 @@ void UMetaPoint::SelectParticles(const FVector& InputPosition)
 		AccumulatedTime = -10000;
 
 		FScopeLock Lock(&FloodFillingGuard); 
-		AsyncTask(ENamedThreads::AnyBackgroundThreadNormalTask, [this]()
+		AsyncTask(ENamedThreads::Type::AnyBackgroundHiPriTask, [this]()
 		{
 			const FVector SelectionWorldPosition = SelectionObject->GetComponentLocation();
 			
@@ -133,13 +144,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);
+			const float Dist = FVector::Distance(SelectStartPosition, SelectionLocalPosition);
+			UE_LOG(LogTemp, Warning, TEXT("Dist: %f"), Dist);
+			float InitThreshold = FUtilities::InterpolateDensityAtPosition(MyDensityField, SelectStartPosition);
+			float Thr = MyDensityField->MaxDensity - (Dist / 100.0f) * (MyDensityField->MaxDensity - MyDensityField->MinDensity) - InitThreshold;
 
-			const float MaxDistance = FVector::Distance(SelectStartPosition, SelectionLocalPosition);
-			FloodedIndices = FUtilities::FloodFilling_2(MyDensityField, LocalMaximumIndex, MetaPointThreshold, MaxDistance);
+			MetaPointThreshold = FUtilities::InterpolateDensityAtPosition(MyDensityField, SelectionLocalPosition);
 			
+			const float MaxDistance = FVector::Distance(SelectStartPosition, SelectionLocalPosition);
+			//FloodedIndices = FUtilities::FloodFilling_2(MyDensityField, LocalMaximumIndex, MetaPointThreshold, MaxDistance);
+			FloodedIndices = FUtilities::FloodFilling(MyDensityField, LocalMaximumIndex, MetaPointThreshold);
+
 			GenerateVoxelMeshSmooth(FloodedIndices);
 		});
 	}
@@ -197,83 +214,106 @@ void UMetaPoint::GenerateVoxelMeshSmooth(const TArray<int32> Voxels)
 	if(IsMarchingCubesRunning)
 		return;
 
+	if(Voxels.IsEmpty())
+	{
+		AccumulatedTime = 0;
+		ProceduralMesh->ClearAllMeshSections();
+		ProceduralMesh->SetVisibility(false);
+		return;
+	}
+
 	IsMarchingCubesRunning = true;
 	
-	FScopeLock Lock(&ProceduralMeshGuard);
-	AsyncTask(ENamedThreads::AnyBackgroundThreadNormalTask, [this, Voxels]()
+	FVector Min(FLT_MAX, FLT_MAX, FLT_MAX);
+	FVector Max(-FLT_MAX, -FLT_MAX, -FLT_MAX);
+	
+	for(const int FloodedIndex : Voxels)
 	{
-		FVector Min(FLT_MAX, FLT_MAX, FLT_MAX);
-		FVector Max(-FLT_MAX, -FLT_MAX, -FLT_MAX);
-		
-		for(const int FloodedIndex : Voxels)
-		{
-			Min = FVector::Min(Min, MyDensityField->IndexToVoxelPosition(FloodedIndex));
-			Max = FVector::Max(Max, MyDensityField->IndexToVoxelPosition(FloodedIndex));
-		}
+		Min = FVector::Min(Min, MyDensityField->IndexToVoxelPosition(FloodedIndex));
+		Max = FVector::Max(Max, MyDensityField->IndexToVoxelPosition(FloodedIndex));
+	}
+	
+	UE::Geometry::FMarchingCubes MarchingCubes;
+	MarchingCubes.Bounds = UE::Geometry::TAxisAlignedBox3(Min, Max);
+	//MarchingCubes.CubeSize = MarchingCubeSize == 0 ? MyDensityField->GetStep().X : MarchingCubeSize;
 
-		UE::Geometry::FMarchingCubes MarchingCubes;
-		MarchingCubes.Bounds = UE::Geometry::TAxisAlignedBox3(Min, Max);
-		MarchingCubes.CubeSize = MarchingCubeSize == 0 ? MyDensityField->GetStep().X : MarchingCubeSize;
-		MarchingCubes.CancelF = [this]() {
-			return AbortMarchingCubes;
-		};
-		MarchingCubes.IsoValue = 0;
-
-		AbortMarchingCubes = false;
-		// 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)
-		{
-			Vertices.Add(FVector(Vertex.X, Vertex.Y, Vertex.Z));
-		}
+	const float Volume = MarchingCubes.Bounds.Volume();
 
-		// 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);
-		}
+	// Logarithmic scaling of the cube size based on volume
+	if(Volume > 50000) {
+		MarchingCubes.CubeSize = (3.0f * log10((Volume + 500) / 1000) - 3);
+	}else if(Volume > 500000)
+	{
+		MarchingCubes.CubeSize = 9;
+	}else {
+		MarchingCubes.CubeSize = MyDensityField->GetStep().X;
+	}
 
-		// 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));
-		}
+	MarchingCubes.CubeSize = FMath::Clamp(MarchingCubes.CubeSize, MyDensityField->GetStep().X, 10.0f);
 
-		TArray<FColor> VertexColors;
-		VertexColors.Reserve(Vertices.Num());
-		for (int32 i = 0; i < Vertices.Num(); ++i)
-		{
-			VertexColors.Add(FColor::Green);
-		}
+	MarchingCubes.CancelF = [this]() {
+		return AbortMarchingCubes;
+	};
+	MarchingCubes.IsoValue = 0;
+	MarchingCubes.bEnableValueCaching = true;
+	AbortMarchingCubes = false;
 
+	// 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;
+	};
 
-		AsyncTask(ENamedThreads::Type::GameThread, [this, Vertices, Indices, Normals, VertexColors]()
-		{
-			FScopeLock MeshLock(&ProceduralMeshGuard);
-			ProceduralMesh->CreateMeshSection(0, Vertices, Indices, Normals, TArray<FVector2D>(), VertexColors, TArray<FProcMeshTangent>(), false);
-			ProceduralMesh->SetVisibility(true);
-			AccumulatedTime = 0.0f;
-			IsMarchingCubesRunning = false;
-		});
+	MarchingCubes.bParallelCompute = true;
+	MarchingCubes.RootMode = UE::Geometry::ERootfindingModes::Bisection;
+
+	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)
+	{
+		Vertices.Add(FVector(Vertex.X, Vertex.Y, Vertex.Z));
+	}
+
+	// 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]()
+	{
+		FScopeLock MeshLock(&ProceduralMeshGuard);
+		ProceduralMesh->CreateMeshSection(0, Vertices, Indices, Normals, TArray<FVector2D>(), VertexColors, TArray<FProcMeshTangent>(), false);
+		ProceduralMesh->SetVisibility(true);
+		AccumulatedTime = 0.0f;
+		IsMarchingCubesRunning = false;
 	});
 }
\ No newline at end of file
diff --git a/Source/MetaCastBachelor/MetaPoint.h b/Source/MetaCastBachelor/MetaPoint.h
index 20cd30ed39349630e349e66a2d68736b6683695b..8a8591ce298e915c093ae347ef8a4651ffc50dcd 100644
--- a/Source/MetaCastBachelor/MetaPoint.h
+++ b/Source/MetaCastBachelor/MetaPoint.h
@@ -58,6 +58,7 @@ public:
 
 	virtual void HandleMetaSelectReleased(const FInputActionInstance& Instance) override;
 	virtual void HandleMetaSelectPressed(const FInputActionInstance& Instance) override;
+	virtual void HandleMetaEraseReleased(const FInputActionInstance& Instance) override;
 	void InitMetaPointSelection();
 
 	void GenerateVoxelMeshWithCubes(const TArray<int32> Voxels) const;
diff --git a/Source/MetaCastBachelor/PointCloud.cpp b/Source/MetaCastBachelor/PointCloud.cpp
index aad4992d3ac917b71487bc944061f228b45301f8..d319ce91f39e29ecdb9443e723229bac0bd796cc 100644
--- a/Source/MetaCastBachelor/PointCloud.cpp
+++ b/Source/MetaCastBachelor/PointCloud.cpp
@@ -134,7 +134,7 @@ void APointCloud::DrawVoxel(const int Index, const float Time) const
 
 }
 
-void APointCloud::ColorPointsInVoxels(const TArray<int32>& VoxelIndices)
+void APointCloud::ColorPointsInVoxels(const TArray<int32> VoxelIndices)
 {
 	FScopeLock Lock(&DataGuard);
 	if (!MyDensityField)
diff --git a/Source/MetaCastBachelor/PointCloud.h b/Source/MetaCastBachelor/PointCloud.h
index 90e89d217240ff7301b4135764d0527b2c1df935..9de2a628025c9b51483411f91b94d78166a0cca7 100644
--- a/Source/MetaCastBachelor/PointCloud.h
+++ b/Source/MetaCastBachelor/PointCloud.h
@@ -67,7 +67,7 @@ public:
 
 	void UpdateSelection();
 	void DrawVoxel(const int Index, float Time) const;
-	void ColorPointsInVoxels(const TArray<int32>& VoxelIndices);
+	void ColorPointsInVoxels(const TArray<int32> VoxelIndices);
 
 	UFUNCTION(BlueprintCallable)
 	void ReadPointCloudFromFile(FFilePath FileNamePoints, FFilePath FileNameFlags);
diff --git a/Source/MetaCastBachelor/Utilities.cpp b/Source/MetaCastBachelor/Utilities.cpp
index 76b42e35cda7a3a9d9cb502a09d6f3188a2160dd..c3581b95a017e689cfede3e732ec0fff2995f108 100644
--- a/Source/MetaCastBachelor/Utilities.cpp
+++ b/Source/MetaCastBachelor/Utilities.cpp
@@ -106,7 +106,7 @@ float FUtilities::GaussianKernel(const float Distance, const float 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;
+	return (Value * Value) / 5;
 }
 
 FVector FUtilities::FollowGradientToMaximum(const FDensityField* DensityField, const FVector& StartPosition)
@@ -116,7 +116,7 @@ FVector FUtilities::FollowGradientToMaximum(const FDensityField* DensityField, c
 
 	int Iterations = 0;
 	constexpr int MaxIterations = 100;  // Prevent infinite loops
-	constexpr float StepSize = 0.05f;  // Define a reasonable step size
+	constexpr float StepSize = 0.05f;  // Define step size
 
 	while (Iterations < MaxIterations)
 	{