Skip to content
Snippets Groups Projects
Select Git revision
  • 0fdf6ec4994be6d9b5f78a06508c1747d79a584e
  • main default protected
  • release
3 results

pandas.ipynb

Blame
    • Utilities.cpp 9.07 KiB 
    #include "Utilities.h"
#include "DensityField.h" 

void FUtilities::QuickSort(TArray<float>& Array, const int32 Begin, const int32 End)
{
	if (Begin < End)
	{
		const float Pivot = Array[Begin];
		int32 Left = Begin;
		int32 Right = End;

		while (Left < Right)
		{
			while (Left < Right && Array[Right] > Pivot)
				Right--;
			Array[Left] = Array[Right];

			while (Left < Right && Array[Left] <= Pivot)
				Left++;
			Array[Right] = Array[Left];
		}

		Array[Left] = Pivot;
		QuickSort(Array, Begin, Left - 1);
		QuickSort(Array, Left + 1, End);
	}
}

double FUtilities::InterpolateDensityAtPosition(const FDensityField* DensityField, const FVector& Position)
{
    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);

    // Trilinear interpolation
    const double DX00 = FMath::Lerp(D000, D100, XFraction);
    const double DX01 = FMath::Lerp(D001, D101, XFraction);
    const double DX10 = FMath::Lerp(D010, D110, XFraction);
    const double DX11 = FMath::Lerp(D011, D111, XFraction);

    const double DY0 = FMath::Lerp(DX00, DX10, YFraction);
    const double DY1 = FMath::Lerp(DX01, DX11, YFraction);

    const double DZ = FMath::Lerp(DY0, DY1, ZFraction);

    return DZ;
}

FVector FUtilities::InterpolateGradientAtPosition(const FDensityField* DensityField, const FVector& Position)
{
	if (!DensityField)
		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);

	// Fractional parts to interpolate
	float XFraction = GridPos.X - XBin;
	float YFraction = GridPos.Y - YBin;
	float ZFraction = GridPos.Z - ZBin;

	// Fetch gradients from corners of the cell
	FVector G000 = DensityField->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);

	// Interpolate along x for each of the yz pairs
	FVector G00 = FMath::Lerp(G000, G100, XFraction);
	FVector G01 = FMath::Lerp(G001, G101, XFraction);
	FVector G10 = FMath::Lerp(G010, G110, XFraction);
	FVector G11 = FMath::Lerp(G011, G111, XFraction);

	// Interpolate along y
	FVector G0 = FMath::Lerp(G00, G10, YFraction);
	FVector G1 = FMath::Lerp(G01, G11, YFraction);

	// Final interpolation along z
	FVector Gradient = FMath::Lerp(G0, G1, ZFraction);

	return Gradient;
}

float FUtilities::GaussianKernel(const float Distance, const float Sigma) {
	return FMath::Exp(-0.5 * FMath::Square(Distance / Sigma));
}

FVector FUtilities::FollowGradientToMaximum(const FDensityField* DensityField, const FVector& StartPosition)
{
	FVector CurrentPosition = StartPosition;
	double CurrentDensity = InterpolateDensityAtPosition(DensityField, CurrentPosition);

	int Iterations = 0;
	constexpr int MaxIterations = 100;  // Prevent infinite loops
	constexpr float StepSize = 0.1f;  // 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());
			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());
			break;  // Next position is out of valid bounds
		}

		const double NewDensity = InterpolateDensityAtPosition(DensityField, NextPosition);
		//UE_LOG(LogTemp, Log, TEXT("Iteration %d: CurrentPos = %s, Gradient = %s, NextPos = %s, NewDensity = %f, CurrentDensity = %f"), Iterations, *CurrentPosition.ToString(), *Gradient.ToString(), *NextPosition.ToString(), NewDensity, CurrentDensity);

		if (NewDensity <= CurrentDensity) {  // Check if density has increased
    			//UE_LOG(LogTemp, Log, TEXT("Density did not increase; stopping at position %s with density %f"), *CurrentPosition.ToString(), CurrentDensity);
			break;  // No improvement in density, stop iteration
		}

		CurrentPosition = NextPosition;  // Update position
		CurrentDensity = NewDensity;  // Update current density
		Iterations++;
	}

	return CurrentPosition;
}

