1、Snake模型簡介(jin ji)及其編程實現Snake模型(mxng)也稱為主動(zhdng)輪廓線模型，最初由Kass等人在1987年第一屆計算機國際視覺會議上提出，一經提出就成為計算機視覺領域研究的熱點。Snake的基本思想是通過人的識別能力，在圖像中目標邊界附近確定初始輪廓線，然后對曲線進行能量最小化變形，使其鎖定在待分割目標的邊界上。Snake模型之所以能得到如此重視，是因為它將圖像目標的先驗知識（如大小、位置、形狀等）與圖像特征（灰度、梯度、紋理等）結合起來，克服了傳統圖像分割方法將二者分離的缺陷。近年來，許多文章從傳統Snake模型的能量函數構造和求解算法方面進行改進，在其基礎上

2、衍生出了許多新的Snake模型。1、Snake模型的基本原理其基本思想是依據圖像信息進行曲線（曲面）演化，使其最終找到目標物體的邊界。這種方法將分割問題轉化為最優化問題，利用閉合曲線（或曲面）形變的特定規律，定義度量閉合曲線（曲面）形變的能量函數，通過最小化能量函數使曲線（曲面）逐漸逼近圖像中目標物體的邊緣。先提供待分割圖像的一個初始輪廓的位置，并對其定義個能量函數，是輪廓沿能量降低的方向靠近。當能量函數達到最小的時候，提供的初始輪廓收斂到圖形中目標的真實輪廓。Snake能量函數是有內部能量函數和外部能量函數組成，內部能量控制輪廓的平滑性和連續性，外部能量由圖像能量和約束能量組成，控制輪廓向著

