從零到一實(shí)現(xiàn)snake算法

• 1、Snake算法原理
• 2、基于曲線演化的實(shí)現(xiàn)方法
• • 2.1演化方程推導(dǎo)
  • 2.2離散化過程
  • 2.3 代碼實(shí)現(xiàn)
• 3、基于水平集的實(shí)現(xiàn)方法
• 4、討論與分析
• 源碼地址[snake](https://github.com/woshimami/snake)

1、Snake算法原理

Kass等人¹最早于1988年提出了主動輪廓模型，該方法通過構(gòu)造能量函數(shù)，從而將圖像分割轉(zhuǎn)化為求解能量泛函極值的問題，其核心思想在于，作者認(rèn)為和之前的分割方法相比，在不同的圖像解讀(image interpretation)任務(wù)中，表觀或者淺層次(low level)的圖像判讀任務(wù)應(yīng)該需要一些深層次(high level)的圖像信息，并論文[1]中進(jìn)行了舉例分析，有興趣的同學(xué)可以自行翻閱論文查看，這也和主動輪廓模型的兩大特點(diǎn)不謀而合：全局性和可拓展性。全局性也就是說在圖像處理任務(wù)進(jìn)行中，不僅要關(guān)注局部區(qū)域或者邊緣信息，同樣也要關(guān)注圖像或曲線的整體信息如曲線整體的長度和曲率等。再者針對不同的圖像特點(diǎn)，能夠設(shè)計或者添加個性化的分割策略，接下來我將根據(jù)公式對此進(jìn)一步分析。
對于圖像$I(x,y),C(s)=C(x(s),y(s))$為演化曲線，定義snake模型的能量函數(shù)為：$$E_{snake}=E_{int}+E_{ext}$$ 這里我們區(qū)分內(nèi)部能量和外部能量的依據(jù)就是看它是否和演化曲線本身的特性相關(guān)。比如說下式中，定義內(nèi)部能量項的兩部分內(nèi)容分別代表的就是曲線的長度和曲率： $$E_{int}={\int }_0^1\frac{1}{2}?(\alpha (s)|c_s(s){|}^2+\beta (s)|c_{ss}(s){|}^2)\mathrm{d}s$$ 其中公式中的兩項內(nèi)容分別表示的就是曲線的長度項和曲線的光滑程度, $\alpha (s)$和$\beta (s)$可以根據(jù)曲線上某一點(diǎn)的特征進(jìn)行個性化處理，但是大多數(shù)情況下我們一般定義它們?yōu)槌?shù),用于調(diào)節(jié)曲線長度項和光滑程度占比。
而外部能量項通常由兩部分組成構(gòu)成：$$E_{ext}=E_{img}+E_{cons}$$ $E_{img}$通常表示和圖像相關(guān)的外部能量項目，主要包含如下三種信息，
(1)

圖像的像素信息(Line functional): $$E_{line}={\int }_0^{1}I(x(s),y(s))\mathrm{d}s$$

圖像的梯度信息(edge functional)(注意負(fù)號):$$E_{edge}={\int }_0^1?|?I(x(s),y(s))\mathrm{d}s$$，

或者曲線周圍小范圍空間的梯度信息(scale space):$$E_{line}={\int }_0^1?g_{\delta }?{?}^{2}I(x(s),y(s))\mathrm{d}s$$在實(shí)際使用中圖像的梯度信息$E_{edge}$使用的較為廣泛，以及后面的實(shí)現(xiàn)過程中，我們也會用$E_{edge}$作為圖像力，負(fù)號的作用就在于,在邊緣梯度較大的情況下,整體的能量泛函越小,曲線也就更加趨向于演化到圖像的邊緣位置。

最后還有一個限制項$E_{con}$，這一項就是對曲線演化添加一定的限制作用,比如限制整體曲線到某一點(diǎn)或者某一區(qū)域$R$的距離:$$E_{c}on={\int }_0^1||C(x(s),y(s))?R||\mathrm{d}s$$具體的情況大家可以根據(jù)圖片特點(diǎn)或者需求靈活定義。

2、基于曲線演化的實(shí)現(xiàn)方法

好了，終于到了激動人心的實(shí)現(xiàn)環(huán)節(jié)了，大家是不是跟我一樣激動的直搓手，哈哈哈哈哈哈哈
開始之前我先說明兩點(diǎn)內(nèi)容：1)本次實(shí)現(xiàn)過程基于較為簡單基礎(chǔ)的能量泛函，大家可以先實(shí)現(xiàn)體驗(yàn)一邊，然后再觸類旁通，舉一反三；2)以方面需要對一些數(shù)學(xué)知識進(jìn)行回顧分析，另一方面希望大家仔細(xì)品一品從連續(xù)到離散的變化過程。

