C++-柱面拟合FitCylinder

Posted 2021-09-01 翟天保Steven

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了C++-柱面拟合FitCylinder相关的知识，希望对你有一定的参考价值。

场景需求

在各大领域的图像处理中，经常会有拟合面的需求，最常见的就是拟合斜平面，也有拟合曲面、球面、柱面的场景，本文介绍的就是拟合柱面。

柱面拟合公式：

基于该公式我们可以得到所要拟合面的各系数，其中C3就是powerx，C4就是powery，这两个参数是柱面分析中常见的评价参数。

在拟合面的过程中，有一个非常非常重要的注意事项，就是x和y的取值。在一个场内，需要先把x和y进行归一化处理。建立x和y的网格数据，类似于matlab的meshgrid，比如1000*1000的图像中，x的第一行数据都是-1，最后一行数据都是1，中间等间隔递增，y的第一列数据都是-1，最后一列数据都是1，中间同样等间隔递增。这样得到的系数数值与光学行业标准一致，比如ZYGO公司的拟合方法就是如此。

话不多说，下方为具体实现函数和测试代码。

功能函数代码

// 拟合柱面
cv::Mat FitCylinder(const cv::Mat&phase)
{
	cv::Mat cyc = ImageCropping(phase);
	cv::Mat mask = cv::Mat::zeros(cyc.size(), CV_8UC1);
	mask.setTo(255, cyc == cyc);
	cv::Mat x, y, ang, mag;
	UnitCart(cyc.cols, cyc.rows, x, y);
	UnitPolar(x, y, mag, ang);
	int samplingInterval = 1;
	vector<cv::Point> points;
	cv::findNonZero(mask, points);
	int pointnumber = static_cast<int>(points.size());
	samplingInterval = Max(samplingInterval, static_cast<int>(sqrt(pointnumber) / 100));

	// 抽样，提升计算速度,
	cv::Mat sampling_roi = GridSampling(mask.size(), samplingInterval, samplingInterval);
	sampling_roi.setTo(0, ~mask);
	// 得到抽样点的坐标
	std::vector<cv::Point> samplingidx_roi;
	cv::findNonZero(sampling_roi, samplingidx_roi);

	int len_sam = (int)samplingidx_roi.size();
	cv::Mat ang_sampling(len_sam, 1, CV_32FC1, 0.0f);
	cv::Mat mag_sampling(len_sam, 1, CV_32FC1, 0.0f);
	auto tmpSam = samplingidx_roi.begin();
	for (int i = 0; i < len_sam; ++i, ++tmpSam) {
		int x = tmpSam->x;
		int y = tmpSam->y;
		ang_sampling.ptr<float>(i)[0] = ang.ptr<float>(y)[x];
		mag_sampling.ptr<float>(i)[0] = mag.ptr<float>(y)[x];
	}
	cv::Mat phase_roi = get<float>(cyc, samplingidx_roi);

	cv::Mat dst(len_sam, 6, CV_32FC1, 0.0f);
	cv::Mat pow1 = mag_sampling;
	cv::Mat X = pow1.mul(cosf(ang_sampling));
	cv::Mat Y = pow1.mul(sinf(ang_sampling));
	cv::Mat X2 = pow(X, 2.0);
	cv::Mat Y2 = pow(Y, 2.0);
	cv::Mat XY = X.mul(Y);
	for (int i = 0; i < 6; ++i) {
		switch (i) {
		case 0: dst.col(i).setTo(1.0f); break; // 0
		case 1: X.copyTo(dst.col(i)); break; // 1
		case 2: Y.copyTo(dst.col(i)); break; // 2
		case 3: X2.copyTo(dst.col(i)); break; // 3
		case 4: Y2.copyTo(dst.col(i)); break; // 4
		case 5: XY.copyTo(dst.col(i)); break; // 5
		default: break;
		}
	}
	cv::Mat result;
	cv::solve(dst, phase_roi, result, cv::DECOMP_NORMAL);    //求解方程A’A * dst = A'B中的dst,
	cv::Mat temp = result.clone();
	return result;
}

// 读取excel图像数据
cv::Mat ReadPicFromExcel(string name)
{
	cv::Mat result, pic;
	
	ifstream infile(name);
	string str;
	int col = 0;
	while (getline(infile, str))
	{
		string temp;
		stringstream input(str);
		col = 0;
		while (input >> temp)
		{
			if (temp == "nan")
			{
				pic.push_back((nan("")));
			}
			else {
				pic.push_back((atof(temp.c_str())));
			}
			col++;
		}
	}
	int number = result.rows;
	int row = number / col;
	result = pic.reshape(row, col);
	infile.close();
	return result;
}

