snake是一種主動輪廓模型，笨妞對主動輪廓模型的理解：你先給它一個初始輪廓，模型以初始輪廓為基準逐步迭代，來改進圖像的輪廓，使其更加精確。主動輪廓模型目前用到了2種：CV和snake。前者沒有看算法內部的原理。而snake，以最原始的論文《Snakes: Active Contour Models》為出發點。

1. snake原理

snake在逐步迭代優化過程的目標是能量函數最小化，這個能量函數指的是輪廓能量和圖像能量的總和（為什么要最小化這個能量總和，還不太清楚，論文也沒有具體說）。snake的目標不像sobel、canny等找到整張圖的輪廓。它只搜索你給出的初始輪廓附近，達到輪廓更精確的目標，至少原版的snake只能達到局部優化的目標。

能量函數：

其中指當前輪廓本身的能量，稱為內部能量，而指圖像上輪廓對應點的能量，稱為外部能量，應該是方差相關的項。

而內部能量由兩部分構成：一階導數的模（稱為彈性能量）和二階導數的模（彎曲能量）

為什么是這樣的呢？據說是因為曲線曲率的關系，閉合的輪廓曲線中，凸曲線按照法向量的方向，具有向內的作用力；凹曲線法向量向外，具有向外的力。而曲率計算就是跟一階導數、二階導數相關的。很復雜，不甚理解。

在迭代過程中，彈性能量能快速的把輪廓壓縮成光滑的圓；彎曲能量將輪廓拉成光滑的曲線或直線，他們的作用是保持輪廓的光滑和連續性。通常alpha越大，輪廓收斂越快；beta越大，輪廓越光滑。

外部圖像能量作者分了三種：線性能量，通常更亮度相關；邊緣能量，由圖像的邊緣組成，而邊緣可以通過sobel算子計算；終端能量。

線性能量：

邊緣能量：

終端（角點）能量：

通常可以根據更期望輪廓趨向于哪方面來選擇以上三種能量。在迭代優化過程中，外部能量會使輪廓朝（灰度）高梯度位置靠近。而通常梯度高的位置都是圖像中前景與背景的界限或者物體與物體之間、物體內部不同部分的界限，適合用于分割。

對于優化，優化的目標是總能量函數局部極小，通過能量函數極小或者迭代次數來控制迭代的終止。極小化能量函數通過歐拉方程計算解，作者在附錄中用了數值方法進行推到，將歐拉方程推到為：

其中

引入外部能量：

再轉化為每一步迭代演進過程：

A+ rI為五對角條帶矩陣。

2. GVF snake

關于圖像能量中line、edge、termatation計算其實都挺復雜的。反倒是計算梯度向量場簡單一些。將圖像能量由線、邊緣、角點的能量替換為梯度向量場，就是GVF snake。

3. 程序

對于snake，skimage里面active_contour函數就是典型的snake算法，借用里面的實現程序：