TArray<int32> FUtilities::FloodFilling(const FDensityField* DensityField, const int32 StartIndex, const float Threshold)
{
	//UE_LOG(LogTemp, Warning, TEXT("Threshold for FloodFilling: %f"), Threshold);
	
	TArray<int32> VisitedIndices;
	if (!DensityField) {
		UE_LOG(LogTemp, Warning, TEXT("DensityField is null."));
		return VisitedIndices;
	}

	if (!DensityField->IsValidIndex(StartIndex)) {
		UE_LOG(LogTemp, Warning, TEXT("Start index %d is invalid."), StartIndex);
		return VisitedIndices;
	}

	TArray<int32> Stack;
	Stack.Push(StartIndex);
	TArray<bool> Visited;
	Visited.Init(false, DensityField->GetVoxelNumber());

	while (Stack.Num() > 0)
	{
		int32 CurrentIndex = Stack.Pop();
		if (!Visited[CurrentIndex])
		{
			Visited[CurrentIndex] = true;
			VisitedIndices.Add(CurrentIndex);

			TArray<int32> Neighbors = DensityField->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("Lookign at Neighbor %d at Position %s with density: %f."), NeighborIndex, *DensityField->IndexToGridPosition(NeighborIndex).ToString(), NeighborDensity);
				
				if (NeighborDensity > Threshold && !Visited[NeighborIndex]) {
					Stack.Push(NeighborIndex);
					//UE_LOG(LogTemp, Log, TEXT("Pushing Neighbor %d to stack."), NeighborIndex);
				}
			}
		}
	}

	// Logging final details
	//UE_LOG(LogTemp, Log, TEXT("Flood filling completed from start index %d. Total voxels selected: %d"), StartIndex, VisitedIndices.Num());
	for (const int32 idx : VisitedIndices) {
		//UE_LOG(LogTemp, Log, TEXT("Selected Voxel Index: %d at Grid Position: %s"), idx, *DensityField->IndexToGridPosition(idx).ToString());
	}

	return VisitedIndices;
}

void FUtilities::FloodFilling_2(const FDensityField* DensityField, const float Threshold, TArray<bool>& VisitedIndices, TArray<int32>& IndicesToVisit, TArray<int32>& FloodedIndices, TArray<int32>& PotentialRevisit)
{
	if (!DensityField) {
		UE_LOG(LogTemp, Warning, TEXT("DensityField is null."));
		return;
    	}

	TArray<int32> CurrentPotentialRevisit;

	
	constexpr int MaxIterations = 2000;
	int IterationCount = 0;
	
	while (IndicesToVisit.Num() > 0)
	{
		
		IterationCount++;
		if(IterationCount > MaxIterations)
		{
			while (IndicesToVisit.Num() > 0) {
				CurrentPotentialRevisit.Add(IndicesToVisit.Pop());
			}
			break;
		}
		
		int32 CurrentIndex = IndicesToVisit.Pop();

		if(DensityField->IsValidIndex(CurrentIndex))
		{
			VisitedIndices[CurrentIndex] = true;
			const double CurrentDensity = DensityField->IndexToVoxelDensity(CurrentIndex);

			if (CurrentDensity >= Threshold) {
				FloodedIndices.Add(CurrentIndex);

				TArray<int32> Neighbors = DensityField->IndexToVoxelNeighbors(CurrentIndex);
				for (int32 NeighborIndex : Neighbors)
				{
					if (!VisitedIndices[NeighborIndex])
					{
						IndicesToVisit.Push(NeighborIndex);
					}
				}
			} else {
				CurrentPotentialRevisit.Add(CurrentIndex);
			}
		}
	}

	// Update PotentialRevisit with the new list from this run
	PotentialRevisit = CurrentPotentialRevisit;
}