// 图像裁剪
cv::Mat ImageCropping(const cv::Mat &phase) {
	// 非测量区一般都进行了NaN处理，所以掩膜绘制只需要判断是否为NaN值即可
	cv::Mat mask = cv::Mat::zeros(phase.size(), CV_8UC1);
	mask.setTo(255, phase == phase);
	int roi_up = 10000;
	int roi_down = 0;
	int roi_left = 10000;
	int roi_right = 0;
	int row = phase.rows;
	int col = phase.cols;
	for (int i = 0; i < row; i++)
	{
		uchar *m = mask.ptr<uchar>(i);
		for (int j = 0; j < col; j++)
		{
			if (m[j] != 0)
			{
				if (j < roi_left)roi_left = j;
				if (j > roi_right)roi_right = j;
				if (i < roi_up)roi_up = i;
				if (i > roi_down)roi_down = i;
			}
		}
	}
	int w = roi_right - roi_left;
	int h = roi_down - roi_up;
	// 一般提取奇数尺寸，方便计算
	if (w % 2 == 0)w++;
	if (h % 2 == 0)h++;
	cv::Mat crop_phase = phase(cv::Rect(roi_left, roi_up, w, h)).clone();
	return crop_phase;
}

// meshgrid
void UnitCart(int squaresizex, int squaresizey, cv::Mat& x, cv::Mat& y) {
	CV_Assert(squaresizex % 2 == 1);
	CV_Assert(squaresizey % 2 == 1);
	x.create(squaresizey, squaresizex, CV_32FC1);
	y.create(squaresizey, squaresizex, CV_32FC1);
	//设置边界
	x.col(0).setTo(-1.0);
	x.col(squaresizex - 1).setTo(1.0f);
	y.row(0).setTo(1.0);
	y.row(squaresizey - 1).setTo(-1.0f);

	float deltax = 2.0f / (squaresizex - 1.0f);  //两个元素的间隔
	float deltay = 2.0f / (squaresizey - 1.0f);  //两个元素的间隔
	//计算其他位置的值
	for (int i = 1; i < squaresizex - 1; ++i) {
		x.col(i) = -1.0f + i * deltax;
	}
	for (int i = 1; i < squaresizey - 1; ++i) {
		y.row(i) = 1.0f - i * deltay;
	}
}

// 极坐标转化
void UnitPolar(cv::Mat& x, cv::Mat& y, cv::Mat& mag, cv::Mat& ang) {
	//cv::cartToPolar(x, y, mag, ang, indegree);   //直角坐标转换为极坐标
	mag = cv::Mat(x.size(), x.type());
	ang = cv::Mat(x.size(), x.type());
	int row = mag.rows;
	int col = mag.cols;
	for (int i = 0; i < row; ++i)
	{
		float*m = mag.ptr<float>(i);
		float*a = ang.ptr<float>(i);
		float*xx = x.ptr<float>(i);
		float*yy = y.ptr<float>(i);
		for (int j = 0; j < col; ++j)
		{
			m[j] = sqrt(xx[j] * xx[j] + yy[j] * yy[j]);
			a[j] = atan2(yy[j], xx[j]);
		}
	}
}

// 采样提速
cv::Mat GridSampling(const cv::Size& size, int rowinterval, int colinterval) {
	cv::Mat dst(size, CV_8UC1, cv::Scalar(0));
	//设置采样的位置点
	int Row = dst.rows;
	int Col = dst.cols;
	for (int row = 0; row < Row; row += rowinterval) {
		for (int col = 0; col < Col; col += colinterval) {
			dst.at<uchar>(row, col) = 255;
		}
	}
	return dst;
}

// 获取计算数据
template <typename T>
cv::Mat get(const cv::Mat& src,const std::vector<cv::Point>& idx) 
{
	int num = (int)idx.size();
	cv::Mat dst(num, 1, src.type());

	/* pragma omp parallel for 是OpenMP中的一个指令，
	表示接下来的for循环将被多线程执行，另外每次循环之间不能有关系 */
#pragma omp parallel for
	for (int i = 0; i < num; ++i) {
		dst.at<T>(i, 0) = src.at<T>(idx[i]);
	}
	return dst;
}

