jellyfish_Hamid_Naderi_Yeganeh with numpy and jax

jellyfish_jax.py

# # https://zhuanlan.zhihu.com/p/1920686890295734431?share_code=iC4x5GhaOYvt&utm_psn=1921130722581776046
import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt

def Cs(C):
    return jnp.sum(jnp.stack(C, axis=-1), axis=-1)

def Cp(C):
    return jnp.prod(jnp.stack(C, axis=-1), axis=-1)

def L(x, y, v):
    return ((7/10) + (3/2) * jnp.exp(-jnp.exp(400 * (y - 27/100))) * (y + 1/5)**2) * x + \
           jnp.exp(-jnp.exp(400 * (y - 3/10))) * ((3/10)*(y - 3/10)**2 + jnp.cos(5*y + 6*x + v)/20)

def R(x, y, v):
    l = L(x, y, v)
    return jnp.exp(
        -jnp.exp(
            -50 * (jnp.sin(v + 30 * jnp.arccos((9/5)*l*(0.5 + (1/jnp.pi)*jnp.arctan(-180*jnp.abs(l) + 100)))
            ) + 1/50 - y**6/10 + (2/jnp.pi)*jnp.arctan(200*(y - 1/4)))
        ) - jnp.exp(1000*((9/5)*jnp.abs(l) - 9/10)) - jnp.exp(100*y - 25)
    )

def Q(x, y, s):
    return y - (10 - s)/30 + 1/5 - x**2 + jnp.cos(8*x + 5*s)/10

def P(x, y, s):
    q = Q(x, y, s)
    terms = [(50 - s)/1000 * (3/5)**u * jnp.arctan(
        3 * jnp.cos(
            2**u * (1 + s) * jnp.arctan(1000*x / (250*(2*y - 1)**2 + 1)) +
            0.5*(7/5)**u * jnp.cos(2**u * (5 + s)*q) + 2*u*s
        )
    ) for u in range(8)]
    return q + Cs(terms)

def J(x, y, s):
    return jnp.exp(
        -jnp.exp(500 * (jnp.abs(x + ((40 - s)/40 * jnp.cos(8*s) + 1/2)*(y - 1/2)**2) - jnp.pi/15
                - (3/20)*jnp.arctan(-8*(y - (10 - s)/30))))
        - jnp.exp(-200 * P(x, y, s))
        - jnp.exp(1000 * (y - 7/10))
    )

def T(x, y, s):
    A = (23**s * 20**(-s) * 10 * (1 + jnp.cos(10*s)))
    B1 = A * (jnp.cos(7*s)*x + jnp.sin(7*s)*y + 2*jnp.cos(17*s))
    B2 = 4 * jnp.cos((23**s * 20**(-s) * 10) * (jnp.cos(7*s)*x + jnp.sin(7*s)*y))
    B3 = 2 * jnp.cos(5*s)
    C1 = A * (jnp.cos(7*s)*y - jnp.sin(7*s)*x + 2*jnp.cos(15*s))
    C2 = 4 * jnp.cos((23**s * 20**(-s) * 10) * (jnp.cos(8*s)*x + jnp.sin(8*s)*y))
    C3 = 2 * jnp.cos(7*s)
    return jnp.cos(B1 + B2 + B3) * jnp.cos(C1 + C2 + C3)

def E(x, y):
    return Cs([3*(25**s / 26**s) * T(x, y, s)/10 for s in range(1, 51)])

def C_fn(x, y):
    def inner(s):
        prod = Cp([1 - (3/4) * jnp.exp(-jnp.exp(-100*(u - 1/2))) * J(x, y, u) for u in range(s)])
        term = jnp.exp(-jnp.exp(3 * Q(x, y, s))) + \
               (2/5) * jnp.exp(-jnp.exp(20*jnp.abs(P(x, y, s)) - 1 + E(x, y)/5)) + \
               (2/5) * jnp.exp(-jnp.exp(100*jnp.abs(P(x, y, s)) - 3 + E(x, y)))
        return prod * J(x, y, s) * term
    return Cs([inner(s) for s in range(1, 36)])

def U(x, y):
    return 9/20 + 2*jnp.exp(-jnp.exp(14 * (jnp.abs(jnp.abs(x)**(137/20 - 5*y) + 2*(y - 57/100)**2 - 27/100) - 1/25 + E(x, y)/100))) + \
           (3/5)*jnp.exp(-jnp.exp(jnp.abs(20*jnp.abs(x)**(13/2 - 5*y) + 2*(y - 1/2)**2 - 3/20) - 4/5 + E(x, y)/5)) + \
           (1/5)*jnp.exp(-jnp.exp(100*(jnp.abs(x)**(13/2 - 5*y) + 2*(y - 1/2)**2 - 3/20)))

