【主动轮廓模型（一）】《Snakes: Active Contour Models》算法原理与OpenCV实现

傷城~ 2022-12-11 11:28 605阅读 0赞

### 文章目录 ###

*  1 概述
 *  2 算法原理
 *   *  2.1 内部能量  E i n t E\_\{int\} Eint
     *  2.2 图像能量  E i m a g e E\_\{image\} Eimage
     *   *  （1）线函数 E l i n e E\_\{line\} Eline
         *  （2）边函数 E e d g e E\_\{edge\} Eedge
         *  （3）末端函数 E t e r m E\_\{term\} Eterm
     *  2.3 外部能量  E c o n E\_\{con\} Econ
 *  3 模型求解
 *  4 算法实现（OpenCV3）

--------------------

# 1 概述 #

主动轮廓模型（也称Active Contour Model、Snake）是Kass等人在1988年提出的，该算法将图像分割问题转换为求解能量泛函最小值的问题。主要思路是通过构造能量泛函，经过算法迭代，轮廓曲线由初始位置逐渐向使能量函数最小（或局部极小）的图像边缘逼近，最终分割出目标。

# 2 算法原理 #

首先需要人为地在图像上给出初始轮廓曲线，确切的说是一组用于控制曲线形状的控制点： v ( s ) = \[ x ( s ) , y ( s ) \] s ∈ \[ 0 , 1 \] v(s)=\[x(s),y(s)\] s\\in\[0,1\] v(s)=\[x(s),y(s)\]s∈\[0,1\]，这些点收尾相连构成一个封闭的轮廓线。其中 x ( s ) x(s) x(s)和 y ( s ) y(s) y(s)分别表示每个控制点在图像中的坐标位置， s s s是以傅立叶变换形式描述边界的自变量，也可以理解为弧长。则Snake曲线的能量函数表示为：

E s n a k e ∗ = ∫ 0 1 E s n a k e ( v ( s ) ) d s = ∫ 0 1 E i n t ( v ( s ) ) + E i m a g e ( v ( s ) ) + E c o n ( v ( s ) ) d s = ∫ 0 1 E i n t ( v ( s ) ) + E e x t ( v ( s ) ) d s (1) \\begin\{aligned\} E\_\{snake\}^\* &= \\int\_0^1 E\_\{snake\}(v(s))ds \\\\ &= \\int\_0^1 E\_\{int\}(v(s))+E\_\{image\}(v(s))+E\_\{con\}(v(s))ds \\\\ &= \\int\_0^1 E\_\{int\}(v(s))+E\_\{ext\}(v(s))ds \\end\{aligned\} \\tag\{1\} Esnake∗=∫01Esnake(v(s))ds=∫01Eint(v(s))\+Eimage(v(s))\+Econ(v(s))ds=∫01Eint(v(s))\+Eext(v(s))ds(1)

其中， E i n t E\_\{int\} Eint为内部能量， E i m a g e E\_\{image\} Eimage为图像能量， E c o n E\_\{con\} Econ为外部约束能量。其中图像能量和外部约束能量统称为外部能量，即： E e x t = E i m a g e + E c o n E\_\{ext\}=E\_\{image\}+E\_\{con\} Eext=Eimage\+Econ。

## 2.1 内部能量  E i n t E\_\{int\} Eint ##

内部能量 E i n t E\_\{int\} Eint由保证曲线连续性的一阶导和保证曲线平滑的二阶导组成，表示为：

E i n t = 1 2 ( α ( s ) ∣ v s ( s ) ∣ 2 + β ( s ) ∣ v s s ( s ) ∣ 2 ) (2) E\_\{int\}=\\frac\{1\}\{2\} (\\alpha(s)|v\_s(s)|^2+\\beta(s)|v\_\{ss\}(s)|^2) \\tag\{2\} Eint=21(α(s)∣vs(s)∣2\+β(s)∣vss(s)∣2)(2)

