请选择 进入手机版 | 继续访问电脑版
查看: 169|回复: 0

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

[复制链接]

3

主题

8

帖子

24

积分

新手上路

Rank: 1

积分
24
发表于 2019-1-4 18:43:44 | 显示全部楼层 |阅读模式
下面就直接上代码吧,这个是目前自己写的一个克里金插值程序的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++版本的克里金插值,如果觉得不错,欢迎评分支持一下!





上一篇:【原创资料】2019最新Unity3D实战打造多款爆款三维游戏开发全套视频教程(3D游戏)
下一篇:【Windows软件】2019最新非常好用的抖音短视频在线无水印播放和黑科技下载神器
发帖求助前要善用【论坛搜索】功能,那里可能会有你要找的答案; 如果你在论坛求助问题,并且已经从坛友或者管理的回复中解决了问题,请把帖子分类或者标题加上【已解决】。
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

微信扫一扫

我爱科技论坛(www.52tech.tech)旨在打造全网最大的免费资源共享平台。目前论坛包括考研资料、编程学习、黑科技/科学上网、开源软件等资源模块,竭力服务于正在学习道路上的每一个人。

QQ|Archiver|我爱科技论坛 快乐学习交流

(请勿发布违反中华人民共和国法律法规的言论,会员观点不代表我爱科技论坛的官方立场)

Powered by 我爱科技论坛 X3.4© 2001-2020.

返回顶部 返回列表