2.1演化方程推導(dǎo)

首先我們定義能量函數(shù):$$E_{snake}={\int }_0^1\frac{1}{2}?(\alpha |c_s(s){|}^2+\beta |c_{ss}(s){|}^2?|?I(x(s),y(s)){|}^2)\mathrm{d}s$$,當(dāng)能量函數(shù)最小時，就可以獲得$C(x(s),y(s))$的曲線方程，本質(zhì)上這也是一個求能量泛函極值的問題，直接根據(jù)歐拉-拉格朗日方程(變分原理)，當(dāng)能量泛函取極值時，應(yīng)滿足：$$?\alpha c^{{}^{\prime \prime }}(s)+\beta c^{{}^{\prime \prime \prime \prime }}(ss)+?E_{ext}=0$$ 可是就算推到這種程度我們以然沒辦法直接求出曲線$c(x(s),y(s))$,尤其是我們希望公式能夠以偏微分方程的形式出現(xiàn)，這樣我們就可以利用計算機(jī)進(jìn)行編程序迭代計算了，于是我們在曲線中引入了一個輔組變量$t$，則解可以表示為$c(t,s)$,隨時間變化的解$c(t,s)$應(yīng)該使得能量函數(shù)逐漸變小(推導(dǎo)過程后續(xù)會補(bǔ)充),從而可以得出:$$\frac{?c(s)}{?t}=\alpha c^{{}^{\prime \prime }}(s)?\beta c^{{}^{\prime \prime \prime \prime }}(ss)??E_{ext}$$

2.2離散化過程

接下來就是整個離散化過程了,首先,我們先明確一下概念，上式中的$s$表示什么意思，因?yàn)槲覀兛吹降膱D片都是離散的點(diǎn)，所以我們剛開定義初始化的snake曲線時，實(shí)際上定義的也就是一系列離散的點(diǎn)的集合，而曲線的演化過程實(shí)際上也就是這些固定數(shù)目的點(diǎn)的演化過程，離散化后，我們用$i$表示$s$,那假設(shè)初始化的snake曲線有$N$個點(diǎn)，那$i\in [0,N?1]$,且$x[N]=x[0],x[N+1]=x[2]....$,那么對于曲線上任一點(diǎn)$c(x(i),y(i))$：$$x^{{}^{\prime \prime }}[i]=x[i?1]?2x[i]+x[i+1]$$ $$x^{{}^{\prime \prime \prime \prime }}[i]=x[i?2]?4x[i?1]+6x[i]?4x[i+1]+x[i+2]$$對y也是同理，假設(shè)圖像$x$方向和$y$方向的梯度為$G(x,y)$(均為正值),那外部力則分別為$G_x[x,y],G_y[x,y]$,假定每次迭代后,曲線變化較小,我們用$G_{t?1}(x,y)$代替$G_t(x,y)$,那么對任一點(diǎn)$c(x(i),y(i))$就可以推導(dǎo)出離散后的公式為:
$$\frac{?x_t[i]}{?t}=\alpha (x_t[i?1]?2xt[i]+x_t[i+1])+\beta (?x_t[i?2]+4x_t[i?1]?6xt[i]+4x_t[i+1]?x_t[i+2])+G_x(x_t[i],y_t[i])$$用$\lambda$表示步長,則可以進(jìn)一步表示為：$$\frac{x_t[i]?x_{t?1}[i]}{\lambda }=\alpha (x_t[i?1]?2xt[i]+x_t[i+1])+\beta (?x_t[i?2]+4x_t[i?1]?6xt[i]+4x_t[i+1]?x_t[i+2])+G_x(x_{t?1}[i],y_{t?1}[i])$$同理， $$\frac{y_t[i]?y_{t?1}[i]}{\lambda }=\alpha (y_t[i?1]?2yt[i]+y_t[i+1])+\beta (?y_t[i?2]+4y_t[i?1]?6yt[i]+4y_t[i+1]?y_t[i+2])+G_y(y_{t?1}[i],y_{t?1}[i])$$ ,如果一條曲線上有$N$個點(diǎn)，那就對$x$和$y$分別列出$N$個式子，如果用矩陣表示的話，那就是:$$\frac{X_t?X_{t?1}}{\lambda }=A{X}_t+G_x(x_{t?1},y_{t?1})$$矩陣A則為：$$?\begin{array}{cccccccccc}?2\alpha +6\beta & \alpha +4\beta & ?\beta & 0& 0& 0& ?& ?\beta & \alpha +4\beta & \\ \alpha +4\beta & ?2\alpha +6\beta & \alpha +4\beta & ?\beta & 0& 0& ?& 0& ?\beta & \\ ?\beta & \alpha +4\beta & ?2\alpha +6\beta & \alpha +4\beta & ?\beta & 0& ?& 0& 0& \\ 0& ?\beta & \alpha +4\beta & ?2\alpha +6\beta & \alpha +4\beta & ?\beta & ?& 0& 0& \\ 0& 0& ?\beta & \alpha +4\beta & ?2\alpha +6\beta & \alpha +4\beta & ?& 0& 0& \\ ?& ?& ?& ?& ?& & & & & \\ ?\beta & 0& 0& 0& 0& 0& ?& ?2\alpha +6\beta & \alpha +4\beta & \\ \alpha +4\beta & ?\beta & 0& 0& 0& 0& ?& \alpha +4\beta & ?2\alpha +6\beta & \end{array}?$$，求解矩陣即可得出：$$X_t=(1?\lambda A{)}^{?1}{X}_{t?1}+\lambda G_x(X_{t?1},Y_{t?1})$$ $$Y_t=(1?\lambda A{)}^{?1}{Y}_{t?1}+\lambda G_y(X_{t?1},Y_{t?1})$$至此，需要的核心演化公式到手。

