点云数据。

首先呢我们得了解点云数据的格式：

第一排为点的数量，其下为三位点数据（X Y Z）
此为输出格式，基于规则格网的grd文件，以及txt文件
DSAA

nx ny nx is the integer number of grid lines along the X axis (columns)
ny is the integer number of grid lines along the Y axis (rows)
xlo xhi xlo is the minimum X value of the grid
xhi is the maximum X value of the grid
ylo yhi ylo is the minimum Y value of the grid
yhi is the maximum Y value of the grid
zlo zhi zlo is the minimum Z value of the grid
zhi is the maximum Z value of the grid
grid row 1
grid row 2
grid row 3
� these are the rows of Z values of the grid, organized in row order. Each row has a constant Y coordinate. Grid row 1 corresponds to ylo and the last grid row corresponds to yhi. Within each row, the Z values are arranged from xlo to xhi.

DSAA

10 10

0.0 9.0

0.0 7.0

25.00 97.19

91.03 77.21 60.55 46.67 52.73 64.05 41.19 54.99 44.30 25.00

96.04 81.10 62.38 48.74 57.50 63.27 48.67 60.81 51.78 33.63

92.10 85.05 65.09 53.01 64.44 65.64 52.53 66.54 59.29 41.33

94.04 85.63 65.56 55.32 73.18 70.88 55.35 76.27 67.20 45.78

97.19 82.00 64.21 61.97 82.99 80.34 58.55 86.28 75.02 48.75

91.36 78.73 64.05 65.60 82.58 81.37 61.16 89.09 81.36 54.87

86.31 77.58 67.71 68.50 73.37 74.84 65.35 95.55 85.92 55.76

80.88 75.56 74.35 72.47 66.93 75.49 86.39 92.10 84.41 55.00

74.77 66.02 70.29 75.16 60.56 65.56 85.07 89.81 74.53 51.69

70.00 54.19 62.27 74.51 55.95 55.42 71.21 74.63 63.14 44.99

第一步

根据输入的数据，设置格网大小。
建立符合精度要求的规则格网

Like this:

第二步

在每个各网点的附近寻找符合要求的已知点，用于高程拟合。
在这我用的是欧式距离。也可以用画圆的方法。

第三步就是最小二乘拟合啦

误差方程式：

求解什么的就不放出来了
下面直接放代码咯。

拟合部分的代码

lob为
struct pointrgb
{
double X, Y, Z;
int R, G, B;
};
的一个vector，存放三位点数据以及对应RGB值。
此为一些比较的定义：

bool compX(const pointrgb& a, const pointrgb& b)
{return a.X < b.X;
}
bool compY(const pointrgb& a, const pointrgb& b)
{return a.Y < b.Y;
}
bool compZ(const pointrgb& a, const pointrgb& b)
{return a.Z < b.Z;
}
bool compx(const Point3d& a, const Point3d& b)
{return a.x < b.x;
}

主函数代码：

