Tweaks params for Density map

Source/MetaCastBachelor/DensityField.cpp

@@ -14,14 +14,12 @@ void FDensityField::InitializeDensityField(const APointCloud* PointCloud, const
 {
 	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);
+	
+	CalculateVoxelDensities_2(PointCloud, PointCloud->InfluenceRadius);
+	//CalculateVoxelDensitiesByNumber(PointCloud, 4);
+	CalculateAndStoreGradients();
 }
 
 void FDensityField::InitializeDataStructures(const FVector& MinBounds, const FVector& MaxBounds, const int32 XAxisNum, const int32 YAxisNum, const int32 ZAxisNum)
@@ -185,6 +183,196 @@ void FDensityField::CalculateVoxelDensitiesByNumber(const APointCloud* PointClou
 
 	UE_LOG(LogTemp, Log, TEXT("Density calculation completed in %f seconds."), ElapsedTime);
 }
+/*
+void FDensityField::CalculateVoxelDensities_2(const APointCloud* PointCloud, const float InfluenceRadius)
+{
+    if (!PointCloud)
+        return;
+
+    const double StartTime = FPlatformTime::Seconds();
+    UE_LOG(LogTemp, Log, TEXT("Starting density calculation."));
+
+    TArray<float> SmoothingLengths;
+    SmoothingLengths.Init(InfluenceRadius, PointCloud->PositionVectors.Num());
+
+    // Clear existing densities
+    for (auto& Node : VoxelList)
+    {
+        Node.SetVoxelDensity(0.0);
+        Node.SetClosePointsNumber(0.0);
+    }
+
+    // Iterate over each particle to update smoothing lengths based on local densities
+    UpdateSmoothingLengths(PointCloud, SmoothingLengths, InfluenceRadius);
+
+    // Iterate over each particle for density calculation
+    for (int i = 0; i < PointCloud->PositionVectors.Num(); i++)
+    {
+        const FVector& ParticlePosition = PointCloud->PositionVectors[i];
+        const float ParticleSmoothingLength = SmoothingLengths[i];
+        const float ParticleGaussianDenominator = 2 * ParticleSmoothingLength * ParticleSmoothingLength;
+
+    	const int32 StartX = FMath::Max(0, FMath::FloorToInt((ParticlePosition.X - ParticleSmoothingLength - PointCloud->MinBounds.X) / XStep));
+    	const int32 EndX = FMath::Min(XNum - 1, FMath::FloorToInt((ParticlePosition.X + ParticleSmoothingLength - PointCloud->MinBounds.X) / XStep));
+    	const int32 StartY = FMath::Max(0, FMath::FloorToInt((ParticlePosition.Y - ParticleSmoothingLength - PointCloud->MinBounds.Y) / YStep));
+    	const int32 EndY = FMath::Min(YNum - 1, FMath::FloorToInt((ParticlePosition.Y + ParticleSmoothingLength - PointCloud->MinBounds.Y) / YStep));
+    	const int32 StartZ = FMath::Max(0, FMath::FloorToInt((ParticlePosition.Z - ParticleSmoothingLength - PointCloud->MinBounds.Z) / ZStep));
+    	const int32 EndZ = FMath::Min(ZNum - 1, FMath::FloorToInt((ParticlePosition.Z + ParticleSmoothingLength - PointCloud->MinBounds.Z) / ZStep));
+
+        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 (VoxelList.IsValidIndex(Index))
+                    {
+                        FVector NodePosition = VoxelList[Index].GetVoxelPosition();
+                        const float Distance = FVector::Dist(NodePosition, ParticlePosition);
+                        if (Distance <= ParticleSmoothingLength)
+                        {
+	                        const float Weight = exp(-(Distance * Distance) / ParticleGaussianDenominator);
+                            VoxelList[Index].AddToVoxelDensity(Weight);
+                            VoxelList[Index].SetClosePointsNumber(VoxelList[Index].GetClosePointsNumber() + 1.0);
+                        }
+                    }
+                }
+    }
+
+    double MaxDensity = 0.0;
+    for (const auto& Node : VoxelList)
+        if (Node.GetVoxelDensity() > MaxDensity)
+            MaxDensity = Node.GetVoxelDensity();
+
+    const double EndTime = FPlatformTime::Seconds();
+    const double ElapsedTime = EndTime - StartTime;
+    UE_LOG(LogTemp, Log, TEXT("Density calculation completed in %f seconds."), ElapsedTime);
+}*/
+
+void FDensityField::SetInitialDensityEstimates(const APointCloud* PointCloud, const float InfluenceRadius)
+{
+   // Iterate over each particle
+	for (const FVector& ParticlePosition : PointCloud->PositionVectors)
+	{
+		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)
+	    {
+	        for (int32 Y = StartY; Y <= EndY; ++Y)
+	        {
+	            for (int32 X = StartX; X <= EndX; ++X)
+	            {
+	                const int32 Index = GridPositionToIndex(X, Y, Z);
+	                if (Index >= 0 && Index < VoxelList.Num())
+	                {
+	                    FVector NodePosition = VoxelList[Index].GetVoxelPosition();
+	                    const float Distance = FVector::Dist(NodePosition, ParticlePosition);
+
+	                    if (Distance <= InfluenceRadius)
+	                    {
+		                    const float NormalizedDistance = Distance / InfluenceRadius;
+	                    	
+	                        double Weight =  1- (NormalizedDistance * NormalizedDistance);
+
+	                    	Weight *= 0.5968310365946 / (InfluenceRadius * InfluenceRadius * InfluenceRadius);
+	                        VoxelList[Index].AddToVoxelDensity(Weight);
+	                    }
+	                }
+	            }
+	        }
+	    }
+	}
+}
+
+void FDensityField::CalculateVoxelDensities_2(const APointCloud* PointCloud, const float InfluenceRadius)
+{
+	if (!PointCloud)
+	    return;
+
+	// Clear existing densities and counts
+	for (auto& Node : VoxelList)
+	{
+	    Node.SetVoxelDensity(0.0);
+	    Node.SetClosePointsNumber(0.0);
+	}
+	
+	SetInitialDensityEstimates(PointCloud, InfluenceRadius);
+	
+	TArray<float> SmoothingLengths;
+	SmoothingLengths.Init(InfluenceRadius, PointCloud->PositionVectors.Num());
+
+	// Iterate over each particle
+	int Count = 0;
+	for (const FVector& ParticlePosition : PointCloud->PositionVectors)
+	{
+		const int VoxelIndex = WorldPositionToIndex(ParticlePosition);
+		const float TotalVoxelDensity = IndexToVoxelDensity(VoxelIndex);
+		const int NumberInVoxel = VoxelPointLookupTable->GetPointsInVoxel(VoxelIndex).Num();
+		const float InterpolatedDensity = FUtilities::InterpolateDensityAtPosition(this, ParticlePosition);
+
+		const float SmoothingLength = FMath::Min(FMathf::Pow(TotalVoxelDensity / NumberInVoxel / InterpolatedDensity,(1.0 / 3.0)), 5 * GetStep().X);
+		SmoothingLengths[Count] = SmoothingLength;
+
+		Count++;
+	}
+
+	for (int i = 0; i < 100; i++)
+	{
+		UE_LOG(LogTemp, Log, TEXT("Smoothing Length: %f"), SmoothingLengths[i]);
+	}
+
+	for (auto& Node : VoxelList)
+	{
+		Node.SetVoxelDensity(0.0);
+		Node.SetClosePointsNumber(0.0);
+	}
+	
+	// Iterate over each particle
+	Count = 0;
+	for (const FVector& ParticlePosition : PointCloud->PositionVectors)
+	{
+		const float NewInfluenceRadius = SmoothingLengths[Count];
+		const int32 StartX = FMath::Max(0, FMath::FloorToInt((ParticlePosition.X - NewInfluenceRadius - PointCloud->MinBounds.X) / XStep));
+		const int32 EndX = FMath::Min(XNum - 1, FMath::FloorToInt((ParticlePosition.X + NewInfluenceRadius - PointCloud->MinBounds.X) / XStep));
+		const int32 StartY = FMath::Max(0, FMath::FloorToInt((ParticlePosition.Y - NewInfluenceRadius - PointCloud->MinBounds.Y) / YStep));
+		const int32 EndY = FMath::Min(YNum - 1, FMath::FloorToInt((ParticlePosition.Y + NewInfluenceRadius - PointCloud->MinBounds.Y) / YStep));
+		const int32 StartZ = FMath::Max(0, FMath::FloorToInt((ParticlePosition.Z - NewInfluenceRadius - PointCloud->MinBounds.Z) / ZStep));
+		const int32 EndZ = FMath::Min(ZNum - 1, FMath::FloorToInt((ParticlePosition.Z + NewInfluenceRadius - PointCloud->MinBounds.Z) / ZStep));
+
+		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 (Index >= 0 && Index < VoxelList.Num())
+					{
+						FVector NodePosition = VoxelList[Index].GetVoxelPosition();
+						const float Distance = FVector::Dist(NodePosition, ParticlePosition);
+
+						if (Distance <= NewInfluenceRadius)
+						{
+							const float NormalizedDistance = Distance / NewInfluenceRadius;
+	                    	
+							double Weight =  1 - (NormalizedDistance * NormalizedDistance);
+
+							Weight *= (NewInfluenceRadius * NewInfluenceRadius * NewInfluenceRadius);
+							VoxelList[Index].AddToVoxelDensity(FUtilities::SmoothingKernel(Distance, NewInfluenceRadius));
+						}
+					}
+				}
+			}
+		}
+
+		Count++;
+	}
+}
+
 
 void FDensityField::CalculateAndStoreGradients()
 {

Source/MetaCastBachelor/DensityField.h

@@ -55,6 +55,8 @@ public:
 
 	// INITIALIZATION FUNCTIONS
 	void InitializeDensityField(const APointCloud* PointCloud, const FVector& MinBounds, const FVector& MaxBounds, const FIntVector AxisNumbers);
+	void CalculateVoxelDensities_2(const APointCloud* PointCloud, const float InfluenceRadius);
+	void SetInitialDensityEstimates(const APointCloud* PointCloud, const float InfluenceRadius);
 
 	//	DEBUG FUNCTIONS
     void DrawDebugVoxelDensity(const UWorld* World, const AActor* Actor, float Duration, float Scale, float DensityThreshold) const;

Source/MetaCastBachelor/MetaPoint.cpp

@@ -55,9 +55,11 @@ void UMetaPoint::HandleMetaSelectReleased(const FInputActionInstance& Instance)
 	//ProceduralMesh->ClearAllMeshSections();
 	//ProceduralMesh->SetVisibility(false);
 
+	MyPointCloud->ColorPointsInVoxels(FloodedIndices);
+
 	AsyncTask(ENamedThreads::Type::AnyBackgroundHiPriTask, [this]()
 	{
-		MyPointCloud->ColorPointsInVoxels(FloodedIndices);
+		
 	});
 }
 