2.3 代碼實(shí)現(xiàn)

看完上面的離散化過程，是不是大家現(xiàn)在都胸有成竹，躍躍欲試，蠢蠢欲動啦，哈哈，接下來就讓我們一起開始愉快的coding吧。
聲明使用的依賴庫主要有：

import numpy as np import cv2 import matplotlib.pyplot as plt from pylab import* import skimage.filters as filt

代碼實(shí)現(xiàn)總共分為如下幾個步驟：

讀入圖片：
隨手用windows畫圖工具畫了個圓(簡約派)：

然后讀入圖片，注意

Image = cv2.imread('xxx.bmp',1) image = cv2.cvtColor(Image,cv2.COLOR_BGR2GRAY) img = np.array(image,dtype = np.float) plt.imshow(img,cmap = 'gray') # 讀入灰度圖 print(shape(img))

定義初始化snake曲線,注意 注意 注意，此處我們需要首先明確x，y分別表示矩陣的列和行，同時也就意味著繪圖的時候x，y分別表示x-y坐標(biāo)系的橫坐標(biāo)軸和縱坐標(biāo)軸呀，一定要品明白哈，不然會出大問題的哈(說多了都是淚)

t = np.linspace(0,2*np.pi,60,endpoint = True) y = 100+30*np.cos(t) x = 100+30*np.sin(t)

然后定義關(guān)鍵參數(shù)，部分參數(shù)值參考原論文：

alpha = 0.001 beta = 0.4 gamma = 100 # 迭代步長 sigma = 20 # 濾波用的高斯分布參數(shù) iterations = 300 # 迭代次數(shù)

計算矩陣，有沒有發(fā)現(xiàn)整個矩陣只由參數(shù)控制，而和圖片無關(guān)：

N = np.size(x) a = gamma*(2*alpha+6*beta)+1 b = gamma*(-alpha-4*beta) c = gamma*betap = np.zeros((N,N),dtype=np.float) p[0] = np.c_[a,b,c,np.zeros((1,N-5)),c,b] print(p[0].shape) for i in range(N):p[i] = np.roll(p[0],i) p = np.linalg.inv(p)

計算外部力，因?yàn)槭掷L圖像噪聲較小，很難直接通過梯度產(chǎn)生外部力，因此先對圖像進(jìn)行高斯濾波處理，注意此處的卷積核和高斯分布參數(shù)都選擇的較大，大家可以品一品為什么要這樣做：

smoothed = cv2.GaussianBlur((img-img.min()) / (img.max()-img.min()),(89,89),sigma) giy,gix = np.gradient(smoothed) gmi = (gix**2+giy**2)**0.5 gmi = (gmi - gmi.min()) / (gmi.max() - gmi.min()) Iy,Ix = np.gradient(gmi)

定義函數(shù)，確保演化曲線不超出圖像邊界，并返回整數(shù)形位置信息：

def fmax(x,y):x[x < 0] = 0y[y < 0] = 0x[x > img.shape[1]-1] = img.shape[1]-1y[y > img.shape[0]-1] = img.shape[0]-1return y.round().astype(int),x.round().astype(int)

