自主三维GIS引擎笔记-实现三维球045
时间：2023-12-05
导读：最小GIS迷你地球实现（实现一套最小的三维GIS球体） V1.0.0.0版本.数据加代码比较大（主要是数据，数据有1G多，代码约5000行),无法上传，如需要微信联系(17381925156).效果图：.相机推进后：.1 . 功能目标.1.1 实现基本的卫片数据浏览.1.2 实现
最小GIS迷你地球实现（实现一套最小的三维GIS球体） V1.0.0.0版本
数据加代码比较大（主要是数据，数据有1G多，代码约5000行),无法上传，如需要微信联系(17381925156)
效果图：
相机推进后：
1 . 功能目标
1.1 实现基本的卫片数据浏览
1.2 实现高程数据的浏览
1.3 实现基本的三维操作，平移，旋转缩放功能
2 必备知识
2.1 图形学基础（OpenGL基础知识）,采用OpenGL3.0标准，纹理数组，shader,Instance绘制
2.2 C++编程基础，vector,list,string,map,function,sort,c++
17,匿名函数(lamda表达式)，智能指针
2.3 多线程支持，多线程采用生产者与消费者模型，信号量,互斥体，队列
2.4 地理信息基本支持，卫片，高程，投影，转换，坐标系
3 软件结构图
3.1 场景管理
3.1.1 相机实现
实现第一人称相机，主要是实现功能，
定点旋转
virtual void rotateViewXByCenter(real angle,real3 pos)
{
if (!(_flag & FLAG_ROT_X))
{
return;
}
//! 计算眼睛到鼠标点的方向
real3 vDir = pos - _eye;
real len1 = length(vDir);
real len = 0;
vDir = normalize(vDir);
matrix4r mat(1);
mat.rotate(angle, _right);
vDir = vDir * mat;
_eye = pos - vDir * len1;
_dir = _dir * mat;
_up = _up * mat;
#ifndef COORD_LH
_right = CELL::normalize(cross(_dir,_up));
_up = CELL::normalize(cross(_right,_dir));
#else
_right = CELL::normalize(cross(_up,_dir));
_up = CELL::normalize(cross(_dir,_right));
#endif
len = CELL::length(_eye - _target);
_target = _eye + _dir * len;
update();
}
定点缩放
/**
* 指定点推进摄像机
*/
virtual void scaleCameraByPos(const real3& pos,real persent)
{
real3 dir = CELL::normalize(pos - _eye);
real dis = CELL::length(pos - _eye) * persent;
real disCam = CELL::length(_target - _eye) * persent;
real3 dirCam = CELL::normalize(_target - _eye);
_eye = pos - dir * dis;
_target = _eye + dirCam * disCam;
update();
}
精确移动
virtual void moveLeft(real fElapsed)
{
_eye -= normalize(_right) * _speed * fElapsed;
_target -= normalize(_right) * _speed * fElapsed;
}
virtual void moveRight(real fElapsed)
{
_eye += normalize(_right) * _speed * fElapsed;
_target += normalize(_right) * _speed * fElapsed;
}
virtual void moveFront(real fElapsed)
{
_eye += _dir * _speed * fElapsed;
_target += _dir * _speed * fElapsed;
}
virtual void moveBack(real fElapsed)
{
_eye -= _dir * _speed * fElapsed;
_target -= _dir * _speed * fElapsed;
}
/**
* 往上看，则向下移动摄像机
*/
virtual void moveUp(real fElapsed)
{
_eye += _up * _speed * fElapsed;
_target += _up * _speed * fElapsed;
}
/**
* 向下看，则摄像机向上移动
*/
virtual void moveDown(real fElapsed)
{
_eye -= _up * _speed * fElapsed;
_target -= _up * _speed * fElapsed;
}
/**
* 根据给定的方向移动摄像机
*/
virtual void moveDir(real3 dir,real fElapsed)
{
_eye += dir * _speed * fElapsed;
_target += dir * _speed * fElapsed;
}
实现屏幕到世界坐标投影转换
/**
* 窗口坐标转化为世界坐标
*/
real3 screenToWorld(const real2& screen)
{
real4 screens(screen.x,screen.y,0,1);
real4 world;
unProject(screens,world);
return real3(world.x,world.y,world.z);
}
世界坐标到屏幕坐标投影转换
/**
* 世界坐标转化为窗口坐标
*/
real2 worldToScreen( const real3& world)
{
real4 worlds(world.x,world.y,world.z,1);
real4 screens;
project(worlds,screens);
return real2(screens.x,screens.y);
}
根据屏幕坐标点创建射线
Ray createRayFromScreen(int x,int y)
{
real4 minWorld;
real4 maxWorld;
real4 screen(real(x),real(y),0,1);
real4 screen1(real(x),real(y),1,1);
unProject(screen,minWorld);
unProject(screen1,maxWorld);
Ray ray;
ray.setOrigin(real3(minWorld.x,minWorld.y,minWorld.z));
real3 dir(maxWorld.x - minWorld.x,maxWorld.y - minWorld.y, maxWorld.z - minWorld.z);
ray.setDirection(normalize(dir));
return ray;
}
3.1.2 裁剪实现
根据投影矩阵，观察矩阵构造视锥
1. 支持判断点是否在视锥内
2. 支持判断是否球体在视锥内
3. 支持判断包围盒是否在视锥内。
template<class T>
class tfrustum
{
public:
enum
{
FRUSTUM_LEFT = 0,
FRUSTUM_RIGHT = 1,
FRUSTUM_TOP = 2,
FRUSTUM_BOTTOM = 3,
FRUSTUM_FAR = 4,
FRUSTUM_NEAR = 5,
};
public:
/**
* project * modleview
*/
void loadFrustum(const tmat4x4<T> &mvp)
{
const T* dataPtr = mvp.data();
_planes[FRUSTUM_LEFT ] = Plane<T>(dataPtr[12] - dataPtr[0], dataPtr[13] - dataPtr[1], dataPtr[14] - dataPtr[2], dataPtr[15] - dataPtr[3]);
_planes[FRUSTUM_RIGHT ] = Plane<T>(dataPtr[12] + dataPtr[0], dataPtr[13] + dataPtr[1], dataPtr[14] + dataPtr[2], dataPtr[15] + dataPtr[3]);
_planes[FRUSTUM_TOP ] = Plane<T>(dataPtr[12] - dataPtr[4], dataPtr[13] - dataPtr[5], dataPtr[14] - dataPtr[6], dataPtr[15] - dataPtr[7]);
_planes[FRUSTUM_BOTTOM] = Plane<T>(dataPtr[12] + dataPtr[4], dataPtr[13] + dataPtr[5], dataPtr[14] + dataPtr[6], dataPtr[15] + dataPtr[7]);
_planes[FRUSTUM_FAR ] = Plane<T>(dataPtr[12] - dataPtr[8], dataPtr[13] - dataPtr[9], dataPtr[14] - dataPtr[10], dataPtr[15] - dataPtr[11]);
_planes[FRUSTUM_NEAR ] = Plane<T>(dataPtr[12] + dataPtr[8], dataPtr[13] + dataPtr[9], dataPtr[14] + dataPtr[10], dataPtr[15] + dataPtr[11]);
}
bool pointInFrustum(const tvec3<T> &pos) const
{
for (int i = 0; i < 6; i++)
{
if (_planes[i].distance(pos) <= 0)
return false;
}
return true;
}
bool sphereInFrustum(const tvec3<T> &pos, const float radius) const
{
for (int i = 0; i < 6; i++)
{
if (_planes[i].distance(pos) <= -radius)
return false;
}
return true;
}
bool cubeInFrustum(T minX,T maxX,T minY,T maxY,T minZ,T maxZ) const
{
for (int i = 0; i < 6; i++)
{
if (_planes[i].distance(tvec3<T>(minX, minY, minZ)) > 0) continue;
if (_planes[i].distance(tvec3<T>(minX, minY, maxZ)) > 0) continue;
if (_planes[i].distance(tvec3<T>(minX, maxY, minZ)) > 0) continue;
if (_planes[i].distance(tvec3<T>(minX, maxY, maxZ)) > 0) continue;
if (_planes[i].distance(tvec3<T>(maxX, minY, minZ)) > 0) continue;
if (_planes[i].distance(tvec3<T>(maxX, minY, maxZ)) > 0) continue;
if (_planes[i].distance(tvec3<T>(maxX, maxY, minZ)) > 0) continue;
if (_planes[i].distance(tvec3<T>(maxX, maxY, maxZ)) > 0) continue;
return false;
}
return true;
}
const Plane<T> &getPlane(const int plane) const
{
return _planes[plane];
}
public:
Plane<T> _planes[6];
};
3.1.2 浏览实现
输入输出接口定义，实现对鼠标键盘事件的抽象，后续可以实现跨平台处理
class CELLInput
{
public:
/// <summary>
/// 鼠标左键按下
/// </summary>
virtual void onLButtonDown(int x,int y) = 0;
/// <summary>
/// 鼠标左键提起
/// </summary>
virtual void onLButtonUp(int x, int y) = 0;
virtual void onLButtonDbClick(int x, int y) = 0;
virtual void onRButtonDbClick(int x, int y) = 0;
virtual void onMButtonDbClick(int x, int y) = 0;
/// <summary>
/// 鼠标移动事件
/// </summary>
virtual void onRButtonDown(int x, int y) = 0;
/// <summary>
/// 鼠标移动事件
/// </summary>
virtual void onRButtonUp(int x, int y) = 0;
/// <summary>
/// 鼠标移动事件
/// </summary>
virtual void onMButtonDown(int x, int y) = 0;
/// <summary>
/// 鼠标移动事件
/// </summary>
virtual void onMButtonUp(int x, int y) = 0;
/// <summary>
/// 鼠标移动事件
/// </summary>
virtual void onMouseMove(int x,int y) = 0;
/// <summary>
/// 鼠标移动事件
/// </summary>
virtual void onMouseWheel(int delta) = 0;
/// <summary>
/// 键盘事件
/// </summary>
virtual void onKeyDown(int key) = 0;
/// <summary>
/// 键盘事件
/// </summary>
virtual void onKeyUp(int key) = 0;
};
具体业务代码实现：处理windows事件
LRESULT eventProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
switch (message)
{
case WM_LBUTTONDOWN:
_frame->onLButtonDown(LOWORD(lParam),HIWORD(lParam));
break;
case WM_LBUTTONUP:
_frame->onLButtonUp(LOWORD(lParam),HIWORD(lParam));
break;
case WM_RBUTTONDOWN:
_frame->onRButtonDown(LOWORD(lParam),HIWORD(lParam));
break;
case WM_RBUTTONUP:
_frame->onRButtonUp(LOWORD(lParam), HIWORD(lParam));
break;
case WM_MBUTTONDOWN:
_frame->onMButtonDown(LOWORD(lParam), HIWORD(lParam));
break;
case WM_MBUTTONUP:
_frame->onMButtonUp(LOWORD(lParam), HIWORD(lParam));
break;
case WM_MOUSEMOVE:
_context._mouseX = LOWORD(lParam);
_context._mouseY = HIWORD(lParam);
_frame->onMouseMove(LOWORD(lParam), HIWORD(lParam));
break;
case WM_MOUSEWHEEL:
_frame->onMouseWheel(GET_WHEEL_DELTA_WPARAM(wParam));
break;
case WM_KEYDOWN:
_frame->onKeyDown((int)wParam);
break;
case WM_KEYUP:
_frame->onKeyUp((int)wParam);
break;
case WM_LBUTTONDBLCLK:
_frame->onLButtonDbClick(LOWORD(lParam), HIWORD(lParam));
break;
case WM_RBUTTONDBLCLK:
_frame->onLButtonDbClick(LOWORD(lParam), HIWORD(lParam));
break;
case WM_MBUTTONDBLCLK:
_frame->onMButtonDbClick(LOWORD(lParam), HIWORD(lParam));
break;
case WM_PAINT:
{
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
EndPaint(hWnd, &ps);
}
break;
case WM_SIZE:
{
if(IsWindow(_hWnd))
{
RECT rect;
GetClientRect(_hWnd,&rect);
_context._width = rect.right - rect.left;
_context._height = rect.bottom - rect.top;
}
}
break;
case WM_DESTROY:
_threadRun = false;
CELLThread::join();
PostQuitMessage(0);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
return S_OK;
具体动作实现：
void CELLFrameBigMap::update(CELLContext& )
{
_context._device->setViewPort(0,0,_context._width,_context._height);
_context._screenPrj = CELL::ortho<real>(0.0f,(real)_context._width,(real)_context._height,0,-1000.0f,1000.0f);
_context._camera.setViewSize(real2(_context._width,_context._height));
_context._camera.perspective(45.0,real(_context._width)/real(_context._height),10.0,FSIZE * 10);
_context._camera.update();
_context._mvp = _context._camera._matProj * _context._camera._matView * _context._camera._matWorld ;
_context._vp = _context._camera._matProj * _context._camera._matView;
_context._timePerFrame = (float)_timeStamp.getElapsedSecond();
_timeStamp.update();
matrix4r matVP = _context._vp.transpose();
_context._frustum.loadFrustum(matVP);
if (_context._keyState[VK_UP])
{
_context._camera.moveFront(_context._timePerFrame);
}
if (_context._keyState[VK_DOWN])
{
_context._camera.moveBack(_context._timePerFrame);
}
if (_context._keyState[VK_LEFT])
{
_context._camera.moveLeft(_context._timePerFrame);
}
if (_context._keyState[VK_RIGHT])
{
_context._camera.moveRight(_context._timePerFrame);
}
_terrain.update(_context);
}
void CELLFrameBigMap::onFrameStart(CELLContext& context)
{
context._device->clear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
context._device->clearColor(0, 0, 0, 1);
context._device->disableRenderState(GL_CULL_FACE);
context._device->enableRenderState(GL_DEPTH_TEST);
_terrain.render(context);
}
void CELLFrameBigMap::onFrameEnd(CELLContext& context)
{
}
void CELLFrameBigMap::onLButtonDown(int x, int y)
{
getPointsFromScreen(x,y,_basePoint);
_bLbuttonDown = true;
_lbuttonDown = int2(x,y);
}
void CELLFrameBigMap::onLButtonUp(int x, int y)
{
_lbuttonDown = int2(x,y);
_bLbuttonDown = false;
}
void CELLFrameBigMap::onLButtonDbClick(int x, int y)
{
/// _basePoint
real3 vDir = CELL::normalize(-_basePoint);
real dist = CELL::length(_context._camera.getEye() - _context._camera.getTarget());
_context._camera.setTarget(_basePoint);
real3 vEye = _basePoint - dist * vDir;
_context._camera.setEye(vEye);
_context._camera.calcDir();
_context._camera.rotateViewXByCenter(0,_basePoint);
_context._camera.update();
}
void CELLFrameBigMap::onRButtonDbClick(int x, int y)
{
}
void CELLFrameBigMap::onMButtonDbClick(int x, int y)
{
}
void CELLFrameBigMap::onRButtonDown(int , int )
{
}
void CELLFrameBigMap::onRButtonUp(int , int )
{
}
void CELLFrameBigMap::onMButtonDown(int x, int y)
{
_mButtonDown = int2(x,y);
_bMButtonDown = true;
}
void CELLFrameBigMap::onMButtonUp(int x, int y)
{
_mButtonDown = int2(x,y);
_bMButtonDown = false;
}
void CELLFrameBigMap::onMouseMove(int x, int y)
{
if (_bLbuttonDown)
{
int2 curPoint(x, y);
int2 offset = curPoint - _lbuttonDown;
_lbuttonDown = curPoint;
_context._camera.rotateViewYByCenter(offset.x, _basePoint);
_context._camera.rotateViewXByCenter(offset.y, _basePoint);
}
if (_bMButtonDown)
{
int2 ofScreen = int2(x,y) - _mButtonDown;
_mButtonDown = int2(x,y);
moveScene(_basePoint,ofScreen);
}
}
void CELLFrameBigMap::onMouseWheel(int delta)
{
double persent = 1;
if (delta > 0)
{
persent = 0.9;
}
else
{
persent = 1.1;
}
_context._camera.scaleCameraByPos(_basePoint,persent);
}
3.2 投影管理
3.2.1 墨卡托投影实现
定义
Web Mercator 坐标系使用的投影方法不是严格意义的墨卡托投影，而是一个被 EPSG（European Petroleum Survey Group）
称为伪墨卡托的投影方法，这个伪墨卡托投影方法的大名是 Popular Visualization Pseudo Mercator，PVPM。因为这个坐标
系统是 Google Map 最先使用的，或者更确切地说，是Google 最先发明的。在投影过程中，将表示地球的参考椭球体近似的
作为正球体处理（正球体半径 R = 椭球体半长轴 a）。后来，Web Mercator 在 Web 地图领域被广泛使用，这个坐标系就名声
大噪。尽管这个坐标系由于精度问题一度不被GIS专业人士接受，但最终 EPSG 还是给了 WKID:3857。以整个世界范围，
赤道作为标准纬线，本初子午线作为中央经线，两者交点为坐标原点，向东向北为正，向西向南为负。
X轴：由于赤道半径为6378137米，则赤道周长为2PIr = 20037508.3427892，因此X轴的取值范围：[-20037508.3427892,20037508.3427892]。
Y轴：由墨卡托投影的公式可知，当纬度φ接近90°时，y值趋向于无穷。
为了方便实现一级地图在长度上一分为二，也就是在高一级精度上用4张切片地图显示，最好将xy轴显示的取值范围统一，
为此Y轴的取值范围也限定在[-20037508.3427892,20037508.3427892]之间。
通过计算纬度最大值就是85.0511287798066°
另外值得一提的是，现在网上大部分地图瓦片瓦片都是基于Web墨卡托的，就是第0层1一张图片，第1层2*2张图片…,如果使用OpenLayers可以很好的调用公共服务器
编号
EPSG:3857
EPSG:102100
EPSG:900913 (谷歌自己起的，没有被EPSG接受，900913=google)
转换实现
1 class GlobalMercator
2 {
3 protected:
4 int tileSize ;
5 real initialResolution;
6 real originShift;
7
8 public:
9 GlobalMercator()
10 {
11 tileSize = 256;
12 initialResolution = 2 * PI * 6378137 / tileSize;
13 originShift = 2 * PI * 6378137 / 2.0;
14 }
15
16 real2 LatLonToMeters(real lat, real lon )
17 {
18 //"Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913"
19 real mx = lon * originShift / 180.0 ;
20 real my = log( tan((90 + lat) *PI / 360.0 )) / (PI / 180.0);
21 my = my * originShift / 180.0;
22 return real2(mx, my);
23 }
24
25 real2 MetersToLatLon(int mx,int my )
26 {
27 real lon = (mx / originShift) * 180.0;
28 real lat = (my / originShift) * 180.0;
29 lat = real(180 / PI * (2 * atan( exp( lat * PI / 180.0)) - PI / 2.0));
30 return real2(lon,lat);
31 }
32
33 real resolution(int zoom )
34 {
35 // return (2 * math.pi * 6378137) / (self.tileSize * 2**zoom)
36 return initialResolution / (pow(2,real(zoom)));
37 }
38 real2 PixelsToMeters(int px, int py, int zoom)
39 {
40 real res = resolution( zoom );
41 real mx = px * res - originShift;
42 real my = py * res - originShift;
43 return real2(mx,my);
44 }
45 int2 MetersToPixels(real mx, real my, int zoom)
46 {
47 real res = resolution( zoom );
48 int px = int((mx + originShift) / res) ;
49 int py = int((my + originShift) / res);
50 return int2(px,py);
51 }
52
53 int2 PixelsToTile(int px, int py)
54 {
55 int tx = int( ceil( px / real(tileSize) ) - 1 );
56 int ty = int( ceil( py / real(tileSize) ) - 1 );
57 return int2(tx,ty);
58 }
59
60 int2 MetersToTile(real mx, real my, int zoom)
61 {
62 int2 vPt = MetersToPixels( mx, my, zoom);
63 return PixelsToTile(vPt.x,vPt.y) ;
64 }
65 int2 LongLatToTile(real dLon,real dLat,int zoom)
66 {
67 real2 vMeter = LatLonToMeters(dLat,dLon);
68
69 return MetersToTile(vMeter.x,vMeter.y,zoom);
70 }
71
72 };
1 class CELLMercator :public CELLSpRef
2 {
3 protected:
4 GlobalMercator _proj;
5 public:
6
7 CELLMercator(void)
8 {
9
10 }
11
12 ~CELLMercator(void)
13 {
14 }
15 /**
16 * 将经度转化为x-tile key
17 */
18 static int long2tilex(real lon, int z)
19 {
20 return (int)(floor((lon + 180.0) / 360.0 * pow(2.0, z)));
21 }
22 /**
23 * 将纬度转化为y-tile key
24 */
25 static int lat2tiley(real lat, int z)
26 {
27 return (int)(floor((1.0 - log( tan(lat * PI/180.0) + 1.0 / cos(lat * PI/180.0)) / PI) / 2.0 * pow(2.0, z)));
28 }
29 /**
30 * 给定 x-tile 获取经度
31 */
32 static real tilex2long(int x, int z)
33 {
34 return x / pow(2.0, z) * 360.0 - 180;
35 }
36 /**
37 * 给定 y-tile 获取纬度
38 */
39 static real tiley2lat(int y, int z)
40 {
41 real n = PI - 2.0 * PI * y / pow(2.0, z);
42 return 180.0 / PI * atan(0.5 * (exp(n) - exp(-n)));
43 }
44
45 static real2 tile2lonlat(int x,int y, int z)
46 {
47 real p = pow(2.0, z);
48 real n = PI - 2.0 * PI * y / p;
49 real lat = 180.0 / PI * atan(0.5 * (exp(n) - exp(-n)));
50 real lon = x / p * 360.0 - 180;
51 return real2(lon,lat);
52 }
53 real2 tile2World(int x,int y, int z)
54 {
55 real2 lonLat = tile2lonlat(x,y,z);
56 return _proj.LatLonToMeters(lonLat.x,lonLat.y);
57 }
58 /**
59 * 根据经纬度与级别返回瓦片Key
60 */
61 virtual int3 getKey(unsigned level, real rLong,real rLat)
62 {
63 //! 注意，这里不做显示，是为了处理地图滚动，即地图的显示范围可以不进行限制
64 #if 1
65 if (rLong <= -180)
66 {
67 rLong = -179.9;
68 }
69 if (rLong >= 180)
70 {
71 rLong = 179.9;
72 }
73 if (rLat < -85)
74 {
75 rLat = -85;
76 }
77 if (rLat > 85)
78 {
79 rLat = 85;
80 }
81 #endif
82 int levelTileNumber = 0;
83 levelTileNumber = 1<<level;
84
85 int xTile = long2tilex(rLong,level);
86
87 int yTile = lat2tiley(rLat,level);
88
89 return int3(xTile,yTile,level);
90 }
91 public:
92 /**
93 * 经纬度转化为世界坐标
94 */
95 virtual real3 longLatToWorld(const real3& longLatx) override
96 {
97 //return longLatx;
98 real3 longLat = longLatx;
99 longLat.x = tmin<real>(179.9999f,longLatx.x);
100 longLat.x = tmax<real>(-179.9999f,longLat.x);
101
102 longLat.y = tmin<real>(85.0f,longLatx.y);
103 longLat.y = tmax<real>(-85.0f,longLat.y);
104 auto res = _proj.LatLonToMeters(longLatx.y,longLatx.x);
105 return real3(res.x,0,res.y);
106 }
107
108 /**
109 * 世界坐标转化为经纬度
110 */
111 virtual real3 worldToLongLat(const real3& world) override
112 {
113 const real worldMin = -20037504.0f * 1000;
114 const real worldMax = 20037504.0f * 1000;
115 real dWordX = (real)tmin<real>(world.x, worldMax);
116 dWordX = (real)tmax<real>((real)dWordX, worldMin);
117
118 real dWordY = (real)tmin<real>(world.y, worldMax);
119 dWordY = (real)tmax<real>((real)dWordY, worldMin);
120 auto res = _proj.MetersToLatLon((int)dWordX, (int)dWordY);
121 return real3(res.x,res.y,0);
122 }
123
124 /**
125 * 得到当前级别下tile的个数
126 */
127 int getTileNumber(int lev)
128 {
129 return (int)pow(2,real(lev));
130 }
131 };
3.2.2 wgs84投影实现
定义
WGS的最新版本为WGS 84（也称作WGS 1984、EPSG:4326），1984年定义、最后修订于2004年。
投影过程如下图所示：
球体被切成一个个小圆弧，一共60个投影带，分别为01,02………60
分为60个投影带，每一个投影带投影通用编码分别为：32601,32650,32660，从上图我们可以看到从南到北，
同一个投影带又被划分为分别为X,W,V等小格子：在同一个投影带内，比如50投影带内，通用编码为32650，依次又形成了，
50X,50W,50V,50S等小格子
代码实现
1 class CELLWgs842d :public CELLSpRef
2 {
3 public:
4
5 /**
6 * 经纬度转化为世界坐标
7 */
8 virtual real3 longLatToWorld(const real3& longLatx) override
9 {
10 real3 world;
11 world.x = longLatx.x * WGS_84_RADIUS_EQUATOR;
12 world.y = longLatx.z;
13 world.z = longLatx.y * WGS_84_RADIUS_EQUATOR;
14 return world;
15 }
16
17 /**
18 * 世界坐标转化为经纬度
19 */
20 virtual real3 worldToLongLat(const real3& world) override
21 {
22 real3 lonlat;
23 lonlat.x = world.x / WGS_84_RADIUS_EQUATOR;
24 lonlat.y = world.y / WGS_84_RADIUS_EQUATOR;
25 lonlat.y = world.z;
26 return lonlat;
27 }
28
29 /**
30 * 根据经纬度与级别返回瓦片Key
31 */
32 virtual int3 getKey(unsigned level, real rLong,real rLat) override
33 {
34 unsigned xTiles = 2<<level;
35 unsigned yTiles = 1<<level;
36 double xWidth = TWO_PI / xTiles;
37 double yHeight = PI / yTiles;
38 double xCoord = (rLong + PI) / xWidth;
39 if (xCoord >= xTiles)
40 {
41 xCoord = xTiles - 1;
42 }
43
44 double yCoord = (HALF_PI - rLat) / yHeight;
45 if (yCoord >= yTiles)
46 {
47 yCoord = yTiles - 1;
48 }
49
50 return int3(int(xCoord),int(yCoord),level);
51 }
52 };
3.3.2 wgs84球体投影实现
定义：
实现
1 class CELLWgs84Sphere :public CELLSpRef
2 {
3 protected:
4 tellipsoidModelIn<double> _spr;
5 public:
6
7 /**
8 * 经纬度转化为世界坐标
9 */
10 virtual real3 longLatToWorld(const real3& longLatx) override
11 {
12 return _spr.longlatAltToWorld(longLatx);
13 }
14
15 /**
16 * 世界坐标转化为经纬度
17 */
18 virtual real3 worldToLongLat(const real3& world) override
19 {
20 return _spr.worldToLonglatAlt(world);
21 }
22
23 /**
24 * 根据经纬度与级别返回瓦片Key
25 */
26 virtual int3 getKey(unsigned level, real rLong,real rLat) override
27 {
28 unsigned xTiles = 2<<level;
29 unsigned yTiles = 1<<level;
30 double xWidth = TWO_PI / xTiles;
31 double yHeight = PI / yTiles;
32 double xCoord = (rLong + PI) / xWidth;
33 if (xCoord >= xTiles)
34 {
35 xCoord = xTiles - 1;
36 }
37
38 double yCoord = (HALF_PI - rLat) / yHeight;
39 if (yCoord >= yTiles)
40 {
41 yCoord = yTiles - 1;
42 }
43
44 return int3(int(xCoord),int(yCoord),level);
45 }
3.3 四叉树
结构定义
1 class CELLQuadTree :public CELLObject
2 {
3 public:
4 enum ChildId
5 {
6 CHILD_LT,
7 CHILD_RT,
8 CHILD_LB,
9 CHILD_RB,
10 };
11 enum
12 {
13 FLAG_HAS_IMAGE = 1<<0,
14 FLAG_HAS_DEM = 1<<1,
15 FLAG_HAS_CULL = 1<<2,
16 FLAG_RENDER = 1<<3,
17 };
18 public:
19 typedef std::vector<CELLQuadTree*> ArrayNode;
20 typedef CELLTerrainInterface TerrainApi;
21 using NodePtr = std::shared_ptr<CELLQuadTree>;
22 public:
23 TerrainApi* _terrain;
24 /// 数据标志
25 uint _flag;
26 /// 对应瓦片id
27 TileId _tileId;
28 real2 _vStart;
29 real2 _vEnd;
30 matrix4r _local = matrix4r(1);
31
32 float2 _uvStart;
33 float2 _uvEnd;
34 /// 对应瓦片的范围(世界坐标)
35 aabb3dr _aabb;
36 /// 位置
37 ChildId _cornerId;
38 /// 当前瓦片的父节点
39 CELLQuadTree* _parent;
40 /// 瓦片的孩子节点
41 NodePtr _childs[4];
42 /// 纹理id
43 uint _textureId;
44 /// dem数据
45 DemData _demData;
46 VertexBufferId _vbo;
47 IndexBufferId _ibo;
48 };
3.3.1 对象管理
GIS 系统中会频繁大量加载卸载地图数据（动态加卸载），时间长了会产生大量的内存碎片，影像加载绘制效率，为了提升性能
系统使用内存池以及对象池来优化瓦片数据。
引入基类定义，方便扩展新的特性
class CELLPool :public CELLObject
{
public:
};
using PoolPtr = std::shared_ptr<CELLPool>;
内存池定义如下：
使用模板实现，采用单例设计模式对外提供内存管理能力
template<size_t blackSize>
class CELLPoolMemory :public CELLPool
{
public:
struct Data
{
char _data[blackSize];
};
using ObjectList = std::list<void*>;
protected:
ObjectList _objests;
public:
~CELLPoolMemory()
{
destroy();
}
void init(size_t nCnt)
{
for (size_t i = 0; i < nCnt; i++)
{
auto obj = new Data();
_objests.push_back(obj);
}
}
void* alloc()
{
if (_objests.empty())
{
return new Data();
}
else
{
void* ptr = _objests.back();
_objests.pop_back();
return ptr;
}
}
/// <summary>
/// 回收函数
/// </summary>
void free(void* obj)
{
_objests.push_back(obj);
}
void destroy()
{
for (auto var : _objests)
{
delete var;
}
_objests.clear();
}
public:
static CELLPoolMemory<blackSize>& instance()
{
static CELLPoolMemory<blackSize> sInstance;
return sInstance;
}
};
对象池定义如下：
使用模板实现，采用单例设计模式对外提供内存管理能力
template<class T>
class CELLPoolObject :public CELLPool
{
public:
using TPtr = T*;
using ObjectList = std::list<TPtr>;
using SPtr = std::shared_ptr<T>;
protected:
ObjectList _objests;
public:
/// <summary>
/// 初始化
/// </summary>
template<class ...Args>
void init(size_t nCnt,Args ...args)
{
for (size_t i = 0; i < nCnt; i++)
{
TPtr ptr = new T(args...);
_objests.push_back(ptr);
}
}
/// <summary>
/// 分配函数
/// </summary>
template<class ...Args>
SPtr alloc(Args ...args)
{
if (_objests.empty())
{
return SPtr(new T(args...),deleteObject);
}
else
{
TPtr ptr = _objests.back();
_objests.pop_back();
return SPtr(ptr,deleteObject);
}
}
/// <summary>
/// 回收函数
/// </summary>
void free(TPtr obj)
{
_objests.push_back(obj);
}
void destroy()
{
_objests.clear();
}
public:
static inline void deleteObject(T* obj)
{
instance().free(obj);
}
static CELLPoolObject<T>& instance()
{
static CELLPoolObject<T> sInstance;
return sInstance;
}
};
3.3.2 网络数据生成
主要是生成可以绘制的三角网络数据，为了统一接口，实现如下；
struct V3U3N4
{
float x, y, z;
float h;
float u, v;
};
template<ushort row,ushort col>
class CELLDemMesh :public CELLObject
{
protected:
V3U3N4 _data[row * col];
ushort3 _faces[(row -1) * (col - 1) * 2];
int _col = 0;
int _row = 0;
aabb3dr _aabb;
public:
CELLDemMesh()
{
_row = row;
_col = col;
}
void redefine(int c,int r)
{
_col = c;
_row = r;
}
int getRow() const
{
return _row;
}
int getCol() const
{
return _col;
}
V3U3N4* vertData()
{
return _data;
}
const V3U3N4* vertData() const
{
return _data;
}
int vertCount() const
{
return _row * _col;
}
int vertByte() const
{
return vertCount() * sizeof(V3U3N4);
}
ushort3* faceData()
{
return _faces;
}
const ushort3* faceData() const
{
return _faces;
}
int faceCount() const
{
return (_row -1) * (_col - 1) * 2;
}
int faceByte() const
{
return faceCount() * sizeof(ushort3);
}
void fillFace()
{
uint fSize = (_row -1) * (_col - 1) * 2;
uint index = 0;
for (ushort r = 0; r < _row - 1; ++r)
{
for (ushort c = 0; c < _col - 1; ++c)
{
_faces[index + 0][0] = (_col * r) + c;
_faces[index + 0][1] = (_col * r) + c + 1;
_faces[index + 0][2] = (_col * (r + 1)) + c;
_faces[index + 1][0] = (_col * r) + c + 1;
_faces[index + 1][1] = _col * (r + 1) + c + 1;
_faces[index + 1][2] = _col * (r + 1) + c;
index += 2;
}
}
}
void fillVertex(const real2& vStart,const real2& vEnd,bool bInitY,CELLSpRef* spRef)
{
real2 vSize = vEnd - vStart;
real2 vGrid = real2(vSize.x/(_col -1),vSize.y/(_row - 1));
_aabb.setNull();
real3 tmp[row * col];
if (bInitY)
{
for (ushort r = 0; r < _row; ++r)
{
for (ushort c = 0; c < _col; ++c)
{
int idx = r * _col + c;
real3 vWorld = spRef->longLatToWorld(real3(vStart.x + c * vGrid.x,vStart.y + r * vGrid.y,0));
tmp[idx] = vWorld;
/// _data[idx].x = vWorld.x;
/// _data[idx].y = vWorld.y;
/// _data[idx].z = vWorld.z;
_aabb.merge(vWorld);
}
}
}
else
{
for (ushort r = 0; r < _row; ++r)
{
for (ushort c = 0; c < _col; ++c)
{
int idx = r * _col + c;
real3 vWorld = spRef->longLatToWorld(real3(vStart.x + c * vGrid.x,vStart.y + r * vGrid.y,_data[idx].h));
tmp[idx] = vWorld;
/// _data[idx].x = vWorld.x;
/// _data[idx].y = vWorld.y;
/// _data[idx].z = vWorld.z;
_aabb.merge(vWorld);
}
}
}
real3 vCenter = _aabb.getCenter();
for (size_t i = 0; i < row * col; i++)
{
real3 vLocal = tmp[i] - vCenter;
_data[i].x = vLocal.x;
_data[i].y = vLocal.y;
_data[i].z = vLocal.z;
}
}
void fillUV(const real2& vStart,const real2& vEnd,int layer)
{
real2 uvSize = vEnd - vStart;
real2 uvStep = real2(uvSize.x/(_col -1),uvSize.y/(_row - 1));
for (ushort r = 0; r < _row; ++r)
{
for (ushort c = 0; c < _col; ++c)
{
int idx = r * _col + c;
_data[idx].u = (vStart.x + real(c) * uvStep.x);
_data[idx].v = (vStart.y + real(r) * uvStep.y);
}
}
}
void fillHeight(CELLFiledHeight* pHeight)
{
float* pData = pHeight->data();
for (size_t i = 0 ;i < _row * _col ; ++ i)
{
_data[i].h = pData[i];
}
}
aabb3dr getAabb() const
{
return _aabb;
}
/// <summary>
/// 瓦片与射线相交测试
/// </summary>
bool intersect(const CELL::Ray& ray,const matrix4r& mat,real& dist)
{
int faceCnt = faceCount();
real time = FLT_MAX;
bool bPick = false;
auto result = ray.intersects(_aabb);
if (!result.first)
return false;
for (int i = 0; i < faceCnt; i++)
{
int i1 = _faces[i].x;
int i2 = _faces[i].y;
int i3 = _faces[i].z;
real3 v1 = mat * real3(_data[i1].x,_data[i1].y,_data[i1].z);
real3 v2 = mat * real3(_data[i2].x,_data[i2].y,_data[i2].z);
real3 v3 = mat * real3(_data[i3].x,_data[i3].y,_data[i3].z);
real cur = FLT_MAX;
real u = 0;
real v = 0;
bool res = CELL::intersectTriangle<real>(ray.getOrigin(),ray.getDirection(),v1,v2,v3,&cur,&u,&v);
if (res)
{
time = CELL::tmin(time,cur);
bPick = true;
}
}
if (bPick)
{
dist = time;
}
return bPick;
}
};
using DemData = CELLDemMesh<DEM_W,DEM_H>;
3.3.3 裂分算法实现（核心部分实现)
CELLCamera& camera = context._camera;
real3 vWCenter = _aabb.getCenter();
real3 vWSize = _aabb.getSize();
real fSize = CELL::length(vWSize) * 0.5;
real dis = CELL::length(vWCenter - camera._eye);
if (context._frustum.cubeInFrustum(
_aabb._minimum.x
,_aabb._maximum.x
,_aabb._minimum.y
,_aabb._maximum.y
,_aabb._minimum.z
,_aabb._maximum.z))
{
_flag &= ~FLAG_HAS_CULL;
}
else
{
_flag |= FLAG_HAS_CULL;
}
if (dis/fSize < 3 && hasNoFlag(FLAG_HAS_CULL))
{
if(!hasChild() && hasImage())
{
real2 vLlCenter = lonLatCenter();
real2 vLLHalf = lonLatRange() * 0.5;
real3 vCenter = real3(vLlCenter.x, 0, vLlCenter.y);
real3 vSize = real3(vLLHalf.x, 0, vLLHalf.y);
if (_terrain->getTaskCount() > 20)
{
return ;
}
_childs[CHILD_LT] = NodePtr(new CELLQuadTree(
_terrain
, this
,real2(vCenter.x - vSize.x,vCenter.z)
,real2(vCenter.x,vCenter.z + vSize.z)
,(int)_tileId._lev + 1
,CHILD_LT
));
_childs[CHILD_RT] = NodePtr(new CELLQuadTree(
_terrain
,this
, real2(vCenter.x, vCenter.z)
, real2(vCenter.x + vSize.x, vCenter.z + vSize.z)
, (int)_tileId._lev + 1
, CHILD_RT
));
_childs[CHILD_LB] = NodePtr(new CELLQuadTree(
_terrain
, this
, real2(vCenter.x - vSize.x, vCenter.z - vSize.z)
, real2(vCenter.x, vCenter.z)
, (int)_tileId._lev + 1
, CHILD_LB
));
_childs[CHILD_RB] = NodePtr(new CELLQuadTree(
_terrain
, this
, real2(vCenter.x, vCenter.z - vSize.z)
, real2(vCenter.x + vSize.x, vCenter.z)
, (int)_tileId._lev + 1
, CHILD_RB
));
}
else
{
for (int i = 0; i < 4; ++i)
{
if (_childs[i] && hasNoFlag(FLAG_HAS_CULL))
{
_childs[i]->update(context);
}
else
{
_flag &= FLAG_RENDER;
}
}
}
}
else
{
for (int i = 0 ;i < 4 ; ++ i)
{
_childs[i] = nullptr;
}
}
3.4 图层管理
##
3.5 文件格式管理
为了统一接口，提炼一个通用的基类如下代码所示：
class ExtDesc
{
public:
std::string _ext;
std::string _version;
std::string _desc;
};
using ExtDescs = std::vector<ExtDesc>;
class CELLFormat :public CELLObject
{
public:
friend class CELLFormatMgr;
protected:
ExtDescs _descs;
public:
/// <summary>
/// 判断是否支持对应扩展名的文件
/// </summary>
bool support(const std::string& ext,const std::string& version = "1.0.0.0")
{
for (auto& var : _descs)
{
if (var._ext != ext )
continue;
if (var._version != version )
continue;
return true;
}
return false;
}
/// <summary>
/// 判断是否有子节点
/// </summary>
virtual ObjectPtr read(const char* fileName) = 0;
};
using FormatPtr = std::shared_ptr<CELLFormat>;
3.5.1 dem文件读取实现
图1为5米格网的DEM，图2为ALOS 12.5米分辨率的DEM，图3为ASTER GDEM V3 30米分辨率的DEM，图4为SRTM3 90米分辨率的DEM。
大多数项目上使用的是30米的数据源，可以从网上上下载，高程数据可以理解成一个二维数据，每一个元素存储的是高度信息（海拔)
GIS中加载考虑到速度，也会把瓦片切分成标准的瓦片，本节是读取cesium高程数据代码：
3.5.2 terrain文件读取实现
FILE* pFile = fopen(fileName,"rb");
if (pFile == nullptr)
return nullptr;
uint16_t data[65 * 65] = {0};
CELL::CELLFieldHeight* pHeight = dynamic_cast<CELL::CELLFieldHeight*>(user);
if (pHeight == nullptr)
{
pHeight = new CELL::FEFieldHeight();
}
pHeight->create(65,65,0);
fread(data,1,sizeof(data),pFile);
fclose(pFile);
float* pData = pHeight->data();
for (size_t i = 0 ;i < 65 * 65;++i)
{
pData[i] = float(data[i]) * 0.2f;
}
for (int j = 0, k = 65; j < 65 / 2; j++)
{
for (int i = 0; i < 65; i++)
{
auto temp = pData[j * k + i];
pData[j * k + i] = pData[(65 - j - 1) * k + i];
pData[(65 - j - 1) * k + i] = temp;
}
}
return pHeight;
3.6 任务系统实现
采用生产者与消费者模式：
一对一模式：
多对多模式：
系统中瓦片加载的典型应用如下图描述,这里注意，创建纹理必须在主线程中完成。
主要功能，支持多线程中执行任务
支持取消任务
支持任务执行失败通知
支持任务执行成功通知
class CELLTaskObserver
{
public:
/// <summary>
/// 任务取消通知
/// </summary>
virtual void onTaskCancel(CELLTask* task) = 0;
/// <summary>
/// 任务完成通知
/// </summary>
virtual void onTaskExe(CELLTask* task) = 0;
/// <summary>
/// 任务完成通知
/// </summary>
virtual void onTaskFinish(CELLTask* task) = 0;
};
class CELLTaskSystem
{
public:
class TaskThread :public CELLThread
{
public:
bool _exitFlag;
CELLTaskSystem* _system;
public:
TaskThread(CELLTaskSystem* pSystem)
{
_system = pSystem;
_exitFlag = true;
}
public:
virtual void join()
{
_exitFlag = true;
CELLThread::join();
}
virtual bool onCreate()
{
_exitFlag = false;
return false;
}
virtual bool onRun()
{
while (!_exitFlag)
{
_system->run();
}
return false;
}
virtual bool onDestroy()
{
return false;
}
};
typedef std::vector<CELLThread*> ArrayThread;
typedef std::list<TaskPtr> ArrayTask;
public:
CELLTaskObserver* _observer;
ArrayThread _threads;
ArrayTask _tasks;
CELLSemaphore _semphore;
mutable CELLMutex _mutex;
public:
CELLTaskSystem()
:_observer(0)
{}
/// <summary>
///
/// </summary>
virtual ~CELLTaskSystem()
{}
/// <summary>
/// 设置观察者指针
/// </summary>
virtual void setObserver(CELLTaskObserver* observer)
{
_observer = observer;
}
/// <summary>
/// 启动任务处理
/// </summary>
virtual void start(int threadNum = 4)
{
destroy();
for (int i = 0 ;i < threadNum ; ++ i)
{
TaskThread* pThread = new TaskThread(this);
pThread->start();
_threads.push_back(pThread);
}
}
/// <summary>
/// 销毁
/// </summary>
virtual void destroy()
{
for (size_t i = 0 ;i < _threads.size() ; ++ i)
{
_threads[i]->join();
delete _threads[i];
}
_threads.clear();
}
/// <summary>
/// 添加任务接口
/// </summary>
virtual void addTask(CELLTask* task)
{
{
CELLMutex::ScopeLock lk(_mutex);
_tasks.push_back(task->toPtr<CELLTask>());
}
_semphore.set(1);
}
/// <summary>
/// 获取任务队列中的任务个数
/// </summary>
virtual size_t getTaskCount() const
{
CELLMutex::ScopeLock lk(_mutex);
return _tasks.size();
}
public:
virtual void run()
{
if(!_semphore.wait())
{
return;
}
TaskPtr pTask = nullptr;
{
CELLMutex::ScopeLock lk(_mutex);
if (_tasks.empty())
{
return;
}
/// 1. 取数据
pTask = _tasks.front();
_tasks.pop_front();
}
/// 2. 执行过程
/// 3. 通知过程
if (_observer && pTask)
{
_observer->onTaskExe(pTask.get());
_observer->onTaskFinish(pTask.get());
}
}
};
上一篇：重磅更新！Sermant 1.2
下一篇：界面重建——Marching c
内容