3、實際輪廓收斂，其中約束能量可根據具體的對象形態定義，使得snake具有很大的靈活性。Snake模型發展10多年來，許多學者對于經典的snake模型做了改進，提出各種改進的snake模型，其中梯度矢量流（Gradient Vector Flow, GVF）模型擴大了經典snake的外力作用范圍，加強了對目標凹輪廓邊緣的吸引力，提高了傳統的snake模型。2. 基本的Snake模型數學上，將活動輪廓表示成一條參數曲線V(s,t)=(x(s,t),y(s,t)，其中，V是曲線點的二維坐標，t是時間參數，s是弧長參數。輪廓的總弧長歸一化到1。改曲線的能量可以用能量泛函表示為E(V)=Eint(V)+E

4、ext(V)，Eint(V)是內部能量泛函， Eext(V)是外部能量泛函。曲線V在圖像的空間域運行使得E(V)最小。其中內部能量泛函定義為：EintV=1201sVs+(s)|2V2s|2ds兩個物理參量表明(biomng)曲線的物理特征： (s)是張緊系數，值越大，輪廓曲線(qxin)收縮越快；(s)是強度(qingd)系數，值越大，輪廓曲線變得僵硬而不易彎曲。外部能量泛函定義為：EextV=01e(V(s)ds其中Eext(V)是定義在圖像平面上的標量勢函數，設I(x,y)是一灰度圖像，一個可能使活動輪廓想邊界運動的外部能量函數可以設計為Eext(x,y)=-c|G(x,y)*I(x,y

5、)|2,其中，c是一個正的權因子，它控制勢的大小，是梯度算子，G(x,y)是標準差為的二維高斯函數，*是二維圖像卷積算子。G*I表明圖像進過一個高斯平滑濾波，用它代替I可降低計算梯度的噪聲。根據變分原理使能量泛函E(V)最小化的曲線V(S)滿足下述歐拉方程： 該方程可解釋為力平衡方程，它表示當輪廓達到平衡點時內力和外力的平衡。其中，前兩項參數表示內部的拉伸彈性力和彎曲剛性力，第三項參數表示曲線所受到的圖像信息外力。上式表明曲線能量最小化過程就是在圖像信息外力和曲線本身的內力作用下運動，達到平衡狀態的過程。3. 改進的Snake模型傳統的snake模型存在不足之處：首先，初始輪廓線必須接近真實的

6、圖像邊緣。其次，Snake模型有可能收斂到局部極值點，無法逼近物體內部的“凹形”輪廓。針對這些缺陷，研究者主要從能量函數構造和求解方面進行改進。（1）. 能量函數構造方面的改進。能量函數構造方面的改進具有代表性的主要有Cohen等人提出的氣球Snake模型和Xu提出的GVF Snake模型。氣球Snake模型在Snake模型中增加了氣球力，大小為常數，方向沿輪廓點的法線外方向。在氣球力的作用下，輪廓線作為一個整體膨脹或收縮，當輪廓線進入圖像能量場的作用范圍后，被吸引向感興趣區域的邊界。Cohen的氣球力模型力場吸引范圍大，克服了初始化缺陷，但為了使克服噪聲和不越過弱邊界兩者達到平衡，在怎樣選取

7、合適的氣球力權值方面仍很難把握。Cohen等人還提出一種距離勢能模型，雖然改進了初始化問題，但Snake曲線無法收斂到凹陷區域。Xu采用梯度向量場作為模型的外部能量場，產生了GVF Snake模型該模型擴大了Snake的捕獲區，初始化輪廓可以不必充分靠近真實邊緣，并能將Snake模型輪廓曲線拖向物體的深度凹陷區。但是該算法計算量太大，實時性較差。（2）. 數值計算方法的改進。Kass采用變分法解能量(nngling)函數方程，一般只能得到局部最優解，還可能出現數值解不穩定的現象，并且要求外力必須是可微的，使算法的應用受到限制。Amini采用動態規劃算法來求解局部最優解，改方法數值穩定，還可以增

8、加硬約束，放寬了對外力可微的限制，但運算速度慢，時間復雜度高。Williams在此基礎上應用貪婪算法，在保持解的穩定性情況下，將時間復雜度明顯降低。貪婪算法實現簡單，速度快，容易加入新約束，最終結點(ji din)分布比較均勻。貪婪算法雖然有上述優點，但在應用中仍然存在一些問題：沒有明確硬強制力對Snake模型的作用方法，也沒有給出具體的硬強制力模型；代表輪廓線的控制點點數固定不變，沒有給出初始輪廓點的選取方法。目前，許多有關Snake模型的文章主要是在上述幾種模型的基礎上進行改進，采用(ciyng)較多的是Xu提出的梯度向量流Snake模型和Williams的貪婪算法。4. Snake的進化

9、模型(1). McInerney 提出一種拓撲自適應snake模型（Topology Adaptive Snake,T-Snake）該算法基于仿射細胞圖像分解（Affine Cell Image Decomposition，ACID）先在待分割圖像上加以三角網格，然后在圖像區域的適當位置做一條初始曲線，最后取曲線與網格的交點作為snake的初始離散點，其第i個snake的離散點的坐標為其中，相鄰兩點，之間由一條彈性樣條連接而成由于T-Snake模型可借助三角形網格和網格點的特征函數來確定邊界三角形，可促使snake模型演化過程中的分裂和合并，從而保證了其具有能夠處理拓撲結果復雜圖像的能力，因此

10、能夠很好的滿足醫學圖像拓撲結果復雜的特點。此算法用于腦部MR切片有良好的性能。(2). 雙T-Snake模型雙T-Snake模型（Dual-T-Snakes）是在T-Snake模型的基礎上產生的，其主要思想是采用內外兩個初始輪廓，其中一個輪廓從目標外向內收縮和分裂，另一個輪廓從目標內部向外膨脹，兩個初始輪廓可以(ky)離目標邊界較遠，迭代的過程中對能量較大的輪廓增加驅動力，使其靠近與之相對應的輪廓，直到連個輪廓收斂到同一個為止(3). Loop Snake 模型(mxng)Loop Snake模型是一種加強了拓撲控制的T-Snake模型，這種方法的關鍵(gunjin)集中在曲線的每一步進化中都

11、要形成循環，其基本思想是，確保圖像輪廓曲線精確地線性地映射到適當的分類中，然后在額外的記過loop-Tree的幫助下，盡可能少的時間內運用已經被snake探究的循環來決定是否進行區域劃分，這種模型的實質是對T-Snake模型的一種改進。由于加強了拓撲控制，使得Loop Snake模型既可以忽略背景中強噪聲又可以在演化過程中進行多次分裂。(4). 連續snake模型在Snake模型中，輪廓曲線由一條給定容許誤差范圍的光滑曲線組成，相對于離散snake來說，連續snake模型所需要的控制點少，比離散的更具優越性。(5). B-Snake模型B-Snake模型是通過B樣條曲線來定義的，其輪廓曲線由各

12、曲線段光滑相連而成，每一個曲線段都是由一個給定次數多項式表示，這種多項式是B樣條曲度函數的一種線性組合，并以控制點為系數。在有些B-Snake模型中并沒有明確應用內部能量，這是因為B樣條本身就含有內部能量，snake輪廓曲線只受外力影響著圖像邊緣移動。可用于對圖像切片分割區域的描述與跟蹤而用于器官的三維重建。5. 基本Snake模型的程序實現使用matlab實現基本的snake模型，主程序如下（主程序中調用的部分函數將在附錄中給出）：cd .; s = cd; s = s, /snake; path(s, path); cd examples; help tradition_ex; % = E

13、xample 1: U-shape object = % 讀取64x64 U形圖像(t xin) I,map = rawread(./images/U64.pgm); %邊緣(binyun)檢測(jin c) disp( 邊緣檢測 .); f = 1 - I/255; %圖像歸一化 f0 = gaussianBlur(f,1); %高斯函數模糊處理 % note: snake potential is the negative of edge map disp( 計算外部力 .); px,py = gradient(f0); %其中px為其水平方向上的梯度，py為其垂直方向上的梯度 % 顯示結

14、果 figure(1); subplot(121); imdisp(-f); title(snake potential); subplot(122); quiver(px,py); %在x-y平面上繪制向量（px,py） axis(image, off, ij); % 確定圖像顯示坐標 title(traditional force); % snake變換 disp( ); disp( 請按任意鍵開始變換：); pause; figure(1); subplot(121); %將圖一分為二，并在左側顯示一張。 colormap(gray(64); %得出一個64色的灰度圖image(1-f)

15、+1)*40); %產生正方形坐標系 axis(square, off); disp(); disp(. 捕獲范圍正在縮小.) t = 0:0.05:6.28; x = 32 + 17*cos(t); y = 32 + 17*sin(t); x,y = snakeinterp(x,y,3,1); %插入或移動節點處理 snakedisp(x,y,r) pause(1); disp(. 按-C 終止(zhngzh)程序！); for i=1:100, x,y = snakedeform(x,y,0.05,0,1,4,px,py,5); x,y = snakeinterp(x,y,3,1); sn

16、akedisp(x,y,r) title(變換中., 迭代(di di)次數 = num2str(i*5) pause(0.5); end disp( ); disp( 請按任意鍵顯示(xinsh)最終結果：); pause; subplot(121); colormap(gray(64); image(1-f)+1)*40); axis(square, off); snakedisp(x,y,r) title(最終結果： 迭代次數 = num2str(i*5);程序運行結果如下：圖1. 原始圖像及圖像的傳統力場圖2. Snake模型(mxng)變換中圖3. Snake模型變換(binhun)

17、的最終結果由實驗結果可以看出，傳統高斯力Snake模型(mxng)在迭代求解過程中，很好的使初始輪廓線向U形圖的邊緣靠近，最終得到U形圖的邊緣輪廓；同時，從圖中可以看出，可變輪廓線無法捕捉U形的底部，輪廓沒有收縮到U形圖的凹陷部分。附錄：rawread.mfunction X,map = rawread(filename,n,m)% RAWREAD 讀入可移植的點陣圖文件, 或者原始文件.%RAWREAD(imagefile.raw, xsize, ysize) 讀入raw格式的圖像文件%RAWREAD(imagefile.pgm) 讀入一 pgm 格式的圖像文件%X,map = RAWREA

18、D(imagefile.raw) 同時返回一個灰度圖和一個彩色圖%X,map = rawread(imagefile.raw,sx,sy);%or X,map = rawread(imagefile.pgm);% image(X)% colormap(map) 顯示出正確(zhngqu)的色彩。dot = max(find(filename = .);suffix = filename(dot+1:dot+3);if strcmp(suffix,pgm) | strcmp(suffix,raw)%strcmp比較(bjio)字符 disp(sprintf(nopens %s filen,fil

19、ename); fp = fopen(filename,rb,b); % Big-endian byte order. if (fp0) error(Cannot open filename .); end if strcmp(suffix,pgm) % Read and crack the header(打開(d ki)標題） head = fread(fp,2,uchar); % pgm magic number : P5 if strcmp(head,P5), fprintf(1,n Magic Number : %sn,head); else fprintf(1,n Bad Magic

20、 Number : %sn,head); error(cannot continue this way, good bye cruel world); end c = fread(fp,1,uchar); %reads the carriage return（回車符） separating P5 from the creator precreator = fread(fp,1,uchar); % look for a # character preceeding a creator signature if precreator = #, c = setstr(20); % any chara

21、cter except carriage return cr = setstr(10); % defines a carriage return while c = cr, c = fread(fp,1,uchar); creator = creator,c; end; fprintf(1,n creator : %sn,creator); else fprintf(n No creator signaturen); fseek(fp,-1,cof); % return one char before end;end if nargin dmax, 則要在i與i+1之間添加一個點。% 如果(r

22、gu)d(i,i+1)dmin, 則要移動i或者i+1。x = x(:); y = y(:); %分別將x,y以列向量的形式輸出N = length(x);d = abs(x(2:N 1)- x(:) + abs(y(2:N 1)- y(:);% remove the points which distance to neighbor points is shorter than dminIDX = (ddmax);z = snakeindex(IDX);%本例中z=1:N,輸出z=1:Np = 1:N+1;xi = interp1(p,x;x(1),z);%interp1(庫函數）；這里x;

23、x(1)的表述是對的yi = interp1(p,y;y(1),z);% 其表示的是一個列向量N = length(xi);d = abs(xi(2:N 1)- xi(:) + abs(yi(2:N 1)- yi(:);while (max(d)dmax), IDX = (ddmax); z = snakeindex(IDX); p = 1:N+1; xi = interp1(p,xi;xi(1),z); yi = interp1(p,yi;yi(1),z); N = length(xi); d = abs(xi(2:N 1)- xi(:) + abs(yi(2:N 1)- yi(:);end

24、snakedisp.mfunction snakedisp(x,y,style)hold on% 轉換(zhunhun)為列向量(xingling)x = x(:); y = y(:);if nargin = 3 plot(x;x(1,1),y;y(1,1),style); hold offelse disp(snakedisp.m: The input parameter is not correct!); end snakedeform.mfunction x,y = snakedeform(x,y,alpha,beta,gamma,kappa,fx,fy,ITER)generates the parameters for snake% alpha=0.05 beta=0.1 gamma=0.5 ITER=5N = length(x);alpha = alpha* ones(1,N); beta = beta*ones(1,N);%產生對角(du jio)向量alpham1 = alpha(2:N) alpha(1);alphap1 = alpha(N) alpha(1:N-1);betam1