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

cv::Vec4d cubicFit(const std::vector<double>& x, const std::vector<double>& y) {  
    // 使用最小二乘法进行三次多项式拟合  
    int n = x.size();  
    cv::Mat A(n, 4, CV_64F);  
    cv::Mat b(n, 1, CV_64F);  
    
    for (int i = 0; i < n; ++i) {  
        A.at<double>(i, 0) = x[i] * x[i] * x[i]; // x^3  
        A.at<double>(i, 1) = x[i] * x[i];       // x^2  
        A.at<double>(i, 2) = x[i];               // x  
        A.at<double>(i, 3) = 1.0;                // constant term  
        b.at<double>(i, 0) = y[i];  
    }  

    cv::Mat coeffs;  
    cv::solve(A, b, coeffs); // 计算回归系数  
    return cv::Vec4d(coeffs.at<double>(0), coeffs.at<double>(1), coeffs.at<double>(2), coeffs.at<double>(3));  
}  

int main() {  
    // 示例输入数据  
    std::vector<double> x = {0, 1, 2, 3, 4, 5};  
    std::vector<double> y = {2, 3, 5, 7, 11, 13}; // 示例数据点  

    // 拟合多项式  
    cv::Vec4d polyParams = cubicFit(x, y);  
    
    // 创建绘图窗口  
    int imgWidth = 400;  
    int imgHeight = 400;  
    cv::Mat img = cv::Mat::zeros(imgHeight, imgWidth, CV_8UC3);  

    // 绘制原始数据点  
    for (size_t i = 0; i < x.size(); ++i) {  
        int plotX = static_cast<int>(x[i] * imgWidth / 5);  
        int plotY = imgHeight - static_cast<int>(y[i] * imgHeight / 15); // 底部对齐  
        cv::circle(img, cv::Point(plotX, plotY), 5, cv::Scalar(0, 0, 255), -1); // 红色点  
    }  

    // 绘制拟合曲线  
    for (double xi = 0; xi <= 5; xi += 0.01) {  
        double yi = polyParams[0] * xi * xi * xi + polyParams[1] * xi * xi + polyParams[2] * xi + polyParams[3];  
        int plotX = static_cast<int>(xi * imgWidth / 5);  
        int plotY = imgHeight - static_cast<int>(yi * imgHeight / 15); // 底部对齐  
        cv::circle(img, cv::Point(plotX, plotY), 1, cv::Scalar(255, 0, 0), -1); // 蓝色曲线  
    }  
    
    // 显示图像  
    cv::imshow("Fitted Curve", img);  
    cv::waitKey(0); // 等待按键按下  
    return 0;  
}