通过调整权值 α ( s ) \\alpha(s) α(s)和 β ( s ) \\beta(s) β(s)可以控制曲线的形状。例如将 β ( s ) \\beta(s) β(s)置为0可以让曲线最终出现拐角，即曲线二阶不连续。

## 2.2 图像能量  E i m a g e E\_\{image\} Eimage ##

在低层次的计算机视觉中，需要能够将轮廓吸引到特定图像特征的能量函数。原始的Snake轮廓模型中提出了3种不同的能量函数，分别将snake轮廓吸引到线、边和末端。完整的图像能量  E i m a g e E\_\{image\} Eimage可以表示为这3个能量函数的权值组合，通过调整这3个权值，可以形成不同的轮廓形状。

E i m a g e = ω l i n e E l i n e + ω e d g e E e d g e + ω t e r m E t e r m (3) E\_\{image\}=\\omega\_\{line\}E\_\{line\}+\\omega\_\{edge\}E\_\{edge\}+\\omega\_\{term\}E\_\{term\} \\tag\{3\} Eimage=ωlineEline\+ωedgeEedge\+ωtermEterm(3)

### （1）线函数 E l i n e E\_\{line\} Eline ###

最简单直接且有用的图像能量函数是图像本身，即图像本身的灰度。令：

E l i n e = I ( x , y ) (4) E\_\{line\}=I(x,y) \\tag\{4\} Eline=I(x,y)(4)

控制 ω l i n e \\omega\_\{line\} ωline的正负号可以控制轮廓被吸引到较暗的线或是较亮的线，也就是使轮廓试图靠近轮廓的最暗或最亮处。

然而，如果snake轮廓的一部分到达了一个低能量的图像特征位置，这个样条项将推动snake轮廓临近的部分朝着这个特征可能的延续方向移动，这会使其在一个最优的局部最小位置引入一个较大的能量。一种解决方案是允许snake轮廓和模糊能量函数平衡，然后慢慢降低模糊程度。

Marr和Hildreth在他们的论文里证明了“*灰度的突然变化会在一阶导数中引起波峰或波谷，或在二阶导数中等效地引起零交点*”。为了显示图像尺度空间连续性和Marr-Hildreth边缘检测理论的关系，snake模型中采用了模糊的边能量函数：

E l i n e = − ( G σ ∗ ∇ 2 I ) 2 (6) E\_\{line\}=-(G\_\{\\sigma\}\*\\nabla^2I)^2 \\tag\{6\} Eline=−(Gσ∗∇2I)2(6)

其中 G σ G\_\{\\sigma\} Gσ是标准差为 σ \\sigma σ的高斯函数，该函数的最小值位于在Marr-Hildreth理论中被定义的 G σ ∗ ∇ 2 I G\_\{\\sigma\}\*\\nabla^2I Gσ∗∇2I的零交点处。在能量函数中加入这一项意味着snake轮廓在被吸引到零交点的同时，仍然受到它自己的平滑限制。

### （2）边函数 E e d g e E\_\{edge\} Eedge ###

在图像上找边缘可以通过梯度来实现，令：

E e d g e = − ∣ ∇ I ( x , y ) ∣ 2 (5) E\_\{edge\}=-|\\nabla I(x,y)|^2 \\tag\{5\} Eedge=−∣∇I(x,y)∣2(5)

在某个点上，梯度越大上式的能量越小，则snake轮廓将被吸引到梯度较大的区域。

### （3）末端函数 E t e r m E\_\{term\} Eterm ###

为了找到轮廓的终止位置，作者将平滑过图像中等高线的曲率加入到能量函数中。令 C ( x , y ) = G σ ( x , y ) ∗ I ( x , y ) C(x,y)=G\_\{\\sigma\}(x,y)\*I(x,y) C(x,y)=Gσ(x,y)∗I(x,y)是模糊后的图像， θ = t a n − 1 ( C y C x ) \\theta=tan^\{-1\}(\\frac\{C\_y\}\{C\_x\}) θ=tan−1(CxCy)是梯度角， n = ( c o s θ , s i n θ ) n=(cos\\theta,sin\\theta) n=(cosθ,sinθ)、 n ⊥ = ( − s i n θ , c o s θ ) n\_\{\\perp\}=(-sin\\theta,cos\\theta) n⊥=(−sinθ,cosθ)是沿着和垂直梯度方向的单位向量。则 C ( x , y ) C(x,y) C(x,y)中等高线的曲率可以表示为：

