diff --git a/Content/MetaPointMap.umap b/Content/MetaPointMap.umap
index 651dc1767d85ce66c1085a001440e676d242a400..f44a6fa3bd8788f0813c8f4c4351e9c5bf22240f 100644
Binary files a/Content/MetaPointMap.umap and b/Content/MetaPointMap.umap differ
diff --git a/Content/SelectionVolume.uasset b/Content/SelectionVolume.uasset
index 73d06d7c6830ac1fc883baf71e935610588d84e9..b4a7c48878aadac755046e2d70262c2c42fa0168 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 3aa6e9bcc1397e4cdfae5b1ab8fb2388fe6895b7..f95293968bea14ae7527061b0f7c92618d168060 100644
--- a/Source/MetaCastBachelor/DensityField.cpp
+++ b/Source/MetaCastBachelor/DensityField.cpp
@@ -104,6 +104,72 @@ void FDensityField::CalculateVoxelDensities(const APointCloud* PointCloud, const
 
 }
 
+void FDensityField::CalculateVoxelDensitiesByNumber(const APointCloud* PointCloud, const float InfluenceRadius)
+{
+    if (!PointCloud)
+        return;
+
+    const double StartTime = FPlatformTime::Seconds();  // Start timing
+    UE_LOG(LogTemp, Log, TEXT("Starting density calculation."));
+
+    // Clear existing densities
+    for (auto& Node : VoxelList)
+    {
+        Node.SetVoxelDensity(0.0);  // Reset the density counter
+    }
+
+    UE_LOG(LogTemp, Log, TEXT("Cleared previous densities."));
+
+    // Iterate over each particle
+    for (const FVector& ParticlePosition : PointCloud->PositionVectors)
+    {
+        // Calculate which voxels the particle influences
+        const int32 StartX = FMath::Max(0, static_cast<int32>((ParticlePosition.X - InfluenceRadius - GridOrigin.X) / XStep));
+        const int32 EndX = FMath::Min(XNum - 1, static_cast<int32>((ParticlePosition.X + InfluenceRadius - GridOrigin.X) / XStep));
+        const int32 StartY = FMath::Max(0, static_cast<int32>((ParticlePosition.Y - InfluenceRadius - GridOrigin.Y) / YStep));
+        const int32 EndY = FMath::Min(YNum - 1, static_cast<int32>((ParticlePosition.Y + InfluenceRadius - GridOrigin.Y) / YStep));
+        const int32 StartZ = FMath::Max(0, static_cast<int32>((ParticlePosition.Z - InfluenceRadius - GridOrigin.Z) / ZStep));
+        const int32 EndZ = FMath::Min(ZNum - 1, static_cast<int32>((ParticlePosition.Z + InfluenceRadius - GridOrigin.Z) / ZStep));
+
+        // Update densities within the influence radius
+        for (int32 Z = StartZ; Z <= EndZ; ++Z)
+        {
+            for (int32 Y = StartY; Y <= EndY; ++Y)
+            {
+                for (int32 X = StartX; X <= EndX; ++X)
+                {
+                    const int32 Index = GridPositionToIndex(X, Y, Z);
+                    if (IsValidIndex(Index))
+                    {
+                        FVector VoxelPosition = IndexToVoxelPosition(Index);
+                        const float Distance = FVector::Dist(VoxelPosition, ParticlePosition);
+
+                        if (Distance < InfluenceRadius)
+                        {
+                            VoxelList[Index].AddToVoxelDensity(1.0);  // Simply increment the density count
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    double MaxDensity = 0.0;
+    for (const auto& Node : VoxelList)
+    {
+        if (Node.GetVoxelDensity() > MaxDensity)
+        {
+            MaxDensity = Node.GetVoxelDensity();
+        }
+    }
+    UE_LOG(LogTemp, Log, TEXT("Maximum density found: %f"), MaxDensity);
+
+	const double EndTime = FPlatformTime::Seconds();  // End timing
+	const double ElapsedTime = EndTime - StartTime;  // Calculate elapsed time
+
+	UE_LOG(LogTemp, Log, TEXT("Density calculation completed in %f seconds."), ElapsedTime);
+}
+
 void FDensityField::CalculateAndStoreGradients()
 {
 	UE_LOG(LogTemp, Log, TEXT("Starting gradient calculation."));
@@ -182,6 +248,8 @@ void FDensityField::DrawDebugVoxelDensity(const UWorld* World, const AActor* Act
 	UE_LOG(LogTemp, Log, TEXT("DensityField drawn with gradient colors using logarithmic scaling!"));
 }
 
+
+
 //	WORLD POSITION CONVERSION FUNCTIONS
 
 int32 FDensityField::WorldPositionToIndex(const FVector &Position) const {
@@ -354,6 +422,21 @@ double FDensityField::IndexToVoxelDensity(const int32 Index) const
 	return 0.0;
 }
 
+void FDensityField::IndexToVoxelBounds(const int32 Index, FVector& OutMinBounds, FVector& OutMaxBounds) const
+{
+	if (!IsValidIndex(Index)) {
+		UE_LOG(LogTemp, Warning, TEXT("Invalid voxel index provided: %d"), Index);
+		return;
+	}
+
+	// Calculate the bounds
+	const FVector VoxelCenter = IndexToVoxelPosition(Index);
+	OutMinBounds = VoxelCenter - GetStep() * 0.5f;
+	OutMaxBounds = VoxelCenter + GetStep() * 0.5f;
+
+	//UE_LOG(LogTemp, Log, TEXT("Voxel %d Bounds: Min(%s), Max(%s)"), Index, *OutMinBounds.ToString(), *OutMaxBounds.ToString());
+}
+
 //	VALIDITY FUNCTIONS
 
 bool FDensityField::IsValidIndex(const int32 Index) const {
diff --git a/Source/MetaCastBachelor/DensityField.h b/Source/MetaCastBachelor/DensityField.h
index 7477fd5c318dc75b9d8a4a11a527eaaaf6160fb5..3f4ff1008b2ee812e351e5e8972346fda050215f 100644
--- a/Source/MetaCastBachelor/DensityField.h
+++ b/Source/MetaCastBachelor/DensityField.h
@@ -60,12 +60,14 @@ class FDensityField
 	FVector GridOrigin;
 	
 public:
+	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();
 
 	//	DEBUG FUNCTIONS
@@ -85,6 +87,7 @@ public:
 	FIntVector3 IndexToGridPosition(int32 Index) const;
 	int32 GetVoxelNumber() const;
     double IndexToVoxelDensity(int32 Index) const;
+	void IndexToVoxelBounds(int32 Index, FVector& OutMinBounds, FVector& OutMaxBounds) const;
 	double GridPositionToVoxelDensity(int32 XBin, int32 YBin, int32 ZBin) const;
 	FVector GridPositionToVoxelGradient(int32 XBin, int32 YBin, int32 ZBin) const;
 	FVector GetStep() const;
@@ -94,6 +97,7 @@ 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/MetaCastBaseline.cpp b/Source/MetaCastBachelor/MetaCastBaseline.cpp
index 7d19389ca4d2712b00351abeee29e81d02e0d1d3..f6bd224b4e055da59cc4b63d22e578552b94d58a 100644
--- a/Source/MetaCastBachelor/MetaCastBaseline.cpp
+++ b/Source/MetaCastBachelor/MetaCastBaseline.cpp
@@ -99,6 +99,7 @@ void UMetaCastBaseline::InitInputBindings()
 		
 		EnhancedInputComponent->BindAction(MetaSelectAction, ETriggerEvent::Started, this, &UMetaCastBaseline::HandleMetaSelectPressed);
 		EnhancedInputComponent->BindAction(MetaSelectAction, ETriggerEvent::Completed, this, &UMetaCastBaseline::HandleMetaSelectReleased);
+		EnhancedInputComponent->BindAction(MetaSelectAction, ETriggerEvent::Canceled, this, &UMetaCastBaseline::HandleMetaSelectReleased);
 		
 		EnhancedInputComponent->BindAction(MetaEraseAction, ETriggerEvent::Started, this, &UMetaCastBaseline::HandleMetaErasePressed);
 		EnhancedInputComponent->BindAction(MetaEraseAction, ETriggerEvent::Completed, this, &UMetaCastBaseline::HandleMetaEraseReleased);
diff --git a/Source/MetaCastBachelor/MetaPoint.cpp b/Source/MetaCastBachelor/MetaPoint.cpp
index dd70f5d1ad2ca193c7bfb2668b6b6e6dbc197ff4..4a71e325d0e0f39d79a93d784a5eba1efab578e7 100644
--- a/Source/MetaCastBachelor/MetaPoint.cpp
+++ b/Source/MetaCastBachelor/MetaPoint.cpp
@@ -1,7 +1,7 @@
 #include "MetaPoint.h"
 #include "Utilities.h"
 
-UMetaPoint::UMetaPoint() : Index(0), MyDensityField(nullptr), Threshold(0), World(nullptr), MyProceduralMesh(nullptr)
+UMetaPoint::UMetaPoint() : LocalMaximumIndex(0), MyDensityField(nullptr), MetaPointThreshold(0), World(nullptr)
 {
 	// Initialize the Marching Cubes algorithm
 	MyMarchingCubes = new MarchingCubes();
@@ -26,6 +26,13 @@ void UMetaPoint::BeginPlay()
 	ProceduralMesh->SetMaterial(2, SelectionVolumeMat);
 }
 
+void UMetaPoint::TickComponent(float DeltaTime, ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction)
+{
+	Super::TickComponent(DeltaTime, TickType, ThisTickFunction);
+
+	AccumulatedTime += DeltaTime;
+}
+
 void UMetaPoint::EraseParticles(const FVector& InputPosition)
 {
 	
@@ -34,92 +41,83 @@ void UMetaPoint::EraseParticles(const FVector& InputPosition)
 void UMetaPoint::HandleMetaSelectReleased(const FInputActionInstance& Instance)
 {
 	Super::HandleMetaSelectReleased(Instance);
+
+	ProceduralMesh->ClearAllMeshSections();
+	ProceduralMesh->SetVisibility(false);
+
+	AsyncTask(ENamedThreads::Type::AnyBackgroundHiPriTask, [this]()
+	{
+		MyPointCloud->ColorPointsInVoxels(FloodedIndices);
+	});
 }
 
 void UMetaPoint::HandleMetaSelectPressed(const FInputActionInstance& Instance)
 {
 	Super::HandleMetaSelectPressed(Instance);
 
-	FindAndMarkLocalMaximum();
+	InitMetaPointSelection();
 }
 
-// Method to perform gradient ascent to find the local maximum density starting from a given position.
-void UMetaPoint::FindAndMarkLocalMaximum()
+void UMetaPoint::InitMetaPointSelection()
 {	
-	const FVector WorldPosition = SelectionObject->GetComponentLocation();
+	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 StartPosition = MyPointCloud->PointCloudVisualizer->GetComponentTransform().InverseTransformPosition(WorldPosition);
+	const FVector SelectionLocalPosition = MyPointCloud->PointCloudVisualizer->GetComponentTransform().InverseTransformPosition(SelectionWorldPosition);
 
 	// Perform gradient ascent to find the local maximum starting from the converted local position
-	LocalMaximum = FUtilities::FollowGradientToMaximum(MyDensityField, StartPosition);
+	LocalMaximum = FUtilities::FollowGradientToMaximum(MyDensityField, SelectionLocalPosition);
 
 	// Convert the local maximum back to world coordinates
-	Index = MyDensityField->WorldPositionToIndex(LocalMaximum);
+	LocalMaximumIndex = MyDensityField->WorldPositionToIndex(LocalMaximum);
 
-	Threshold = FUtilities::InterpolateDensityAtPosition(MyDensityField, LocalMaximum);
-	TestingThresholdFactor = 1;
+	MetaPointThreshold = FUtilities::InterpolateDensityAtPosition(MyDensityField, SelectionLocalPosition);
 
 	// Initialize Visited array to false for all voxels
 	Visited.Init(false, MyDensityField->GetVoxelNumber());
 
 	// Initialize IndicesToVisit with the starting index
 	IndicesToVisit.Empty();
-	IndicesToVisit.Add(Index);
+	IndicesToVisit.Add(LocalMaximumIndex);
 
 	// Clear FloodedIndices
 	FloodedIndices.Empty();
-	
-	World = GetWorld();
-}
-
-void UMetaPoint::TickComponent(float DeltaTime, ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction)
-{
-	Super::TickComponent(DeltaTime, TickType, ThisTickFunction);
+	RevisitIndices.Empty();
 
-	//TestingThresholdFactor -= DeltaTime * 0.1f;
-
-	AccumulatedTime += DeltaTime;
-
-	constexpr float DecayRate = 0.1f; // Adjust this value to control the rate of decay
-	TestingThresholdFactor *= FMath::Exp(-DecayRate * DeltaTime);
-
-	// Ensure TestingThresholdFactor never goes below a small threshold to prevent it from becoming zero
-	constexpr float MinThreshold = 0.001f; // Adjust as needed
-	if (TestingThresholdFactor < MinThreshold)
-	{
-		TestingThresholdFactor = MinThreshold;
-	}
+	ProceduralMesh->ClearAllMeshSections();
+	World = GetWorld();
 }
 
 void UMetaPoint::SelectParticles(const FVector& InputPosition)
 {
-	//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, 0);
-
-	if (AccumulatedTime >=  0.1f)
+	if (AccumulatedTime >=  1 / MetaCastPerSecond)
 	{
-		AccumulatedTime = -10.0f;
-
+		AccumulatedTime = -1000;
+		
 		AsyncTask(ENamedThreads::AnyBackgroundThreadNormalTask, [this]()
 		{
-			//TArray<int32> Voxels = FUtilities::FloodFilling(MyDensityField, Index, Threshold * TestingThresholdFactor);
-			FUtilities::FloodFilling_2(MyDensityField, Threshold * TestingThresholdFactor, Visited, IndicesToVisit, FloodedIndices, IndicesToVisit);
-
-			auto Voxels = FloodedIndices;
-			AsyncTask(ENamedThreads::GameThread, [this, Voxels]()
-			{
-				GenerateVoxelMesh(FloodedIndices);
-
-				AccumulatedTime = 0.0f;
-			});
+			const FVector SelectionWorldPosition = SelectionObject->GetComponentLocation();
+			// 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);
+
+			// Perform gradient ascent to find the local maximum starting from the converted local position
+			LocalMaximum = FUtilities::FollowGradientToMaximum(MyDensityField, SelectionLocalPosition);
+
+			// Convert the local maximum back to world coordinates
+			LocalMaximumIndex = MyDensityField->WorldPositionToIndex(LocalMaximum);
+
+			MetaPointThreshold = FUtilities::InterpolateDensityAtPosition(MyDensityField, SelectionLocalPosition);
+			FloodedIndices = FUtilities::FloodFilling(MyDensityField, LocalMaximumIndex, MetaPointThreshold);
+			
+			GenerateVoxelMesh(FloodedIndices);
+			
+			AccumulatedTime = 0.0f;
 		});
 	}
 }
 
-void UMetaPoint::GenerateVoxelMesh(const TArray<int32>& Voxels) const
+void UMetaPoint::GenerateVoxelMesh(const TArray<int32> Voxels) const
 {
 	TArray<FVector> Vertices;
 	TArray<int32> Triangles;
@@ -135,12 +133,6 @@ void UMetaPoint::GenerateVoxelMesh(const TArray<int32>& Voxels) const
 
 		TArray<FVector> VoxelCorners = MyDensityField->IndexToVoxelCornersWorld(VoxelIndex);
 
-		// Convert corners to world space
-		for (FVector& Corner : VoxelCorners)
-		{
-			//Corner = MyPointCloud->PointCloudVisualizer->GetComponentTransform().TransformPosition(Corner);
-		}
-
 		// Add vertices for the voxel
 		const int32 BaseIndex = Vertices.Num();
 		Vertices.Append(VoxelCorners);
@@ -159,10 +151,48 @@ void UMetaPoint::GenerateVoxelMesh(const TArray<int32>& Voxels) const
 		// Add red color for each corner vertex
 		for (int32 i = 0; i < 8; ++i)
 		{
-			VertexColors.Add(FColor::Red);
+			VertexColors.Add(FColor::Green);
 		}
 	}
 
 	// Create the mesh section
-	ProceduralMesh->CreateMeshSection(0, Vertices, Triangles, TArray<FVector>(), TArray<FVector2D>(), VertexColors, TArray<FProcMeshTangent>(), true);
+	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 UMetaPoint::CreateMeshSections(const TArray<FVector>& Vertices, const TArray<int32>& Triangles, const TArray<FColor>& VertexColors, const int32 MaxVerticesPerSection) const
+{
+	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)
+	{
+		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)
+		{
+			SectionVertices.Add(Vertices[VertexIndex + i]);
+			SectionVertexColors.Add(VertexColors[VertexIndex + i]);
+		}
+
+		TArray<int32> SectionTriangles;
+		SectionTriangles.Reserve(NumTrianglesInSection * 3);
+		for (int32 i = 0; i < NumTrianglesInSection * 3; ++i)
+		{
+			SectionTriangles.Add(Triangles[TriangleIndex * 3 + i] - VertexIndex);
+		}
+
+		// Create or update the mesh section
+		ProceduralMesh->CreateMeshSection(SectionIndex, SectionVertices, SectionTriangles, TArray<FVector>(), TArray<FVector2D>(), SectionVertexColors, TArray<FProcMeshTangent>(), false);
+	}
+}
\ No newline at end of file
diff --git a/Source/MetaCastBachelor/MetaPoint.h b/Source/MetaCastBachelor/MetaPoint.h
index ce3e88b1466bceed87a1c4f53bb1e05394e40320..e221d596fa68eed248fa363f616f08c7d2d4eb3f 100644
--- a/Source/MetaCastBachelor/MetaPoint.h
+++ b/Source/MetaCastBachelor/MetaPoint.h
@@ -14,10 +14,10 @@ class UMetaPoint : public UMetaCastBaseline
 	
 	MarchingCubes* MyMarchingCubes;
 	FVector WorldMaximum;
-	int32 Index;
+	int32 LocalMaximumIndex;
 	FVector LocalMaximum;
 	FDensityField* MyDensityField;
-	float Threshold;
+	float MetaPointThreshold;
 	UPROPERTY()
 	UWorld* World;
 	float AccumulatedTime = 0.0;
@@ -32,15 +32,11 @@ class UMetaPoint : public UMetaCastBaseline
 	TArray<int32> IndicesToVisit;
 	TArray<int32> FloodedIndices;
 	TArray<int32> RevisitIndices;
-	
-public:
-
-	UPROPERTY(EditAnywhere, BlueprintReadWrite)
-	UProceduralMeshComponent* MyProceduralMesh;
 
 	UPROPERTY(EditAnywhere)
-	float TestingThresholdFactor = 2.0;
+	int MetaCastPerSecond = 10;
 	
+public:
 	UMetaPoint();
 	
 	virtual void BeginPlay() override;
@@ -50,9 +46,10 @@ public:
 
 	virtual void HandleMetaSelectReleased(const FInputActionInstance& Instance) override;
 	virtual void HandleMetaSelectPressed(const FInputActionInstance& Instance) override;
-	void FindAndMarkLocalMaximum();
+	void InitMetaPointSelection();
 
-	void GenerateVoxelMesh(const TArray<int32>& Voxels) const;
+	void GenerateVoxelMesh(const TArray<int32> Voxels) const;
+	void CreateMeshSections(const TArray<FVector>& Vertices, const TArray<int32>& Triangles, const TArray<FColor>& VertexColors, int32 MaxVerticesPerSection) const;
 
 	virtual void TickComponent(float DeltaTime, ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction) override;
 };
diff --git a/Source/MetaCastBachelor/PointCloud.cpp b/Source/MetaCastBachelor/PointCloud.cpp
index d6482df6e1cee64a437fc509e340452d6b98d930..747d186f28458c275465ed59f67bcbad0c023d37 100644
--- a/Source/MetaCastBachelor/PointCloud.cpp
+++ b/Source/MetaCastBachelor/PointCloud.cpp
@@ -22,7 +22,7 @@ void APointCloud::BeginPlay()
 
 void APointCloud::ReadPointCloudFromFile(const FFilePath FileNamePoints, const FFilePath FileNameFlags)
 {
-	TArray<FVector> LoadedPointCloud = PointCloudDataReader::LoadPointCloudData(FileNamePoints.FilePath);
+	const TArray<FVector> LoadedPointCloud = PointCloudDataReader::LoadPointCloudData(FileNamePoints.FilePath);
 	TArray<int32> LoadedPointCloudFlags = PointCloudDataReader::LoadFlags(FileNameFlags.FilePath);
 	
 	// Initialize the boolean array for flags with all elements set to false
@@ -30,7 +30,7 @@ void APointCloud::ReadPointCloudFromFile(const FFilePath FileNamePoints, const F
 	BooleanFlags.Init(false, LoadedPointCloud.Num());
 
 	// Set true for indices that are included in the flags array
-	for (int32 FlagIndex : LoadedPointCloudFlags)
+	for (const int32 FlagIndex : LoadedPointCloudFlags)
 	{
 		if (FlagIndex < BooleanFlags.Num())
 		{
@@ -91,10 +91,10 @@ 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 Sigma = 3.0f;  // Standard deviation for the Gaussian kernel
-	MyDensityField->CalculateVoxelDensities(this, InfluenceRadius, Sigma);
-
+	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())
@@ -107,12 +107,24 @@ void APointCloud::UpdateSelection()
 {
 	for (int32 i = 0; i < PositionVectors.Num(); i++)
 	{
-		if (SelectionFlags[i])
+		if(DefaultFlags[i])
 		{
-			PointColors[i] = FColor::Green;
+			if (SelectionFlags[i])
+			{
+				PointColors[i] = FColor::Red;
+			}else
+			{
+				PointColors[i] = FColor::Orange;
+			}
 		}else
 		{
-			PointColors[i] = FColor::Blue;
+			if (SelectionFlags[i])
+			{
+				PointColors[i] = FColor::Green;
+			}else
+			{
+				PointColors[i] = FColor::Blue;
+			}
 		}
 	}
 	
@@ -135,10 +147,59 @@ void APointCloud::DrawVoxel(const int Index, const float Time) const
 
 }
 
+void APointCloud::ColorPointsInVoxels(const TArray<int32> VoxelIndices)
+{
+	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)
+		{
+			if(!MyPositionVectors.IsValidIndex(i))
+			{
+				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;
+			}
+		}
+	}
+
+	//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 e22fd462812578614a97b6d4284ea46c5acc4348..4770f1b4e6560b2f146abcf5a251f69d33c80ec4 100644
--- a/Source/MetaCastBachelor/PointCloud.h
+++ b/Source/MetaCastBachelor/PointCloud.h
@@ -17,6 +17,7 @@ public:
 	TArray<bool> DefaultFlags;
 	TArray<bool> SelectionFlags;
 	TArray<FColor> PointColors;
+	mutable FCriticalSection DataGuard;
 
 	UPROPERTY(VisibleAnywhere, BlueprintReadOnly, Category = "Bounds")
 	FVector MinBounds;
@@ -53,12 +54,13 @@ protected:
 	virtual void BeginPlay() override;
 	
 
-public:	
+public:
 	// Called every frame
 	virtual void Tick(float DeltaTime) override;
 
 	void UpdateSelection();
 	void DrawVoxel(const int Index, float Time) const;
+	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 343a2ee4f57a50551b4b9088c17277e65e9153ac..236264834ff8020d4e791220713fb65bcd11b198 100644
--- a/Source/MetaCastBachelor/Utilities.cpp
+++ b/Source/MetaCastBachelor/Utilities.cpp
@@ -138,7 +138,7 @@ FVector FUtilities::FollowGradientToMaximum(const FDensityField* DensityField, c
 		//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);
+			//UE_LOG(LogTemp, Log, TEXT("Density did not increase; stopping at position %s with density %f"), *CurrentPosition.ToString(), CurrentDensity);
 			break;  // No improvement in density, stop iteration
 		}
 
@@ -203,38 +203,6 @@ TArray<int32> FUtilities::FloodFilling(const FDensityField* DensityField, const
 	return VisitedIndices;
 }
 
-/*
-TArray<int32> FUtilities::FloodFilling_2(const FDensityField* DensityField, const float Threshold, TArray<bool>& Visited, TArray<int32>& IndicesToVisit, TArray<int32>& FloodedIndices)
-{	
-	if (!DensityField) {
-		UE_LOG(LogTemp, Warning, TEXT("DensityField is null."));
-		return FloodedIndices;
-	}
-
-	while (IndicesToVisit.Num() > 0)
-	{
-		int32 CurrentIndex = IndicesToVisit.Pop();
-		if (!Visited[CurrentIndex])
-		{
-			Visited[CurrentIndex] = true;
-			FloodedIndices.Add(CurrentIndex);
-
-			TArray<int32> Neighbors = DensityField->IndexToVoxelNeighbors(CurrentIndex);
-            
-			for (int32 NeighborIndex : Neighbors)
-			{
-				const double NeighborDensity = DensityField->IndexToVoxelDensity(NeighborIndex);
-				
-				if (NeighborDensity > Threshold && !Visited[NeighborIndex]) {
-					IndicesToVisit.Push(NeighborIndex);
-				}
-			}
-		}
-	}
-
-	return FloodedIndices;
-}*/
-
 void FUtilities::FloodFilling_2(const FDensityField* DensityField, const float Threshold, TArray<bool>& VisitedIndices, TArray<int32>& IndicesToVisit, TArray<int32>& FloodedIndices, TArray<int32>& PotentialRevisit)
 {
 	if (!DensityField) {
@@ -244,10 +212,13 @@ void FUtilities::FloodFilling_2(const FDensityField* DensityField, const float T
 
 	TArray<int32> CurrentPotentialRevisit;
 
-	constexpr int MaxIterations = 100;
+	
+	constexpr int MaxIterations = 2000;
 	int IterationCount = 0;
+	
 	while (IndicesToVisit.Num() > 0)
 	{
+		
 		IterationCount++;
 		if(IterationCount > MaxIterations)
 		{
@@ -258,28 +229,29 @@ void FUtilities::FloodFilling_2(const FDensityField* DensityField, const float T
 		}
 		
 		int32 CurrentIndex = IndicesToVisit.Pop();
-		VisitedIndices[CurrentIndex] = true;
-		const double CurrentDensity = DensityField->IndexToVoxelDensity(CurrentIndex);
 
-		if (CurrentDensity >= Threshold) {
-			FloodedIndices.Add(CurrentIndex);
+		if(DensityField->IsValidIndex(CurrentIndex))
+		{
+			VisitedIndices[CurrentIndex] = true;
+			const double CurrentDensity = DensityField->IndexToVoxelDensity(CurrentIndex);
 
-			TArray<int32> Neighbors = DensityField->IndexToVoxelNeighbors(CurrentIndex);
-			for (int32 NeighborIndex : Neighbors)
-			{
-				if (!VisitedIndices[NeighborIndex])
+			if (CurrentDensity >= Threshold) {
+				FloodedIndices.Add(CurrentIndex);
+
+				TArray<int32> Neighbors = DensityField->IndexToVoxelNeighbors(CurrentIndex);
+				for (int32 NeighborIndex : Neighbors)
 				{
-					IndicesToVisit.Push(NeighborIndex);
+					if (!VisitedIndices[NeighborIndex])
+					{
+						IndicesToVisit.Push(NeighborIndex);
+					}
 				}
+			} else {
+				CurrentPotentialRevisit.Add(CurrentIndex);
 			}
-		} else {
-			CurrentPotentialRevisit.Add(CurrentIndex);
 		}
 	}
 
 	// Update PotentialRevisit with the new list from this run
 	PotentialRevisit = CurrentPotentialRevisit;
 }
-
-
-