场景需求
在各大领域的图像处理中，经常会有拟合面的需求，最常见的就是拟合斜平面，也有拟合曲面、球面、柱面的场景，本文介绍的就是拟合柱面。
柱面拟合公式：
基于该公式我们可以得到所要拟合面的各系数，其中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 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 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 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 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; }
测试效果
图1 柱面灰度图像
图2 拟合结果
图3 3维直观图
图4 数据对比
在测试案例中，我加载了一个柱面数据，因为干涉测量的结果数值一般都很低，所以我将数值都乘了100，这样得到的灰度图像能明显看出来黑白区别，如图1所示，在图4中也可以看出那个倾斜度，白色的区域就是数值大于0的部分；因为我在VS中截取的区域和在我们软件中截取的区域并不完全一致，所以powerx和powery的数值有所差异，只要量级一致就说明没有错误；如图2所示，VS中输出的结果除以100再对比，coef3就是powerx，coef4就是powery。
如果函数有什么可以改进完善的地方，非常欢迎大家指出，一同进步何乐而不为呢~
如果文章帮助到你了，可以点个赞让我知道，我会很快乐~加油！