迭代計算，并顯示圖像：

plt.plot(y,x,'.') for i in range(iterations):fex = Ix[fmax(x,y)]fey = Iy[fmax(x,y)]x = np.dot(p,x+gamma*fex)y = np.dot(p,y+gamma*fey)if i%10 ==0:plt.plot(x,y,'.') plt.show()

計算結(jié)果，撒花，撒花，撒花！

3、基于水平集的實(shí)現(xiàn)方法

本來這一段打算拖一拖的，但是看到評論區(qū)有哥們兒問到了，我就仔細(xì)找了找文獻(xiàn)，看看怎么給處理下，結(jié)果發(fā)現(xiàn)自己給自己挖了個大坑，對于原始的snake算法來講，使用level set并不是一個很可行的方案或者說需要對snake算法做一些改變，那么為什么將level set方法直接引用到snake模型上不合適呢？
1）首先，大家需要明白常見的主動輪廓模型的常見處理流程：1、構(gòu)建能量泛函; 2、極小化能量泛函(歐拉-拉格朗日方程+梯度下降法)獲得用于演化的偏微分方程; 3、離散化偏微分方程; 4、代碼實(shí)現(xiàn)與調(diào)試，一般都是從步驟1或者步驟2中引入level set，可以參見我的另一篇博客(水平集方法引入主動輪廓模型算法中的兩種方法)；2）但是，對于原始的snake算法來講，只存在從步驟1即從能量函數(shù)中引入水平集函數(shù)的可能性，不過因?yàn)槟菢涌赡芴幚砥饋硖闊?#xff0c;所以后來Caselles 等人就在1993年直接對原始snake算法就行了調(diào)整，提出了改進(jìn)的幾何活動輪廓模型，不僅引入了水平集方法還擺脫了原始snake模型對參數(shù)的依賴問題。
因此，總的來說，將level set方法引入snake是理論上的可行性的，只不過比較麻煩，這個博主最近時間比較緊張，等過一段時間閑了再想辦法給實(shí)現(xiàn)下吧(又立flag…)。

4、討論與分析

接下來就讓我們愉快的討(tu)論(cao)下snake算法吧！ 有一說一，不得不承認(rèn)，Kass等人¹當(dāng)初將圖像分割問題轉(zhuǎn)化為求解能量泛函極值的問題確實(shí)很有想象力，這個控制領(lǐng)域的李雅普諾夫函數(shù)簡直有異曲同工之妙之妙，很好奇當(dāng)時作者哪來的這么大的腦洞哈。

先說優(yōu)點(diǎn)，因?yàn)槟芰亢瘮?shù)的構(gòu)造方式，就決定了這種方法具有很好的擴(kuò)展性，這也是主動輪廓模型能夠針對具體應(yīng)用問題不斷改進(jìn)的重要原因；

但是缺點(diǎn)也很明顯：1）snake模型太過于依賴參數(shù)，針對不同的圖像圖像都需要找出"合適的參數(shù)"才能完美的分割輪廓，大家可以嘗試在代碼中調(diào)試參數(shù)值就會有所體會，甚至減少初始演化點(diǎn)的個數(shù)都會導(dǎo)致演化失敗，如下兩圖分別是相同參數(shù)，初始時設(shè)定60個演化點(diǎn)和30個演化點(diǎn)的演化結(jié)果區(qū)別，，在不改變圖像及參數(shù)的情況下，僅僅因?yàn)檠莼c(diǎn)個數(shù)減少，演化結(jié)果就迥然不同，這就凸顯了snake模型對參數(shù)的依賴性以及算法較差的魯棒性，這顯然不能符合實(shí)際應(yīng)用需求。因此后來Caselles等人²提出了幾何活動輪廓模型對snake進(jìn)行了改進(jìn)：1)曲線能夠自動地根據(jù)圖像特征如梯度而不是依賴人工設(shè)定參數(shù)演化到目標(biāo)邊緣；2)引入了水平集方法，用面的演化代替了點(diǎn)的演化，并且使得圖像能夠處理拓?fù)浣Y(jié)構(gòu)的變化。

源碼地址snake

[1] (KASS,M. snake Active contour models[J]. Proc.first Intl Conf.computer Vision, 1988, 4.) [2] (Caselles V . Geometric models for active contours[C]// International Conference on Image Processing. IEEE Computer Society, 1995.)

總結(jié)

以上是生活随笔為你收集整理的基于边缘的主动轮廓模型——从零到一用python实现snake的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。