insert();CString tmpstr;/*if (Txt.IsEmpty())AfxMessageBox(_T("AD"));*/sort(lob.begin(), lob.end(), compX);xlo = lob.at(0).X;xhi = lob[count-1].X;double dx = xhi - xlo;cols = ceil(dx);sort(lob.begin(), lob.end(), compY);ylo = lob.at(0).Y;yhi = lob[count-1].Y;double dy = yhi - ylo;rows = ceil(dy);//losX.resize(cols);losY.resize(rows);vector<pointrgb>  ptz;ptz.resize(size_t(rows*cols));for (int i = 0; i < rows; i++){losY.at(i).resize(cols);}for (int i = 0; i < (int)lob.size(); i++){int locateY1 = int(int(lob.at(i).Y - ylo+0.5));int locateX1 = int(int(lob.at(i).X - xlo+0.5));pointrgb YP1;YP1.X = lob.at(i).X;YP1.Y = lob.at(i).Y;YP1.Z = lob.at(i).Z;YP1.R = lob.at(i).R;YP1.G = lob.at(i).G;YP1.B = lob.at(i).B;if (locateY1 < int(rows)&&locateX1 < int(cols)){losY.at(locateY1).at(locateX1).push_back(YP1);ptz.at(size_t(locateY1) * cols + locateX1).R = YP1.R;ptz.at(size_t(locateY1) * cols + locateX1).G = YP1.G;ptz.at(size_t(locateY1) * cols + locateX1).B = YP1.B;}}//四舍五入预处理点位数据//存放计算出的高程Mat Hmatrix,Hmatrixres;Hmatrix.create(rows, cols, CV_64F);Hmatrix.copyTo(Hmatrixres);tmpstr.Format(_T("%s\n"),_T("DSAA"));Txt += tmpstr;//Txt += _T("\n");tmpstr.Format(_T("%d %d\n"), cols, rows);Txt += tmpstr;tmpstr.Format(_T("%lf %lf\n"), xlo, xhi);Txt += tmpstr;tmpstr.Format(_T("%lf %lf\n"), ylo, yhi);Txt += tmpstr;sort(lob.begin(), lob.end(), compZ);double zlo, zhi;zlo = lob[0].Z;zhi=lob[count-1].Z;//sort(ptz.begin(), ptz.end());tmpstr.Format(_T("%lf %lf\n"), lob[0].Z, lob[count - 1].Z);Txt += tmpstr;tmpstr.Format(_T("%d\n"), rows * cols);CCdata += tmpstr;//文件头for (int i = 0; i < rows; i++){for (int j = 0; j < cols; j++){int num1=0;int size = 4;vector<Point3d>  Ptvec;while (num1 <= 6){num1 = getnum(i, j, size, Ptvec);size = size + 2;}Point2d pt2dtmp;double x = xlo + j;double y = ylo + i;pt2dtmp.x = x;pt2dtmp.y = y;//double meanX = 0, meanY = 0;//for (int i = 0; i < (int)Ptvec.size(); i++)//{//    meanX += Ptvec.at(i).x;// meanY += Ptvec.at(i).y;//}//meanX /= Ptvec.size();//meanY /= Ptvec.size();//double erry = i - meanY + ylo;//double errx = j - meanX + xlo;计算选点重心与格网点的偏移量Hmatrix.at<double>(i, j) = interplation2(i,j,Ptvec);if ((i<rows / 5 || i>(rows * 4) / 5) || ((j<cols / 5 || j>(cols * 4) / 5))||size>4)Hmatrix.at<double>(i, j) = interplation(Ptvec);//拟合ptz.at(size_t(i)* cols + j).Z = interplation(Ptvec);ptz.at(size_t(i)* cols + j).X = xlo + j;ptz.at(size_t(i)* cols + j).Y = ylo + i;tmpstr.Format(_T("%lf %lf %lf %d %d %d\n"), ptz.at(size_t(i)* cols + j).X, ptz.at(size_t(i)* cols + j).Y, ptz.at(size_t(i)* cols + j).Z, ptz.at(size_t(i)* cols + j).R, ptz.at(size_t(i)* cols + j).G, ptz.at(size_t(i)* cols + j).B);CCdata += tmpstr;/*if ((size_t(i) * cols + j)){if ((fabs(errx) < 4 && fabs(erry) <= 4) || (ptz.at(size_t(i) * cols + j) > ptz.at(size_t(i) * cols + j - 1) + 5 || ptz.at(size_t(i) * cols + j) < ptz.at(size_t(i) * cols + j - 1) - 5)||i==0){ptz.at(size_t(i) * cols + j) = interplation(Ptvec);}if (i > 0){if (ptz.at(size_t(i) * cols + j) > ptz.at(size_t(i - 1) * cols + j) + 5 || ptz.at(size_t(i - 1) * cols + j) < ptz.at(size_t(i - 1) * cols + j) - 5){ptz.at(size_t(i) * cols + j) = interplation(Ptvec);}}}*///Hmatrix.at<double>(i, j) = ptz[size_t(i) * cols + j];Ptvec.clear();}}gaussTest(Hmatrix, Hmatrixres);for (int i = 0; i < rows; i++){for (int j = 0; j < cols; j++){tmpstr.Format(_T("%.4f"), Hmatrix.at<double>(i,j));Txt += tmpstr;Txt += _T(" ");if (j == cols - 1)Txt += _T("\n");}}AfxMessageBox(_T("计算完成！请选择文件保存路径！"));fileout();

一些功能函数，包括重心化，最小二乘解算，以及测试代码：

