Handling singularities in squircle parametric equations

In order to plot a squircle using its parametric equations:

The equations can be found in page 9 of this this article:

Fong, Chamberlain. "Squircular Calculations." arXiv preprint arXiv:1604.02174 (2016).

Unfortunately the following singularities are not handled by those equations:

• For t = 0, t = pi/2, t = pi, t = 3pi/2, and t=2*pi

Here is the code to plot these equations in Python

"""Plot squircle shape using its parametric equations"""
import numpy as np
from numpy import pi, sqrt, cos, sin, sign
import matplotlib.pyplot as plt

# FG stands for Fernandez-Guasti who introduced the algebric equations.
def FG_param(r, s, t):
x = np.ones_like(t)
y = np.ones_like(t)
for i, _t in enumerate(t):
if np.isclose(_t, 0.0):
x[i] = r
y[i] = 0.0
elif np.isclose(_t, pi/2.0):
x[i] = 0.0
y[i] = r
elif np.isclose(_t, pi):
x[i] = -r
y[i] = 0.0
elif np.isclose(_t, 3*pi/2.0):
x[i] = 0.0
y[i] = -r
elif np.isclose(_t, 2*pi):
x[i] = r
y[i] = 0.0
else:
x[i] = r*sign(cos(t[i]))/(s*sqrt(2)*np.abs(sin(t[i])))*sqrt(1 - sqrt(1-s**2*sin(2*t[i])**2))
y[i] = r*sign(sin(t[i]))/(s*sqrt(2)*np.abs(cos(t[i])))*sqrt(1 - sqrt(1-s**2*sin(2*t[i])**2))
return x, y

s = 0.7
r = 1.0
NTHETA = 90
t = np.linspace(0, 2*pi, NTHETA)

x, y = FG_param(r, s, t)

plt.gca().set_aspect('equal')
plt.scatter(x, y, s=5)
plt.show()


The function FG_param looks very clumsy. I appreciate any help to make it more compact and neat.

Thank you

• It is possible to simplify via knowing that boolean statements return 0 or 1 if they are coerced to; for example, np.isclose(_t, pi/2.0) returns 1 if true, so that can be inverted (display zero if matched) and then multiplied by the RHS to remove it; the equation can be restructured to be r + xsin(...), where x is the condition. This can be repeated for all remaining conditions and structured accordingly. I may provide a code sample Dec 6, 2019 at 4:39
• Something sort of like this, but this is not 100% correct: x[i] = r * \ (np.isclose(_t, 0.0) or np.isclose(_t, 2*pi)) * \ (not np.isclose(_t, 0.0) or np.isclose(_t, 2*pi)) * \ (1 - (not np.isclose(_t, pi))) * \ (not np.isclose(_t, 3*pi/2.0)) Dec 6, 2019 at 4:40

Here's a different solution. The singularities occur when sin(_t) or cos(_t) is zero. So check for those conditions.

def FG_param(r, s, t):
x = np.ones_like(t)
y = np.ones_like(t)

for i, _t in enumerate(t):
sin_t = sin(_t)
cos_t = cos(_t)

if np.isclose(sin_t, 0.0):
x[i] = -r if np.isclose(_t, pi) else r
y[i] = 0.0

elif np.isclose(cos_t, 0.0):
x[i] = 0.0
y[i] = r if _t < pi else -r

else:
radical = sqrt(1 - sqrt(1 - s**2 * sin(2*_t)**2))

x[i] = r*sign(cos_t) / (s*sqrt(2)*np.abs(sin_t)) * radical
y[i] = r*sign(sin_t) / (s*sqrt(2)*np.abs(cos_t)) * radical

return x, y

• The singularities occur also when s=0 Dec 6, 2019 at 12:13
• @Navaro: Exactly. In that case (i.e: s=0) the map will be a circle with radius r. Dec 7, 2019 at 0:31
• @Navaro s==0 wasn't handled in the original code, so I didn't handle it either. It would be an easy case to handle. Dec 7, 2019 at 1:08
• @RootTwo: Yes, you could just add something like this at the top of the function: if s == 0: return r*cos(t), r*sin(t) Dec 7, 2019 at 1:12

This solution is an order of magnitude faster than the other ones. It uses numpy functions, so there aren't any explicit loops. np.where chooses values that keep the equations from blowing up where the denominator would go to zero. Case for s == 0 is handled separately.

def FG_param1(r, s, t):

if s == 0:
return r*cos(t), r*sin(t)

sin_t = sin(t)
cos_t = cos(t)

radical = sqrt(1 - sqrt(1 - s**2 * sin(2*t)**2))

x_denom = np.where(np.isclose(sin_t, 0.0), 1.0, s*sqrt(2)*np.abs(sin_t))
x = r * sign(cos_t) * np.where(np.isclose(sin_t, 0.0, 1e-12), 1.0, radical / x_denom)

y_denom = np.where(np.isclose(cos_t, 0.0), 1.0, s*sqrt(2)*np.abs(cos_t))
y = r * sign(sin_t) * np.where(np.isclose(cos_t, 0.0, 1e-12), 1.0, radical / y_denom)

return x, y


def FG_param(r, s, t):
x = np.ones_like(t)
y = np.ones_like(t)

for i, _t in enumerate(t):
radical = sqrt(1 - sqrt(1 - s**2 * sin(2*_t)**2))

try:
x[i] = r*sign(cos(_t)) / (s*sqrt(2)*np.abs(sin(_t))) * radical

except ZeroDivisionError:
x[i] = -r if np.isclose(_t, pi) else r

try:
y[i] = r*sign(sin(_t)) / (s*sqrt(2)*np.abs(cos(_t))) * radical

except ZeroDivisionError:
y[i] = r if _t < pi/2.0 else -r

return x, y