目錄
- 概述
- 詳論
- 資料準備
- 轉換演算法
- TIN構建
- 具體實作
- 實驗結果
- 參考
??概述
在GIS（地理資訊科學）中，地形有兩種表達方式，一種是格網DEM，一種是不規則三角網TIN，一般情況下規則格網DEM用的比較多，因為可以將高程當作像素，將其存盤為圖片型別的資料（例如.tif），但是規則格網存盤的資料量大，按規則取點，并不能最大程度的保證地形特征，所以很多情況下需要將其表達為不規則三角網，也就是TIN，
??詳論
1??資料準備
下載SRTM30的DEM資料，找到美國大峽谷附近的地形，通過UTM投影，將其轉換成30米的平面坐標的DEM（.tif格式），通過Global Mapper打開，顯示的效果如下：
2??轉換演算法
格網DEM本身也可以看作是一個三角網，每個方格由兩個三角形組成，N個方格據組成了一個地形格網，所以在參考文獻一中提到了一種保留重要點法，將格網DEM中認為不重要的點去除掉，剩下的點構建成不規則三角網即可，那么怎么直到有的點重要，有的點不重要呢？參考文獻一中提到了一種約束：
可以看到這類似于影像處理中的濾波操作，通過比較每個高程點與周圍的平均高差，如果大于一個閾值，則為重要點，否則為不重要點，其中的關鍵點就是求空間點與直線的距離，具體演算法可參看這篇文章《空間點與直線距離演算法》，
3??TIN構建
經過保留重要點法過濾之后，剩下的點就要進行構網了，一般來說最好構建成Delaunay三角網（因為Delaunay三角網具有很多最優特性），Delaunay三角網的構建演算法也挺復雜，不過可以通過計算幾何演算法庫CGAL來構建，
查閱CGAL的檔案，發現CGAL居然已經有了GIS專題，里面有許多與地形處理相關的示例，其中一個示例就是通過點集生成了Delaunay三角網，并且生成了.ply檔案，.ply檔案正好是一種三維資料格式，能夠被很多三維軟體打開，
4??具體實作
解決了兩個關鍵演算法，具體實作就很簡單了：引入GDAL資料來處理地形資料（.tif），遍歷每個像素點（高程點）做濾波操作，通過CGAL來構建TIN：
#include <iostream>
#include <string>
#include <Vec3.hpp>
#include <threeCGAL.h>
#include <gdal_priv.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Projection_traits_xy_3.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/boost/graph/graph_traits_Delaunay_triangulation_2.h>
#include <CGAL/boost/graph/copy_face_graph.h>
#include <CGAL/Point_set_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
using Projection_traits = CGAL::Projection_traits_xy_3<Kernel>;
using Point_2 = Kernel::Point_2;
using Point_3 = Kernel::Point_3;
using Segment_3 = Kernel::Segment_3;
// Triangulated Irregular Network
using TIN = CGAL::Delaunay_triangulation_2<Projection_traits>;
using namespace std;
int main(int argc, char *argv[])
{
GDALAllRegister();
string demPath = "D:/Work/DEM2TIN/DEM.tif";
string tinPath = "D:/Work/DEM2TIN/Tin.ply";
GDALDataset* img = (GDALDataset *)GDALOpen(demPath.c_str(), GA_ReadOnly);
if (!img)
{
cout << "Can't Open Image!" << endl;
return 1;
}
int imgWidth = img->GetRasterXSize(); //影像寬度
int imgHeight = img->GetRasterYSize(); //影像高度
int bandNum = img->GetRasterCount(); //波段數
//int depth = GDALGetDataTypeSize(img->GetRasterBand(1)->GetRasterDataType()) / 8; //影像深度
int depth = sizeof(float); //影像深度
double padfTransform[6];
img->GetGeoTransform(padfTransform);
double dx = padfTransform[1];
double startx = padfTransform[0] + 0.5 * dx;
double dy = -padfTransform[5];
double starty = padfTransform[3] - imgHeight * dy + 0.5 * dy;
//申請buf
int bufWidth = imgWidth;
int bufHeight = imgHeight;
size_t imgBufNum = (size_t)bufWidth * bufHeight * bandNum;
size_t imgBufOffset = (size_t)bufWidth * (bufHeight - 1) * bandNum;
float *pblock = new float[imgBufNum];
//讀取
img->RasterIO(GF_Read, 0, 0, bufWidth, bufHeight, pblock + imgBufOffset, bufWidth, bufHeight,
GDT_Float32, bandNum, nullptr, bandNum*depth, -bufWidth * bandNum*depth, depth);
CGAL::Point_set_3<Point_3> points;
double zThreshold = 5;
//
for (int yi = 0; yi < imgHeight; yi++)
{
for (int xi = 0; xi < imgWidth; xi++)
{
//將四個角點的約束加入，保證與DEM范圍一致
if ((xi == 0 && yi == 0) || (xi == imgWidth - 1 && yi == 0) ||
(xi == imgWidth - 1 && yi == imgHeight - 1) || (xi == 0 && yi == imgHeight - 1))
{
double gx1 = startx + dx * xi;
double gy1 = starty + dy * yi;
size_t m11 = (size_t)(imgWidth)* yi + xi;
tinyCG::Vec3d P(gx1, gy1, pblock[m11]);
points.insert(Point_3(P.x(), P.y(), P.z()));
}
else
{
double gx0 = startx + dx * (xi - 1);
double gy0 = starty + dy * (yi - 1);
double gx1 = startx + dx * xi;
double gy1 = starty + dy * yi;
double gx2 = startx + dx * (xi + 1);
double gy2 = starty + dy * (yi + 1);
size_t m00 = (size_t)imgWidth * (yi - 1) + xi - 1;
size_t m01 = (size_t)imgWidth * (yi - 1) + xi;
size_t m02 = (size_t)imgWidth * (yi - 1) + xi + 1;
size_t m10 = (size_t)imgWidth* yi + xi - 1;
size_t m11 = (size_t)imgWidth* yi + xi;
size_t m12 = (size_t)imgWidth* yi + xi + 1;
size_t m20 = (size_t)imgWidth * (yi + 1) + xi - 1;
size_t m21 = (size_t)imgWidth * (yi + 1) + xi;
size_t m22 = (size_t)imgWidth * (yi + 1) + xi + 1;
tinyCG::Vec3d P(gx1, gy1, pblock[m11]);
double zMeanDistance = 0;
int counter = 0;
if(m00 < imgBufNum && m22 < imgBufNum)
{
tinyCG::Vec3d A(gx0, gy0, pblock[m00]);
tinyCG::Vec3d E(gx2, gy2, pblock[m22]);
zMeanDistance = zMeanDistance + tinyCG::threeCGAL::CalDistancePointAndLine(P, A, E);
counter++;
}
if (m02 < imgBufNum && m20 < imgBufNum)
{
tinyCG::Vec3d C(gx2, gy0, pblock[m02]);
tinyCG::Vec3d G(gx0, gy2, pblock[m20]);
zMeanDistance = zMeanDistance + tinyCG::threeCGAL::CalDistancePointAndLine(P, C, G);
counter++;
}
if (m01 < imgBufNum && m21 < imgBufNum)
{
tinyCG::Vec3d B(gx1, gy0, pblock[m01]);
tinyCG::Vec3d F(gx1, gy2, pblock[m21]);
zMeanDistance = zMeanDistance + tinyCG::threeCGAL::CalDistancePointAndLine(P, B, F);
counter++;
}
if (m12 < imgBufNum && m10 < imgBufNum)
{
tinyCG::Vec3d D(gx2, gy1, pblock[m12]);
tinyCG::Vec3d H(gx0, gy1, pblock[m10]);
zMeanDistance = zMeanDistance + tinyCG::threeCGAL::CalDistancePointAndLine(P, D, H);
counter++;
}
zMeanDistance = zMeanDistance / counter;
if (zMeanDistance > zThreshold)
{
points.insert(Point_3(P.x(), P.y(), P.z()));
}
}
}
}
delete[] pblock;
pblock = nullptr;
GDALClose(img);
// Create DSM
TIN dsm (points.points().begin(), points.points().end());
using Mesh = CGAL::Surface_mesh<Point_3>;
Mesh dsm_mesh;
CGAL::copy_face_graph (dsm, dsm_mesh);
std::ofstream dsm_ofile (tinPath, std::ios_base::binary);
CGAL::set_binary_mode (dsm_ofile);
CGAL::write_ply (dsm_ofile, dsm_mesh);
dsm_ofile.close();
return 0;
}
5??實驗結果
將最終生成的三維模型檔案.ply通過MeshLab打開，渲染效果如下：
通過Global Mapper還可以看到具體的三角構網效果：
??參考
- DEM模型之間的相互轉換
代碼地址1
代碼地址2 提取碼:x0wt
轉載請註明出處，本文鏈接：https://www.uj5u.com/qiye/282232.html
標籤：其他
上一篇：格網DEM生成不規則三角網TIN