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