等角写像による画像の変換〜Schwarz-Christoffel 変換〜part 2 - IMACEL Academy -人工知能・画像解析の技術応用に向けて-| エルピクセル株式会社

等角写像による画像の変換〜Schwarz-Christoffel 変換〜 part1 - IMACEL Academy -人工知能・画像解析の技術応用に向けて-|LPixel(エルピクセル)

等角写像の一つであるSchwarz Christoffel 変換を用いて，画像の変換をしてみます．python によるコードも記載しております．画像はhttps://uk.mathworks.com/help/images/examples/exploring-a-conformal-mapping_ja_JP.htmlより．

1. Schwarz- Christoffel 変換

Schwarz - Christoffel 変換とは，半平面を等角写像に基づいて三角形，四角形などの多角形に変換する写像です．具体的な変換式は以下のようになります．

Schwarz-Christoffel 変換

$w$平面上の多角形に関して，その各頂点の内角を順番に $\alpha_1 \pi, \alpha_2 \pi,\cdots, \alpha_n \pi$とするとき \begin{align} \frac{dw}{dz} = \gamma(z-a_1)^{\alpha_1 - 1}(z-a_2)^{\alpha_2 - 1} \cdots
(z-a_n)^{\alpha_n - 1} \end{align} によって，半平面から$n$角形に写される．

今回の場合，四角形に写すことを考えていますので，$\alpha_1 = \alpha_2 = \alpha_3 = \alpha_4 = 1/2$を代入し，さらに$a_1=-1/k, a_2 =-1, a_3 = 1, a_4 = 1/k$とすると， \begin{align} \frac{dw}{dz} =\frac{\tilde{\gamma}}{\sqrt{(1-z^2)(1-k^2z^2)}} \end{align}

となります．これは，ヤコビの楕円関数と深く関係しています．詳しくは以下のページをご覧ください．

https://ja.wikipedia.org/wiki/ヤコビの楕円関数

2. 今回のフローチャート

今回は，正方形から円への写像を実装します．具体的には，

Schwartz- Christoffel の逆変換により，正方形を半平面に
半平面を円に

の方法によって行います．

四角形から円に写す等角写像

Schwarz-Christoffel 変換と半平面から円に写す写像を組み合わせることで，四角形から円への写像を実現する．

3. 実装

それでは実装をしていきます．まずは，以下のように格子点を定義しておきます．

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

パッケージのインストール.py

# Schwarz Christoffel Mapping

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

for x_i in x_m:
    plt.plot(x_i*np.ones(m),x,color='black')
    plt.plot(x,x_i*np.ones(m),color='black', 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()

格子点定義.py

次に，必要となる楕円関数などを定義します．python のscipy の関数は引数として複素数値を受け入れないので，級数を用いて定義し直します．

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):
        t = t + np.exp(pi*1j* n*n*tau + 2*pi*n*1j*z)
    return t

def teta01(z,tau):
    t = teta(z+1.0/2,tau,0)
    return t


def teta11(z,tau):
    t = teta(z+1.0/2,tau,1)
    return t

def teta10(z,tau):
    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)
    l = 1.0/2 * (1.0 - np.sqrt(k_dash))/(1.0 + np.sqrt(k_dash)) 
    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)
    # 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)
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))

関数定義.py

なお，ここでのテータ関数は，以下のように定義したものです．

\begin{align} \vartheta(z,\tau) = \sum_{n=-\infty}^{\infty} \exp\left(\pi i n^2 \tau + 2\pi i n z \right) \end{align}

$q = \exp(\pi i \tau)$として

\begin{align} \vartheta_0(z,\tau) &:= \vartheta_{01}(z, \tau) = \sum_{n=-\infty}^{\infty} e^{\pi i \tau n^{2} + 2 \pi i n \left(z + \frac{1}{2}\right)} =1 + 2 \sum_{n=1}^{\infty} (-1)^{n} q^{n^{2}} \cos 2 n \pi z \end{align} \begin{align} \vartheta_{1}(z, \tau) &:= - \vartheta_{11}(z, \tau) = - \sum_{n=-\infty}^{\infty} e^{\pi i \tau \left(n + \frac{1}{2}\right)^{2} + 2 \pi i \left(n + \frac{1}{2}\right) \left(z + \frac{1}{2}\right)} =2\sum_{n=0}^{\infty}{(-1)^{n} q^{{\left(n + \frac{1}{2}\right)}^2} \sin(2 n + 1) \pi z} \end{align} \begin{align} \vartheta_{2}(z, \tau) &:= \vartheta_{10}(z, \tau) = \sum_{n=-\infty}^{\infty} e^{\pi i \tau \left(n + \frac{1}{2}\right)^2 + 2 \pi i \left(n + \frac{1}{2}\right)z} = 2 \sum_{n=0}^{\infty} q^{{\left(n+\frac{1}{2}\right)}^2} \cos (2 n + 1) \pi z \end{align} \begin{align} \vartheta_3(z, \tau) &:= \vartheta_{00}(z, \tau) = \sum_{n=-\infty}^{\infty} e^{i \pi \tau n^{2} + 2 \pi i n z} = 1 + 2 \sum_{n=1}^{\infty} q^{n^{2}} \cos 2 n \pi z. \end{align}

上の関数を用いて変換してみます．

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='black')
    plt.plot(np.real(w_v),np.imag(w_v),color='black', 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()

変換.py

見事，円に変換されましたね！

4. 画像への適用

それでは，次は実際の画像に関して，この変換を適用してみます．
今回用いるのは，このような将棋の画像です．

# load image 
img = cv2.imread('picture_c/base200.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()

画像のロード.py

これに関して，変換してみます．

# 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.show()
plt.savefig("output/base.png")

Schwarz-Christoffel.py

なかなか面白い結果が得られましたね！

5. まとめ

今回のシリーズでは，等角写像を実際に実装して，確認してみました．そして，等角写像の最重要ポイントである，Schwarz-Christoffel 変換に関して実装し，どのように変換されるのか，確認してみました．

実は，この等角写像は数学的・物理的・また医学的にも重要です．例えば脳の位置合わせや，はたまた透明マントの作り方など，実は奥が深いのです．
古典的な方法ではありますが，興味を持つ人が一人でも増えることを願って，この記事を終わりにしたいと思います．

参考文献

Analytical Methods for Squaring the Disc https://arxiv.org/pdf/1509.06344.pdf