[推荐]：kringing插值源代码

zzz8662

最优秀的网格插值方法贡献给大家！
另外：哪位老兄能把这段C代码转化为VB代码？偶看不懂C代码：（  谢谢了。*-*6 *-*6 *-*6
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 = 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 * 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 = 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 * 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;
sum = B[ll];
B[ll] = B;
if(ii != 0)
for(int j=ii; j sum -= A(i, j) * B[j];
else if(sum != 0.0)
ii = i;
B = sum;
}

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