@@ -101,14 +103,17 @@ void UMetaPoint::InitMetaPointSelection()
 	World = GetWorld();
 
 	LastHandInput = SelectionWorldPosition;
+	AbortMarchingCubes = true;
+	AccumulatedTime = 0.0f;
 }
 
 void UMetaPoint::SelectParticles(const FVector& InputPosition)
 {
 	if (AccumulatedTime >=  1 / MetaCastPerSecond)
 	{
-		AccumulatedTime = -1000;
+		AccumulatedTime = -10000;
 
+		FScopeLock Lock(&FloodFillingGuard); 
 		AsyncTask(ENamedThreads::AnyBackgroundThreadNormalTask, [this]()
 		{
 			const FVector SelectionWorldPosition = SelectionObject->GetComponentLocation();
@@ -187,8 +192,14 @@ void UMetaPoint::GenerateVoxelMeshWithCubes(const TArray<int32> Voxels) const
 	});
 }
 
-void UMetaPoint::GenerateVoxelMeshSmooth(TArray<int32> Voxels)
+void UMetaPoint::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);
@@ -202,9 +213,13 @@ void UMetaPoint::GenerateVoxelMeshSmooth(TArray<int32> Voxels)
 
 		UE::Geometry::FMarchingCubes MarchingCubes;
 		MarchingCubes.Bounds = UE::Geometry::TAxisAlignedBox3(Min, Max);
