From e3068ee50a0689dc4e50a70bbdb42f075b06b93f Mon Sep 17 00:00:00 2001 From: Francesco Date: Fri, 27 Mar 2026 11:36:42 +0100 Subject: [PATCH] Implemented the superformula and the membrane generator --- src/general/utils/Superformula.cs | 180 +++++++++++ src/general/utils/Superformula.cs.uid | 1 + src/microbe_stage/MembraneShapeGenerator.cs | 287 +++++++++++++++--- src/microbe_stage/editor/SuperformulaPlot.cs | 8 + .../editor/SuperformulaPlot.cs.uid | 1 + 5 files changed, 427 insertions(+), 50 deletions(-) create mode 100644 src/general/utils/Superformula.cs create mode 100644 src/general/utils/Superformula.cs.uid create mode 100644 src/microbe_stage/editor/SuperformulaPlot.cs create mode 100644 src/microbe_stage/editor/SuperformulaPlot.cs.uid diff --git a/src/general/utils/Superformula.cs b/src/general/utils/Superformula.cs new file mode 100644 index 00000000000..75d398e2816 --- /dev/null +++ b/src/general/utils/Superformula.cs @@ -0,0 +1,180 @@ +using System; +using System.Collections.Generic; +using System.Runtime.CompilerServices; +using System.Runtime.InteropServices; + +/// +/// A simplified instance of a cached Superformula evaluator. +/// +public readonly struct Superformula +{ + private const double FiniteDifferenceEpsilon = 10e-4; + private const double FiniteDifferenceEpsilonSquared = FiniteDifferenceEpsilon * FiniteDifferenceEpsilon; + private const double MaxStepSize = double.Pi / 32.0; + private const double MinStepSize = double.Pi / 256.0; + + // The threshold to make the algorithm suspect we have skipped an infinite (or very high) curvature point. + // Maybe we should use a more rigorous value, as this is purely empirical. + private const double SuspiciousStepThreshold = 0.1; + + private readonly (double Rho, double Theta)[] polarCache; + private readonly List<(double X, double Y)> cartesianCache; + + public Superformula(int n1, int n2, int n3, int m) + { + var samples = new List<(double, double)>(); + + double theta = -double.Pi; + double previousTheta = double.NaN; + double previousRho = double.NaN; + while (theta < double.Pi) + { + double rho = EvalDirect(theta, n1, n2, n3, m); + double currentTheta = theta; + + // A bisection method to increment the precision around a high curvature or non-differentiable point. + // This won't yield a good result if in the current interval there are multiple such points, but this case + // should be too rare to even be considered. + while (Math.Abs(rho - previousRho) >= SuspiciousStepThreshold) + { + // We have probably skipped a discontinuity in the first derivative, aka a high curvature point, so we + // want to increase the precision around here. + + double newTheta = (theta + previousTheta) / 2.0; + + if (Math.Abs(theta - previousTheta) < MinStepSize) + { + // We are below the required tolerance, so let's break out of this loop. + // First we need to revert to the theta value before the loop. + theta = currentTheta; + + break; + } + + double newRho = EvalDirect(newTheta, n1, n2, n3, m); + + if (Math.Abs(newRho - previousRho) >= SuspiciousStepThreshold) + { + // We are still past the suspected non-differentiable point, so let's make this the actual current + // theta. + theta = newTheta; + rho = newRho; + + samples.Add((rho, theta)); + + continue; + } + + // The last case is if we are now behind the suspected non-differentiable point, so we make this the + // previous theta. + previousTheta = theta; + previousRho = rho; + + samples.Add((rho, theta)); + } + + samples.Add((rho, theta)); + + // Use the curvature radius as a dynamic step. + double kappa = Curvature(theta, n1, n2, n3, m); + double step = double.IsNaN(kappa) || kappa <= 0.0001 ? MaxStepSize : 1.0 / kappa; + theta += Math.Clamp(step, MinStepSize, MaxStepSize); + + previousRho = rho; + previousTheta = theta; + + if (theta > double.Pi && previousTheta < double.Pi) + { + theta = double.Pi; + } + } + + // Ensure all the samples are ordered wrt theta. + samples.Sort((a, b) => a.Item2.CompareTo(b.Item2)); + + polarCache = samples.ToArray(); + + // Convert to cartesian in-place and accumulate curve length. + double dx, dy; + for (int i = 0; i < samples.Count; ++i) + { + (double Rho, double Theta) polar = samples[i]; + samples[i] = PolarToCartesian(polar.Rho, polar.Theta); + + if (i <= 0) + continue; + + (double X, double Y) currentPoint = samples[i]; + (double X, double Y) previousPoint = samples[i - 1]; + + dx = currentPoint.X - previousPoint.X; + dy = currentPoint.Y - previousPoint.Y; + + Length += Math.Sqrt(dx * dx + dy * dy); + } + + (double X, double Y) firstPoint = samples[0]; + (double X, double Y) lastPoint = samples[^1]; + + dx = firstPoint.X - lastPoint.X; + dy = firstPoint.Y - lastPoint.Y; + + Length += Math.Sqrt(dx * dx + dy * dy); + + cartesianCache = samples; + } + + public ReadOnlySpan<(double X, double Y)> CartesianData => CollectionsMarshal.AsSpan(cartesianCache); + public ReadOnlySpan<(double Rho, double Theta)> PolarData => polarCache; + + public double Length { get; } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static double EvalDirect(double theta, int n1, int n2, int n3, int m) + { + double cos = Math.Cos(theta * m / 4.0); + double sin = Math.Sin(theta * m / 4.0); + + cos = Math.Pow(Math.Abs(cos), n2); + sin = Math.Pow(Math.Abs(sin), n3); + + double value = Math.Pow(cos + sin, 1 / (double)n1); + + return 1 / value; + } + + /// + /// Finite difference method to numerically approximate curvature. + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static double Curvature(double theta, int n1, int n2, int n3, int m) + { + double radiusCenter = EvalDirect(theta, n1, n2, n3, m); + double radiusPlus = EvalDirect(theta + FiniteDifferenceEpsilon, n1, n2, n3, m); + double radiusMinus = EvalDirect(theta - FiniteDifferenceEpsilon, n1, n2, n3, m); + + double firstDerivative = (radiusPlus - radiusMinus) / (2 * FiniteDifferenceEpsilon); + double secondDerivative = (radiusPlus - 2 * radiusCenter + radiusMinus) / FiniteDifferenceEpsilonSquared; + + double radiusCenterSquared = radiusCenter * radiusCenter; + double firstDerivativeSquared = firstDerivative * firstDerivative; + + double curvatureNumerator = Math.Abs(radiusCenterSquared + 2 * firstDerivativeSquared - + radiusCenter * secondDerivative); + double curvatureDenominator = Math.Pow(radiusCenterSquared + firstDerivativeSquared, 1.5); + + if (curvatureDenominator < FiniteDifferenceEpsilon) + return double.PositiveInfinity; + + return curvatureNumerator / curvatureDenominator; + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private static (float X, float Y) PolarToCartesian(double rho, double theta) + { + double x = rho * Math.Cos(theta); + double y = rho * Math.Sin(theta); + + return ((float)x, (float)y); + } +} diff --git a/src/general/utils/Superformula.cs.uid b/src/general/utils/Superformula.cs.uid new file mode 100644 index 00000000000..c39c048d30b --- /dev/null +++ b/src/general/utils/Superformula.cs.uid @@ -0,0 +1 @@ +uid://cddqq431ehc5g diff --git a/src/microbe_stage/MembraneShapeGenerator.cs b/src/microbe_stage/MembraneShapeGenerator.cs index 925a22869af..420281dbc03 100644 --- a/src/microbe_stage/MembraneShapeGenerator.cs +++ b/src/microbe_stage/MembraneShapeGenerator.cs @@ -1,5 +1,7 @@ using System; +using System.Buffers; using System.Collections.Generic; +using System.Runtime.CompilerServices; using System.Threading; using Godot; using Array = Godot.Collections.Array; @@ -36,6 +38,19 @@ public class MembraneShapeGenerator /// private readonly List startingBuffer = new(); + /// + /// How sharply the gate for the superformula blending is. It is the reciprocal of the smooth union operation + /// factor. + /// + private float superformulaEnvelopeSharpness = 5.0f; + + /// + /// Scale factor applied to every superformula radius sample. + /// + private float superformulaEnvelopeScale = 4.0f; + + public Superformula? SuperformulaEnvelope { get; private set; } = new Superformula(1, 4, 3, 5); + /// /// Gets a generator for the current thread. This is required to be used as the generators are not thread safe. /// @@ -45,6 +60,25 @@ public static MembraneShapeGenerator GetThreadSpecificGenerator() return ThreadLocalGenerator.Value!; } + /// + /// Configures superformula parameters. + /// + public void ConfigureSuperformula(int n1, int n2, int n3, int m, float sharpness = 5.0f, float scale = 4.0f) + { + SuperformulaEnvelope = new Superformula(n1, n2, n3, m); + + superformulaEnvelopeSharpness = sharpness; + superformulaEnvelopeScale = scale; + } + + /// + /// Disables the superformula envelope. + /// + public void ClearSuperformula() + { + SuperformulaEnvelope = null; + } + /// /// Generates the 2D points for a membrane given the parameters that affect the shape /// @@ -54,61 +88,18 @@ public static MembraneShapeGenerator GetThreadSpecificGenerator() /// public MembranePointData GenerateShape(Vector2[] hexPositions, int hexCount, MembraneType membraneType) { - // The length in pixels (probably not accurate?) of a side of the square that bounds the membrane. - // Half the side length of the original square that is compressed to make the membrane. - int cellDimensions = 10; + bool useSuperformula = UsesSuperformulaEnvelope(membraneType); - for (int i = 0; i < hexCount; ++i) + if (useSuperformula) { - var pos = hexPositions[i]; - if (MathF.Abs(pos.X) + 1 > cellDimensions) - { - cellDimensions = (int)MathF.Abs(pos.X) + 1; - } - - if (MathF.Abs(pos.Y) + 1 > cellDimensions) - { - cellDimensions = (int)MathF.Abs(pos.Y) + 1; - } + GenerateSuperformulaEnvelope(hexPositions, hexCount); } - - // Make the length longer to guarantee that everything fits easily inside the square - cellDimensions *= 100; - - startingBuffer.Clear(); - - // Integer divides are intentional here - // ReSharper disable PossibleLossOfFraction - - for (int i = MembraneResolution; i > 0; i--) + else { - startingBuffer.Add(new Vector2(-cellDimensions, - cellDimensions - 2 * cellDimensions / MembraneResolution * i)); + // This uses the original generation pipeline + GenerateBlobEnvelope(hexPositions, hexCount, membraneType); } - for (int i = MembraneResolution; i > 0; i--) - { - startingBuffer.Add(new Vector2(cellDimensions - 2 * cellDimensions / MembraneResolution * i, - cellDimensions)); - } - - for (int i = MembraneResolution; i > 0; i--) - { - startingBuffer.Add(new Vector2(cellDimensions, - -cellDimensions + 2 * cellDimensions / MembraneResolution * i)); - } - - for (int i = MembraneResolution; i > 0; i--) - { - startingBuffer.Add(new Vector2(-cellDimensions + 2 * cellDimensions / MembraneResolution * i, - -cellDimensions)); - } - - // ReSharper restore PossibleLossOfFraction - - // Get new membrane points for vertices2D - GenerateMembranePoints(hexPositions, hexCount, membraneType); - // This makes a copy of the vertices so the data is safe to modify in further calls to this method return new MembranePointData(hexPositions, hexCount, membraneType, vertices2D); } @@ -462,6 +453,127 @@ private static ArrayMesh BuildEngulfMesh(Vector2[] vertices2D, int vertexCount, return generatedMesh; } + private static float SampleRadius(ReadOnlySpan<(double Rho, double Theta)> polarData, double theta) + { + int n = polarData.Length; + + switch (n) + { + case 0: + return 1.0f; + case 1: + (double r, double _) = polarData[0]; + return (float)r; + } + + int low = 0, high = n - 1; + + while (low < high) + { + int mid = low + high >> 1; + double midT = polarData[mid].Theta; + + if (midT < theta) + low = mid + 1; + else + high = mid; + } + + // low is the first index with angle >= theta. + int iB = low; + int iA = (iB - 1 + n) % n; + + (double rA, double tA) = polarData[iA]; + (double rB, double tB) = polarData[iB]; + + double dT = tB - tA; + + if (Math.Abs(dT) < 1e-9) + return (float)rA; + + if (dT > Math.PI) + dT -= 2.0 * Math.PI; + + if (dT < -Math.PI) + dT += 2.0 * Math.PI; + + double offset = theta - tA; + + if (offset > Math.PI) + offset -= 2.0 * Math.PI; + + if (offset < -Math.PI) + offset += 2.0 * Math.PI; + + double t = Math.Clamp(offset / dT, 0.0, 1.0); + + return (float)(rA + t * (rB - rA)); + } + + [MethodImpl(MethodImplOptions.AggressiveInlining)] + private bool UsesSuperformulaEnvelope(MembraneType membraneType) + { + return true; //SuperformulaEnvelope.HasValue && membraneType.InternalName == "silica"; + } + + private void GenerateBlobEnvelope(Vector2[] hexPositions, int hexCount, MembraneType membraneType) + { + // The length in pixels (probably not accurate?) of a side of the square that bounds the membrane. + // Half the side length of the original square that is compressed to make the membrane. + int cellDimensions = 10; + + for (int i = 0; i < hexCount; ++i) + { + var pos = hexPositions[i]; + if (MathF.Abs(pos.X) + 1 > cellDimensions) + { + cellDimensions = (int)MathF.Abs(pos.X) + 1; + } + + if (MathF.Abs(pos.Y) + 1 > cellDimensions) + { + cellDimensions = (int)MathF.Abs(pos.Y) + 1; + } + } + + // Make the length longer to guarantee that everything fits easily inside the square + cellDimensions *= 100; + + startingBuffer.Clear(); + + // Integer divides are intentional here + // ReSharper disable PossibleLossOfFraction + + for (int i = MembraneResolution; i > 0; i--) + { + startingBuffer.Add(new Vector2(-cellDimensions, + cellDimensions - 2 * cellDimensions / MembraneResolution * i)); + } + + for (int i = MembraneResolution; i > 0; i--) + { + startingBuffer.Add(new Vector2(cellDimensions - 2 * cellDimensions / MembraneResolution * i, + cellDimensions)); + } + + for (int i = MembraneResolution; i > 0; i--) + { + startingBuffer.Add(new Vector2(cellDimensions, + -cellDimensions + 2 * cellDimensions / MembraneResolution * i)); + } + + for (int i = MembraneResolution; i > 0; i--) + { + startingBuffer.Add(new Vector2(-cellDimensions + 2 * cellDimensions / MembraneResolution * i, + -cellDimensions)); + } + + // ReSharper restore PossibleLossOfFraction + + // Get new membrane points for vertices2D + GenerateMembranePoints(hexPositions, hexCount, membraneType); + } + private void GenerateMembranePoints(Vector2[] hexPositions, int hexCount, MembraneType membraneType) { // Move all the points in the source buffer close to organelles @@ -474,7 +586,9 @@ private void GenerateMembranePoints(Vector2[] hexPositions, int hexCount, Membra var direction = (startingBuffer[i] - closestOrganelle).Normalized(); var movement = direction * Constants.MEMBRANE_ROOM_FOR_ORGANELLES; - startingBuffer[i] = closestOrganelle + movement; + var newPosition = closestOrganelle + movement; + + startingBuffer[i] = newPosition; } float circumference = 0.0f; @@ -549,4 +663,77 @@ private void GenerateMembranePoints(Vector2[] hexPositions, int hexCount, Membra vertices2D[i] = point + movement; } } + + private void GenerateSuperformulaEnvelope(Vector2[] hexPositions, int hexCount) + { + var superformula = SuperformulaEnvelope!.Value; + var polarData = superformula.PolarData; + + // Find extreme values and the center of the hex collection. + float maxX = float.MinValue, maxY = float.MinValue; + float minX = float.MaxValue, minY = float.MaxValue; + foreach (var (organellePositionX, organellePositionY) in hexPositions) + { + if (organellePositionX > maxX) + maxX = organellePositionX; + + if (organellePositionY > maxY) + maxY = organellePositionY; + + if (organellePositionX < minX) + minX = organellePositionX; + + if (organellePositionY < minY) + minY = organellePositionY; + } + + var center = new Vector2(maxX + minX, maxY + minY) / 2.0f; + + var hexPositionsTranslated = ArrayPool.Shared.Rent(hexCount); + + // Translate hexes to local space + for (int i = 0; i < hexCount; ++i) + { + hexPositionsTranslated[i] = hexPositions[i] - center; + } + + // Map everything + var amplitudeX = MathF.Max(MathF.Abs(maxX - minX) / 2.0f, 1.0f); + var amplitudeY = MathF.Max(MathF.Abs(maxY - minY) / 2.0f, 1.0f); + var factor = 1.0f / (superformulaEnvelopeSharpness * superformulaEnvelopeSharpness); + + vertices2D.Clear(); + vertices2D.Capacity = polarData.Length; + + for (int i = 0, end = polarData.Length; i < end; ++i) + { + double theta = polarData[i].Theta; + + Vector2 rayDirection = new Vector2((float)Math.Cos(theta), (float)Math.Sin(theta)); + + float signedRadius = SampleRadius(polarData, theta) * superformulaEnvelopeScale; + + float ellipseRadius = amplitudeX * amplitudeY / MathF.Sqrt(MathF.Pow(amplitudeY * rayDirection.X, 2) + + MathF.Pow(amplitudeX * rayDirection.Y, 2)); + + float stretchedSignedRadius = signedRadius * ellipseRadius; + + Vector2 pointAtInfinity = rayDirection * 100000.0f; + Vector2 closestOrganelle = FindClosestOrganelleHex(hexPositionsTranslated, hexCount, pointAtInfinity); + + float organelleProjection = closestOrganelle.Dot(rayDirection) + Constants.MEMBRANE_ROOM_FOR_ORGANELLES; + float difference = stretchedSignedRadius - organelleProjection; + + // Custom math magic that acts as a smooth max function. It uses a modified hyperbola to join the operand + // values. + float union = 0.5f * (MathF.Sqrt(factor + difference * difference) + stretchedSignedRadius + + organelleProjection); + + Vector2 point = rayDirection * union; + + vertices2D.Add(point + center); + } + + ArrayPool.Shared.Return(hexPositionsTranslated); + } } diff --git a/src/microbe_stage/editor/SuperformulaPlot.cs b/src/microbe_stage/editor/SuperformulaPlot.cs new file mode 100644 index 00000000000..e135a824386 --- /dev/null +++ b/src/microbe_stage/editor/SuperformulaPlot.cs @@ -0,0 +1,8 @@ +using Godot; + +/// +/// A plotter for superformula instances. +/// +public partial class SuperformulaPlot : Control +{ +} diff --git a/src/microbe_stage/editor/SuperformulaPlot.cs.uid b/src/microbe_stage/editor/SuperformulaPlot.cs.uid new file mode 100644 index 00000000000..bfdc044bdee --- /dev/null +++ b/src/microbe_stage/editor/SuperformulaPlot.cs.uid @@ -0,0 +1 @@ +uid://c8wxu6gluhb3y