動的輪郭法(snake)をpythonで書いてみた

画像の領域分割方法には様々な手法がある
そのなかでも簡単なアルゴリズムである動的輪郭モデルをライブラリ１行でポチー！
ではなく実際に書いてみることにした

動的輪郭モデルとは

輪郭を形成する頂点集合$A$を考える。その集合$A$を用いて定義されるエネルギー$E(A)$を最小にする問題を解くことで得られる集合$A_{ans}$を正解の輪郭にしましょうねっていう考え方らしい。
千葉大のスライド
http://www.cfme.chiba-u.jp/~haneishi/class/7_segmentation.pdf
あるまじろさんのスライド
https://www.slideshare.net/Arumaziro/ss-37035661
栄藤稔さんのpdf
https://www.jstage.jst.go.jp/article/mit/12/1/12_9/_pdf

そんなわけで定義されるエネルギーを考えることになるが、どうやらこのエネルギーを定義する時に用いる関数は流派がある。（というか目的輪郭画像の形状にある程度あわせていろいろ種類があるみたい。）
というわけで基本的な定義から、今回の抽出に向けて適当にかんがえたエネルギーをpython3で実装してみた。

エネルギーの定義

栄藤稔さんの資料を元にざっくりと作成した(https://www.jstage.jst.go.jp/article/mit/12/1/12_9/_pdf)
$$E_{snake}^{*} = \int_{0}^{1}E_{In}(\boldsymbol{v}(s)) + E_{Ex}(\boldsymbol{v}(s))+ E_{Con}(\boldsymbol{v}(s))ds$$
$$E_{in}(\boldsymbol{v}(s)) = \frac{\alpha|\boldsymbol{v_s}(s)|^2 + \beta |\boldsymbol{v_{ss}}(s)|^2}{2}$$
$$E_{Ex}(\boldsymbol{v}(s)) = -\gamma |\nabla I(\boldsymbol{v_{ss}}(s))|^2$$
$$E_{Con}(\boldsymbol{v}(s)) = \kappa |\boldsymbol{v}(s)-\boldsymbol{v_g}(s)|^2$$

今回は以上のように定義した$E_{snake}^{*}$を最小にする頂点を計算している。

実装結果

+日の丸をinputとして用意した

• 領域分割した結果
  
  これはこれは300ループで頂点数は200でやっている

• にゃんちゅうをinputとして用意した
  

• 領域分割した結果
  

うーん微妙な精度で領域を見つけている
これは500ループで頂点数は２０００でやっている

実装

Snake.py

#coding:utf-8
import numpy as np
import cv2
import math
import copy


#画像の読み込み
test = cv2.imread("./Img/input.jpg", cv2.IMREAD_COLOR)#BGRなので気をつける

gray_test = cv2.imread("./Img/input.jpg",cv2.IMREAD_GRAYSCALE)
height = test.shape[0]
width = test.shape[1] 
#画像の書き出し
cv2.imwrite('./Img/test.jpg', test)
cv2.imwrite('./Img/gray_test.jpg',gray_test)

N = 2000

v = np.zeros((N,2))
start_v = np.zeros((N,2))
vec_g = np.zeros(2)
for i in range(0,N):
    if(i<N/4):
        v[i] = [height/(N/4)*i,0]
    elif(i<2*N/4):
        v[i] = [height-1,width/(N/4)*(i-N/4)]
    elif(i<3*N/4):
        v[i] = [height-1-height/(N/4)*(i-2*N/4),width-1]
    else:
        v[i] = [0,width-1 - width/(N/4)*(i-3*N/4)]
       
        
#初期輪郭点が円形の場合はこっち
# for i in range(0,N):
#     v[i] = [ height/2*math.sin(2*math.pi*i/N)+height/2, width/2*math.cos(2*math.pi*i/N)+width/2]
#     start_v[i] = [ height/2*math.sin(2*math.pi*i/N)+height/2, width/2*math.cos(2*math.pi*i/N)+width/2]


start_v = copy.deepcopy(v)
display = copy.deepcopy(test)

#パラメータ
alpha =1
beta = 1
gamma = 10
kappa =1

def EpsIn(vec0,vec1,vec2):#test
    value = 0
    value += alpha*np.linalg.norm(vec1-vec0)**2+beta*np.linalg.norm(vec2-2*vec1+vec0)**2
    value /= 2
#     print("In:"+str(value))
    return value

def EpsEx(vec0,pix):#gray
    value = 0
    x = int(vec0[0])
    y = int(vec0[1])
    
    if(x+1 >= height or y+1 >= width):
        return float('inf') 
    else:
        I = [abs(int(pix[x+1,y]) - int(pix[x,y])) ,abs(int(pix[x,y+1])-int(pix[x,y]))]
        value = -gamma*np.linalg.norm(I)**2
#         print("Ex:"+str(value))
        return value

def EpsCon(vec0,vec_g):#test
    value = 0
    value += kappa*np.linalg.norm((vec0[0] - vec_g[0],vec0[1]-vec_g[1]))**2
#     print("Con:"+str(value))
    return value
    
def Energy(vec0,vec1,vec2,vec_g,pix):
    value = 0
    value = EpsIn(vec0,vec1,vec2)+EpsEx(vec0,pix)+EpsCon(vec0,vec_g)
#     print("Result:"+str(value))
    return value
#探索
n = 500
dx = [-1, -1, -1, 0, 0, 0, 1, 1, 1]
dy = [1, 0, -1, 1, 0, -1, 1, 0, -1]
# dx = [1,1,1,0,0,0,-1,-1,-1]
# dy = [1,-1,0,1,-1,0,1,-1,0]
#210
#543
#876

flag = 4
for loop in range(0,n):
    for i in range(0,N):
        flag = 4
        eps_min = float('inf') 
        vec_g = [0,0]
        
        #重心中心にするならこれ
#         for j in range(0,N):
#             vec_g += [v[j,0],v[j,1]]
            
        for j in range(0,9):            
            move  = [v[i,0]+dx[j], v[i,1]+dy[j]]
            if(move[0] < 0 or move[1] < 0 or move[0] >= height  or move[1] >= width):
                continue #はみ出し処理
                
            #重心中心にするならこれ
            #vec_g += [dx[j],dy[j]]
            #vec_g =[vec_g[0]/N, vec_g[1]/N]
            #画像中心を基準に
            vec_g = [int(height/2),int(width/2)]
        
            energy = Energy(move,v[(i+1)%N],v[(i+2)%N],vec_g,gray_test)
            if(eps_min>energy):
                eps_min = energy
                flag = j
        v[i] += [dx[flag],dy[flag]]
        
        #逐次書き出し
    if(loop%10==0):
        cv2.imwrite('./Img/result'+str(loop)+'.jpg', display)
        display = copy.deepcopy(test)
        for i in range(0,N):
            cv2.line(display, (int(v[i,1]),int(v[i,0])), (int(v[(i+1)%N,1]),int(v[(i+1)%N,0])), (0, 255, 0), 2)
            


for i in range(0,N):
    cv2.line(display, (int(v[i,1]),int(v[i,0])), (int(v[(i+1)%N,1]),int(v[(i+1)%N,0])), (0, 255, 0), 2)

for i in range(0,N):
    cv2.line(display, (int(start_v[i,1]),int(start_v[i,0])), (int(start_v[(i+1)%N,1]),int(start_v[(i+1)%N,0])), (255, 0, 0), 2)

cv2.imwrite('./Img/result.jpg', display)

実装のコメント

８近傍を探索しまくるスクリプト。
このエネルギー関数では拡大するような項をかんがえていないので、
いったん考えている境界の内部に頂点が入ってしまった場合ずんずん中に染み込んでいってしまう可能性がある。
（実際ループ４５０あたりからにゃんちゅうのヒゲの部分が侵食されてしまった）
エネルギー関数をもっと適切にすることでいい感じになっていくと思われるが追求はまた今度。。。