【演算法+OpenCV】基於三次Bezier原理的曲線擬合演算法C++與OpenCV實現
阿新 • • 發佈：2019-02-04
近期，因為要實現經過多個控制點的曲線擬合，研究起了曲線擬合演算法。綜合搜尋到的資料，發現Bezier曲線擬合演算法是一種相對較容易實現、且擬合的效果較好的演算法。關於Bezier曲線原理，請參照（Bezier曲線原理），這裡就不再做具體介紹了，我們使用的是Besier三次曲線擬合原理。下面主要介紹演算法的實現過程。
如下圖中，P0、P1、P2、P3四個點，我們最終是想獲取過這四個點的封閉平滑曲線。
根據Bezier三次曲線擬合的原理，我們可以分別擬合P0P1、P1P2、P2P3、P3P0四段曲線，進而連線成一個封閉的曲線。但是，Bezier三次曲線擬合需要在兩點之間找到兩個控制點。每個點的控制點可以根據其前後相鄰的兩點獲得，具體實現如下：
void get_control_points(double x0, double y0, double x1, double y1, double x2, double y2, double& p1x, double& p1y, double& p2x, double& p2y, double t) { double d01 = sqrt(pow(x1 - x0, 2) + pow(y1 - y0, 2)); double d12 = sqrt(pow(x2 - x1, 2) + pow(y2 - y1, 2)); double fa = t * d01 / (d01 + d12); double fb = t * d12 / (d01 + d12); p1x = x1 - fa * (x2 - x0); p1y = y1 - fa * (y2 - y0); p2x = x1 + fb * (x2 - x0); p2y = y1 + fb * (y2 - y0); return; }
其中，（x0，y0）、（x1，y1）、（x2，y2）分別為P0、P1、P2三點的座標；t為曲率因子（取值範圍0-1.0），影響的是擬合曲線的曲率，後面將對其作進一步介紹。根據三點的座標和t，即可求得P1點的兩個控制點C10（p1x，p1y）、C11(p2x，p2y)。
以此類推，我們可分別求得P2、P3、P0等各點的控制點，如下圖所示。
接下來，我們我們逐段繪製Besier曲線。通過兩個頂點P0、P1和兩個控制點C01、C10，根據Bezier曲線擬合原理，即可獲得連線P0、P1兩點的曲線。
void get_bezier(double x1, double y1, double x2, double y2, double p12x, double p12y, double p21x, double p21y, std::vector<int>& vec_x, std::vector<int> & vec_y) { int prev_x = (int)round(x1); int prev_y = (int)round(y1); int last_x = (int)round(x2); int last_y = (int)round(y2); for (double s = 0.0; s < (1.0 + 0.00001); s += DELTA_S) { double J0 = pow(1 - s, 3); double J1 = pow(1 - s, 2) * s * 3; double J2 = pow(s, 2) * (1 - s) * 3; double J3 = pow(s, 3); double ptx = x1 * J0 + p12x * J1 + p21x * J2 + x2 * J3; double pty = y1 * J0 + p12y * J1 + p21y * J2 + y2 * J3; int iptx = (int)round(ptx); int ipty = (int)round(pty); vec_x.push_back(iptx); vec_y.push_back(ipty); } return; }
其中，DELTA_S是擬合的步長。再通過輪廓查詢演算法和插值演算法，即可得到一段完整的Besier曲線，如下圖所示。
同樣，以此類推，我們可以分別得到P1P2、P2P3、P3P0之間的曲線，將這幾段曲線連線在一起，即可得到一條完整的封閉曲線。
最後，我們再看一看前文提到的曲率因子t對擬合出來的曲線的影響。分別令 t = 0.0、0.2、0.5、0.8、1.0，得到的曲線分別如下圖所示。
下面是基於OpenCV的完整實現程式碼：
#include "spline_curve.h"
#include <vector>
#define DELTA_S 0.01
#define T 0.5
void get_control_point(cv::Point2d& point0, cv::Point2d& point1, cv::Point2d& point2,
cv::Point2d& c01, cv::Point2d& c12, double t)
{
double d01 = sqrt(pow(point1.x - point0.x, 2) + pow(point1.y - point0.y, 2));
double d12 = sqrt(pow(point2.x - point1.x, 2) + pow(point2.y - point1.y, 2));
double fa = t * d01 / (d01 + d12);
double fb = t * d12 / (d01 + d12);
c01.x = point1.x - fa * (point2.x - point0.x);
c01.y = point1.y - fa * (point2.y - point0.y);
c12.x = point1.x + fb * (point2.x - point0.x);
c12.y = point1.y + fb * (point2.y - point0.y);
return;
}
void get_control_points_array(std::vector<cv::Point2d>& key_points, std::vector<cv::Point2d>& vec_c01,
std::vector<cv::Point2d>& vec_c02, double t)
{
int N = key_points.size();
for (int i = 0; i < N; i++)
{
cv::Point2d c01, c02;
if (i == 0)
{
get_control_point(key_points[N - 1], key_points[i], key_points[i + 1], c01, c02, t);
vec_c01.push_back(c01);
vec_c02.push_back(c02);
}
else if (i < (N - 1))
{
get_control_point(key_points[i - 1], key_points[i], key_points[i + 1], c01, c02, t);
vec_c01.push_back(c01);
vec_c02.push_back(c02);
}
else
{
get_control_point(key_points[i - 1], key_points[i], key_points[0], c01, c02, t);
vec_c01.push_back(c01);
vec_c02.push_back(c02);
}
}
}
bool is_adjcent_point(cv::Point2i& point1, cv::Point2i& point2)
{
if (((point1.x == point2.x) && (point1.y == point2.y)) ||
(std::abs(point1.x - point2.x) > 1) || (std::abs(point1.y - point2.y) > 1))
{
return false;
}
return true;
}
bool is_same_point(cv::Point2i& point1, cv::Point2i& point2)
{
if ((point1.x == point2.x) && (point1.y == point2.y))
{
return true;
}
return false;
}
// interpolation between not adjacent points
void get_line_points(cv::Point2i& point1, cv::Point2i& point2, std::vector<cv::Point2i>& line_points)
{
line_points.push_back(point1);
int dx = abs(point1.x - point2.x);
int dy = abs(point1.y - point2.y);
if (dx == 0 && dy == 0)
{
return;
}
if (dx > dy)
{
if (point1.x < point2.x)
{
for (int i = point1.x + 1; i < point2.x; i++)
{
int y = (int)(((point1.y - point2.y + 0.0) / (point1.x - point2.x)) * (i - point1.x) + point1.y);
line_points.push_back(cv::Point2i(i, y));
}
}
else
{
for (int i = point1.x - 1; i > point2.x; i--)
{
int y = (int)(((point1.y - point2.y + 0.0) / (point1.x - point2.x)) * (i - point1.x) + point1.y);
line_points.push_back(cv::Point2i(i, y));
}
}
}
else
{
if (point1.y < point2.y)
{
for (int i = point1.y + 1; i < point2.y; i++)
{
int x = (int)(((point1.x - point2.x + 0.0) / (point1.y - point2.y)) * (i - point1.y) + point1.x);
line_points.push_back(cv::Point2i(x, i));
}
}
else
{
for (int i = point1.y - 1; i > point2.y; i--)
{
int x = (int)(((point1.x - point2.x + 0.0) / (point1.y - point2.y)) * (i - point1.y) + point1.x);
line_points.push_back(cv::Point2i(x, i));
}
}
}
line_points.push_back(point2);
return;
}
bool get_spline(cv::Point2d& point1, cv::Point2d& point2, cv::Point2d& c01,
cv::Point2d& c12, std::vector<cv::Point2i>& spline_points, double delta_s)
{
cv::Point2i point_prev = (cv::Point2i)point1;
cv::Point2i point_last = (cv::Point2i)point2;
spline_points.push_back(point_prev);
for (double s = 0.0; s < (1.0 + 0.0001); s += delta_s)
{
double J0 = pow(1 - s, 3);
double J1 = pow(1 - s, 2) * s * 3;
double J2 = pow(s, 2) * (1 - s) * 3;
double J3 = pow(s, 3);
double ptx = point1.x * J0 + c01.x * J1 + c12.x * J2 + point2.x * J3;
double pty = point1.y * J0 + c01.y * J1 + c12.y * J2 + point2.y * J3;
cv::Point2i ipoint;
ipoint.x = (int)round(ptx);
ipoint.y = (int)round(pty);
if (is_same_point(ipoint, point_last))
{
get_line_points(point_prev, point_last, spline_points);
break;
}
if (is_adjcent_point(point_prev, ipoint))
{
spline_points.push_back(ipoint);
point_prev = ipoint;
}
else if (is_same_point(point_prev, ipoint))
{
continue;
}
else
{
get_line_points(point_prev, ipoint, spline_points);
point_prev = ipoint;
}
}
return true;
}
void smooth_curve(std::vector<cv::Point2i>& curve_in, std::vector<cv::Point2i>& curve_out, bool is_closed)
{
int vec_size = curve_in.size();
for (int i = 0; i < (vec_size - 2); i += 2)
{
if (i == 0 && is_closed)
{
if (is_adjcent_point(curve_in[vec_size - 1], curve_in[1]))
{
curve_out.push_back(curve_in[1]);
}
else
{
curve_out.push_back(curve_in[0]);
curve_out.push_back(curve_in[1]);
}
}
if (is_adjcent_point(curve_in[i], curve_in[i + 2]))
{
curve_out.push_back(curve_in[i + 2]);
}
else
{
curve_out.push_back(curve_in[i + 1]);
curve_out.push_back(curve_in[i + 2]);
}
}
return;
}
bool get_spline_curve(std::vector<cv::Point2d>& key_points, std::vector<cv::Point2i>& spline_curve, double t, bool is_closed)
{
if (key_points.size() < 2)
{
std::cout << "Key points is less than two!!!" << std::endl;
return false;
}
if (key_points.size() == 2)
{
cv::Point2i point1 = (cv::Point2i)key_points[0];
cv::Point2i point2 = (cv::Point2i)key_points[1];
get_line_points(point1, point2, spline_curve);
return true;
}
std::vector<cv::Point2d> vec_c01, vec_c12;
get_control_points_array(key_points, vec_c01, vec_c12, t);
std::vector<cv::Point2i> temp_spline;
for (int i = 0; i < key_points.size(); i++)
{
if (i < (key_points.size() - 1))
{
get_spline(key_points[i], key_points[i + 1], vec_c12[i], vec_c01[i + 1], temp_spline, DELTA_S);
continue;
}
if (is_closed)
{
get_spline(key_points[i], key_points[0], vec_c12[i], vec_c01[0], temp_spline, DELTA_S);
}
}
smooth_curve(temp_spline, spline_curve, is_closed);
return true;
}2017.03.09完成初稿