文件名称:
Kriging 算法实现 2维和3维地图等高线
开发工具:
文件大小: 163kb
下载次数: 0
上传时间: 2006-02-23
详细说明：
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 Kr iging 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:templatedouble 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));}templatedouble GetDistance(double xpos, double ypos, const ForwardIterator start, int i){ return ::sqrt(::pow(((*(start+i)).x - xpos), 2) + ::pow(((*(start+i)).y - ypos), 2));}templateclass 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 vecB(m_nSize); for(int i=0; i 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:templatevoid 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) 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= amax) { imax = i; amax = dum; } } if(j != imax) { for(int k=0; kvoid 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=0; i--) { sum = B[i]; if(i< n) { for(int j=i+1; j 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. ...展开收缩
(系统自动生成,下载前可以参看下载内容)
下载文件列表
相关说明
- 本站资源为会员上传分享交流与学习,如有侵犯您的权益,请联系我们删除.
- 本站是交换下载平台，提供交流渠道，下载内容来自于网络，除下载问题外，其它问题请自行百度。
- 本站已设置防盗链，请勿用迅雷、QQ旋风等多线程下载软件下载资源，下载后用WinRAR最新版进行解压.
- 如果您发现内容无法下载，请稍后再次尝试；或者到消费记录里找到下载记录反馈给我们.
- 下载后发现下载的内容跟说明不相乎，请到消费记录里找到下载记录反馈给我们，经确认后退回积分.
- 如下载前有疑问，可以通过点击"提供者"的名字，查看对方的联系方式，联系对方咨询.
