最近，在做一個(gè)關(guān)于遙感信息提取的工作，發(fā)現(xiàn)一個(gè)Snake模型可以應(yīng)用于圖像分割，初學(xué)Python和Matlab，因此使用這兩門(mén)語(yǔ)言進(jìn)行了一定的了解和深入學(xué)習(xí)，這也是我的第一篇博客，歡迎一起討論。

Snake算法原理

這個(gè)模型是基于“分割線是有能量的”這一基礎(chǔ)進(jìn)行的研究，當(dāng)輪廓線穩(wěn)定到能量最小處時(shí)，就是分割的終點(diǎn)，因此需要輸入一個(gè)初始輪廓線用以計(jì)算能量，在運(yùn)算過(guò)程中，能量是朝著減小的方向進(jìn)行的。

相關(guān)公式如下：

Snake的邊界由一組控制點(diǎn)表示，s是以傅立葉變換形式描述邊界的自變量

Snake能量函數(shù)是有內(nèi)部能量函數(shù)和外部能量函數(shù)組成。

• 內(nèi)部能量控制輪廓的平滑性和連續(xù)性，由彈性勢(shì)能和彎曲勢(shì)能組成，其中彈性勢(shì)能為v一階導(dǎo)數(shù)的模，彎曲勢(shì)能為v二階導(dǎo)數(shù)的模。

• 外部能量（又稱圖像力）由線性能量Eline、邊緣能量Eedge和終端能量Eterminal三者加權(quán)相加而成。

• 設(shè)定Eline為圖像強(qiáng)度，Eedge為亮度的梯度變化

• 定義C(x,y)為高斯濾波（可否不進(jìn)行濾波操作？）后的圖像，θ是(x,y)處的梯度角度

可以得到Eterminal

關(guān)于優(yōu)化：作者用數(shù)值方法進(jìn)行推導(dǎo)，得到了一個(gè)迭代的方程式（這個(gè)我不太清楚），主要來(lái)源于以下文章

老笨妞 snake算法總結(jié)

田問(wèn)渠Carlnait 主動(dòng)輪廓模型——Snake分割算法 matlab源碼

在此也借鑒了許多的博客內(nèi)容，在這一部分不過(guò)多贅述。

?

由于輪廓點(diǎn)信息都是投影坐標(biāo)，因此需要將投影坐標(biāo)先轉(zhuǎn)換為行列式坐標(biāo)。在轉(zhuǎn)換時(shí)，matlab數(shù)組序號(hào)起始是1，而python數(shù)組起始是0，因此在計(jì)算x，y的行列式坐標(biāo)時(shí)，應(yīng)相應(yīng)加減0.5

Matlab 進(jìn)行計(jì)算（代碼主要來(lái)源于以上所列博客及matlab官網(wǎng)，并根據(jù)實(shí)際情況進(jìn)行修改）

1.打開(kāi)影像與分割線

Former=shaperead("修改前線/6_former.shp"); [Igs,R]=geotiffread("校正背景/6.tif"); Igs=double(Igs); x=(Former.X-R.XWorldLimits(1))/R.CellExtentInWorldX + 0.5; y=(Former.Y-R.YWorldLimits(2))/-R.CellExtentInWorldY + 0.5;%第y行 第x列 if(length(x)==length(y))c=length(x);for i=1:c-1x_temporary(i)=x(i);y_temporary(i)=y(i);endx=x_temporary;y=y_temporary; % x(:,c)=x(:,1);% shp文件中已經(jīng)是環(huán)狀點(diǎn)集合，因此不需要把第一個(gè)點(diǎn)添加到最后 % y(:,c)=y(:,1);%xy = [x;y]; end

2.樣條曲線差值（不過(guò)我沒(méi)有進(jìn)行插值）

% 樣條曲線差值 t=1:c; ts = 1:0.1:c; xys = spline(t,xy,ts); xs = xys(1,:); ys = xys(2,:); % 樣條差值效果 hold on temp=plot(x,y,'ro',xs,ys,'b.'); legend(temp,'原點(diǎn)','插值點(diǎn)');

?3.Snake主體部分（這一部分除參數(shù)外基本一樣）

******需要改動(dòng)的一個(gè)地方是：柵格圖像的行列號(hào)坐標(biāo)與實(shí)際坐標(biāo)是有差異的，差異在0.5個(gè)空間分辨率上，例如：行列坐標(biāo)計(jì)算為2的坐標(biāo)點(diǎn)事實(shí)上在柵格中是位于第二個(gè)像元和第三個(gè)像元之間的共享邊界上，因此不能直接使用計(jì)算得到的距離坐標(biāo)來(lái)代入到矩陣坐標(biāo)中進(jìn)行計(jì)算。