double CDem::getnum(int row, int col, int size, vector<Point3d>& Ptvec)
{int countpt = 0;for (int i = row - size; i < row + size + 1; i++){for (int j = col - size; j < col + size + 1; j++){if (i >= 0 && j >= 0 && i < rows && j < cols){countpt += losY.at(i).at(j).size();for (int k = 0; k < (int)losY.at(i).at(j).size(); k++){Point3d tmp;tmp.x = losY.at(i).at(j).at(k).X;tmp.y = losY.at(i).at(j).at(k).Y;tmp.z = losY.at(i).at(j).at(k).Z;Ptvec.push_back(tmp);}}}}return countpt;
}
void CDem::fileout()
{CFileDialog dlgFile(TRUE, _T("grd"), NULL, OFN_ALLOWMULTISELECT | OFN_EXPLORER, _T("(文本文件)|*.grd"));if (dlgFile.DoModal() == IDCANCEL)return;CString strFileName = dlgFile.GetPathName();setlocale(LC_ALL, "");CStdioFile sf;if (!sf.Open(strFileName, CFile::modeCreate | CFile::modeWrite))return;sf.WriteString(Txt);sf.Close();
}
double CDem::interplation(int row, int col, vector<Point3d>ptvec)
{//1 重心化Mat M;Mat P = Mat::zeros(Size(ptvec.size(), ptvec.size()), CV_64F);Mat Z;M.create(ptvec.size(), 6, CV_64F);Z.create(ptvec.size(), 1, CV_64F);point3d* pt3d = new point3d[ptvec.size()];double* d = new double[ptvec.size()];for (int i = 0; i < (int)ptvec.size(); i++){pt3d[i].x = (double)(ptvec.at(i).x) - col-xlo;pt3d[i].y = (double)(ptvec.at(i).y) - row-ylo;d[i] = double(pt3d[i].x * pt3d[i].x + pt3d[i].y * pt3d[i].y);d[i] = 1 / d[i];M.at<double>(i, 0) = pt3d[i].x * pt3d[i].x;M.at<double>(i, 1) = pt3d[i].x * pt3d[i].y;M.at<double>(i, 2) = pt3d[i].y * pt3d[i].y;M.at<double>(i, 3) = pt3d[i].x;M.at<double>(i, 4) = pt3d[i].y;M.at<double>(i, 5) = 1;Z.at<double>(i, 0) = ptvec.at(i).z;P.at<double>(i, i) = d[i];}Mat M_ = M.t();Mat mpm;mpm.create(6, 6, CV_64F);mpm = (M_ * P) * M;Mat mpm1;mpm1.create(6, 6, CV_64F);mpm1 = mpm.inv();Mat X;X.create(6, 1, CV_64F);X = mpm1 * M_ * P * Z;delete[] pt3d;delete[] d;return X.at<double>(5, 0);
}

整个程序的源代码较繁琐，在这就不全部展示。
本人博客里面提供完整下载。链接: link.

数字高程(移动曲面）拟合（C++）相关推荐

DEM数据知识介绍-数字高程模型
1. DEM数据源特征 2. DEM数据采样理论 3. DEM数据采样方法 4. DEM数据采样质量控制 5. DEM数据共享与利 ...
数字高程模型DEM和构建学习1
DEM DEM,一般指数字高程模型.Digital Elevation Model,简称DEM. 是通过有限的地形高程数据实现对地面地形的数字化模拟,它是用一组有序数值阵列形式表示地面高程的一种实体地 ...
[GIS原理] 9 数字地形分析DTA、数字地形模型DTM、数字高程模型DEM、数字地表模型DSM、不规则三角网TIN
在知识传播途中,向涉及到的相关著作权人谨致谢意! 文章目录 1 数字地形分析(DTA) 1.1 数字地形模型(DTM) 1.1.1 DSM与DEM对比 1.2 数字地形分析研究与应用进展 1.2.1 ...
「干货」12.5米数字高程DEM专题图制作教程
[干货]12.5米数字高程DEM专题图制作教程概述数字高程模型(Digital Elevation Model),简称DEM,是表达地面高程起伏形态的实体地面模型. 可用于绘制等高线.坡度图.坡向 ...
【QGIS入门实战精品教程】4.8：QGIS如何下载SRTM数字高程模型DEM？
本文讲解QGIS中下载SRTM数字高程模型DEM,以黑龙江省塔河县为例. 图幅效果: 最终效果: 文章目录 1. 下载安装STRM Download插件 2. 加载矢量数据,读取范围 3. 下载STR ...
高精度数字高程数据1m的dem
最近一直在做一个三维场景的项目,老板一直让我下载高精度的数据,从而构建三维场景,一开始在各种网站上下载如中国科学院镜像站点下载,德国航天局的DLR数字高程数据,美国航天局NASA网站上等一些地方下载, ...
全球数字高程数据：ASTER GDEM
简介 ASTER GDEM,即先进星载热发射和反射辐射仪全球数字高程模型,与SRTM一样为数字高程DEM,其全球空间分辨率为30米.该数据是根据 NASA的新一代对地观测卫星Terra的详尽观测结果制 ...
12米数字高程DEM现已上线！附DEM专题图制作教程
概述: 数字高程模型(Digital Elevation Model),简称DEM,是表达地面高程起伏形态的实体地面模型.可用于绘制等高线.坡度图.坡向图.立体透视图.立体景观图,并应用于制作正射影像 ...
2022年台湾省矢量数据（点线面）及数字高程数据下载
鉴于最近台湾形势,百度高德地图开放道路和poi数据,各大媒体争相关注.因此将台湾省截至目前可获取的最新矢量数据,及srtm格式30米分辨的dem数据分享,需要的人可自行下载. 一.矢量数据(shp格式 ...

数字高程(移动曲面）拟合（C++）

基于点云数据的数字高程最小二乘拟合

点云数据。

相关原理

第一步

第二步

第三步就是最小二乘拟合啦

拟合部分的代码

数字高程(移动曲面）拟合（C++）相关推荐

最新文章

热门文章