def W(x, y, v):
    return jnp.exp(-jnp.exp(v * (jnp.abs(x)**(7 - 5*y) + 2*(y - 3/5)**2 - 1/4)))

def K(x, y):
    return jnp.exp(-jnp.exp(jnp.cos(20*x * ((1/1000)*(jnp.abs(250 - 20*(10*y - 6)**2) + 1))**(-1000/(1000*jnp.abs(7 - 5*y) + 1)))) + E(x, y)/3)

def H(x, y, v):
    part1 = ((19 - 21*v + 7*v**2)/10) * (R(x, y, 0) + R(x, y, 1)*(1 - R(x, y, 0)) * Cp([1 - (3/4)*J(x, y, s) for s in range(1, 36)]))
    part2 = (9*v/40)*(1 + 2*y/3 - x**2/20 + E(x, y)/20)
    part3 = ((5 - 3*v + v**2)/5) * C_fn(x, y) * (1 - R(x, y, 0)) * (1 - (3/5)*W(x, y, 500))
    part4 = ((5 - 3*v + v**2)/5) * (W(x, y, 5) + y/4 - 3/20 + (3/10 - v/20 + (2 - v)/5 * (y - 2/5)) * K(x, y) + \
            (8 - 7*v)/20 * jnp.exp(-jnp.exp(10*y - 5 + E(x, y)))) * W(x, y, 500) * U(x, y)
    return part1 + part2 + part3 + part4

def F(x):
    return jnp.floor(255 * jnp.exp(-jnp.exp(-1000 * x)) * jnp.abs(x)**jnp.exp(-jnp.exp(1000 * (x - 1))))

# 创建网格
step=2
m, n = jnp.meshgrid(jnp.arange(1, 2001, step), jnp.arange(1, 1201, step))
x = (m - 1000) / 600
y = (651 - n) / 600

# 计算图像
Rp = F(H(x, y, 0)).astype(jnp.uint8)
Gp = F(H(x, y, 1)).astype(jnp.uint8)
Bp = F(H(x, y, 2)).astype(jnp.uint8)

# 合并为RGB图像
Pp = jnp.stack([Rp, Gp, Bp], axis=-1)

# 显示图像
plt.imshow(Pp)
plt.axis('off')                                        # 去掉坐标轴
plt.subplots_adjust(left=0, right=1, top=1, bottom=0)  # 图像填满整个Figure
plt.margins(0, 0)                                      # 去掉图像边距
plt.gca().set_position([0, 0, 1, 1])                   # 图像铺满Axes区域
plt.show()

jellyfish_numpy.py

# https://zhuanlan.zhihu.com/p/1920686890295734431?share_code=iC4x5GhaOYvt&utm_psn=1921130722581776046
import numpy as np
import matplotlib.pyplot as plt

def my_exp(x):
    if isinstance(x, float):
        x = min(x, 400)
    elif isinstance(x, np.ndarray):
        x_greater_400 = x > 400
        x[x_greater_400] = 400
    return np.exp(x, dtype=np.float64)

# 定义函数
def Cs(C):
    return np.sum(np.stack(C, axis=2), axis=2)

def Cp(C):
    return np.prod(np.stack(C, axis=2), axis=2)

def L(x, y, v):
    return ((7/10) + (3/2)*my_exp(-my_exp(400*(y - 27/100)))*(y + 1/5)**2)*x + my_exp(-my_exp(400*(y - 3/10)))*((3/10)*(y - 3/10)**2 + np.cos(5*y + 6*x + v)/20)

def R(x, y, v):
    return my_exp(-my_exp(-50*(np.sin(v + 30*np.arccos((9/5)*L(x, y, v)*(0.5 + (1/np.pi)*np.arctan(-180*np.abs(L(x, y, v)) + 100)))) + 1/50 - y**6/10 + (2/np.pi)*np.arctan(200*(y - 1/4)))) - my_exp(1000*((9/5)*np.abs(L(x, y, v)) - 9/10)) - my_exp(100*y - 25))

def Q(x, y, s):
    return y - (10 - s)/30 + 1/5 - x**2 + np.cos(8*x + 5*s)/10

def P(x, y, s):
    return Q(x, y, s) + Cs([((50 - s)/1000)*(3/5)**u*np.arctan(3*np.cos(2**u*(1 + s)*np.arctan(1000*x/(250*(2*y - 1)**2 + 1)) + 0.5*(7/5)**u*np.cos(2**u*(5 + s)*Q(x, y, s)) + 2*u*s)) for u in range(8)])

