底层算法引用自 https://www.iteye.com/blog/baitai-1037388

我做一层封装和扩展。

#pragma once

/*
    给定n个point,拟合一条曲线,
    并可以获取这个曲线上任意点
    的坐标
*/
namespace GZH
{

void gauss_solve(int n, double A[], double x[], double b[]);
void polyfit(int n, double x[], double y[], int poly_n, double a[]);

// Order 阶 
// 推荐 小于五个点,3阶,多于五个点,5阶
template<class T, int Order>
class PolyFit
{
public:
    using ContainerType = T;
    using PointType = typename T::value_type;

public:
    PolyFit& append(const PointType& point) 
    {
        isCalCoef = true;
        points_.push_back(point);
        return *this;
    }

    PolyFit& append(const ContainerType& points) 
    {
        isCalCoef = true;
        for(const auto &it : points)
            points_.push_back(it);
        return *this;
    }

    PolyFit& operator<<(const PointType& point) 
    {
        return append(point);
    }

    PolyFit& operator<<(const ContainerType& points) 
    {
        return append(points);
    }

    double y(double x) 
    {
        // 计算系数
        calCoef();

        double sum = 0.0;
        for (int i = 0; i < coefficient_.size(); i++)
            sum += coefficient_[i] * pow(x, i);
        return sum;
    }

    PointType point(double x) 
    {
        return PointType(x, y(x));
    }

    ContainerType range(double beginX, double endX, double step) 
    {
        ContainerType points;
        for (double x = beginX; x < endX; x += step)
            points.push_back(point(x));
        points.push_back(point(endX));
        return points;
    }

private:
    // 计算系数
    void calCoef() 
    {
        if (!isCalCoef)
            return;

        if (points_.size() == 0)
            return;

        double coef[Order + 1] = { 0 };
        double* x = new double[points_.size()];
        double* y = new double[points_.size()];
        memset(x, 0, points_.size());
        memset(y, 0, points_.size());

        for (int i = 0; i < points_.size(); i++)
        {
            x[i] = points_[i].x();
            y[i] = points_[i].y();
        }

        auto n1 = x[0];
        auto n2 = x[1];
        auto n3 = x[2];

        polyfit(points_.size(), x, y, Order, coef);

        for (int i = 0; i < Order + 1; i++)
        {
            coefficient_.push_back(coef[i]);
        }

        isCalCoef = false;
    }

private:
    // 系数
    std::vector<double> coefficient_;
    // 原始点
    ContainerType points_;
    // 是否需要计算系数
    bool isCalCoef = true;
};

/*==================polyfit(n,x,y,poly_n,a)===================*/
/*=======拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n========*/
/*=====n是数据个数 xy是数据值 poly_n是多项式的项数======*/
/*===返回a0,a1,a2,……a[poly_n],系数比项数多一(常数项)=====*/
void polyfit(int n, double x[], double y[], int poly_n, double a[])
{
    int i, j;
    double* tempx, * tempy, * sumxx, * sumxy, * ata;
    tempx = (double*)calloc(n, sizeof(double));
    sumxx = (double*)calloc(poly_n * 2 + 1, sizeof(double));
    tempy = (double*)calloc(n, sizeof(double));
    sumxy = (double*)calloc(poly_n + 1, sizeof(double));
    ata = (double*)calloc((poly_n + 1) * (poly_n + 1), sizeof(double));
    for (i = 0; i < n; i++)
    {
        tempx[i] = 1;
        tempy[i] = y[i];
    }
    for (i = 0; i < 2 * poly_n + 1; i++)
        for (sumxx[i] = 0, j = 0; j < n; j++)
        {
            sumxx[i] += tempx[j];
            tempx[j] *= x[j];
        }
    for (i = 0; i < poly_n + 1; i++)
        for (sumxy[i] = 0, j = 0; j < n; j++)
        {
            sumxy[i] += tempy[j];
            tempy[j] *= x[j];
        }
    for (i = 0; i < poly_n + 1; i++)
        for (j = 0; j < poly_n + 1; j++)
            ata[i * (poly_n + 1) + j] = sumxx[i + j];
    gauss_solve(poly_n + 1, ata, a, sumxy);

    free(tempx);
    free(sumxx);
    free(tempy);
    free(sumxy);
    free(ata);
}

void gauss_solve(int n, double A[], double x[], double b[])
{
    int i, j, k, r;
    double max;
    for (k = 0; k < n - 1; k++)
    {
        max = fabs(A[k * n + k]); /*find maxmum*/
        r = k;
        for (i = k + 1; i < n - 1; i++)
            if (max < fabs(A[i * n + i]))
            {
                max = fabs(A[i * n + i]);
                r = i;
            }
        if (r != k)
            for (i = 0; i < n; i++)         /*change array:A[k]&A[r] */
            {
                max = A[k * n + i];
                A[k * n + i] = A[r * n + i];
                A[r * n + i] = max;
            }
        max = b[k];                    /*change array:b[k]&b[r]     */
        b[k] = b[r];
        b[r] = max;
        for (i = k + 1; i < n; i++)
        {
            for (j = k + 1; j < n; j++)
                A[i * n + j] -= A[i * n + k] * A[k * n + j] / A[k * n + k];
            b[i] -= A[i * n + k] * b[k] / A[k * n + k];
        }
    }

    for (i = n - 1; i >= 0; x[i] /= A[i * n + i], i--)
    {
        for (j = i + 1, x[i] = b[i]; j < n; j++)
        {
            x[i] -= A[i * n + j] * x[j];
        }
    }
}

} // ! namespace GZH

第一个模板参数是容器类型,可以为T 为容器必须有push_back函数,T::value_type必须有x(), y() 函数,设计为模板是因为想同时兼容Qt和STL容器, std::vector<MyPoint>等等都可以。

使用方式:

QList<QPointF> points = {
    { 6, 70 }, { 7, 110 }, { 9, 160 }, { 10, 140 },
    { 11, 130 }, { 12, 105 }, { 13, 95 },{ 14, 80 },
    { 15, 80 }, { 16, 72 }, { 17, 66 }
};
GZH::PolyFit<QList<QPointF>, 5> polyfit;
polyfit << points;
auto rangePoints = polyfit.range(6, 17, 0.1);

5阶多项式拟合,range 返回 x从6 到17 的所有坐标,间隔是0.1,例如,(6, xx), (6.1, xxx), (6.2, xxx) … (16.9, xxx), (17, xxx)

0

Publication author

offline 2年

lengyuewuxin

0
Comments: 0Publics: 25Registration: 15-07-2018
Authorization
*
*
Registration
*
*
*
Password generation