% ========================================================================= % Snakes算法實(shí)現(xiàn)部分 % ========================================================================= NIter =200; % 迭代次數(shù) alpha=0.2; beta=0.2; gamma = 1; kappa = 0.1; wl = 0.06; we=0.16; wt=0.78; [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:rowfor j= 1:colEterm(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% 外部力 Eext = Eimage + Econ Eext = wl*Eline + we*Eedge + wt*Eterm; % 計(jì)算梯度 [fx,fy]=gradient(Eext);xs=x';%數(shù)組序列與坐標(biāo)之間的差異 ys=y'; [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); figure(1) for i=1:NIter%限制輪廓點(diǎn)在圖像內(nèi)for p=1:mif(xs(p)<1)xs(p)=1;endif(ys(p)<1)ys(p)=1;endif(xs(p)>col)xs(p)=col;endif(ys(p)>row)ys(p)=row;endendssx = gamma*xs - kappa*interp2(fx,xs,ys);ssy = gamma*ys - kappa*interp2(fy,xs,ys);% 計(jì)算snake的新位置xs = Ainv * ssx;ys = Ainv * ssy; end

4.顯示圖像

imshow(Igs,[],'InitialMagnification','fit'); %'[]'是灰度顯示，"'InitialMagnification','fit'"是把圖像放大顯示以匹配對(duì)話框大小 hold on; plot(x',y','b-'); plot([xs;xs(1)], [ys; ys(1)], 'r-'); hold off;

imshow中的參數(shù)是為了放大圖像可移至matlab小技巧：imshow放大圖像顯示位置(以合適大小顯示)

Python

1.導(dǎo)入相關(guān)的包

import gdal import numpy as np import osr import ogr import os import csv from scipy import signal,ndimage,linalg import matplotlib.pyplot as plt

2.定義附屬函數(shù)?

cirshift（）? #用于做矩陣的上下循環(huán)移位操作，我沒(méi)有找見(jiàn)python中這個(gè)函數(shù)，但我認(rèn)為一定有，有知道的朋友麻煩告知一聲哦，謝謝

def circshift(A,K,m=1):B=np.array(A)x,y=A.shapefor i in range(0,x):for j in range(0,y):B[(i+K)%x,j]=A[i,j]return B

3.調(diào)用部分

Imagefilename="D:\\6.tif" Linefilename="D:\\6_former.shp" Igs,xy,xys=Snake_Python(Imagefilename,Linefilename,wl=0.06,we=0.14,wt=0.8)plt.imshow(Igs,cmap='gray') plt.plot(xy[0],xy[1]) plt.plot(xys[0],xys[1])#按列疊加在一起 plt.show()

4.Snake函數(shù)（主體）

使用python編寫(xiě)的Snake算法基本上是參照Matlab那個(gè)算法的流程寫(xiě)下來(lái)的，只在部分區(qū)域有差異。

def Snake_Python(Imagefilename,Linefilename,alpha=0.2,beta=0.2,gamma = 1,kappa = 0.1,wl = 0.06,we=0.16,wt=0.78,NIter =200):#打開(kāi)遙感影像dataset=gdal.Open(Imagefilename)im_width=dataset.RasterXSize #列數(shù)im_height=dataset.RasterYSize #行數(shù)im_bands=dataset.RasterCount #波段數(shù)im_geotrans=dataset.GetGeoTransform() #仿射矩陣im_proj=dataset.GetProjection() #地圖投影信息im_band=dataset.GetRasterBand(1)Igs=im_band.ReadAsArray(0,0,im_width,im_height)del dataset#關(guān)閉圖像進(jìn)程Igs=np.double(Igs)######################### 讀點(diǎn)數(shù)據(jù)(存在excel里處理過(guò))#with open("D:\\Point_6.csv") as pointprefile:# pointdata=csv.reader(pointprefile)# for pointrow in pointdata:# point.append([pointrow[3],pointrow[4]])########################## #打開(kāi)shapefile文件得到坐標(biāo)ds=ogr.Open(Linefilename,False) #打開(kāi)Shape文件（False - read only, True - read/write）layer=ds.GetLayer(0)#獲取圖層geomlist=[]#SF數(shù)據(jù)記錄– 幾何對(duì)象feature = layer.GetNextFeature() #獲得第一個(gè)SFwhile feature is not None:geom = feature.GetGeometryRef()#創(chuàng)建Geometry引用geomgeomlist += [geom.ExportToWkt()]feature = layer.GetNextFeature()ds.Destroy() #關(guān)閉shapefile進(jìn)程point=ogr.CreateGeometryFromWkt(geomlist[0]).GetPoints()#導(dǎo)出為string再生成Geometry，其實(shí)相當(dāng)于geom，可以使用geom在循環(huán)中使用GetPointspoint_change=np.array(point)x=point_change[:,0]y=point_change[:,1]x=(x-im_geotrans[0])/im_geotrans[1] -0.5 #投影坐標(biāo)轉(zhuǎn)換為行列號(hào)y=(y-im_geotrans[3])/im_geotrans[5] -0.5# =========================================================================# Snakes算法實(shí)現(xiàn)部分# =========================================================================#NIter =200; # 迭代次數(shù)#alpha=0.2; beta=0.2; gamma = 1; kappa = 0.1;#wl = 0.06; we=0.16; wt=0.78;row,col = Igs.shape#數(shù)組與矩陣關(guān)于方向的說(shuō)法正好相反，矩陣的x向是同行運(yùn)算，二維數(shù)組的x向是同列運(yùn)算（第一維度，axis=0）#圖像力 線函數(shù)Eline=Igs#圖像力 邊函數(shù)gy,gx=np.gradient(Igs)#############################################################################The gradient is computed using second order accurate central differences ##in the interior points and either first or second order accurate one-sides##(forward or backwards) differences at the boundaries. #############################################################################Eedge=-1*np.sqrt(gx*gx+gy*gy) #python中 *，dot是點(diǎn)乘；@，multiply是叉乘#圖像力-終點(diǎn)函數(shù)#卷積是為了求解偏導(dǎo)數(shù)，而離散點(diǎn)的偏導(dǎo)即差分求解m1 = [-1,1]m2 = [1,-2,1]m3 = [[1,-1],[-1,1]]#卷積 矩陣操作 axis=1 第二維度 行操作；axis=0第一緯度 列操作cx = ndimage.filters.convolve1d(Igs,m1,mode='constant',axis=1)cy = ndimage.filters.convolve1d(Igs,m1,mode='constant',axis=0)cxx = ndimage.filters.convolve1d(Igs,m2,mode='constant',axis=1)cyy = ndimage.filters.convolve1d(Igs,m2,mode='constant',axis=0)cxy = ndimage.filters.convolve(Igs,m3,mode='constant')Eterm=np.zeros([row,col])#初始化數(shù)組大小for i in range(0,row):for j in range(0,col):Eterm[i,j] =(cyy[i,j]*cx[i,j]**2 + cxx[i,j]*cy[i,j]**2 -2*cxy[i,j]*cx[i,j]*cy[i,j])/((1+cx[i,j]**2+cy[i,j]**2)**1.5)#計(jì)算外部力Eext = wl*Eline + we*Eedge + wt*Eterm;fy, fx = np.gradient(Eext);#計(jì)算五對(duì)角狀矩陣xs = x.copy()ys = y.copy()m= xs.shape[0]mm,nn = fx.shapeb=np.zeros([5,])b[0]=betab[1]=-(alpha + 4*beta)b[2]=(2*alpha + 6*beta) #b(i) 表示v(i)系數(shù)，從(i-2)到(i+2)b[3]=b[1];b[4]=b[0];A = b[0]*circshift(np.eye(m),2);A = A + b[1]*circshift(np.eye(m),1);A = A + b[2]*circshift(np.eye(m),0);A = A + b[3]*circshift(np.eye(m),-1);A = A + b[4]*circshift(np.eye(m),-2);##計(jì)算矩陣A+γI逆LU=linalg.lu(A+gama*np.eye(m),permute_l=True)L=LU[0]U=LU[1]Ainv=linalg.inv(U)@linalg.inv(L)for i in range(0,NIter):#限制邊界條件，約束點(diǎn)一直在圖像里xs[xs<0] = 0ys[ys<0] = 0xs[xs>(col-1)] = col - 1ys[ys>(row-1)] = row - 1#迭代計(jì)算xt=(A+γI)inv@(x(t-1)-fx(x(t-1),y(t-1)));yt同理ssx = gama*xs - kappa*ndimage.map_coordinates(fx,[ys,xs]);#fx中插值計(jì)算出[xs,ys]的點(diǎn)ssy = gama*ys - kappa*ndimage.map_coordinates(fy,[ys,xs]);xs = Ainv@ssx;ys = Ainv@ssy;return Igs,[x,y],[np.concatenate((xs,[xs[0]])),np.concatenate((ys,[ys[0]]))]

總結(jié)

此算法核心思想是一樣的，通過(guò)不同方法實(shí)現(xiàn)，在這個(gè)過(guò)程中也意識(shí)到了很多narray和matlab矩陣的區(qū)別，以及一些函數(shù)的區(qū)別，比如ndimage.map_coordinates和interp2的插值區(qū)別，當(dāng)然numpy中也有Matrix矩陣，在這個(gè)地方?jīng)]有用是因?yàn)橛X(jué)得numpy功能還是更強(qiáng)大一些哈哈哈，雖然可能這個(gè)強(qiáng)大并沒(méi)有用到，初學(xué)python和matlab，在這個(gè)改寫(xiě)的過(guò)程中學(xué)到了很多，另外，后來(lái)我發(fā)現(xiàn)了一個(gè)Python的skimage和Matlab的Active Contour Segmentation可以直接調(diào)用，但還沒(méi)有進(jìn)行深入研究，待研究后更新

總結(jié)

以上是生活随笔為你收集整理的Snake算法与遥感影像应用，python matlab对比的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

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