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<class ForwardIterator>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<class ForwardIterator>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 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<T> 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;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<class T>void LUDecompose(TMatrix<T>& A, std::vector<int>& Permutation, int& d) throw(NumericException){ int n = A.GetHeight(); vector<T> vv(n); Permutation.resize(n); d=1; T amax; for(int i=0; i<n; i++) { amax = 0.0; for(int j=0; j<n; 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<n; j++) { for (i=0; i<j; i++) { sum = A(i, j); for(int k=0; k<i; k++) sum -= A(i, k) * A(k, j); A(i, j) = sum; } amax = 0.0; for(i=j; i<n; i++) { sum = A(i, j); for(int k=0; k<j; 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<n; 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<n; i++) A(i, j) *= dum; } }}template<class T>void LUBackSub(TMatrix<T>& A, std::vector<int>& Permutation, std::vector<T>& B) throw(NumericException){ int n = A.GetHeight(); T sum; int ii = 0; int ll; for(int i=0; i<n; i++) { ll = Permutation[i]; sum = B[ll]; B[ll] = B[i]; if(ii != 0) for(int j=ii; j<i; 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<n; 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<Point3D> input // assume this vector has KNOWN 3D points Interpolater* pInterpolater = new Kriging(input.begin(), input.end(), 4); vector<double> 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.
- 1
- 粉丝: 882
- 资源: 2万+
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
最新资源
- C语言基础-C语言编程基础之Leetcode编程题解之第30题串联所有单词的子串.zip
- C语言基础-C语言编程基础之Leetcode编程题解之第29题两数相除.zip
- C语言基础-C语言编程基础之Leetcode编程题解之第28题找出字符串中第一个匹配项的下标.zip
- 实验报告模板(1).docx
- C语言基础-C语言编程基础之Leetcode编程题解之第26题删除有序数组中的重复项.zip
- C语言基础-C语言编程基础之Leetcode编程题解之第25题K个一组翻转链表.zip
- hnu计算机系统作业-计算机系统基础课程大作业.zip
- 树莓派app.apk
- C++的基于同态加密技术的匿名电子投票系统源码.zip
- SW建模格式图.zip
- 1
- 2
- 3
- 4
- 5
- 6
前往页