-		MarchingCubes.CubeSize = MyDensityField->GetStep().X;
+		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);
@@ -254,12 +269,11 @@ void UMetaPoint::GenerateVoxelMeshSmooth(TArray<int32> Voxels)
 
 		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

Source/MetaCastBachelor/MetaPoint.h

@@ -24,6 +24,12 @@ class UMetaPoint : public UMetaCastBaseline
 	FVector LastHandInput;
 	UPROPERTY(EditAnywhere)
 	float MinimalHandMovement = 0.01;
+	bool IsMarchingCubesRunning = false;
+	bool AbortMarchingCubes = false;
+
+	// Critical section to protect shared variables
+	mutable FCriticalSection FloodFillingGuard;
+	mutable FCriticalSection ProceduralMeshGuard;
 
 	UPROPERTY()
 	UProceduralMeshComponent* ProceduralMesh;
@@ -36,8 +42,11 @@ class UMetaPoint : public UMetaCastBaseline
 	TArray<int32> FloodedIndices;
 	TArray<int32> RevisitIndices;
 
+
 	UPROPERTY(EditAnywhere)
-	int MetaCastPerSecond = 1;
+	float MarchingCubeSize = 0;
+	UPROPERTY(EditAnywhere)
+	int MetaCastPerSecond = 10;
 	
 public:
 	UMetaPoint();
@@ -56,4 +65,3 @@ public:
 
 	virtual void TickComponent(float DeltaTime, ELevelTick TickType, FActorComponentTickFunction* ThisTickFunction) override;
 };
\ No newline at end of file
-

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)

Source/MetaCastBachelor/PointCloud.h

@@ -12,6 +12,9 @@ class APointCloud : public AActor
 	GENERATED_BODY()
 	
 public:
+	UPROPERTY(EditAnywhere)
+	float InfluenceRadius = 3.0f;
+
 	FDensityField* MyDensityField;
 	UPROPERTY()
 	TArray<FVector> PositionVectors;
@@ -64,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);

Source/MetaCastBachelor/Utilities.cpp

@@ -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.005f;  // Define a reasonable step size
+	constexpr float StepSize = 0.05f;  // Define a reasonable step size
 
 	while (Iterations < MaxIterations)
 	{
@@ -250,7 +250,7 @@ TArray<int32> FUtilities::FloodFilling_2(const FDensityField* DensityField, cons
 
 					}else
 					{
-						const float Multiplier = 1 + FMath::Pow(DistanceToStart - MaxDistance, 2);
+						const float Multiplier = 1 + FMath::Pow(DistanceToStart - MaxDistance, 1.1);
 						if(NeighborDensity > Threshold * Multiplier)
 						{
 							Stack.Push(NeighborIndex);