作者：翟天保Steven
版权声明：著作权归作者所有，商业转载请联系作者获得授权，非商业转载请注明出处
场景需求
在各大领域的图像处理中，经常会有拟合面的需求，最常见的就是拟合斜平面，也有拟合曲面、球面、柱面的场景，本文介绍的就是拟合柱面。
柱面拟合公式：
基于该公式我们可以得到所要拟合面的各系数，其中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。
如果函数有什么可以改进完善的地方，非常欢迎大家指出，一同进步何乐而不为呢~
如果文章帮助到你了，可以点个赞让我知道，我会很快乐~加油！
