snake是一種主動(dòng)輪廓模型，笨妞對(duì)主動(dòng)輪廓模型的理解：你先給它一個(gè)初始輪廓，模型以初始輪廓為基準(zhǔn)逐步迭代，來(lái)改進(jìn)圖像的輪廓，使其更加精確。主動(dòng)輪廓模型目前用到了2種：CV和snake。前者沒(méi)有看算法內(nèi)部的原理。而snake，以最原始的論文《Snakes: Active Contour Models》為出發(fā)點(diǎn)。

1. snake原理

snake在逐步迭代優(yōu)化過(guò)程的目標(biāo)是能量函數(shù)最小化，這個(gè)能量函數(shù)指的是輪廓能量和圖像能量的總和（為什么要最小化這個(gè)能量總和，還不太清楚，論文也沒(méi)有具體說(shuō)）。snake的目標(biāo)不像sobel、canny等找到整張圖的輪廓。它只搜索你給出的初始輪廓附近，達(dá)到輪廓更精確的目標(biāo)，至少原版的snake只能達(dá)到局部?jī)?yōu)化的目標(biāo)。

能量函數(shù)：

其中指當(dāng)前輪廓本身的能量，稱(chēng)為內(nèi)部能量，而指圖像上輪廓對(duì)應(yīng)點(diǎn)的能量，稱(chēng)為外部能量，應(yīng)該是方差相關(guān)的項(xiàng)。

而內(nèi)部能量由兩部分構(gòu)成：一階導(dǎo)數(shù)的模（稱(chēng)為彈性能量）和二階導(dǎo)數(shù)的模（彎曲能量）

為什么是這樣的呢？據(jù)說(shuō)是因?yàn)榍€曲率的關(guān)系，閉合的輪廓曲線中，凸曲線按照法向量的方向，具有向內(nèi)的作用力；凹曲線法向量向外，具有向外的力。而曲率計(jì)算就是跟一階導(dǎo)數(shù)、二階導(dǎo)數(shù)相關(guān)的。很復(fù)雜，不甚理解。

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

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

線性能量：

邊緣能量：

終端（角點(diǎn)）能量：

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

對(duì)于優(yōu)化，優(yōu)化的目標(biāo)是總能量函數(shù)局部極小，通過(guò)能量函數(shù)極小或者迭代次數(shù)來(lái)控制迭代的終止。極小化能量函數(shù)通過(guò)歐拉方程計(jì)算解，作者在附錄中用了數(shù)值方法進(jìn)行推到，將歐拉方程推到為：

其中

引入外部能量：

再轉(zhuǎn)化為每一步迭代演進(jìn)過(guò)程：

A+ rI為五對(duì)角條帶矩陣。

2. GVF snake

關(guān)于圖像能量中l(wèi)ine、edge、termatation計(jì)算其實(shí)都挺復(fù)雜的。反倒是計(jì)算梯度向量場(chǎng)簡(jiǎn)單一些。將圖像能量由線、邊緣、角點(diǎn)的能量替換為梯度向量場(chǎng)，就是GVF snake。

3. 程序

對(duì)于snake，skimage里面active_contour函數(shù)就是典型的snake算法，借用里面的實(shí)現(xiàn)程序：

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

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

另外一個(gè)簡(jiǎn)化的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這個(gè)gvf snake基本的骨架是對(duì)的，但是太簡(jiǎn)單了，效果確實(shí)不好。

總結(jié)

以上是生活随笔為你收集整理的snake算法总结的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

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