[kriging](一)网上下载的kriging克里金的C++程序的初步调试

笔者在网上下载了一份克里金的C++程序，水平有限，正在逐步地调试中。

初步

克里金法现在在许多软件都已经有集成了，据笔者所知：

• arcgis ： 看过有的arcgis培训视频里面简略介绍了里面的插值方法，有克里金。
• surfer ：笔者用过这个，效果很好。
• gocad ：听说有
• landmark &jason ：研究生时候用的专业的石油物探软件，里面有克里金法，万恶的开始。
• pykrig : python实现的克里金
• gstl : 一个地质统计学的C++的库，里面有克里金方法。
• gslib : 一个Fortran的地质统计学的库，里面有克里金方法。
• SGeMs ：这个比较出名，应该是斯坦福大学的开源的软件，但是年代很久远了。
• dace : 貌似是各matlab的工具箱，里面有克里金方法。
  以上的库基本可以满足日常的需求了，但是，笔者目前需要从事一些深入地研究，将克里金改用在其他领域，就需要从代码层面上了解了。

网上下的一段代码

从网上下载了一份代码，同样年代久远，使用MFC界面，内核里面插值方法有反距离权及克里金法。感觉上大多数想学克里金的人都会下载到这样一份不知所云的代码，每个文件都不知道在干什么。

如果想要完整编译该代码，VS2017必须完整地添加古老的MFC功能。在visual studio installer中添加如下功能，但是笔者并未这样做。

经过简略的分析
代码中，Kriking.h和InverseDist.h既包含了头又包含了实现的方法，同时继承自Interpolater.h


反距离权插值的思想为克里金的基础，所以先实现基础的反距离权
新建工程，然后在工程中把下载下来的代码中的和反距离权插值法的内容添加进项目工程

主函数main.cpp内容如下：

#include "InputReader.h"
#include "InverseDist.h"
#include <Eigen/Dense>
//#include "Kriging.h"

using namespace Eigen;
int main()
{
	InputReader m_ir;
	int nDiameter = 200;
	vector<double> vecZs;
	Interpolater* pInterpolater = NULL;
	m_ir.Read("./testdata.txt"); 

	vector<Point3D>& input = const_cast<vector<Point3D>&>(m_ir.Get3DPoints());
	
	for (int i = 0; i < input.size(); i++)
	{
		cout <<i<<" "<< input[i].x <<" " << input[i].y << " " << input[i].z << endl;
	}

	pInterpolater = new InverseDist(200, 4);
	//pInterpolater = new Kriging(input.begin(), input.end(), 4);
	int nRadius = nDiameter / 2;
	for (int i = 0; i < nDiameter; i++) {
		for (int j = 0; j < nDiameter; j++) {
			double z = pInterpolater->GetInterpolatedZ(j - nRadius, i - nRadius, &input[0], &input[input.size()-1]);
			vecZs.push_back(z);
		}
	}
	delete pInterpolater;

	/*转换到Eigen矩阵方便在cmd中查看输出的插值后的z**/
	MatrixXd m_matrix = Map<Matrix<double, 1,4000,RowMajor>>(vecZs.data());
	cout << "After Inverse Dist Interpolate:\n"<<m_matrix.transpose() << endl;

	return 0;
}

运行结果

数据txt文件读取结果

原始数据可视化
反距离权插值后的z值，该程序初始化的距离为radius=200，距离的p次方位反比权，p初始值为8，这里使用了下载程序中的示例数据。

按道理，应该还需要实现插值后的每个z值点对应的x和y，就能够查看反距离插值后的三维数据，这个留待以后完成。