希望呢感有哪位大人帮我解释一下这个段代码~
template<class T, class ForwardIterator>
class TKriging : public TInterpolater<ForwardIterator>
{
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<m_nSize; j++) {
for(int i=0; i<m_nSize; 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<double> vecB(m_nSize);
for(int i=0; i<m_nSize; 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<m_nSize-1; i++) {
double inputz = (*(first+i)).z;
z += vecB[i] * inputz;
}
if(z < 0)
z = 0;
return z;
}
private:
TMatrix<T> m_matA;
vector<int> m_Permutation;
int m_nSize;
double m_dSemivariance;
};
typedef TKriging<double, Point3D*> Kriging;