參考博客：

https://www.mathworks.com/matlabcentral/fileexchange/19567-active-contour-segmentation

數(shù)字圖像處理-圖像分割：Snake主動(dòng)輪廓模型 Matlab代碼及運(yùn)行結(jié)果

簡(jiǎn)介

在“圖像分割之（一）概述”中咱們簡(jiǎn)單了解了目前主流的圖像分割方法。下面咱們主要學(xué)習(xí)下基于能量泛函的分割方法。這里學(xué)習(xí)下Snake模型簡(jiǎn)單的知識(shí)，Level Set（水平集）模型會(huì)在后面的博文中說到。

基于能量泛函的分割方法：

該類方法主要指的是活動(dòng)輪廓模型（active contour model）以及在其基礎(chǔ)上發(fā)展出來的算法，其基本思想是使用連續(xù)曲線來表達(dá)目標(biāo)邊緣，并定義一個(gè)能量泛函使得其自變量包括邊緣曲線，因此分割過程就轉(zhuǎn)變?yōu)榍蠼饽芰糠汉淖钚≈档倪^程，一般可通過求解函數(shù)對(duì)應(yīng)的歐拉(Euler．Lagrange)方程來實(shí)現(xiàn)，能量達(dá)到最小時(shí)的曲線位置就是目標(biāo)的輪廓所在。

主動(dòng)輪廓線模型是一個(gè)自頂向下定位圖像特征的機(jī)制，用戶或其他自動(dòng)處理過程通過事先在感興趣目標(biāo)附近放置一個(gè)初始輪廓線，在內(nèi)部能量（內(nèi)力）和外部能量（外力）的作用下變形外部能量吸引活動(dòng)輪廓朝物體邊緣運(yùn)動(dòng)，而內(nèi)部能量保持活動(dòng)輪廓的光滑性和拓?fù)湫?#xff0c;當(dāng)能量達(dá)到最小時(shí)，活動(dòng)輪廓收斂到所要檢測(cè)的物體邊緣。

一、曲線演化理論

曲線演化理論在水平集中運(yùn)用到，但我感覺在主動(dòng)輪廓線模型的分割方法中，這個(gè)知識(shí)是公用的，所以這里我們簡(jiǎn)單了解下。

曲線可以簡(jiǎn)單的分為幾種：

曲線存在曲率，曲率有正有負(fù)，于是在法向曲率力的推動(dòng)下，曲線的運(yùn)動(dòng)方向之間有所不同：有些部分朝外擴(kuò)展，而有些部分則朝內(nèi)運(yùn)動(dòng)。這種情形如下圖所示。圖中藍(lán)色箭頭處的曲率為負(fù)，而綠色箭頭處的曲率為正。

簡(jiǎn)單曲線在曲率力（也就是曲線的二次導(dǎo)數(shù)）的驅(qū)動(dòng)下演化所具有的一種非常特殊的數(shù)學(xué)性質(zhì)是：一切簡(jiǎn)單曲線，無論被扭曲得多么嚴(yán)重，只要還是一種簡(jiǎn)單曲線，那么在曲率力的推動(dòng)下最終將退化成一個(gè)圓，然后消逝（可以想象下，圓的所有點(diǎn)的曲率力都向著圓心，所以它將慢慢縮小，以致最后消逝）。
描述曲線幾何特征的兩個(gè)重要參數(shù)是單位法矢和曲率，單位法矢描述曲線的方向，曲率則表述曲線彎曲的程度。曲線演化理論就是僅利用曲線的單位法矢和曲率等幾何參數(shù)來研究曲線隨時(shí)間的變形。曲線的演變過程可以認(rèn)為是表示曲線在作用力 F 的驅(qū)動(dòng)下，朝法線方向 N 以速度 v 演化。而速度是有正負(fù)之分的，所以就有如果速度 v 的符號(hào)為負(fù)，表示活動(dòng)輪廓演化過程是朝外部方向的，如為正，則表示朝內(nèi)部方向演化，活動(dòng)曲線是單方向演化的，不可能同時(shí)往兩個(gè)方向演化。

所以曲線的演變過程，就是不同力在曲線上的作用過程，力也可以表達(dá)為能量。世界萬物都趨向于能量最小而存在。因?yàn)榇藭r(shí)它是最平衡的，消耗最小的（不知理解對(duì)不？）。那么在圖像分割里面，我們目標(biāo)是把目標(biāo)的輪廓找到，那么在目標(biāo)的輪廓這個(gè)地方，整個(gè)輪廓的能量是最小的，那么曲線在圖像任何一個(gè)地方，都可以因?yàn)榱Τ@個(gè)能量最小的輪廓演變，當(dāng)演變到目標(biāo)的輪廓的時(shí)候，因?yàn)槟芰孔钚?#xff0c;力平衡了，速度為0了，也就不動(dòng)了，這時(shí)候目標(biāo)就被我們分割出來了。

