aya_xiyi 发表于 2019-1-4 18:43:44

【原创源码】2019最新全网亲测可用的克里金插值的C++版本(完整代码,亲测可用)

下面就直接上代码吧,这个是目前自己写的一个克里金插值程序的C++版本,目前在VS2017开发环境的测试下测试可用,就分享给大家吧。至此分享的资料源码如下:

第一个是kringing.h的头文件
#pragma once
#include <vector>

using namespace std;

//二维点,z代表属性值(可以是时间或者其他)
typedef struct _tagKrigingSeedPointInfo
{
        double          x,y,z;
        _tagKrigingSeedPointInfo(double a,double b,double c){x = a;y = b;z = c;};
        _tagKrigingSeedPointInfo();
} SKriSeedPnt;

class CKriging
{
private:
        //二维点
        vector<SKriSeedPnt>      m_vSeedPnt;
        double                   *m_A;

private:
        //求矩阵逆
        voidreverse(double * A,int n);

        //变差函数
        double VarFunction(double h);

        //计算两点之间的距离(用于二维点)
        double Distance(SKriSeedPnt p1,SKriSeedPnt p2);



public:
        //插值
        double Insert(double x,double y);

public:

        CKriging(vector<SKriSeedPnt> &vSeedpnt);
        ~CKriging(void);
};



接下来的这个是kringing.cpp核心文件
#include <math.h>
#include <memory.h>
#include"stdafx.h"

//克里金插值
CKriging::CKriging(vector<SKriSeedPnt> &vSeedpnt)
{

        m_vSeedPnt = vSeedpnt;
        int pntnum = m_vSeedPnt.size();
        m_A = new double [(pntnum + 1) * (pntnum + 1)];

        for (int i = 0; i < pntnum + 1; i++)
        {
                for (int j = i; j < pntnum + 1; j++)
                {
                        if (i == j)
                        {
                                m_A = 0.0;
                        }

                        else if (i == pntnum || j == pntnum)
                        {
                                m_A = 1;
                        }
                        else
                        {
                                m_A = VarFunction(Distance(vSeedpnt, vSeedpnt));
                        }
                        m_A = m_A;
                }

        }

        reverse(m_A, pntnum + 1);

}

CKriging::~CKriging(void)
{
        delete m_A;
}

double CKriging::VarFunction(double h)
{
        //if(fabs(h) < 0.000001) return 0.0;
        //double a = 80.0;
        //double c0 = 100.0;
        //double c = 500.0;

        //if (h > a)
        //{
        //        return c0+c;
        //}
        //else
        //{
        //        return c0 + c*(1.5*(h/a) - 0.5*(h*h*h/a/a/a));
        //}
        return h * 10;
}

//二维插值
double CKriging::Insert(double x, double y)
{
        double * M = new double;//非齐次项

        SKriSeedPnt tp = SKriSeedPnt(x, y, 0);

        for (unsigned int i = 0; i < m_vSeedPnt.size() + 1; i++)//求解M
        {
                if (i == m_vSeedPnt.size())
                {
                        M = 1;
                        continue;
                }
                M = VarFunction(Distance(m_vSeedPnt, tp));
        }

        double z = 0;
        for (unsigned int i = 0; i < m_vSeedPnt.size(); i++)//线性求和,求出插值z
        {
                double lamda = 0.0;
                for (int j = 0; j < m_vSeedPnt.size() + 1; j++)
                {
                        lamda += m_A * M;
                }
                z += lamda * m_vSeedPnt.z;

        }

        delete M;
        return z;
}


//二维点的距离
double CKriging::Distance(SKriSeedPnt p1, SKriSeedPnt p2)
{
        double sqtr = sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y));
        return sqtr;
}

voidCKriging::reverse(double * A, int n)
{
        int *is, *js, i, j, k, l, u, v;
        double d, p;
        double * a;
        a = new double;
        is = new int;
        js = new int;
        for (k = 0; k < n*n; k++)
                a = A;
        for (k = 0; k <= n - 1; k++)
        {
                d = 0.0;
                for (i = k; i <= n - 1; i++)
                        for (j = k; j <= n - 1; j++)
                        {
                                l = i*n + j; p = fabs(a);
                                if (p > d) { d = p; is = i; js = j; }
                        }
                        if (d + 1.0 == 1.0)
                        {
                                delete[]is; delete[]js; delete[]a;
                                //printf("err**not inv\n");
                                return;
                        }
                        if (is != k)
                                for (j = 0; j <= n - 1; j++)
                                {
                                        u = k*n + j; v = is * n + j;
                                        p = a; a = a; a = p;
                                }
                                if (js != k)
                                        for (i = 0; i <= n - 1; i++)
                                        {
                                                u = i*n + k; v = i*n + js;
                                                p = a; a = a; a = p;
                                        }
                                        l = k*n + k;
                                        a = 1.0 / a;
                                        for (j = 0; j <= n - 1; j++)
                                                if (j != k)
                                                {
                                                        u = k*n + j; a = a * a;
                                                }
                                                for (i = 0; i <= n - 1; i++)
                                                        if (i != k)
                                                                for (j = 0; j <= n - 1; j++)
                                                                        if (j != k)
                                                                        {
                                                                                u = i*n + j;
                                                                                a = a - a * a;
                                                                        }
                                                                        for (i = 0; i <= n - 1; i++)
                                                                                if (i != k)
                                                                                {
                                                                                        u = i*n + k; a = -a * a;
                                                                                }
        }
        for (k = n - 1; k >= 0; k--)
        {
                if (js != k)
                        for (j = 0; j <= n - 1; j++)
                        {
                                u = k*n + j; v = js * n + j;
                                p = a; a = a; a = p;
                        }
                        if (is != k)
                                for (i = 0; i <= n - 1; i++)
                                {
                                        u = i*n + k; v = i*n + is;
                                        p = a; a = a; a = p;
                                }
        }
        memcpy(A, a, n*n * sizeof(double));
        delete[]is; delete[]js; delete[]a;
}


基本上直接新建一个项目,把上面代码复制进去,就可以运行了。
以上就是C++版本的克里金插值,如果觉得不错,欢迎评分支持一下!

页: [1]
查看完整版本: 【原创源码】2019最新全网亲测可用的克里金插值的C++版本(完整代码,亲测可用)