天津大學(xué)研究生數(shù)字圖像處理作業(yè)-Snake模型(共13頁)

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

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

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

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

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

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

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

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

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

10、能夠很好的滿足醫(yī)學(xué)圖像拓?fù)浣Y(jié)果復(fù)雜的特點(diǎn)。此算法用于腦部MR切片有良好的性能。(2). 雙T-Snake模型雙T-Snake模型（Dual-T-Snakes）是在T-Snake模型的基礎(chǔ)上產(chǎn)生的，其主要思想是采用內(nèi)外兩個(gè)初始輪廓，其中一個(gè)輪廓從目標(biāo)外向內(nèi)收縮和分裂，另一個(gè)輪廓從目標(biāo)內(nèi)部向外膨脹，兩個(gè)初始輪廓可以(ky)離目標(biāo)邊界較遠(yuǎn)，迭代的過程中對(duì)能量較大的輪廓增加驅(qū)動(dòng)力，使其靠近與之相對(duì)應(yīng)的輪廓，直到連個(gè)輪廓收斂到同一個(gè)為止(3). Loop Snake 模型(mxng)Loop Snake模型是一種加強(qiáng)了拓?fù)淇刂频腡-Snake模型，這種方法的關(guān)鍵(gunjin)集中在曲線的每一步進(jìn)化中都

11、要形成循環(huán)，其基本思想是，確保圖像輪廓曲線精確地線性地映射到適當(dāng)?shù)姆诸愔校缓笤陬~外的記過loop-Tree的幫助下，盡可能少的時(shí)間內(nèi)運(yùn)用已經(jīng)被snake探究的循環(huán)來決定是否進(jìn)行區(qū)域劃分，這種模型的實(shí)質(zhì)是對(duì)T-Snake模型的一種改進(jìn)。由于加強(qiáng)了拓?fù)淇刂?，使得Loop Snake模型既可以忽略背景中強(qiáng)噪聲又可以在演化過程中進(jìn)行多次分裂。(4). 連續(xù)snake模型在Snake模型中，輪廓曲線由一條給定容許誤差范圍的光滑曲線組成，相對(duì)于離散snake來說，連續(xù)snake模型所需要的控制點(diǎn)少，比離散的更具優(yōu)越性。(5). B-Snake模型B-Snake模型是通過B樣條曲線來定義的，其輪廓曲線由各

12、曲線段光滑相連而成，每一個(gè)曲線段都是由一個(gè)給定次數(shù)多項(xiàng)式表示，這種多項(xiàng)式是B樣條曲度函數(shù)的一種線性組合，并以控制點(diǎn)為系數(shù)。在有些B-Snake模型中并沒有明確應(yīng)用內(nèi)部能量，這是因?yàn)锽樣條本身就含有內(nèi)部能量，snake輪廓曲線只受外力影響著圖像邊緣移動(dòng)?？捎糜趯?duì)圖像切片分割區(qū)域的描述與跟蹤而用于器官的三維重建。5. 基本Snake模型的程序?qū)崿F(xiàn)使用matlab實(shí)現(xiàn)基本的snake模型，主程序如下（主程序中調(diào)用的部分函數(shù)將在附錄中給出）：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)檢測(cè)(jin c) disp( 邊緣檢測(cè) .); f = 1 - I/255; %圖像歸一化 f0 = gaussianBlur(f,1); %高斯函數(shù)模糊處理 % note: snake potential is the negative of edge map disp( 計(jì)算外部力 .); px,py = gradient(f0); %其中px為其水平方向上的梯度，py為其垂直方向上的梯度 % 顯示結(jié)

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

15、+1)*40); %產(chǎn)生正方形坐標(biāo)系 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); %插入或移動(dòng)節(jié)點(diǎn)處理 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)次數(shù) = num2str(i*5) pause(0.5); end disp( ); disp( 請(qǐng)按任意鍵顯示(xinsh)最終結(jié)果：); pause; subplot(121); colormap(gray(64); image(1-f)+1)*40); axis(square, off); snakedisp(x,y,r) title(最終結(jié)果： 迭代次數(shù) = num2str(i*5);程序運(yùn)行結(jié)果如下：圖1. 原始圖像及圖像的傳統(tǒng)力場(chǎng)圖2. Snake模型(mxng)變換中圖3. Snake模型變換(binhun)

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

18、D(imagefile.raw) 同時(shí)返回一個(gè)灰度圖和一個(gè)彩色圖%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)標(biāo)題） 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之間添加一個(gè)點(diǎn)。% 如果(r

22、gu)d(i,i+1)dmin, 則要移動(dòng)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(庫函數(shù)）；這里x;

23、x(1)的表述是對(duì)的yi = interp1(p,y;y(1),z);% 其表示的是一個(gè)列向量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% 轉(zhuǎn)換(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);%產(chǎn)生對(duì)角(du jio)向量alpham1 = alpha(2:N) alpha(1);alphap1 = alpha(N) alpha(1:N-1);betam1