從零到一實現snake算法

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

1、Snake算法原理

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

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

圖像的梯度信息(edge functional)(注意負號):$$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$$在實際使用中圖像的梯度信息$E_{edge}$使用的較為廣泛，以及后面的實現過程中，我們也會用$E_{edge}$作為圖像力，負號的作用就在于,在邊緣梯度較大的情況下,整體的能量泛函越小,曲線也就更加趨向于演化到圖像的邊緣位置。

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

2、基于曲線演化的實現方法

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

2.1演化方程推導

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

2.2離散化過程

接下來就是整個離散化過程了,首先,我們先明確一下概念，上式中的$s$表示什么意思，因為我們看到的圖片都是離散的點，所以我們剛開定義初始化的snake曲線時，實際上定義的也就是一系列離散的點的集合，而曲線的演化過程實際上也就是這些固定數目的點的演化過程，離散化后，我們用$i$表示$s$,那假設初始化的snake曲線有$N$個點，那$i\in [0,N?1]$,且$x[N]=x[0],x[N+1]=x[2]....$,那么對于曲線上任一點$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也是同理，假設圖像$x$方向和$y$方向的梯度為$G(x,y)$(均為正值),那外部力則分別為$G_x[x,y],G_y[x,y]$,假定每次迭代后,曲線變化較小,我們用$G_{t?1}(x,y)$代替$G_t(x,y)$,那么對任一點$c(x(i),y(i))$就可以推導出離散后的公式為:
$$\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$表示步長,則可以進一步表示為：$$\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$個點，那就對$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 代碼實現

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

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

代碼實現總共分為如下幾個步驟：

讀入圖片：
隨手用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坐標系的橫坐標軸和縱坐標軸呀，一定要品明白哈，不然會出大問題的哈(說多了都是淚)

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

然后定義關鍵參數，部分參數值參考原論文：

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

計算矩陣，有沒有發現整個矩陣只由參數控制，而和圖片無關：

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)

計算外部力，因為手繪圖像噪聲較小，很難直接通過梯度產生外部力，因此先對圖像進行高斯濾波處理，注意此處的卷積核和高斯分布參數都選擇的較大，大家可以品一品為什么要這樣做：

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)

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

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()

計算結果，撒花，撒花，撒花！

3、基于水平集的實現方法

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

4、討論與分析

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

先說優點，因為能量函數的構造方式，就決定了這種方法具有很好的擴展性，這也是主動輪廓模型能夠針對具體應用問題不斷改進的重要原因；

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

源碼地址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.)

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。