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

Snake算法原理

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

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

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

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

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

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

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

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

可以得到Eterminal

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

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

田問渠Carlnait 主動輪廓模型——Snake分割算法 matlab源碼

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

?

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

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

1.打開影像與分割線

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)狀點集合，因此不需要把第一個點添加到最后 % y(:,c)=y(:,1);%xy = [x;y]; end

2.樣條曲線差值（不過我沒有進行插值）

% 樣條曲線差值 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,'原點','插值點');

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

******需要改動的一個地方是：柵格圖像的行列號坐標與實際坐標是有差異的，差異在0.5個空間分辨率上，例如：行列坐標計算為2的坐標點事實上在柵格中是位于第二個像元和第三個像元之間的共享邊界上，因此不能直接使用計算得到的距離坐標來代入到矩陣坐標中進行計算。

% ========================================================================= % Snakes算法實現(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)); % 圖像力-終點函數(shù) % 卷積是為了求解偏導數(shù)，而離散點的偏導即差分求解 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; % 計算梯度 [fx,fy]=gradient(Eext);xs=x';%數(shù)組序列與坐標之間的差異 ys=y'; [m,n] = size(xs); [mm,nn] = size(fx); % 計算五對角狀矩陣 % 附錄: 公式（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);% 計算矩陣的逆 [L,U] = lu(A + gamma * eye(m)); Ainv = inv(U) * inv(L); figure(1) for i=1:NIter%限制輪廓點在圖像內(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);% 計算snake的新位置xs = Ainv * ssx;ys = Ainv * ssy; end

4.顯示圖像

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

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

Python

1.導入相關(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)移位操作，我沒有找見python中這個函數(shù)，但我認為一定有，有知道的朋友麻煩告知一聲哦，謝謝

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編寫的Snake算法基本上是參照Matlab那個算法的流程寫下來的，只在部分區(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):#打開遙感影像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)閉圖像進程Igs=np.double(Igs)######################### 讀點數(shù)據(jù)(存在excel里處理過)#with open("D:\\Point_6.csv") as pointprefile:# pointdata=csv.reader(pointprefile)# for pointrow in pointdata:# point.append([pointrow[3],pointrow[4]])########################## #打開shapefile文件得到坐標ds=ogr.Open(Linefilename,False) #打開Shape文件（False - read only, True - read/write）layer=ds.GetLayer(0)#獲取圖層geomlist=[]#SF數(shù)據(jù)記錄– 幾何對象feature = layer.GetNextFeature() #獲得第一個SFwhile feature is not None:geom = feature.GetGeometryRef()#創(chuàng)建Geometry引用geomgeomlist += [geom.ExportToWkt()]feature = layer.GetNextFeature()ds.Destroy() #關(guān)閉shapefile進程point=ogr.CreateGeometryFromWkt(geomlist[0]).GetPoints()#導出為string再生成Geometry，其實相當于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 #投影坐標轉(zhuǎn)換為行列號y=(y-im_geotrans[3])/im_geotrans[5] -0.5# =========================================================================# Snakes算法實現(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)于方向的說法正好相反，矩陣的x向是同行運算，二維數(shù)組的x向是同列運算（第一維度，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是點乘；@，multiply是叉乘#圖像力-終點函數(shù)#卷積是為了求解偏導數(shù)，而離散點的偏導即差分求解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)#計算外部力Eext = wl*Eline + we*Eedge + wt*Eterm;fy, fx = np.gradient(Eext);#計算五對角狀矩陣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);##計算矩陣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):#限制邊界條件，約束點一直在圖像里xs[xs<0] = 0ys[ys<0] = 0xs[xs>(col-1)] = col - 1ys[ys>(row-1)] = row - 1#迭代計算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中插值計算出[xs,ys]的點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é)

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

總結(jié)

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

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