// 图像数据cos处理
cv::Mat cosf(const cv::Mat& src) {
	CV_Assert(src.type() == CV_32FC1);
	cv::Mat dst(src.size(), src.type());

	int cols = src.cols;
	int rows = src.rows;

	//返回bool值，判断存储是否连续。
	if (src.isContinuous() && dst.isContinuous())
	{
		cols *= rows;
		rows = 1;
	}
	//计算每个元素的cos()
	for (int i = 0; i < rows; i++)
	{
		const float* srci = src.ptr<float>(i);
		float* dsti = dst.ptr<float>(i);
		for (int j = 0; j < cols; j++) {
			dsti[j] = std::cosf(srci[j]);
		}
	}
	return dst;
}

// 图像数据sin处理
cv::Mat sinf(const cv::Mat& src) {
	CV_Assert(src.type() == CV_32FC1);
	cv::Mat dst(src.size(), src.type());

	int cols = src.cols;
	int rows = src.rows;

	//返回bool值，判断存储是否连续。
	if (src.isContinuous() && dst.isContinuous())
	{
		cols *= rows;
		rows = 1;
	}
	//计算每个元素的sin()
	for (int i = 0; i < rows; i++)
	{
		const float* srci = src.ptr<float>(i);
		float* dsti = dst.ptr<float>(i);
		for (int j = 0; j < cols; j++) {
			dsti[j] = std::sinf(srci[j]);
		}
	}
	return dst;
}

// 图像数据平方处理
cv::Mat pow(cv::InputArray src, double power) {
	cv::Mat dst;
	cv::pow(src, power, dst);
	return dst;
}

C++测试代码

#include<iostream>
#include<opencv2/opencv.hpp>
#include<ctime>
#include<string>
#include<sstream>
#include<fstream>
using namespace std;
using namespace cv;
#define Max(a, b) a > b ? a : b

cv::Mat FitCylinder(const cv::Mat&phase);
cv::Mat ReadPicFromExcel(string name);
cv::Mat ImageCropping(const cv::Mat &phase);
void UnitCart(int squaresizex, int squaresizey, cv::Mat& x, cv::Mat& y);
void UnitPolar(cv::Mat& x, cv::Mat& y, cv::Mat& mag, cv::Mat& ang);
cv::Mat GridSampling(const cv::Size& size, int rowinterval, int colinterval);
template <typename T>
cv::Mat get(const cv::Mat& src, const std::vector<cv::Point>& idx);
cv::Mat cosf(const cv::Mat& src);
cv::Mat sinf(const cv::Mat& src);
cv::Mat pow(cv::InputArray src, double power);

int main(void)
{
	cv::Mat phase = ReadPicFromExcel("test.xls");
	cv::Mat coef = FitCylinder(phase);
	for (int i = 0; i < coef.rows; ++i)
	{
		cout << "coef " << i << " : " << coef.at<float>(i, 0) << endl;
	}
	system("pause");
	return 0;
}

// 拟合柱面
cv::Mat FitCylinder(const cv::Mat&phase)
{
	cv::Mat cyc = ImageCropping(phase);
	cv::Mat mask = cv::Mat::zeros(cyc.size(), CV_8UC1);
	mask.setTo(255, cyc == cyc);
	cv::Mat x, y, ang, mag;
	UnitCart(cyc.cols, cyc.rows, x, y);
	UnitPolar(x, y, mag, ang);
	int samplingInterval = 1;
	vector<cv::Point> points;
	cv::findNonZero(mask, points);
	int pointnumber = static_cast<int>(points.size());
	samplingInterval = Max(samplingInterval, static_cast<int>(sqrt(pointnumber) / 100));

	// 抽样，提升计算速度,
	cv::Mat sampling_roi = GridSampling(mask.size(), samplingInterval, samplingInterval);
	sampling_roi.setTo(0, ~mask);
	// 得到抽样点的坐标
	std::vector<cv::Point> samplingidx_roi;
	cv::findNonZero(sampling_roi, samplingidx_roi);

	int len_sam = (int)samplingidx_roi.size();
	cv::Mat ang_sampling(len_sam, 1, CV_32FC1, 0.0f);
	cv::Mat mag_sampling(len_sam, 1, CV_32FC1, 0.0f);
	auto tmpSam = samplingidx_roi.begin();
	for (int i = 0; i < len_sam; ++i, ++tmpSam) {
		int x = tmpSam->x;
		int y = tmpSam->y;
		ang_sampling.ptr<float>(i)[0] = ang.ptr<float>(y)[x];
		mag_sampling.ptr<float>(i)[0] = mag.ptr<float>(y)[x];
	}
	cv::Mat phase_roi = get<float>(cyc, samplingidx_roi);

	cv::Mat dst(len_sam, 6, CV_32FC1, 0.0f);
	cv::Mat pow1 = mag_sampling;
	cv::Mat X = pow1.mul(cosf(ang_sampling));
	cv::Mat Y = pow1.mul(sinf(ang_sampling));
	cv::Mat X2 = pow(X, 2.0);
	cv::Mat Y2 = pow(Y, 2.0);
	cv::Mat XY = X.mul(Y);
	for (int i = 0; i < 6; ++i) {
		switch (i) {
		case 0: dst.col(i).setTo(1.0f); break; // 0
		case 1: X.copyTo(dst.col(i)); break; // 1
		case 2: Y.copyTo(dst.col(i)); break; // 2
		case 3: X2.copyTo(dst.col(i)); break; // 3
		case 4: Y2.copyTo(dst.col(i)); break; // 4
		case 5: XY.copyTo(dst.col(i)); break; // 5
		default: break;
		}
	}
	cv::Mat result;
	cv::solve(dst, phase_roi, result, cv::DECOMP_NORMAL);    //求解方程A’A * dst = A'B中的dst,
	cv::Mat temp = result.clone();
	return result;
}