那現(xiàn)在關(guān)鍵就在于：1）這個(gè)輪廓我們?cè)趺幢硎?#xff1b;2）這些力怎么構(gòu)造，構(gòu)造哪些力才可以讓目標(biāo)輪廓這個(gè)地方的能量最小？

這兩個(gè)問題的描述和解決就衍生出了很多的基于主動(dòng)輪廓線模型的分割方法。第一個(gè)問題的回答，就形成了兩大流派：如果這個(gè)輪廓是參數(shù)表示的，那么就是參數(shù)活動(dòng)輪廓模型（parametric active contour model），典型為snake模型，如果這個(gè)輪廓是幾何表示的，那么就是幾何活動(dòng)輪廓模型（geometric active contour model），即水平集方法（Level Set），它是把二維的輪廓嵌入到三維的曲面的零水平面來表達(dá)的（可以理解為一座山峰的等高線，某個(gè)等高線把山峰切了，這個(gè)高度山峰的水平形狀就出來了，也就是輪廓了），所以低維的演化曲線或曲面，表達(dá)為高維函數(shù)曲面的零水平集的間接表達(dá)形式（這個(gè)輪廓的變化，直觀上我們就可以調(diào)整山峰的形狀或者調(diào)整登高線的高度來得到）。

那對(duì)于第二個(gè)問題，是兩大流派都遇到的問題，是他們都需要解決的最關(guān)鍵的問題。哪些力才可以達(dá)到分割的目標(biāo)呢？這將在后面聊到。

二、Snakes模型

自1987年Kass提出Snakes模型以來，各種基于主動(dòng)輪廓線的圖像分割理解和識(shí)別方法如雨后春筍般蓬勃發(fā)展起來。Snakes模型的基本思想很簡(jiǎn)單，它以構(gòu)成一定形狀的一些控制點(diǎn)為模板（輪廓線），通過模板自身的彈性形變，與圖像局部特征相匹配達(dá)到調(diào)和，即某種能量函數(shù)極小化，完成對(duì)圖像的分割。再通過對(duì)模板的進(jìn)一步分析而實(shí)現(xiàn)圖像的理解和識(shí)別。

簡(jiǎn)單的來講，SNAKE模型就是一條可變形的參數(shù)曲線及相應(yīng)的能量函數(shù)，以最小化能量目標(biāo)函數(shù)為目標(biāo)，控制參數(shù)曲線變形，具有最小能量的閉合曲線就是目標(biāo)輪廓。

構(gòu)造Snakes模型的目的是為了調(diào)和上層知識(shí)和底層圖像特征這一對(duì)矛盾。無論是亮度、梯度、角點(diǎn)、紋理還是光流，所有的圖像特征都是局部的。所謂局部性就是指圖像上某一點(diǎn)的特征只取決于這一點(diǎn)所在的鄰域，而與物體的形狀無關(guān)。但是人們對(duì)物體的認(rèn)識(shí)主要是來自于其外形輪廓。如何將兩者有效地融合在一起正是Snakes模型的長(zhǎng)處。Snakes模型的輪廓線承載了上層知識(shí)，而輪廓線與圖像的匹配又融合了底層特征。這兩項(xiàng)分別表示為Snakes模型中能量函數(shù)的內(nèi)部力和圖像力。

模型的形變受到同時(shí)作用在模型上的許多不同的力所控制，每一種力所產(chǎn)生一部分能量，這部分能量表示為活動(dòng)輪廓模型的能量函數(shù)的一個(gè)獨(dú)立的能量項(xiàng)。

Snake模型首先需要在感興趣區(qū)域的附近給出一條初始曲線，接下來最小化能量泛函，讓曲線在圖像中發(fā)生變形并不斷逼近目標(biāo)輪廓。