import numpy as np import scipy.linalg from scipy.interpolate import RectBivariateSpline from skimage.util import img_as_float from skimage.filters import sobeldef active_contour(image, snake, alpha=0.01, beta=0.1,w_line=0, w_edge=1, gamma=0.01,bc='periodic', max_px_move=1.0,max_iterations=2500, convergence=0.1):"""Active contour model.Active contours by fitting snakes to features of images. Supports singleand multichannel 2D images. Snakes can be periodic (for segmentation) orhave fixed and/or free ends.The output snake has the same length as the input boundary.As the number of points is constant, make sure that the initial snakehas enough points to capture the details of the final contour.Parameters----------image : (N, M) or (N, M, 3) ndarrayInput image.snake : (N, 2) ndarrayInitialisation coordinates of snake. For periodic snakes, it shouldnot include duplicate endpoints.alpha : float, optionalSnake length shape parameter. Higher values makes snake contractfaster.beta : float, optionalSnake smoothness shape parameter. Higher values makes snake smoother.w_line : float, optionalControls attraction to brightness. Use negative values to attract todark regions.w_edge : float, optionalControls attraction to edges. Use negative values to repel snake fromedges.gamma : float, optionalExplicit time stepping parameter.bc : {'periodic', 'free', 'fixed'}, optionalBoundary conditions for worm. 'periodic' attaches the two ends of thesnake, 'fixed' holds the end-points in place, and'free' allows freemovement of the ends. 'fixed' and 'free' can be combined by parsing'fixed-free', 'free-fixed'. Parsing 'fixed-fixed' or 'free-free'yields same behaviour as 'fixed' and 'free', respectively.max_px_move : float, optionalMaximum pixel distance to move per iteration.max_iterations : int, optionalMaximum iterations to optimize snake shape.convergence: float, optionalConvergence criteria.Returns-------snake : (N, 2) ndarrayOptimised snake, same shape as input parameter.References----------.. [1] Kass, M.; Witkin, A.; Terzopoulos, D. "Snakes: Active contourmodels". International Journal of Computer Vision 1 (4): 321(1988).Examples-------->>> from skimage.draw import circle_perimeter>>> from skimage.filters import gaussianCreate and smooth image:>>> img = np.zeros((100, 100))>>> rr, cc = circle_perimeter(35, 45, 25)>>> img[rr, cc] = 1>>> img = gaussian(img, 2)Initiliaze spline:>>> s = np.linspace(0, 2*np.pi,100)>>> init = 50*np.array([np.cos(s), np.sin(s)]).T+50Fit spline to image:>>> snake = active_contour(img, init, w_edge=0, w_line=1) #doctest: +SKIP>>> dist = np.sqrt((45-snake[:, 0])**2 +(35-snake[:, 1])**2) #doctest: +SKIP>>> int(np.mean(dist)) #doctest: +SKIP25"""max_iterations = int(max_iterations)if max_iterations <= 0:raise ValueError("max_iterations should be >0.")convergence_order = 10valid_bcs = ['periodic', 'free', 'fixed', 'free-fixed','fixed-free', 'fixed-fixed', 'free-free']if bc not in valid_bcs:raise ValueError("Invalid boundary condition.\n" +"Should be one of: "+", ".join(valid_bcs)+'.')img = img_as_float(image) height = img.shape[0]width = img.shape[1]RGB = img.ndim == 3# Find edges using sobel:if w_edge != 0:if RGB:edge = [sobel(img[:, :, 0]), sobel(img[:, :, 1]),sobel(img[:, :, 2])]else:edge = [sobel(img)]for i in range(3 if RGB else 1):edge[i][0, :] = edge[i][1, :]edge[i][-1, :] = edge[i][-2, :]edge[i][:, 0] = edge[i][:, 1]edge[i][:, -1] = edge[i][:, -2]else:edge = [0]# Superimpose intensity and edge images:if RGB:img = w_line*np.sum(img, axis=2) \+ w_edge*sum(edge)else:img = w_line*img + w_edge*edge[0]# Interpolate for smoothness:intp = RectBivariateSpline(np.arange(img.shape[1]),np.arange(img.shape[0]),img.T, kx=2, ky=2, s=0)x, y = snake[:, 0].astype(np.float), snake[:, 1].astype(np.float)xsave = np.empty((convergence_order, len(x)))ysave = np.empty((convergence_order, len(x)))# Build snake shape matrix for Euler equationn = len(x)a = np.roll(np.eye(n), -1, axis=0) + \np.roll(np.eye(n), -1, axis=1) - \2*np.eye(n) # second order derivative, central differenceb = np.roll(np.eye(n), -2, axis=0) + \np.roll(np.eye(n), -2, axis=1) - \4*np.roll(np.eye(n), -1, axis=0) - \4*np.roll(np.eye(n), -1, axis=1) + \6*np.eye(n) # fourth order derivative, central differenceA = -alpha*a + beta*b# Impose boundary conditions different from periodic:sfixed = Falseif bc.startswith('fixed'):A[0, :] = 0A[1, :] = 0A[1, :3] = [1, -2, 1]sfixed = Trueefixed = Falseif bc.endswith('fixed'):A[-1, :] = 0A[-2, :] = 0A[-2, -3:] = [1, -2, 1]efixed = Truesfree = Falseif bc.startswith('free'):A[0, :] = 0A[0, :3] = [1, -2, 1]A[1, :] = 0A[1, :4] = [-1, 3, -3, 1]sfree = Trueefree = Falseif bc.endswith('free'):A[-1, :] = 0A[-1, -3:] = [1, -2, 1]A[-2, :] = 0A[-2, -4:] = [-1, 3, -3, 1]efree = True# Only one inversion is needed for implicit spline energy minimization:inv = scipy.linalg.inv(A+gamma*np.eye(n))# Explicit time stepping for image energy minimization:for i in range(max_iterations):fx = intp(x, y, dx=1, grid=False)fy = intp(x, y, dy=1, grid=False)if sfixed:fx[0] = 0fy[0] = 0if efixed:fx[-1] = 0fy[-1] = 0if sfree:fx[0] *= 2fy[0] *= 2if efree:fx[-1] *= 2fy[-1] *= 2xn = np.dot(inv, gamma*x + fx)yn = np.dot(inv, gamma*y + fy)# Movements are capped to max_px_move per iteration:dx = max_px_move*np.tanh(xn-x)dy = max_px_move*np.tanh(yn-y)if sfixed:dx[0] = 0dy[0] = 0if efixed:dx[-1] = 0dy[-1] = 0x += dxy += dy ? ? ? ? x[x<0] = 0y[y<0] = 0x[x>(height-1)] = height - 1y[y>(width-1)] = width - 1# Convergence criteria needs to compare to a number of previous# configurations since oscillations can occur.j = i % (convergence_order+1)if j < convergence_order:xsave[j, :] = xysave[j, :] = yelse:dist = np.min(np.max(np.abs(xsave-x[None, :]) +np.abs(ysave-y[None, :]), 1))if dist < convergence:breakreturn np.array([x, y]).T

