Cesium 目前只支持wgs84地理坐標系和web魔卡托投影坐標系。其中地形一般是wgs84。
問題1：如何使用web魔卡托投影坐標系的影像正確的貼圖到wgs84地形上？
問題2：是否可以擴展使Cesium支持多種坐標系的影像數據？
假設現有一份wgs84地理坐標系地形數據和一份web魔卡托投影坐標系影像數據，如何完成影像貼圖？
步驟一：根據地形瓦片的四至范圍計算對應的影像瓦片
1 ImageryLayer.prototype._createTileImagerySkeletons = function(tile, terrainProvider, insertionPoint) { 2 //...省略部分代碼 3 //影像圖層的地理范圍 4 var imageryBounds = Rectangle.intersection(imageryProvider.rectangle, this._rectangle, imageryBoundsScratch); 5 //獲取影像范圍與當前地形瓦塊范圍的交集 6 var rectangle = Rectangle.intersection(tile.rectangle, imageryBounds, tileImageryBoundsScratch); 7 8 //...省略部分代碼 9 10 //獲取當前影像的切片規則 11 var imageryTilingScheme = imageryProvider.tilingScheme; 12 13 //根據影像與地形瓦塊的范圍交集計算左上角對應影像的XY行列號 14 var northwestTileCoordinates = imageryTilingScheme.positionToTileXY(Rectangle.northwest(rectangle), imageryLevel); 15 //根據影像與地形瓦塊的范圍交集計算右下角對應影像的XY行列號 16 var southeastTileCoordinates = imageryTilingScheme.positionToTileXY(Rectangle.southeast(rectangle), imageryLevel); 17 18 //...省略部分代碼 19 20 //通過起始行列號獲取該地形瓦塊需要的所有影像瓦片 21 for ( var i = northwestTileCoordinates.x; i <= southeastTileCoordinates.x; i++) { 22 23 //...省略部分代碼 24 25 for ( var j = northwestTileCoordinates.y; j <= southeastTileCoordinates.y; j++) { 26 27 //...省略部分代碼 28 29 var imagery = this.getImageryFromCache(i, j, imageryLevel); 30 surfaceTile.imagery.splice(insertionPoint, 0, new TileImagery(imagery, texCoordsRectangle, useWebMercatorT)); 31 ++insertionPoint; 32 } 33 } 34 35 return true; 36 };
其中用到了imageryTilingScheme.positionToTileXY方法，該方法做了wgs84到web魔卡托的投影變換，通過該變換獲取wgs84坐標系對應的web魔卡托投影坐標系的影像瓦片。
1 WebMercatorTilingScheme.prototype.positionToTileXY = function(position, level, result) { 2 var rectangle = this._rectangle; 3 if (!Rectangle.contains(rectangle, position)) { 4 // outside the bounds of the tiling scheme 5 return undefined; 6 } 7 8 var xTiles = this.getNumberOfXTilesAtLevel(level); 9 var yTiles = this.getNumberOfYTilesAtLevel(level); 10 11 var overallWidth = this._rectangleNortheastInMeters.x - this._rectangleSouthwestInMeters.x; 12 var xTileWidth = overallWidth / xTiles; 13 var overallHeight = this._rectangleNortheastInMeters.y - this._rectangleSouthwestInMeters.y; 14 var yTileHeight = overallHeight / yTiles; 15 16 var projection = this._projection; 17 18 var webMercatorPosition = projection.project(position); 19 var distanceFromWest = webMercatorPosition.x - this._rectangleSouthwestInMeters.x; 20 var distanceFromNorth = this._rectangleNortheastInMeters.y - webMercatorPosition.y; 21 22 var xTileCoordinate = distanceFromWest / xTileWidth | 0; 23 if (xTileCoordinate >= xTiles) { 24 xTileCoordinate = xTiles - 1; 25 } 26 var yTileCoordinate = distanceFromNorth / yTileHeight | 0; 27 if (yTileCoordinate >= yTiles) { 28 yTileCoordinate = yTiles - 1; 29 } 30 31 if (!defined(result)) { 32 return new Cartesian2(xTileCoordinate, yTileCoordinate); 33 } 34 35 result.x = xTileCoordinate; 36 result.y = yTileCoordinate; 37 return result; 38 };
下面是wgs84到web魔卡托和web魔卡托到wgs84的互換代碼：
1 WebMercatorProjection.prototype.project = function(cartographic, result) { 2 var semimajorAxis = this._semimajorAxis; 3 var x = cartographic.longitude * semimajorAxis; 4 var y = WebMercatorProjection.geodeticLatitudeToMercatorAngle(cartographic.latitude) * semimajorAxis; 5 var z = cartographic.height; 6 7 if (!defined(result)) { 8 return new Cartesian3(x, y, z); 9 } 10 11 result.x = x; 12 result.y = y; 13 result.z = z; 14 return result; 15 }; 16 17 /** 18 * Converts Web Mercator X, Y coordinates, expressed in meters, to a {@link Cartographic} 19 * containing geodetic ellipsoid coordinates. The Z coordinate is copied unmodified to the 20 * height. 21 * 22 * @param {Cartesian3} cartesian The web mercator Cartesian position to unrproject with height (z) in meters. 23 * @param {Cartographic} [result] The instance to which to copy the result, or undefined if a 24 * new instance should be created. 25 * @returns {Cartographic} The equivalent cartographic coordinates. 26 */ 27 WebMercatorProjection.prototype.unproject = function(cartesian, result) { 28 //>>includeStart('debug', pragmas.debug); 29 if (!defined(cartesian)) { 30 throw new DeveloperError('cartesian is required'); 31 } 32 //>>includeEnd('debug'); 33 34 var oneOverEarthSemimajorAxis = this._oneOverSemimajorAxis; 35 var longitude = cartesian.x * oneOverEarthSemimajorAxis; 36 var latitude = WebMercatorProjection.mercatorAngleToGeodeticLatitude(cartesian.y * oneOverEarthSemimajorAxis); 37 var height = cartesian.z; 38 39 if (!defined(result)) { 40 return new Cartographic(longitude, latitude, height); 41 } 42 43 result.longitude = longitude; 44 result.latitude = latitude; 45 result.height = height; 46 return result; 47 };
到這里，第一步差不多已經完成了。
步驟二：從web魔卡托投影到wgs84紋理貼圖坐標生成。
最初，有一個簡單的頂點着色器和包含計算web魔卡托紋理坐標的片元着色器，實現動態投影效果很好。但是在一些移動設備上，片元着色器的精度有限。這樣導致在中等縮放級別上影像塗抹問題，因為根據不同的地理坐標計算投影到web魔卡托的紋理坐標值相同。
比如：vs傳給fs頂點參數為：（114.0,30.0），在fs中計算出web魔卡托的紋理坐標為（0.5,0.5），然而當vs傳給fs頂點參數為：（114.1,30.0），在fs中計算出web魔卡托的紋理坐標也為（0.5,0.5），這樣會導致不同的地理位置紋理相同的情況。
后來，其解決方案是在頂點着色器中計算重新投影到web魔卡托的紋理坐標，而不是在片元着色器中計算。這樣就需要更多的頂點數據，使用片元着色器計算紋理坐標只需要一個四邊形。要通過頂點着色器重新投影實現相同的精度，我們需要為每個輸出像素提供一個頂點。因此，cesium使用了256x256頂點的網格，因為大多數圖像圖塊都是256x256。幸運的是，可以只創建一次網格並將其上傳到GPU，然后將其重新用於所有的投影，因此與原始的片元着色器方法相比，性能幾乎沒有變化。及采用RTT的方式，生成一張新的紋理，替代原始的web魔卡托紋理。但是由於一些GPU上的GLSL sin函數在精細范圍內不連續，會出現嚴重偽像問題，才有了Cesium現在的解決方案。
當前，通過基於CORDIC算法實現更可靠的sin函數，Cesium解決了上面問題，即使每個頂點要執行大量代碼，但在大多數GPU上，性能似乎都還不錯。
不幸的是，在某些GPU上，性能很糟糕。
因此，Cesium在CPU上計算Web墨卡托紋理坐標，並將T坐標與每個頂點一起存儲（S坐標在“地理和Web墨卡托”中是相同的）。為了使其更快，Cesium將重投影網格減小為僅2個頂點寬和64個頂點高。高度為64從技術上講，重投影精度比以前有所降低，但速度上提高了四倍。
1 function reprojectToGeographic(command, context, texture, rectangle) { 2 var reproject = context.cache.imageryLayer_reproject; 3 4 if (!defined(reproject)) { 5 reproject = context.cache.imageryLayer_reproject = { 6 vertexArray : undefined, 7 shaderProgram : undefined, 8 sampler : undefined, 9 destroy : function() { 10 if (defined(this.framebuffer)) { 11 this.framebuffer.destroy(); 12 } 13 if (defined(this.vertexArray)) { 14 this.vertexArray.destroy(); 15 } 16 if (defined(this.shaderProgram)) { 17 this.shaderProgram.destroy(); 18 } 19 } 20 }; 21 22 var positions = new Float32Array(2 * 64 * 2); 23 var index = 0; 24 for (var j = 0; j < 64; ++j) { 25 var y = j / 63.0; 26 positions[index++] = 0.0; 27 positions[index++] = y; 28 positions[index++] = 1.0; 29 positions[index++] = y; 30 } 31 32 var reprojectAttributeIndices = { 33 position : 0, 34 webMercatorT : 1 35 }; 36 37 var indices = TerrainProvider.getRegularGridIndices(2, 64); 38 var indexBuffer = Buffer.createIndexBuffer({ 39 context : context, 40 typedArray : indices, 41 usage : BufferUsage.STATIC_DRAW, 42 indexDatatype : IndexDatatype.UNSIGNED_SHORT 43 }); 44 45 reproject.vertexArray = new VertexArray({ 46 context : context, 47 attributes : [{ 48 index : reprojectAttributeIndices.position, 49 vertexBuffer : Buffer.createVertexBuffer({ 50 context : context, 51 typedArray : positions, 52 usage : BufferUsage.STATIC_DRAW 53 }), 54 componentsPerAttribute : 2 55 },{ 56 index : reprojectAttributeIndices.webMercatorT, 57 vertexBuffer : Buffer.createVertexBuffer({ 58 context : context, 59 sizeInBytes : 64 * 2 * 4, 60 usage : BufferUsage.STREAM_DRAW 61 }), 62 componentsPerAttribute : 1 63 }], 64 indexBuffer : indexBuffer 65 }); 66 67 var vs = new ShaderSource({ 68 sources : [ReprojectWebMercatorVS] 69 }); 70 71 reproject.shaderProgram = ShaderProgram.fromCache({ 72 context : context, 73 vertexShaderSource : vs, 74 fragmentShaderSource : ReprojectWebMercatorFS, 75 attributeLocations : reprojectAttributeIndices 76 }); 77 78 reproject.sampler = new Sampler({ 79 wrapS : TextureWrap.CLAMP_TO_EDGE, 80 wrapT : TextureWrap.CLAMP_TO_EDGE, 81 minificationFilter : TextureMinificationFilter.LINEAR, 82 magnificationFilter : TextureMagnificationFilter.LINEAR 83 }); 84 } 85 86 texture.sampler = reproject.sampler; 87 88 var width = texture.width; 89 var height = texture.height; 90 91 uniformMap.textureDimensions.x = width; 92 uniformMap.textureDimensions.y = height; 93 uniformMap.texture = texture; 94 95 var sinLatitude = Math.sin(rectangle.south); 96 var southMercatorY = 0.5 * Math.log((1 + sinLatitude) / (1 - sinLatitude)); 97 98 sinLatitude = Math.sin(rectangle.north); 99 var northMercatorY = 0.5 * Math.log((1 + sinLatitude) / (1 - sinLatitude)); 100 var oneOverMercatorHeight = 1.0 / (northMercatorY - southMercatorY); 101 102 var outputTexture = new Texture({ 103 context : context, 104 width : width, 105 height : height, 106 pixelFormat : texture.pixelFormat, 107 pixelDatatype : texture.pixelDatatype, 108 preMultiplyAlpha : texture.preMultiplyAlpha 109 }); 110 111 // Allocate memory for the mipmaps. Failure to do this before rendering 112 // to the texture via the FBO, and calling generateMipmap later, 113 // will result in the texture appearing blank. I can't pretend to 114 // understand exactly why this is. 115 if (CesiumMath.isPowerOfTwo(width) && CesiumMath.isPowerOfTwo(height)) { 116 outputTexture.generateMipmap(MipmapHint.NICEST); 117 } 118 119 var south = rectangle.south; 120 var north = rectangle.north; 121 122 var webMercatorT = float32ArrayScratch; 123 124 var outputIndex = 0; 125 for (var webMercatorTIndex = 0; webMercatorTIndex < 64; ++webMercatorTIndex) { 126 var fraction = webMercatorTIndex / 63.0; 127 var latitude = CesiumMath.lerp(south, north, fraction); 128 sinLatitude = Math.sin(latitude); 129 var mercatorY = 0.5 * Math.log((1.0 + sinLatitude) / (1.0 - sinLatitude)); 130 var mercatorFraction = (mercatorY - southMercatorY) * oneOverMercatorHeight; 131 webMercatorT[outputIndex++] = mercatorFraction; 132 webMercatorT[outputIndex++] = mercatorFraction; 133 } 134 135 reproject.vertexArray.getAttribute(1).vertexBuffer.copyFromArrayView(webMercatorT); 136 137 command.shaderProgram = reproject.shaderProgram; 138 command.outputTexture = outputTexture; 139 command.uniformMap = uniformMap; 140 command.vertexArray = reproject.vertexArray; 141 }
以上將整個動態投影變換的關鍵步驟已指出，利用Proj4更改相應代碼就可以擴展實現其他坐標系的影像數據動態投影。