Kass等提出的原始Snakes模型由一組控制點(diǎn)：v(s)=[x(s), y(s)] s∈[0, 1] 組成，這些點(diǎn)首尾以直線相連構(gòu)成輪廓線。其中x(s)和y(s)分別表示每個(gè)控制點(diǎn)在圖像中的坐標(biāo)位置。 s 是以傅立葉變換形式描述邊界的自變量。在Snakes的控制點(diǎn)上定義能量函數(shù)（反映能量與輪廓之間的關(guān)系）：
其中第1項(xiàng)稱為彈性能量是v的一階導(dǎo)數(shù)的模，第2項(xiàng)稱為彎曲能量，是v的二階導(dǎo)數(shù)的模，第3項(xiàng)是外部能量（外部力），在基本Snakes模型中一般只取控制點(diǎn)或連線所在位置的圖像局部特征例如梯度：
也稱圖像力。（當(dāng)輪廓C靠近目標(biāo)圖像邊緣，那么C的灰度的梯度將會(huì)增大，那么上式的能量最小，由曲線演變公式知道該點(diǎn)的速度將變?yōu)?，也就是停止運(yùn)動(dòng)了。這樣，C就停在圖像的邊緣位置了，也就完成了分割。那么這個(gè)的前提就是目標(biāo)在圖像中的邊緣比較明顯了，否則很容易就越過邊緣了。）

彈性能量和彎曲能量合稱內(nèi)部能量（內(nèi)部力），用于控制輪廓線的彈性形變，起到保持輪廓連續(xù)性和平滑性的作用。而第三項(xiàng)代表外部能量，也被稱為圖像能量，表示變形曲線與圖像局部特征吻合的情況。內(nèi)部能量?jī)H僅跟snake的形狀有關(guān)，而跟圖像數(shù)據(jù)無關(guān)。而外部能量?jī)H僅跟圖像數(shù)據(jù)有關(guān)。在某一點(diǎn)的α和β的值決定曲線可以在這一點(diǎn)伸展和彎曲的程度。

最終對(duì)圖像的分割轉(zhuǎn)化為求解能量函數(shù)Etotal(v)極小化（最小化輪廓的能量）。在能量函數(shù)極小化過程中，彈性能量迅速把輪廓線壓縮成一個(gè)光滑的圓，彎曲能量驅(qū)使輪廓線成為光滑曲線或直線，而圖像力則使輪廓線向圖像的高梯度位置靠攏。基本Snakes模型就是在這3個(gè)力的聯(lián)合作用下工作的。

因?yàn)閳D像上的點(diǎn)都是離散的，所以我們用來優(yōu)化能量函數(shù)的算法都必須在離散域里定義。所以求解能量函數(shù)Etotal(v)極小化是一個(gè)典型的變分問題（微分運(yùn)算中，自變量一般是坐標(biāo)等變量，因變量是函數(shù)；變分運(yùn)算中，自變量是函數(shù)，因變量是函數(shù)的函數(shù)，即數(shù)學(xué)上所謂的泛函。對(duì)泛函求極值的問題，數(shù)學(xué)上稱之為變分法）。

在離散化條件（數(shù)字圖像）下，由歐拉方程可知最終問題的答案等價(jià)于求解一組差分方程：（歐拉方程是泛函極值條件的微分表達(dá)式，求解泛函的歐拉方程，即可得到使泛函取極值的駐函數(shù)，將變分問題轉(zhuǎn)化為微分問題。）

記外部力 F = ?? P， Kass等將上式離散化后，對(duì)x(s)和y(s)分別構(gòu)造兩個(gè)五對(duì)角陣的線性方程組，通過迭代計(jì)算進(jìn)行求解。在實(shí)際應(yīng)用中一般先在物體周圍手動(dòng)點(diǎn)出控制點(diǎn)作為Snakes模型的起始位置，然后對(duì)能量函數(shù)迭代求解。

Reference：

李天慶等，Snake模型綜述，計(jì)算機(jī)工程，2005，第31卷 第9期

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

