結構難しいです

大学時代習ってきたことの復習かも知れませんが結構難しいです

説明文が多すぎますがお許し下さい

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import scipy.special
import cv2

#四角形
s = 1
m = 1000
x_m = np.linspace(-s, s, m/50)
x = np.linspace(-s, s, m)

plt.figure(figsize=(15,15))
plt.rcParams["font.size"] = 30

for i in x_m:
    plt.plot(i*np.ones(m), x, color='c')   # np.oens(m):要素がm個の１次元配列 array([1.0, 1.0, 1.0 ・・・1.0])
    plt.plot(x, i*np.ones(m), color='c', linestyle='dashed')

plt.gca().set_aspect('equal', adjustable='box')

sc = 1.2
plt.xlim(-s*sc, s*sc)
plt.ylim(-s*sc, s*sc)
plt.show()

def teta(z,tau,c): # c = 0 or 1
    c = np.double(c)/2
    n_max = 100.0 + c
    pi = np.pi
    t = 0.0 + 0.0j;
    for n in np.linspace(-n_max,n_max,2*n_max+1):    # 0を含めて-n_maxからn_maxまで2*n_max+1分割
        t = t + np.exp(pi*1j* n*n*tau + 2*pi*n*1j*z) # t + exp(πi n^2τ + 2πn i z)-----①
    return t 

def teta01(z,tau):
    t = teta(z+1.0/2,tau,0)   # ①のzの部分がz+1/2、n の部分はそのまま
    return t

def teta11(z,tau):
    t = teta(z+1.0/2,tau,1)  # ①のzの部分がz+1/2、n の部分が n+1/2
    return t

def teta10(z,tau):            # ①のzの部分はそのまま、n の部分がn+1/2
    t = teta(z,tau,1)
    return t

def teta00(z,tau):           # ①の式そのまま
    t = teta(z,tau,0)
    return t

def cn(u,k):
    # k -> tau
    # Reference: Wikipedia
    # https://ja.wikipedia.org/wiki/%E3%83%A4%E3%82%B3%E3%83%93%E3%81%AE%E6%A5%95%E5%86%86%E9%96%A2%E6%95%B0
    k_dash = np.sqrt(1.0 - k*k)    # k'(τ)を(1 - k^2)^0.5とする
    l = 1.0/2 * (1.0 - np.sqrt(k_dash))/(1.0 + np.sqrt(k_dash))   # l = 1/2 (1-√k')/(1+√k')と定める
    # q=exp(πiτ)としてlの式をベキ級数展開すると
    # l = (q + q^9 + q^25 + ・・・)/(1 + 2q^4 + 2q^16・・・)
    # これの級数の反転を行うと
    # q = l + 2l^5 + 15l^9 + 150l^13 + 1707l^17 + 20910l^21 + 268616l^25・・・となる
    q = l + 2*np.power(l,5) + 15*np.power(l,9)  + 150*np.power(l,13) + 1707*np.power(l,17)
    tau = np.log(q) / (np.pi * 1j)  # qの定義より当然こうなる
    # calculate
    t10 = teta10(0.0+0.0j , tau)
    t00 = teta00(0.0+0.0j , tau)
    t01 = teta01(0.0+0.0j, tau)
    # check print(np.sqrt(1.0 - np.power(t01/t00,4) )) = k    
    z = u/(np.pi* t00 * t00)
    return t01 * teta10(z,tau) / (t10*teta01(z,tau))

Ke = sp.special.ellipk(1.0/2) # 第１種完全楕円積分 E(k)=∫_0->2/π (1 - k^2 son^2θ)^-0.5 dθ、k^2 < 1
def schw_chrf(z):
    # reference: https://arxiv.org/pdf/1509.06344.pdf
    return (1.0 - 1.0j)/np.sqrt(2) * cn( Ke*(1.0+1.0j)/2 * z - Ke, 1.0/np.sqrt(2))

plt.figure(figsize=(15,15))
plt.rcParams["font.size"] = 30
   
for x_i in x_m:
    z_h = x_i*np.ones(m) + x*1j
    z_v = x+1j*x_i*np.ones(m)
    w_h = schw_chrf(z_h)
    w_v = schw_chrf(z_v)
    plt.plot(np.real(w_h),np.imag(w_h),color='c')
    plt.plot(np.real(w_v),np.imag(w_v),color='c', linestyle='dashed')   

plt.gca().set_aspect('equal', adjustable='box')

sc = 1.1
plt.xlim(-s*sc,s*sc)
plt.ylim(-s*sc,s*sc)
plt.show()   

# load image
img = cv2.imread('img/station.png')
#img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)

# show image with matplotlib
pix = img.shape[0]
x = np.linspace(-1,1,pix+1) # define square

plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.gray()

# text mapping
# schwartz christoffel mapping
import sys

for i in range(pix):
    sys.stdout.write(str(i))
    sys.stdout.write(" ")
    for j in range(pix):
        tx_z = np.array([x[i],x[i+1],x[i+1],x[i]])
        ty_z = np.array([x[j],x[j],x[j+1],x[j+1]])
        zt = tx_z + 1j*ty_z
        w_h = schw_chrf(zt)
        tx = np.real(w_h)
        ty = np.imag(w_h)
        clr_b = np.double(img[ pix - j - 1 , i,0])/256
        clr_g = np.double(img[pix - j - 1 , i,1])/256
        clr_r = np.double(img[pix - j - 1 , i,2])/256       
        plt.fill(tx,ty,color=(clr_r,clr_g,clr_b))

plt.gca().set_aspect('equal', adjustable='box')
sc = 1.1
plt.xlim(-s*sc,s*sc)
plt.ylim(-s*sc,s*sc)
plt.savefig("output/cercle.png")
plt.show()