E t e r m = ∂ θ ∂ n ⊥ = ∂ C 2 / ∂ n ⊥ 2 ∂ C / ∂ n = C y y C x 2 − 2 C x y C x C y + C x x C y 2 ( C x 2 + C y 2 ) 3 / 2 (6) \\begin\{aligned\} E\_\{term\}&=\\frac\{\\partial \\theta\}\{\\partial n\_\{\\perp\}\} \\\\ &=\\frac\{\\partial C^2 / \\partial n\_\{\\perp\}^2\}\{\\partial C / \\partial n\} \\\\ &=\\frac\{C\_\{yy\}C\_x^2-2C\_\{xy\}C\_xC\_y+C\_\{xx\}C\_y^2\}\{(C\_x^2+C\_y^2)^\{3/2\}\} \\end\{aligned\} \\tag\{6\} Eterm=∂n⊥∂θ=∂C/∂n∂C2/∂n⊥2=(Cx2\+Cy2)3/2CyyCx2−2CxyCxCy\+CxxCy2(6)

通过组合 E e d g e E\_\{edge\} Eedge和 E t e r m E\_\{term\} Eterm，我们能够创建一个被边和末端吸引的snake曲线。

## 2.3 外部能量  E c o n E\_\{con\} Econ ##

外部能量 E c o n E\_\{con\} Econ来自于外部的约束力。在论文原文中，作者给了一个外部约束的例子， 包括外部的固定点、连接两条snake曲线的锚点或鼠标拖动的点。例如，为了在 x 1 x\_1 x1和 x 2 x\_2 x2点之间创建一个连接的弹性力，就可以将 − k ( x 1 − x 2 ) 2 -k(x\_1-x\_2)^2 −k(x1−x2)2添加到外部能量 E c o n E\_\{con\} Econ中。

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2NmYW45Mjc_size_16_color_FFFFFF_t_70_pic_center_pic_center]

# 3 模型求解 #

由欧拉方程，求解能量 E s n a k e ∗ E\_\{snake\}^\* Esnake∗的最小值（局部极小值），公式（1）的导数必须满足：

α v s s ( s ) − β v s s s s ( s ) − ∇ E i m a g e ( v ( s ) ) − ∇ E c o n ( v ( s ) ) = 0 (7) \\alpha v\_\{ss\}(s) - \\beta v\_\{ssss\}(s) - \\nabla E\_\{image\}(v(s)) - \\nabla E\_\{con\}(v(s)) = 0 \\tag\{7\} αvss(s)−βvssss(s)−∇Eimage(v(s))−∇Econ(v(s))=0(7)

由 v ( s ) = \[ x ( s ) , y ( s ) \] v(s)=\[x(s),y(s)\] v(s)=\[x(s),y(s)\]，将上式改写成x和y两个方向有：