def J(x, y, s):
    return my_exp(-my_exp(500*(np.abs(x + ((40 - s)/40*np.cos(8*s) + 1/2)*(y - 1/2)**2) - np.pi/15 - (3/20)*np.arctan(-8*(y - (10 - s)/30)))) - my_exp(-200*P(x, y, s)) - my_exp(1000*(y - 7/10)))

def T(x, y, s):
    return np.cos((23**s*20**(-s)*10*(1 + np.cos(10*s))*(np.cos(7*s)*x + np.sin(7*s)*y + 2*np.cos(17*s)) + 4*np.cos(23**s*20**(-s)*10*(np.cos(7*s)*x + np.sin(7*s)*y)) + 2*np.cos(5*s)))*np.cos((23**s*20**(-s)*10*(1 + np.cos(10*s))*(np.cos(7*s)*y - np.sin(7*s)*x + 2*np.cos(15*s)) + 4*np.cos(23**s*20**(-s)*10*(np.cos(8*s)*x + np.sin(8*s)*y)) + 2*np.cos(7*s)))

def E(x, y):
    return Cs([3*(25**s/26**s)*T(x, y, s)/10 for s in range(1, 51)])

def C(x, y):
    return Cs([Cp([1 - (3/4)*my_exp(-my_exp(-100*(u - 1/2)))*J(x, y, u) for u in range(s)])*J(x, y, s)*(my_exp(-my_exp(3*(y - (10 - s)/30 + 1/5 - x**2 + np.cos(8*x + 5*s)/10))) + (2/5)*my_exp(-my_exp(20*np.abs(P(x, y, s)) - 1 + E(x, y)/5)) + (2/5)*my_exp(-my_exp(100*np.abs(P(x, y, s)) - 3 + E(x, y)))) for s in range(1, 36)])

def U(x, y):
    return 9/20 + 2*my_exp(-my_exp(14*(np.abs(np.abs(x)**(137/20 - 5*y) + 2*(y - 57/100)**2 - 27/100) - 1/25 + E(x, y)/100))) + (3/5)*my_exp(-my_exp(np.abs(20*np.abs(x)**(13/2 - 5*y) + 2*(y - 1/2)**2 - 3/20) - 4/5 + E(x, y)/5)) + (1/5)*my_exp(-my_exp(100*(np.abs(x)**(13/2 - 5*y) + 2*(y - 1/2)**2 - 3/20)))

def W(x, y, v):
    return my_exp(-my_exp(v*(np.abs(x)**(7 - 5*y) + 2*(y - 3/5)**2 - 1/4)))

def K(x, y):
    return my_exp(-my_exp(np.cos(20*x*((1/1000)*(np.abs(250 - 20*(10*y - 6)**2) + 1))**(-1000/(1000*np.abs(7 - 5*y) + 1)))) + E(x, y)/3)

def H(x, y, v):
    return ((19 - 21*v + 7*v**2)/10)*(R(x, y, 0) + R(x, y, 1)*(1 - R(x, y, 0))*Cp([1 - (3/4)*J(x, y, s) for s in range(1, 36)])) + (9*v/40)*(1 + 2*y/3 - x**2/20 + E(x, y)/20) + ((5 - 3*v + v**2)/5)*C(x, y)*(1 - R(x, y, 0))*(1 - 3/5*W(x, y, 500)) + ((5 - 3*v + v**2)/5)*(W(x, y, 5) + y/4 - 3/20 + (3/10 - v/20 + (2 - v)/5*(y - 2/5))*K(x, y) + (8 - 7*v)/20*my_exp(-my_exp(10*y - 5 + E(x, y))))*W(x, y, 500)*U(x, y)

def F(x):
    return np.floor(255*my_exp(-my_exp(-1000*x))*np.abs(x)**my_exp(-my_exp(1000*(x - 1))))

# 创建网格
step=2
m, n = np.meshgrid(np.arange(1, 2001, step), np.arange(1, 1201, step))
x = (m - 1000) / 600 
y = (651 - n) / 600

# 计算图像
Rp = F(H(x, y, 0)).astype(np.uint8)
Gp = F(H(x, y, 1)).astype(np.uint8)
Bp = F(H(x, y, 2)).astype(np.uint8)

# 合并为RGB图像
Pp = np.stack((Rp, Gp, Bp), axis=2)

# 显示图像
plt.imshow(Pp)
plt.axis('off')                                        # 去掉坐标轴
plt.subplots_adjust(left=0, right=1, top=1, bottom=0)  # 图像填满整个Figure
plt.margins(0, 0)                                      # 去掉图像边距
plt.gca().set_position([0, 0, 1, 1])                   # 图像铺满Axes区域
plt.show()

numpy=2.x.x
jax=0.6.2