在自動駕駛系統(tǒng)中，通常會用起點、終點和一個三階多項式來表示一條車道線，多項式系數(shù)的求解一般用最小二乘法來實現(xiàn)。

本文首先介紹兩種基于最小二乘法的多項式擬合方法的原理，然后基于OpenCV用c++編寫了這兩種擬合方法的代碼，最后通過一個完整的示例來展示如何通過一個離散點集擬合出一條多項式曲線。

多項式曲線擬合是指基于一系列的觀測點去尋找一個多項式來表示這些點的關(guān)系，最小二乘法通過最小化誤差的平方和去尋找數(shù)據(jù)的最佳匹配函數(shù)。假設(shè)有點集

其中，和的關(guān)系滿足函數(shù)

假設(shè)次多項式函數(shù)為

其中為多項式系數(shù)。如果要用這個次多項式來表示與的關(guān)系，那么多項式值與真實值之間的誤差為

采用最小二乘法進行多項式擬合的目的就是尋找一組最佳的多項式系數(shù)使得擬合后整個點集的總誤差最小，而求總誤差最小的問題可以轉(zhuǎn)化為求誤差平方和最小。整個點集的誤差平方和為

要使最小，可以對()求偏導(dǎo)并令其為零：

可得

寫成矩陣形式可得

把上式寫為，那么只要計算出

和

，再代入上面的式子中構(gòu)建出矩陣和，那么多項式系數(shù)就可以通過求解線性方程組得到。

從上一節(jié)的推導(dǎo)過程可以看到，用代數(shù)法求解多項式系數(shù)的方法非常繁瑣而且計算量比較大，我們其實還可以用矩陣求解的方法來簡化計算過程。

令

那么整個點集的誤差平方和可以表示為

對求導(dǎo)可得

令導(dǎo)數(shù)為零，那么可以得到

解得

學(xué)習(xí)了前面的理論知識，我們用c++來編碼實現(xiàn)一下，畢竟光看一堆數(shù)學(xué)公式還是不夠直觀。從前面的理論知識可以知道，用最小二乘法做多項式曲線擬合其實質(zhì)就是一個矩陣求解的問題，我們可以借助一些常用的數(shù)學(xué)庫比如OpenCV或Eigen來求解。本文將采用OpenCV來實現(xiàn)。

如果理解了前面的理論，寫代碼其實就比較簡單了。對照前面理論部分的公式，構(gòu)建好矩陣A和B，然后調(diào)用OpenCV的solve()函數(shù)來求解即可：

// 代數(shù)方式求解
void PolyFit(const std::vector<cv::Point2f> &points, const int order,
             cv::Mat *coeff) {
  const int n = points.size();
  cv::Mat A = cv::Mat::zeros(order + 1, order + 1, CV_64FC1);
  cv::Mat B = cv::Mat::zeros(order + 1, 1, CV_64FC1);

  // 構(gòu)建A矩陣
  for (int i = 0; i < order + 1; ++i) {
    for (int j = 0; j < order + 1; ++j) {
      for (int k = 0; k < n; ++k) {
        A.at<double>(i, j) += std::pow(points.at(k).x, i + j);
      }
    }
  }

  // 構(gòu)建B矩陣
  for (int i = 0; i < order + 1; ++i) {
    for (int k = 0; k < n; ++k) {
      B.at<double>(i, 0) += std::pow(points.at(k).x, i) * points.at(k).y;
    }
  }

  (*coeff) = cv::Mat::zeros(order + 1, 1, CV_64FC1);
  // 求解
  if (!cv::solve(A, B, *coeff, cv::DECOMP_LU)) {
    std::cout << "Failed to solve !" << std::endl;
  }
}

矩陣方式求解的代碼更加簡單：

// 矩陣方式求解
void PolyFit(const std::vector<cv::Point2f> &points, const int order,
             cv::Mat *coeff) {
  const int n = points.size();
  cv::Mat A = cv::Mat::ones(n, order + 1, CV_64FC1);
  cv::Mat B = cv::Mat::zeros(n, 1, CV_64FC1);

  for (int i = 0; i < n; ++i) {
    const double a = points.at(i).x;
    const double b = points.at(i).y;
    B.at<double>(i, 0) = b;
    if (i > 0) {
      for (int j = 1, v = a; j < order + 1; ++j, v *= a) {
        A.at<double>(i, j) = v;
      }
    }
  }

  (*coeff) = (A.t() * A).inv() * A.t() * B;
}

接下來，我們來看一個簡單的曲線擬合示例：

int main(int argc, char **argv) {
  constexpr int kOrder = 3; // 多項式階數(shù)
  constexpr int kWidth = 1000;
  constexpr int kHeight = 500;

  cv::Mat canvas =
      cv::Mat(cv::Size(kWidth, kHeight), CV_8UC3, cv::Scalar(255, 255, 255));

  cv::RNG rng(0xFFFFFFFF); // 隨機數(shù)

  // 生成點集，y坐標(biāo)添加一些隨機噪聲
  std::vector<cv::Point2f> raw_points;
  for (int i = 10; i < kWidth; i += 10) {
    cv::Point2f p;
    p.x = i;
    const auto noise = rng.uniform(-kHeight / 10, kHeight / 10);
    p.y = kHeight - p.x / kWidth * kHeight + noise;

    cv::circle(canvas, p, 5, cv::Scalar(0, 0, 255), -1);
    raw_points.emplace_back(p);
  }

  // 多項式擬合
  cv::Mat coeff;
  PolyFit(raw_points, kOrder, &coeff);

  // 用擬合后的系數(shù)重新生成點集
  std::vector<cv::Point> poly_points;
  for (const auto &rp : raw_points) {
    cv::Point p;
    p.x = rp.x;
    p.y = PolyValue(coeff, kOrder, rp.x);

    cv::circle(canvas, p, 5, cv::Scalar(0, 255, 0), -1);
    poly_points.emplace_back(p);
  }
  
  cv::polylines(canvas, poly_points, false, cv::Scalar(0, 255, 0), 3,
                cv::LINE_4);

  cv::imshow("PolyFit", canvas);
  cv::waitKey(0);
  cv::destroyAllWindows();

  return 0;
}

代碼的流程如下：

1. 生成一個離散點集，其中每個點的y坐標(biāo)被添加一定的隨機噪聲；
2. 調(diào)用前面實現(xiàn)的多項式擬合函數(shù)求解多項式系數(shù)；
3. 用求解出的多項式系數(shù)重新計算每個點的y坐標(biāo)值；
4. 可視化原始點集和擬合后的點集。

其中用于計算多項式值的函數(shù)PolyValue()實現(xiàn)如下：

float PolyValue(const cv::Mat &coeff, const int order, const float x) {
  float v = 0;
  for (int i = 0; i < order; ++i) {
    v += coeff.at<double>(i, 0) * std::pow(x, i);
  }
  return v;
}

可視化的結(jié)果如下圖所示，其中紅色點為原始點，綠色點為擬合后的點。

本文詳細(xì)介紹了基于最小二乘法的多項式擬合的原理和代碼實現(xiàn)，希望對各位讀者有所幫助。如果對本文內(nèi)容有疑問，歡迎評論區(qū)留言與我交流。