\{ α x s s ( s ) − β x s s s s ( s ) − ∂ E e x t ∂ x = 0 α y s s ( s ) − β y s s s s ( s ) − ∂ E e x t ∂ y = 0 (8) \\begin\{cases\} \\alpha x\_\{ss\}(s) - \\beta x\_\{ssss\}(s) - \\frac\{\\partial\{E\}\_\{ext\}\}\{\\partial x\} = 0 \\\\ \\alpha y\_\{ss\}(s) - \\beta y\_\{ssss\}(s) - \\frac\{\\partial\{E\}\_\{ext\}\}\{\\partial y\} = 0 \\end\{cases\} \\tag\{8\} \{ αxss(s)−βxssss(s)−∂x∂Eext=0αyss(s)−βyssss(s)−∂y∂Eext=0(8)

利用顺序连接的锚点的坐标差分来近似导数：

\{ x s s ( s ) = x ( s + 1 ) + x ( s − 1 ) − 2 x ( s ) x s s s s ( s ) = ( x ( s + 2 ) + x ( s ) − 2 x ( s + 1 ) ) + ( x ( s ) + x ( s − 2 ) − 2 x ( s − 1 ) ) − 2 ( x ( s + 1 ) + x ( s − 1 ) − 2 x ( s ) ) (9) \\begin\{cases\} x\_\{ss\}(s)=x(s+1)+x(s-1)-2x(s) \\\\ \\begin\{aligned\} x\_\{ssss\}(s)=&(x(s+2)+x(s)-2x(s+1)) \\\\ &+(x(s)+x(s-2)-2x(s-1)) \\\\ &-2(x(s+1)+x(s-1)-2x(s)) \\end\{aligned\} \\end\{cases\} \\tag\{9\} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧xss(s)=x(s\+1)\+x(s−1)−2x(s)xssss(s)=(x(s\+2)\+x(s)−2x(s\+1))\+(x(s)\+x(s−2)−2x(s−1))−2(x(s\+1)\+x(s−1)−2x(s))(9)

将（9）带入（8），同时令 f x = ∂ E e x t ∂ x , f y = ∂ E e x t ∂ y f\_x=\\frac\{\\partial E\_\{ext\}\}\{\\partial x\}, f\_y=\\frac\{\\partial E\_\{ext\}\}\{\\partial y\} fx=∂x∂Eext,fy=∂y∂Eext有：

\{ β x ( s − 2 ) − ( α + 4 β ) x ( s − 1 ) + ( 2 α + 6 β ) x ( s ) − ( α + 4 β ) x ( s + 1 ) + β x ( s + 2 ) + f x = 0 β y ( s − 2 ) − ( α + 4 β ) y ( s − 1 ) + ( 2 α + 6 β ) y ( s ) − ( α + 4 β ) y ( s + 1 ) + β y ( s + 2 ) + f y = 0 (10) \\begin\{cases\} \\beta x(s-2)-(\\alpha+4\\beta)x(s-1)+(2\\alpha+6\\beta)x(s)-(\\alpha+4\\beta)x(s+1)+\\beta x(s+2)+f\_x=0 \\\\ \\beta y(s-2)-(\\alpha+4\\beta)y(s-1)+(2\\alpha+6\\beta)y(s)-(\\alpha+4\\beta)y(s+1)+\\beta y(s+2)+f\_y=0 \\end\{cases\} \\tag\{10\} \{ βx(s−2)−(α\+4β)x(s−1)\+(2α\+6β)x(s)−(α\+4β)x(s\+1)\+βx(s\+2)\+fx=0βy(s−2)−(α\+4β)y(s−1)\+(2α\+6β)y(s)−(α\+4β)y(s\+1)\+βy(s\+2)\+fy=0(10)

令：

\{ a = 2 α + 6 β b = − ( α + 4 β ) c = β (11) \\begin\{cases\} a=2\\alpha + 6\\beta \\\\ b=-(\\alpha+4\\beta) \\\\ c=\\beta \\end\{cases\} \\tag\{11\} ⎩⎪⎨⎪⎧a=2α\+6βb=−(α\+4β)c=β(11)

则将（10）写成矩阵形式为：

\{ A x + f x ( x , y ) = 0 A y + f y ( x , y ) = 0 (12) \\begin\{cases\} Ax+f\_x(x,y)=0 \\\\ Ay+f\_y(x,y)=0 \\end\{cases\} \\tag\{12\} \{ Ax\+fx(x,y)=0Ay\+fy(x,y)=0(12)

其中A为五对角带状矩阵：

A = \[ a b c ⋯ c b b a b c ⋯ c c b a b c ⋯ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋯ c b a b c c ⋯ c b a b b c ⋯ c b a \] , x = \[ x 1 x 2 x 3 ⋮ x n − 1 x n x 1 \] , f x = \[ f x 1 f x 2 f x 3 ⋮ f x n − 1 f x n f x 1 \] (13) A= \\begin\{bmatrix\} a & b & c & \\cdots & c & b \\\\ b & a & b & c & \\cdots & c \\\\ c & b & a & b & c & \\cdots \\\\ \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots \\\\ \\cdots & c & b & a & b & c \\\\ c & \\cdots & c & b & a & b \\\\ b & c & \\cdots & c & b & a \\end\{bmatrix\}, x= \\begin\{bmatrix\} x\_1 \\\\ x\_2 \\\\ x\_3 \\\\ \\vdots \\\\ x\_\{n-1\} \\\\ x\_n \\\\ x\_1 \\end\{bmatrix\}, f\_x= \\begin\{bmatrix\} f\_\{x\_1\} \\\\ f\_\{x\_2\} \\\\ f\_\{x\_3\} \\\\ \\vdots \\\\ f\_\{x\_\{n-1\}\} \\\\ f\_\{x\_n\} \\\\ f\_\{x\_1\} \\end\{bmatrix\} \\tag\{13\} A=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡abc⋮⋯cbbab⋮c⋯ccba⋮bc⋯⋯cb⋮abcc⋯c⋮babbc⋯⋮cba⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤,x=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡x1x2x3⋮xn−1xnx1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤,fx=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡fx1fx2fx3⋮fxn−1fxnfx1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤(13)

利用梯度下降法，在第t次迭代有：

\{ A x t + f x ( x t − 1 , y t − 1 ) = − γ ( x t − x t − 1 ) A y t + f y ( y t − 1 , y t − 1 ) = − γ ( y t − y t − 1 ) (14) \\begin\{cases\} Ax\_t+f\_x(x\_\{t-1\},y\_\{t-1\})=-\\gamma(x\_t-x\_\{t-1\}) \\\\ Ay\_t+f\_y(y\_\{t-1\},y\_\{t-1\})=-\\gamma(y\_t-y\_\{t-1\}) \\end\{cases\} \\tag\{14\} \{ Axt\+fx(xt−1,yt−1)=−γ(xt−xt−1)Ayt\+fy(yt−1,yt−1)=−γ(yt−yt−1)(14)

其中 γ \\gamma γ是迭代步长，公式（14）的解为：

\{ x t = ( A + γ I ) − 1 ( x t − 1 − f x ( x t − 1 , y t − 1 ) ) y t = ( A + γ I ) − 1 ( y t − 1 − f y ( x t − 1 , y t − 1 ) ) (15) \\begin\{cases\} x\_t=(A+\\gamma I)^\{-1\}(x\_\{t-1\}-f\_x(x\_\{t-1\},y\_\{t-1\})) \\\\ y\_t=(A+\\gamma I)^\{-1\}(y\_\{t-1\}-f\_y(x\_\{t-1\},y\_\{t-1\})) \\end\{cases\} \\tag\{15\} \{ xt=(A\+γI)−1(xt−1−fx(xt−1,yt−1))yt=(A\+γI)−1(yt−1−fy(xt−1,yt−1))(15)

由于 A A A为五对角带状矩阵，因此 A + γ I A+\\gamma I A\+γI也是一个五对角带状矩阵，原文中作者用LU分解来求其逆矩阵。

# 4 算法实现（OpenCV3） #

网上关于snake算法的实现多是matlab代码。OpenCV2中可用cvSnakeImage来实现，但这个函数在OpenCV3中已经被删除。本文将在OpenCV3.14中实现snake算法代码。

cv::Mat Interate(
    	cv::Mat image,
    	cv::Mat xs,
    	cv::Mat ys,
    	double alpha,
    	double beta,
    	double gamma,
    	double kappa,
    	double wl,
    	double we,
    	double wt,
    	int iterations
    )
    { 
    	// 相关参数
    	int N = iterations;
    	cv::Mat smth = image.clone();
    
    	// 图像大小
    	qDebug() << "Calculating size of image";
    	cv::Size size = image.size();
    	int row = size.height;
    	int col = size.width;
    
    	// 计算外部力（图像力）
    	qDebug() << "Computing external forces";
    	cv::Mat E_line = smth.clone(); // E_line is simply the image intensities
    
    	cv::Mat gradx, grady;
    	cv::Sobel(smth, gradx, smth.depth(), 1, 0, 1, 1, 0, cv::BORDER_CONSTANT);
    	cv::Sobel(smth, grady, smth.depth(), 0, 1, 1, 1, 0, cv::BORDER_CONSTANT);
    
    	qDebug() << "Computing gradx and grady";
    	cv::Mat E_edge(row, col, CV_32FC1);
    	for (int i = 0; i < gradx.rows; i++)
    	{ 
    		for (int j = 0; j < gradx.cols; j++)
    		{ 
    			float v_gradx = gradx.at<float>(i, j);
    			float v_grady = grady.at<float>(i, j);
    			
    			E_edge.at<float>(i, j) = -1 * std::sqrt(v_gradx * v_gradx + v_grady * v_grady); // E_edge is measured by gradient in the image
    		}
    	}
    
    	// 导数mask
    	qDebug() << "masks for taking various derivatives";
    	cv::Mat m1 = (cv::Mat_<float>(1, 2) << -1, 1);
    	cv::Mat m2 = (cv::Mat_<float>(2, 1) << -1, 1);
    	cv::Mat m3 = (cv::Mat_<float>(1, 3) << 1, -2, 1);
    	cv::Mat m4 = (cv::Mat_<float>(3, 1) << 1, -2, 1);
    	cv::Mat m5 = (cv::Mat_<float>(2, 2) << 1, -1, -1, 1);
    
    	cv::Mat cx, cy, cxx, cyy, cxy;
    	filter2D(smth, cx, -1, m1);
    	filter2D(smth, cy, -1, m2);
    	filter2D(smth, cxx, -1, m3);
    	filter2D(smth, cyy, -1, m4);
    	filter2D(smth, cxy, -1, m5);
    
    	// 计算 E_term
    	cv::Mat E_term(row, col, CV_32FC1);
    	for (int i = 0; i < row; i++)
    	{ 
    		for (int j = 0; j < col; j++)
    		{ 
    			int v_cx = cx.at<float>(i, j);
    			int v_cy = cy.at<float>(i, j);
    			int v_cxx = cxx.at<float>(i, j);
    			int v_cyy = cyy.at<float>(i, j);
    			int v_cxy = cxy.at<float>(i, j);
    
    			E_term.at<float>(i, j) = (v_cyy*v_cx*v_cx - 2 * v_cxy*v_cx*v_cy + v_cxx * v_cy*v_cy) / (std::pow((1 + v_cx * v_cx + v_cy * v_cy), 1.5));
    		}
    	}
    
    	// 计算E_ext
    	cv::Mat E_ext = (wl*E_line + we * E_edge - wt * E_term);
    
    	// 计算梯度
    	cv::Mat fx, fy;
    	cv::Sobel(E_ext, fx, E_ext.depth(), 1, 0, 1, 0.5, 0, cv::BORDER_CONSTANT);
    	cv::Sobel(E_ext, fy, E_ext.depth(), 0, 1, 1, 0.5, 0, cv::BORDER_CONSTANT);
    
    	cv::transpose(xs, xs);
    	cv::transpose(ys, ys);
    
    	int m = xs.rows;
    	int n = 1;
    
    	int mm = fx.cols;
    	int nn = fx.rows;
    
    	// 计算五对角状矩阵，b(i)表示vi系数(i = i - 2 到 i + 2)
    	double b[5];
    	b[0] = beta;
    	b[1] = -(alpha + 4 * beta);
    	b[2] = 2 * alpha + 6 * beta;
    	b[3] = b[1];
    	b[4] = b[0];
    
    	cv::Mat A = cv::Mat::eye(m, m, CV_32FC1);
    	cv::Mat eyeMat0 = cv::Mat::eye(m, m, CV_32FC1);
    	circRowShift(eyeMat0, 2);
    	eyeMat0.convertTo(eyeMat0, CV_32FC1);
    	A = b[0] * eyeMat0;
    
    	cv::Mat eyeMat1 = cv::Mat::eye(m, m, CV_32FC1);
    	circRowShift(eyeMat1, 1);
    	eyeMat1.convertTo(eyeMat1, CV_32FC1);
    	A = A + b[1] * eyeMat1;
    
    	cv::Mat eyeMat2 = cv::Mat::eye(m, m, CV_32FC1);
    	circRowShift(eyeMat2, 0);
    	eyeMat2.convertTo(eyeMat2, CV_32FC1);
    	A = A + b[2] * eyeMat2;
    
    	cv::Mat eyeMat3 = cv::Mat::eye(m, m, CV_32FC1);
    	circRowShift(eyeMat3, -1);
    	eyeMat3.convertTo(eyeMat3, CV_32FC1);
    	A = A + b[3] * eyeMat3;
    
    	cv::Mat eyeMat4 = cv::Mat::eye(m, m, CV_32FC1);
    	circRowShift(eyeMat4, -2);
    	eyeMat4.convertTo(eyeMat4, CV_32FC1);
    	A = A + b[4] * eyeMat4;
    
    	// 计算矩阵的逆
    	cv::Mat Ainv(A.size(), CV_32FC1);
    	A = A + gamma * cv::Mat::eye(m, m, CV_32FC1);
    	cv::invert(A, Ainv); // Computing Ainv 
    
    	cv::Mat srcImg = cv::imread("D:/endo_image.jpg");
    	// 迭代更新曲线
    	for (int i = 0; i < N; i++)
    	{ 
    		cv::Mat intFx(fx.size(), CV_32FC1);
    		cv::Mat intFy(fy.size(), CV_32FC1);
    
    		cv::remap(fx, intFx, xs, ys, cv::INTER_LINEAR, cv::BORDER_CONSTANT);
    		cv::remap(fy, intFy, xs, ys, cv::INTER_LINEAR, cv::BORDER_CONSTANT);
    
    		cv::Mat ssx(xs.size(), CV_32FC1);
    		cv::Mat ssy(ys.size(), CV_32FC1);
    		for (int k = 0; k < xs.rows; k++)
    		{ 
    			for (int l = 0; l < xs.cols; l++)
    			{ 
    				ssx.at<float>(k, l) = gamma * xs.at<float>(k, l) - kappa * intFx.at<float>(k, l);
    				ssy.at<float>(k, l) = gamma * ys.at<float>(k, l) - kappa * intFy.at<float>(k, l);
    			}
    		}
    
    		// 更新曲线位置
    		xs = Ainv * ssx;
    		ys = Ainv * ssy;
    
    		cv::Mat resultImg = srcImg.clone();
    		for (int j = 0; j < xs.rows; j++)
    		{ 
    			cv::Point center = cv::Point(xs.at<float>(j, 0), ys.at<float>(j, 0));
    			cv::circle(resultImg, center, 4, cv::Scalar(0, 255, 0), 2);
    		}
    
    		// 显示
    		cv::imshow("result", resultImg);
    		cv::waitKey(30);
    	}
    
    	return image;
    }

运行结果：

![在这里插入图片描述][20201003211832516.gif_pic_center]

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2NmYW45Mjc_size_16_color_FFFFFF_t_70_pic_center_pic_center]: /images/20221123/0036379a0cf3434fa97d1ec4f5e87359.png
[20201003211832516.gif_pic_center]: /images/20221123/58c7649824d9449486aec6c38d5ff44c.png