opencv 曲線擬合

最小二乘法多項式曲線擬合原理與實現 https://blog.csdn.net/jairuschan/article/details/7517773/
算法+OpenCV】基於opencv的直線和曲線擬合與繪制（最小二乘法） https://www.cnblogs.com/fengliu-/p/8031406.html

基於opencv c++代碼如下：

#include <iostream>
#include <opencv.hpp>
#include<opencv2/opencv.hpp>

using namespace std;
using namespace cv;

void FitPolynomialCurve(const std::vector<cv::Point>& points, int n, cv::Mat& A){
    //最小二乘法多項式曲線擬合原理與實現 https://blog.csdn.net/jairuschan/article/details/7517773/
    //https://www.cnblogs.com/fengliu-/p/8031406.html
    int N = points.size();
    cv::Mat X = cv::Mat::zeros(n + 1, n + 1, CV_64FC1);
    for (int i = 0; i < n + 1; i++){
        for (int j = 0; j < n + 1; j++){
            for (int k = 0; k < N; k++){
                X.at<double>(i, j) = X.at<double>(i, j) +
                        std::pow(points[k].x, i + j);
            }
        }
    }
    cv::Mat Y = cv::Mat::zeros(n + 1, 1, CV_64FC1);
    for (int i = 0; i < n + 1; i++){
        for (int k = 0; k < N; k++){
            Y.at<double>(i, 0) = Y.at<double>(i, 0) +
                    std::pow(points[k].x, i) * points[k].y;
        }
    }
    A = cv::Mat::zeros(n + 1, 1, CV_64FC1);
    cv::solve(X, Y, A, cv::DECOMP_LU);
}


int main(int argc, char **argv)
{
    string path = "/data_1/everyday/1224/2.jpeg";
    Mat img = imread(path);
    Mat img_gray,img_bi;
    cvtColor(img,img_gray,CV_BGR2GRAY);
    threshold(img_gray,img_bi,80,255,THRESH_BINARY_INV);

    vector<vector<Point> > contours;
    vector<Vec4i> hierarchy;
    findContours( img_bi, contours, hierarchy,  CV_RETR_EXTERNAL, CV_CHAIN_APPROX_SIMPLE , Point(0, 0) );
    std::cout<<contours[0].size()<<std::endl;

    cv::Mat img_draw = cv::Mat(img.rows,img.cols,CV_8UC3,Scalar(0,0,255));
    drawContours(img_draw,contours,-1,Scalar(255,255,255));


    int n = 3;
    cv::Mat A;
    FitPolynomialCurve(contours[0], n, A);
    std::vector<cv::Point> points_fitted;
    for (int x = 0; x < 800; x++)
    {
        double y = A.at<double>(0, 0) + A.at<double>(1, 0) * x +
                A.at<double>(2, 0)*std::pow(x, 2) + A.at<double>(3, 0)*std::pow(x, 3);
        points_fitted.push_back(cv::Point(x, y));
    }

    cv::polylines(img_draw, points_fitted, false, cv::Scalar(0, 0, 0), 1, 8, 0);

    imshow("img_src",img);
    imshow("img_draw",img_draw);
    imshow("img_bi",img_bi);
    waitKey(0);


    return 0;
}

效果圖如下：

但是我后面又整了個S形狀的圖像，找不到能夠很好擬合的函數階數。

#include <iostream>
#include <opencv.hpp>
#include<opencv2/opencv.hpp>

using namespace std;
using namespace cv;

void FitPolynomialCurve(const std::vector<cv::Point>& points, int n, cv::Mat& A){
    //最小二乘法多項式曲線擬合原理與實現 https://blog.csdn.net/jairuschan/article/details/7517773/
    //https://www.cnblogs.com/fengliu-/p/8031406.html
    int N = points.size();
    cv::Mat X = cv::Mat::zeros(n + 1, n + 1, CV_64FC1);
    for (int i = 0; i < n + 1; i++){
        for (int j = 0; j < n + 1; j++){
            for (int k = 0; k < N; k++){
                X.at<double>(i, j) = X.at<double>(i, j) +
                        std::pow(points[k].x, i + j);
            }
        }
    }
    cv::Mat Y = cv::Mat::zeros(n + 1, 1, CV_64FC1);
    for (int i = 0; i < n + 1; i++){
        for (int k = 0; k < N; k++){
            Y.at<double>(i, 0) = Y.at<double>(i, 0) +
                    std::pow(points[k].x, i) * points[k].y;
        }
    }
    A = cv::Mat::zeros(n + 1, 1, CV_64FC1);
    cv::solve(X, Y, A, cv::DECOMP_LU);
}

int main(int argc, char **argv)
{
    string path = "/data_1/everyday/1224/3.jpeg";
    Mat img = imread(path);
    Mat img_gray,img_bi;
    cvtColor(img,img_gray,CV_BGR2GRAY);
    threshold(img_gray,img_bi,80,255,THRESH_BINARY_INV);

    vector<vector<Point> > contours;
    vector<Vec4i> hierarchy;
    findContours( img_bi, contours, hierarchy,  CV_RETR_EXTERNAL, CV_CHAIN_APPROX_SIMPLE , Point(0, 0) );
    std::cout<<contours[0].size()<<std::endl;

    cv::Mat img_draw = cv::Mat(img.rows,img.cols,CV_8UC3,Scalar(0,0,255));
    drawContours(img_draw,contours,-1,Scalar(255,255,255));


    int n = 9;
    cv::Mat A;
    FitPolynomialCurve(contours[0], n, A);
    std::vector<cv::Point> points_fitted;
    for (int x = 0; x < 800; x++)
    {
        double y = A.at<double>(0, 0) + A.at<double>(1, 0) * x +
                A.at<double>(2, 0)*std::pow(x, 2) + A.at<double>(3, 0)*std::pow(x, 3) + A.at<double>(4, 0)*std::pow(x, 4) + A.at<double>(5, 0)*std::pow(x, 5)
                + A.at<double>(6, 0)*std::pow(x, 6) + A.at<double>(7, 0)*std::pow(x, 7) + A.at<double>(8, 0)*std::pow(x, 8) + A.at<double>(9, 0)*std::pow(x, 9);
                //+ A.at<double>(10, 0)*std::pow(x, 10) + A.at<double>(11, 0)*std::pow(x, 11) + A.at<double>(12, 0)*std::pow(x, 12);
        points_fitted.push_back(cv::Point(x, y));
    }

    cv::polylines(img_draw, points_fitted, false, cv::Scalar(0, 0, 0), 1, 8, 0);


    imshow("img_src",img);
    imshow("img_draw",img_draw);
    imshow("img_bi",img_bi);
    waitKey(0);


    return 0;
}

突然想明白，這個S形狀曲線一個x對應好幾個y，不行。需要一個x唯一對應一個y的曲線才能擬合。然后又順手畫了一個，果真可以擬合。

當然代碼每次根據不同的階數寫好多A.at (6, 0)*std::pow(x, 6)，可以用如下函數自動根據x得到y：

double CurveY(double x, cv::Mat& A){
    double y = 0.0;
    double *a = A.ptr<double>();
    for (int i = 0; i < A.rows; i++){
        y += a[i] * pow(x, i);
    }
    return y;
}