// 读取excel图像数据
cv::Mat ReadPicFromExcel(string name)
{
	cv::Mat result, pic;
	
	ifstream infile(name);
	string str;
	int col = 0;
	while (getline(infile, str))
	{
		string temp;
		stringstream input(str);
		col = 0;
		while (input >> temp)
		{
			if (temp == "nan")
			{
				pic.push_back((nan("")));
			}
			else {
				pic.push_back((atof(temp.c_str())));
			}
			col++;
		}
	}
	int number = result.rows;
	int row = number / col;
	result = pic.reshape(row, col);
	infile.close();
	return result;
}

// 图像裁剪
cv::Mat ImageCropping(const cv::Mat &phase) {
	// 非测量区一般都进行了NaN处理，所以掩膜绘制只需要判断是否为NaN值即可
	cv::Mat mask = cv::Mat::zeros(phase.size(), CV_8UC1);
	mask.setTo(255, phase == phase);
	int roi_up = 10000;
	int roi_down = 0;
	int roi_left = 10000;
	int roi_right = 0;
	int row = phase.rows;
	int col = phase.cols;
	for (int i = 0; i < row; i++)
	{
		uchar *m = mask.ptr<uchar>(i);
		for (int j = 0; j < col; j++)
		{
			if (m[j] != 0)
			{
				if (j < roi_left)roi_left = j;
				if (j > roi_right)roi_right = j;
				if (i < roi_up)roi_up = i;
				if (i > roi_down)roi_down = i;
			}
		}
	}
	int w = roi_right - roi_left;
	int h = roi_down - roi_up;
	// 一般提取奇数尺寸，方便计算
	if (w % 2 == 0)w++;
	if (h % 2 == 0)h++;
	cv::Mat crop_phase = phase(cv::Rect(roi_left, roi_up, w, h)).clone();
	return crop_phase;
}

// meshgrid
void UnitCart(int squaresizex, int squaresizey, cv::Mat& x, cv::Mat& y) {
	CV_Assert(squaresizex % 2 == 1);
	CV_Assert(squaresizey % 2 == 1);
	x.create(squaresizey, squaresizex, CV_32FC1);
	y.create(squaresizey, squaresizex, CV_32FC1);
	//设置边界
	x.col(0).setTo(-1.0);
	x.col(squaresizex - 1).setTo(1.0f);
	y.row(0).setTo(1.0);
	y.row(squaresizey - 1).setTo(-1.0f);

	float deltax = 2.0f / (squaresizex - 1.0f);  //两个元素的间隔
	float deltay = 2.0f / (squaresizey - 1.0f);  //两个元素的间隔
	//计算其他位置的值
	for (int i = 1; i < squaresizex - 1; ++i) {
		x.col(i) = -1.0f + i * deltax;
	}
	for (int i = 1; i < squaresizey - 1; ++i) {
		y.row(i) = 1.0f - i * deltay;
	}
}