這個程序有個bug，沒有對新計算出的輪廓做約束，這樣的話如果初始snake比較靠近圖像的邊，那么輪廓就會溢出到圖像外。紅色部分是個人做的修改。

另外一個簡化的gvf snake程序

file_name = '24_82_11.bmp' #file_path = os.path.join('cell_split/23-2018-05-1708.57.22', file_name) img = skimage.color.rgb2gray(skimage.data.imread(file_name)) f = filters.gaussian_gradient_magnitude(img,2.0) fx = filters.gaussian_filter(f,2.0,(1,0)) fy = filters.gaussian_filter(f,2.0,(0,1)) u = fx.copy() v = fy.copy() mu = 0.1 for i in range(10000):m = (fx**2+fy**2)ut = mu*filters.laplace(u)-(u-fx)*mvt = mu*filters.laplace(v)-(v-fy)*mu += utv += vtt = np.linspace(0,2*np.pi,50)for t in range(5000):x2 = filters.gaussian_filter(path,1.0,(2,0),mode='wrap')x4 = filters.gaussian_filter(x2,1.0,(2,0),mode='wrap')dx = np.array([u[int(x),int(y)] for y,x in path])dy = np.array([v[int(x),int(y)] for y,x in path])delta = np.array([dx,dy]).Tpath += 0.001*x2-0.001*x4+1.0*delta#print(path.shape)path[:,0][path[:,0]>129] = 129path[:,1][path[:,1]>105] = 105這個gvf snake基本的骨架是對的，但是太簡單了，效果確實不好。

總結

以上是生活随笔為你收集整理的snake算法总结的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。