加载中…
个人资料
铁科院_郭辉
铁科院_郭辉
  • 博客等级:
  • 博客积分:0
  • 博客访问:223,012
  • 关注人气:151
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
相关博文
推荐博文
谁看过这篇博文
加载中…
正文 字体大小:

kringing插值源代码

(2012-09-28 15:29:39)
标签:

杂谈

分类: 研究方法

 参考来源:http://www.dixue.com/forum.php?mod=viewthread&tid=16118

           http://motioo.blog.163.com/blog/static/117718291200954102822505/

  A year and a half year ago, I published this article to the Codeguru site and got a number of requests about the Kriging algorithm contour map. Unfortunately, my project was changed shortly after that article and later I quit the company so I couldn‘t find time to finish this Contour business. A week ago, I happened to need a contour map again so I decided to solve the Kriging algorithm. I searched the Internet for a commercial library but they all look ugly and hard to use. So, I made up my mind to make my own algorithm. The Kriging algorithm is easy to find, but this algorithm needs a Matrix and solver (LU-Decomposition). Again, I couldn‘t find suitable code for this. I tried to use GSL first but this made my code too big and was slower. Finally, I went back to "Numerical Recipe in C"—yes, that horrible-looking C code—and changed the code there to my taste.

If you read this article before, the rendering part hasn‘t been changed much. I added the Kriging algorithm and revised the codes a little bit. Following is the Kriging Algorithm:

template
double GetDistance(const ForwardIterator start, int i, int j)
{
return ::sqrt(::pow(((*(start+i)).x - (*(start+j)).x), 2) +
::pow(((*(start+i)).y - (*(start+j)).y), 2));
}

template
double GetDistance(double xpos, double ypos,
const ForwardIterator start, int i)
{
return ::sqrt(::pow(((*(start+i)).x - xpos), 2) +
::pow(((*(start+i)).y - ypos), 2));
}

template
class TKriging : public TInterpolater
{
public:
TKriging(const ForwardIterator first, const ForwardIterator last,
double dSemivariance) : m_dSemivariance(dSemivariance)
{
m_nSize = 0;
ForwardIterator start = first;
while(start != last) {
++m_nSize;
++start;
}

m_matA.SetDimension(m_nSize, m_nSize);

for(int j=0; j for(int i=0; i if(i == m_nSize-1 || j == m_nSize-1) {
m_matA(i, j) = 1;
if(i == m_nSize-1 && j == m_nSize-1)
m_matA(i, j) = 0;
continue;
}
m_matA(i, j) = ::GetDistance(first, i, j) * dSemivariance;
}
}
int nD;
LUDecompose(m_matA, m_Permutation, nD);
}
double GetInterpolatedZ(double xpos, double ypos,
ForwardIterator first,
ForwardIterator last)
throw(InterpolaterException)
{
std::vector vecB(m_nSize);
for(int i=0; i double dist = ::GetDistance(xpos, ypos, first, i);
vecB[i] = dist * m_dSemivariance;
}
vecB[m_nSize-1] = 1;

LUBackSub(m_matA, m_Permutation, vecB);

double z = 0;
for(i=0; i double inputz = (*(first+i)).z;
z += vecB[i] * inputz;
}
if(z < 0)
z = 0;
return z;
}
private:
TMatrix m_matA;
vector m_Permutation;
int m_nSize;
double m_dSemivariance;
};

typedef TKriging Kriging;

Because of the template, this doesn‘t look that clean but you can get the idea if you look at it carefully. The matrix solver is as follows:

template
void LUDecompose(TMatrix& A, std::vector&
Permutation, int& d) throw(NumericException)
{
int n = A.GetHeight();
vector vv(n);
Permutation.resize(n);

d=1;

T amax;
for(int i=0; i amax = 0.0;
for(int j=0; j if(fabs(A(i, j)) > amax)
amax = fabs(A(i, j));

if(amax < TINY_VALUE)
throw NumericException();

vv[i] = 1.0 / amax;
}

T sum, dum;
int imax;
for(int j=0; j for (i=0; i sum = A(i, j);
for(int k=0; k sum -= A(i, k) * A(k, j);
A(i, j) = sum;
}
amax = 0.0;

for(i=j; i sum = A(i, j);
for(int k=0; k sum -= A(i, k) * A(k, j);

A(i, j) = sum;
dum = vv[i] * fabs(sum);

if(dum >= amax) {
imax = i;
amax = dum;
}
}

if(j != imax) {
for(int k=0; k dum = A(imax, k);
A(imax, k) = A(j, k);
A(j, k) = dum;
}
d = -d;
vv[imax] = vv[j];
}
Permutation[j] = imax;

if(fabs(A(j, j)) < TINY_VALUE)
A(j, j) = TINY_VALUE;

if(j != n) {
dum = 1.0 / A(j, j);
for(i=j+1; i A(i, j) *= dum;
}
}
}

template
void LUBackSub(TMatrix& A, std::vector&
Permutation, std::vector& B)
throw(NumericException)
{
int n = A.GetHeight();
T sum;
int ii = 0;
int ll;
for(int i=0; i ll = Permutation[i];
sum = B[ll];
B[ll] = B[i];
if(ii != 0)
for(int j=ii; j sum -= A(i, j) * B[j];
else if(sum != 0.0)
ii = i;
B[i] = sum;
}

for(i=n-1; i>=0; i--) {
sum = B[i];
if(i< n) {
for(int j=i+1; j sum -= A(i, j) * B[j];
}
B[i] = sum / A(i, i);
}
}

By using this algorithm, making a 3D grid is easy. Let‘s assume we‘re making a 200x200 grid and we have some scattered data. Then, what we need to do is this:

vector input // assume this vector has KNOWN 3D points

Interpolater* pInterpolater = new Kriging(input.begin(),
input.end(), 4);

vector vecZs;

for(int j=0; j<200; j++) {
for(int i=0; i<200; i++) {
vecZs.push_back(pInterpolater->GetInterpolatedZ(i, j,
input.begin(),
input.end()));
}
}
// Now, vecZs has 40000 z values

delete pInterpolater;

If you have all the grid points with 3D data, you can make a bitmap file with it, or make a triangle strip to render with OpenGL. If you remember that the old contour map was produced from an InverseDistanced algorithm (you can switch to Inverse Distance in the Option menu), you‘ll find a vast improvement over it. I compared the Kriging generated contour map with some commercial programs, and they were almost identical. I hope this helps programmers who want to make a contour map.

0

阅读 评论 收藏 转载 喜欢 打印举报/Report
  • 评论加载中,请稍候...
发评论

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

      

    新浪BLOG意见反馈留言板 电话:4000520066 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有