自主三维GIS引擎笔记-实现三维球045
最小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()); } } };
