OpenCV 拟合抛物线精简版

直接上代码，原理自行寻找，网上多得很

#include <opencv2/opencv.hpp> // opencv
#include <stdio.h> // 内存操作
#include <stdlib.h>
#include <iostream>

using namespace cv;

// 二维离散点拟合抛物线的三种写法，三种表达式均为 y = a·X^2 + b·X + c
// 方法一：利用 CvMat
void fitParabola_1(const std::vector<cv::Point2f> &vecPoints, float &a, float &b, float &c) {

    // 初始化
    CvMat *matA = cvCreateMat(vecPoints.size(), 3, CV_32FC1);
    CvMat *matB = cvCreateMat(vecPoints.size(), 1, CV_32FC1);
    CvMat *matC = cvCreateMat(3, 1, CV_32FC1);

    // 构造 A
    float *A = (float *)malloc(sizeof(float) * vecPoints.size() * 3);
    for (int i = 0; i < vecPoints.size(); ++i) {
        A[3 * i + 0] = vecPoints[i].x * vecPoints[i].x;
        A[3 * i + 1] = vecPoints[i].x;
        A[3 * i + 2] = 1;
    }
    // 构造 B
    float *B = (float *)malloc(sizeof(float) * vecPoints.size());
    for (int i = 0; i < vecPoints.size(); ++i) {
        B[i] = vecPoints[i].y;
    }
    // 赋值
    cvSetData(matA, A, CV_AUTOSTEP);
    cvSetData(matB, B, CV_AUTOSTEP);
    // 求抛物线参数
    cvSolve(matA, matB, matC, CV_LU);
    // 返回值
    a = matC->data.fl[0];
    b = matC->data.fl[1];
    c = matC->data.fl[2];
    // 释放空间
    if (!A) free(A);
    if (!B) free(B);
}

// 方法二：只能针对三个点进行拟合，多个点要报错，尚且不知道为什么
void fitParabola_2(const std::vector<cv::Point2f> &vecPoints, float &a, float &b, float &c) {

    // 目前只能用三个点进行拟合，后期修复
    if (vecPoints.size() != 3) {
        return;
    }

    // 初始化 Mat
    cv::Mat matA(vecPoints.size(), 3, CV_32FC1);
    cv::Mat matB(vecPoints.size(), 1, CV_32FC1);
    cv::Mat matC(3, 1, CV_32FC1);

    // 构造元素
    for (int i = 0; i < vecPoints.size(); ++i) {
        matA.at<float>(i, 0) = vecPoints[i].x * vecPoints[i].x;
        matA.at<float>(i, 1) = vecPoints[i].x;
        matA.at<float>(i, 2) = 1;
    }
    for (int i = 0; i < vecPoints.size(); ++i) {
        matB.at<float>(i, 0) = vecPoints[i].y;
    }
    cv::solve(matA, matB, matC, CV_LU);

    // 返回值
    a = matC.at<float>(0, 0);
    b = matC.at<float>(0, 1);
    c = matC.at<float>(0, 2);
}

// 方法三
// 参考博客：OpenCV最小二乘法：cv::solve函数
// https://blog.csdn.net/qq_20797273/article/details/83930101
void fitParabola_3(const std::vector<cv::Point2f> &vecPoints, float &a, float &b, float &c) {
    // 代表拟合的是二次方程，取 3 则为三次曲线样条
    int pow_number = 2;
    // 离散点个数
    int N = vecPoints.size();
    // 构造 X
    cv::Mat X = cv::Mat::zeros(pow_number + 1, pow_number + 1, CV_64FC1);
    for (int i = 0; i < pow_number + 1; i++) {
        for (int j = 0; j < pow_number + 1; j++) {
            for (int k = 0; k < N; k++) {
                X.at<double>(i, j) = X.at<double>(i, j) + std::pow(vecPoints[k].x, i + j);
            }
        }
    }
    // 构造 Y
    cv::Mat Y = cv::Mat::zeros(pow_number + 1, 1, CV_64FC1);
    for (int i = 0; i < pow_number + 1; i++) {
        for (int k = 0; k < N; k++) {
            Y.at<double>(i, 0) = Y.at<double>(i, 0) + std::pow(vecPoints[k].x, i) * vecPoints[k].y;
        }
    }

    // 求解 A
    cv::Mat A = cv::Mat::zeros(pow_number + 1, 1, CV_64FC1);
    cv::solve(X, Y, A, CV_LU/*DECOMP_LU*/);

    // 输出结果，如果是三次表达式，则继续倒着往回取值
    a = A.at<double>(2, 0);
    b = A.at<double>(1, 0);
    c = A.at<double>(0, 0);
}


int main() {

    // 随机新建离散点
    std::vector<cv::Point2f> vecPoints;
    vecPoints.emplace_back(1, 1);
    vecPoints.emplace_back(7, 6);
    vecPoints.emplace_back(2, 8);
    vecPoints.emplace_back(10, 36);

    // 分别测试拟合函数
    {
        float a, b, c;
        fitParabola_1(vecPoints, a, b, c);
        cout << a << " " << b << " " << c << endl;
        int x = 2;
        std::cout << "fitParabola_1 - y = " << a *x *x + b *x + c << std::endl;
    }

    {
        // 注意：方法二只能针对三个点，点多要崩溃，正在解决中
        //float a, b, c;
        //fitParabola_2(vecPoints, a, b, c);
        //cout << a << " " << b << " " << c << endl;
        //int x = 2;
        //std::cout << "fitParabola_2 - y = " << a * x * x + b * x + c << std::endl;
    }

    {
        float a, b, c;
        fitParabola_3(vecPoints, a, b, c);
        cout << a << " " << b << " " << c << endl;
        int x = 2;
        std::cout << "fitParabola_3 - y = " << a *x *x + b *x + c << std::endl;
    }

    return 0;
}

运行结果