近期，因為要實現經過多個控制點的曲線擬合，研究起了曲線擬合算法。綜合搜索到的資料，發現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.5void 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完成初稿

總結

以上是生活随笔為你收集整理的【算法+OpenCV】基于三次Bezier原理的曲线拟合算法C++与OpenCV实现的全部內容，希望文章能夠幫你解決所遇到的問題。

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