作者：翟天保Steven
版权声明：著作权归作者所有，商业转载请联系作者获得授权，非商业转载请注明出处

场景需求

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

柱面拟合公式：

$Z=C_{0}+C_{1}x+C_{2}y+C_{3}x^{2}+C_{4}y^{2}+C_{5}xy$

基于该公式我们可以得到所要拟合面的各系数，其中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; // 0case 1: X.copyTo(dst.col(i)); break; // 1case 2: Y.copyTo(dst.col(i)); break; // 2case 3: X2.copyTo(dst.col(i)); break; // 3case 4: Y2.copyTo(dst.col(i)); break; // 4case 5: XY.copyTo(dst.col(i)); break; // 5default: 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 row = pic.rows / col;result = pic.reshape(0, row);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 forfor (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 : bcv::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; // 0case 1: X.copyTo(dst.col(i)); break; // 1case 2: Y.copyTo(dst.col(i)); break; // 2case 3: X2.copyTo(dst.col(i)); break; // 3case 4: Y2.copyTo(dst.col(i)); break; // 4case 5: XY.copyTo(dst.col(i)); break; // 5default: 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 row = pic.rows / col;result = pic.reshape(0, row);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 forfor (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相关推荐

RANSAC与圆柱拟合
RANSAC是"RANdom SAmple Consensus(随机抽样一致)"的缩写.它可以从一组包含"局外点"的观测数据集中,通过迭代方式估计数学模型的参数 ...
点云拟合—圆柱面非线性最小二乘实现
最近项目需要这方面的工作,于是开始研究这个了: 圆柱几何特征:圆柱面上的点到其轴线的距离恒等于半径圆柱的方程: 首先是 PCL库自带的圆柱模型拟合,由于在查找最佳圆柱面的过程中会过滤很多点,因此考 ...
python中使用scipy.optimize.leastsq进行球面、圆柱面、平面拟合
scipy.optimize.leastsq官方文档之前自己用最小二乘进行过球面.柱面和平面拟合,球面和平面效果都很好,但是柱面效果一直不好(对初值十分敏感,迭代过程中经常出现奇异矩阵,也有可能是自 ...
【论文笔记】基于点云柱面投影图的平面特征提取SLAM：　Fast planar surface 3D SLAM using LIDAR 2017
论文首先将三维点云投影在了三个平面上得到无直线畸变的投影视图,然后利用<Fast segmentation of range imagery into planar regions>中的方 ...
机器学习概念 — 监督学习、无监督学习、半监督学习、强化学习、欠拟合、过拟合、后向传播、损失和优化函数、计算图、正向传播、反向传播
1. 监督学习和无监督学习监督学习 ( Supervised Learning ) 和无监督学习 ( Unsupervised Learning ) 是在机器学习中经常被提及的两个重要的学习方法. ...
用python的numpy作线性拟合、多项式拟合、对数拟合
转自:http://blog.itpub.net/12199764/viewspace-1743145/ 项目中有涉及趋势预测的工作,整理一下这3种拟合方法: 1.线性拟合-使用math import ...
RANSAC算法（2）：（拟合平面）本文以地面为基础以及源码分布解读
本章代码是本人根据一个未曾谋面的好人学习的(要怀抱希望,世界上好人真的是很多的,我要做一个去给别人带去正能量积极态度的人,加油喽),如需转载学习请注明.谢谢 ---------------基于rans ...
PCL：拟合平面直线和曲线以及空间曲线的原理到算法实现
使用两种思路进行直线拟合: 1.利用逆矩阵思想 --------------进行下列公式的推导需要理解逆矩阵(求A矩阵的逆矩阵,则A矩阵必须是方阵)的知识: (1)为什么要引入逆矩阵呢? 逆矩阵可以类 ...
【camera】自动驾驶感知系统实现（车道线检测和拟合、目标检测与跟踪、道路可行驶区域分割、深度估计、图像视野到BEV空间映射、像平面到地平面映射）
自动驾驶感知系统实现(车道线检测和拟合.目标检测与跟踪.道路可行驶区域分割.深度估计.图像视野到BEV空间映射.像平面到地平面映射) 项目下载地址:项目下载地址推理引擎下载地址:推理引擎下载地址支 ...
CS131专题-4：拟合（最小二乘、RANSAC、霍夫变换）
本专题目的:了解最小二乘.RANSAC.霍夫变换这3个算法的基本原理,能够做到脱口而出,并从零编程实现. 目录 1 前言 2 最小二乘 2.1 基本原理 2.2 求解方法 3 RANSAC 算法 3. ...

C++-柱面拟合FitCylinder

场景需求

功能函数代码

C++测试代码

测试效果

C++-柱面拟合FitCylinder相关推荐

最新文章

热门文章