% ========================================================================= % Snakes：Active Contour Models % ========================================================================= % By gujinjin 2012/12/10-12 Sunny % 基于KASS等的論文思想 % 參考文獻(xiàn)： % [1] KASS etc.Snakes：Active Contour Models % [2] CSDN 博客 - Author：樂不思蜀Tone % [3] Ritwik Kumar(Harvard University),D.Kroon（Twente University)的程序 % [4] 《數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)》 % ------ clc;clf;clear all; % ========================================================================= % 獲取手動(dòng)取點(diǎn)坐標(biāo) % ========================================================================= % 讀取顯示圖像 %I = imread('Coronary_MPR.jpg'); I = imread('test3.png'); % 轉(zhuǎn)化為雙精度型 %I = im2double(I); % 若為彩色，轉(zhuǎn)化為灰度 if(size(I,3)==3), I=rgb2gray(I); end figure(1),imshow(I); %--------------------------- % 高斯濾波 %--------------------------- sigma=1.0; % 創(chuàng)建特定形式的二維高斯濾波器H H = fspecial('gaussian',ceil(3*sigma), sigma); % 對(duì)圖像進(jìn)行高斯濾波,返回和I等大小矩陣 Igs = filter2(H,I,'same'); %--------------------------- % 獲取Snake的點(diǎn)坐標(biāo) %--------------------------- figure(2),imshow(Igs); x=[];y=[];c=1;N=100; %定義取點(diǎn)個(gè)數(shù)c,上限N % 獲取User手動(dòng)取點(diǎn)的坐標(biāo) % [x,y]=getpts while c<N [xi,yi,button]=ginput(1); % 獲取坐標(biāo)向量 x=[x xi]; y=[y yi]; hold on % text(xi,yi,'o','FontSize',10,'Color','red'); plot(xi,yi,'ro'); % 若為右擊，則停止循環(huán) if(button==3), break; end c=c+1; end % 將第一個(gè)點(diǎn)復(fù)制到矩陣最后，構(gòu)成Snake環(huán) xy = [x;y]; c=c+1; xy(:,c)=xy(:,1); % 樣條曲線差值 t=1:c; ts = 1:0.1:c; xys = spline(t,xy,ts); xs = xys(1,:); ys = xys(2,:); % 樣條差值效果 hold on temp=plot(x(1),y(1),'ro',xs,ys,'b.'); legend(temp,'原點(diǎn)','插值點(diǎn)'); % ========================================================================= % Snakes算法實(shí)現(xiàn)部分 % ========================================================================= NIter =1000; % 迭代次數(shù) alpha=0.2; beta=0.2; gamma = 1; kappa = 0.1; wl = 0; we=0.4; wt=0; [row col] = size(Igs); % 圖像力-線函數(shù) Eline = Igs; % 圖像力-邊函數(shù) [gx,gy]=gradient(Igs); Eedge = -1*sqrt((gx.*gx+gy.*gy)); % 圖像力-終點(diǎn)函數(shù) % 卷積是為了求解偏導(dǎo)數(shù)，而離散點(diǎn)的偏導(dǎo)即差分求解 m1 = [-1 1]; m2 = [-1;1]; m3 = [1 -2 1]; m4 = [1;-2;1]; m5 = [1 -1;-1 1]; cx = conv2(Igs,m1,'same'); cy = conv2(Igs,m2,'same'); cxx = conv2(Igs,m3,'same'); cyy = conv2(Igs,m4,'same'); cxy = conv2(Igs,m5,'same'); for i = 1:row for j= 1:col Eterm(i,j) = (cyy(i,j)*cx(i,j)*cx(i,j) -2 *cxy(i,j)*cx(i,j)*cy(i,j) + cxx(i,j)*cy(i,j)*cy(i,j))/((1+cx(i,j)*cx(i,j) + cy(i,j)*cy(i,j))^1.5); end end %figure(3),imshow(Eterm); %figure(4),imshow(abs(Eedge)); % 外部力 Eext = Eimage + Econ Eext = wl*Eline + we*Eedge + wt*Eterm; % 計(jì)算梯度 [fx,fy]=gradient(Eext); xs=xs'; ys=ys'; [m n] = size(xs); [mm nn] = size(fx); % 計(jì)算五對(duì)角狀矩陣 % 附錄: 公式（14） b(i)表示vi系數(shù)(i=i-2 到 i+2) b(1)=beta; b(2)=-(alpha + 4*beta); b(3)=(2*alpha + 6 *beta); b(4)=b(2); b(5)=b(1); A=b(1)*circshift(eye(m),2); A=A+b(2)*circshift(eye(m),1); A=A+b(3)*circshift(eye(m),0); A=A+b(4)*circshift(eye(m),-1); A=A+b(5)*circshift(eye(m),-2); % 計(jì)算矩陣的逆 [L U] = lu(A + gamma.* eye(m)); Ainv = inv(U) * inv(L); % ========================================================================= % 畫圖部分 % ========================================================================= %text(10,10,'+','FontName','Time','FontSize',20,'Color','red'); % 迭代計(jì)算 figure(3) for i=1:NIter; ssx = gamma*xs - kappa*interp2(fx,xs,ys); ssy = gamma*ys - kappa*interp2(fy,xs,ys); % 計(jì)算snake的新位置 xs = Ainv * ssx; ys = Ainv * ssy; % 顯示snake的新位置 imshow(I); hold on; plot([xs; xs(1)], [ys; ys(1)], 'r-'); hold off; pause(0.001) end

運(yùn)行效果

總結(jié)

以上是生活随笔為你收集整理的Active Contour Models 主动轮廓模型的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

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