// 极坐标转化
void UnitPolar(cv::Mat& x, cv::Mat& y, cv::Mat& mag, cv::Mat& ang) {
	//cv::cartToPolar(x, y, mag, ang, indegree);   //直角坐标转换为极坐标
	mag = cv::Mat(x.size(), x.type());
	ang = cv::Mat(x.size(), x.type());
	int row = mag.rows;
	int col = mag.cols;
	for (int i = 0; i < row; ++i)
	{
		float*m = mag.ptr<float>(i);
		float*a = ang.ptr<float>(i);
		float*xx = x.ptr<float>(i);
		float*yy = y.ptr<float>(i);
		for (int j = 0; j < col; ++j)
		{
			m[j] = sqrt(xx[j] * xx[j] + yy[j] * yy[j]);
			a[j] = atan2(yy[j], xx[j]);
		}
	}
}

// 采样提速
cv::Mat GridSampling(const cv::Size& size, int rowinterval, int colinterval) {
	cv::Mat dst(size, CV_8UC1, cv::Scalar(0));
	//设置采样的位置点
	int Row = dst.rows;
	int Col = dst.cols;
	for (int row = 0; row < Row; row += rowinterval) {
		for (int col = 0; col < Col; col += colinterval) {
			dst.at<uchar>(row, col) = 255;
		}
	}
	return dst;
}

// 获取计算数据
template <typename T>
cv::Mat get(const cv::Mat& src,const std::vector<cv::Point>& idx) 
{
	int num = (int)idx.size();
	cv::Mat dst(num, 1, src.type());

	/* pragma omp parallel for 是OpenMP中的一个指令，
	表示接下来的for循环将被多线程执行，另外每次循环之间不能有关系 */
#pragma omp parallel for
	for (int i = 0; i < num; ++i) {
		dst.at<T>(i, 0) = src.at<T>(idx[i]);
	}
	return dst;
}

// 图像数据cos处理
cv::Mat cosf(const cv::Mat& src) {
	CV_Assert(src.type() == CV_32FC1);
	cv::Mat dst(src.size(), src.type());

	int cols = src.cols;
	int rows = src.rows;

	//返回bool值，判断存储是否连续。
	if (src.isContinuous() && dst.isContinuous())
	{
		cols *= rows;
		rows = 1;
	}
	//计算每个元素的cos()
	for (int i = 0; i < rows; i++)
	{
		const float* srci = src.ptr<float>(i);
		float* dsti = dst.ptr<float>(i);
		for (int j = 0; j < cols; j++) {
			dsti[j] = std::cosf(srci[j]);
		}
	}
	return dst;
}

// 图像数据sin处理
cv::Mat sinf(const cv::Mat& src) {
	CV_Assert(src.type() == CV_32FC1);
	cv::Mat dst(src.size(), src.type());

	int cols = src.cols;
	int rows = src.rows;

	//返回bool值，判断存储是否连续。
	if (src.isContinuous() && dst.isContinuous())
	{
		cols *= rows;
		rows = 1;
	}
	//计算每个元素的sin()
	for (int i = 0; i < rows; i++)
	{
		const float* srci = src.ptr<float>(i);
		float* dsti = dst.ptr<float>(i);
		for (int j = 0; j < cols; j++) {
			dsti[j] = std::sinf(srci[j]);
		}
	}
	return dst;
}

// 图像数据平方处理
cv::Mat pow(cv::InputArray src, double power) {
	cv::Mat dst;
	cv::pow(src, power, dst);
	return dst;
}

测试效果

在测试案例中，我加载了一个柱面数据，因为干涉测量的结果数值一般都很低，所以我将数值都乘了100，这样得到的灰度图像能明显看出来黑白区别，如图1所示，在图4中也可以看出那个倾斜度，白色的区域就是数值大于0的部分；因为我在VS中截取的区域和在我们软件中截取的区域并不完全一致，所以powerx和powery的数值有所差异，只要量级一致就说明没有错误；如图2所示，VS中输出的结果除以100再对比，coef3就是powerx，coef4就是powery。

如果函数有什么可以改进完善的地方，非常欢迎大家指出，一同进步何乐而不为呢~

如果文章帮助到你了，可以点个赞让我知道，我会很快乐~加油！

以上是关于C++-柱面拟合FitCylinder的主要内容，如果未能解决你的问题，请参考以下文章

UG NX二次开发（C#）--建模--识别曲面类型（圆柱面）

在matlab画出三维球面并绕轴旋转一定角度

C++曲线拟合代码

一般将硬盘分为五个部分，( )通常位于硬盘的0磁道1柱面1扇区。

用matlab生成在单位球面上均匀的100个点几种方法？