AI 助理
备案控制台
开发者社区 开发与运维 文章 正文

C++-柱面拟合FitCylinder

2023-10-17 81
版权
版权声明：
本文内容由阿里云实名注册用户自发贡献，版权归原作者所有，阿里云开发者社区不拥有其著作权，亦不承担相应法律责任。具体规则请查看《 阿里云开发者社区用户服务协议》和 《阿里云开发者社区知识产权保护指引》。如果您发现本社区中有涉嫌抄袭的内容，填写 侵权投诉表单进行举报，一经查实，本社区将立刻删除涉嫌侵权内容。
简介： 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 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。


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


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

文章标签：
C++
计算机视觉
翟天保Steven
目录
相关文章
闲谈社
MATALB运用——最小二乘法拟合
MATALB运用——最小二乘法拟合
闲谈社
144 0
aliyun7689123603-22772
|
4月前
|
机器学习/深度学习
避免最小二乘法导致的过拟合问题
8月更文挑战第12天
aliyun7689123603-22772
92 4
米果粒
|
5月前
|
机器学习/深度学习 算法
欠拟合
【7月更文挑战第25天】欠拟合。
米果粒
46 2
拓端数据部落
|
7月前
R方和线性回归拟合优度
R方和线性回归拟合优度
拓端数据部落
51 3
拓端数据部落
|
7月前
R语言Poisson回归的拟合优度检验
R语言Poisson回归的拟合优度检验
拓端数据部落
70 4
Hacoj
|
7月前
|
机器学习/深度学习
欠拟合与过拟合
欠拟合与过拟合
Hacoj
74 0
拓端数据部落
|
7月前
|
C++
拟合R语言中的多项式回归
拟合R语言中的多项式回归
拓端数据部落
32 0
星辰同学wwq
|
API
三、欠拟合和过拟合
三、欠拟合和过拟合
星辰同学wwq
71 0
Deephub
线性回归中的L1与L2正则化
线性回归中的L1与L2正则化
Deephub
197 0
线性回归中的L1与L2正则化
初商
|
算法 机器学习/深度学习
秒懂“线性回归预测”
线性回归是机器学习中的概念，线性回归预测算法一般用以解决“使用已知样本对未知公式参数的估计”类问题。
初商
942 0