Handling singularities in squircle parametric equations

Question

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

Answer 1

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

Answer 2

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

Answer 3

Off the top of my head, and untested, how about:

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