三次B样条曲线插值c++

最新推荐文章于 2022-03-05 16:39:40 发布

原创最新推荐文章于 2022-03-05 16:39:40 发布 · 3.9k 阅读

41 ·

CC 4.0 BY-SA版权

文章标签：

#计算机视觉

本文介绍了如何使用C++实现三次B样条曲线的插值方法。通过离散点反求控制点，然后利用控制点进行拟合。具体原理参考了相关文章，并提供了实现代码，最终展示出插值曲线的效果。

B样条拟合通常把离散点集作为控制点来求拟合曲线，原理和过程可参看https://blog.youkuaiyun.com/cnmgbmsdn/article/details/108141194。
如果想要求得过离散点的插值曲线，需要先用离散点反求控制点，然后用求得的控制点拟合曲线。原理在这篇文章中介绍较详细http://www.doc88.com/p-5714423317458.html
在这里插入图片描述如图，P是控制点，V是离散点，左侧参数矩阵的逆矩阵乘以V矩阵就能得到P矩阵。c++代码如下：

/*
通过离散点集反求控制点
*/
vector<Point2f> BSplineInterpolate(vector<Point2f> discretePoints) {
	int N = discretePoints.size();
	Mat weight = Mat::zeros(Size(N,N), CV_32F);
	for (int i = 0; i < weight.rows; i++) {
			weight.ptr<float>(i)[i] = 1;
			weight.ptr<float>(i)[(i+1)%weight.cols] = 4;
			weight.ptr<float>(i)[(i+2)%weight.cols] = 1;
	}
	
	Mat V = Mat::zeros(Size(2, N), CV_32F);
	float* vdata = (float*)V.data;
	int idx = 0;
	for_each(discretePoints.begin(), discretePoints.end(), [&vdata, &idx](Point2f p2f) {vdata[idx++] = p2f.x; vdata[idx++] = p2f.y; });
	V *= 6;
	Mat P = weight.inv()*V;
	vector<Point2f> res(discretePoints.size());
	float* pdata = (float*)P.data;
	idx = 0;
	for_each(res.begin(), res.end(), [&idx, pdata](Point2f& p2f) {p2f.x = pdata[idx++]; p2f.y=pdata[idx++]; });
	return res;
}

求得控制点矩阵之后使用B样条曲线拟合函数可得到插值曲线图。

/*B样条曲线拟合
@return 返回拟合得到的曲线
@discretePoints 输入的离散点，至少4个点
@closed 是否拟合闭合曲线，true表示闭合，false不闭合
@stride 拟合精度
*/
vector<Point2f> BSplineFit(vector<Point2f> discretePoints, bool closed, double stride) {
	vector<Point2f> fittingPoints;
	for (int i = 0; i < (closed ? discretePoints.size() : discretePoints.size() - 1); i++) {
		Point2f xy[4];
		xy[0] = (discretePoints[i] + 4 * discretePoints[(i + 1) % discretePoints.size()] + discretePoints[(i + 2) % discretePoints.size()]) / 6;
		xy[1] = -(discretePoints[i] - discretePoints[(i + 2) % discretePoints.size()]) / 2;
		xy[2] = (discretePoints[i] - 2 * discretePoints[(i + 1) % discretePoints.size()] + discretePoints[(i + 2) % discretePoints.size()]) / 2;
		xy[3] = -(discretePoints[i] - 3 * discretePoints[(i + 1) % discretePoints.size()] + 3 * discretePoints[(i + 2) % discretePoints.size()] - discretePoints[(i + 3) % discretePoints.size()]) / 6;
		for (double t = 0; t <= 1; t += stride) {
			Point2f totalPoints = Point2f(0, 0);
			for (int j = 0; j < 4; j++) {
				totalPoints += xy[j] * pow(t, j);
			}
			fittingPoints.push_back(totalPoints);
		}
	}
	return fittingPoints;
}