请选择 进入手机版 | 继续访问电脑版

我爱科技论坛

 找回密码
 立即注册

QQ登录

只需一步,快速开始

查看: 263|回复: 1

[原创资源] 普通克里金插值的C++版本(完整代码,亲测可用)源代码免费下载

[复制链接]

3

主题

8

帖子

24

积分

新手上路

Rank: 1

积分
24
发表于 2019-1-4 18:43:44 | 显示全部楼层 |阅读模式
克里金插值主要分为K邻域搜索,构建协方差矩阵,求解方程组等环节,一下就是普通克里金插值的代码。
下面就直接上代码吧,这个是目前自己写的一个克里金插值程序的C++版本,目前在VS2017开发环境的测试下测试可用,就分享给大家吧。至此分享的免费资源源码如下:

第一个是kringing.h的头文件
[C++] 纯文本查看 复制代码
#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:
        //求矩阵逆
        void  reverse(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核心文件
[C++] 纯文本查看 复制代码
#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[i *(pntnum + 1) + j] = 0.0;
                        }

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

        }

        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[m_vSeedPnt.size() + 1];//非齐次项

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

        for (unsigned int i = 0; i < m_vSeedPnt.size() + 1; i++)//求解M
        {
                if (i == m_vSeedPnt.size())
                {
                        M[i] = 1;
                        continue;
                }
                M[i] = VarFunction(Distance(m_vSeedPnt[i], 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[i * (m_vSeedPnt.size() + 1) + j] * M[j];
                }
                z += lamda * m_vSeedPnt[i].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;
}

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



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





上一篇:最新 Unity 5.X从入门到精通PDF(近2G)资源免费下载
下一篇:【Windows软件】2019最新非常好用的抖音短视频在线无水印播放和黑科技下载神器
发帖求助前要善用【论坛搜索】功能,那里可能会有你要找的答案; 如果你在论坛求助问题,并且已经从坛友或者管理的回复中解决了问题,请把帖子分类或者标题加上【已解决】。
回复过本主题
的还回复过:

4

主题

19

帖子

46

积分

新手上路

Rank: 1

积分
46
发表于 2020-3-28 12:23:33 | 显示全部楼层
感谢分享
发帖求助前要善用【论坛搜索】功能,那里可能会有你要找的答案; 如果你在论坛求助问题,并且已经从坛友或者管理的回复中解决了问题,请把帖子分类或者标题加上【已解决】。
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

微信扫一扫

快速回复 返回顶部 返回列表