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更改相应代码就可以扩展实现其他坐标系的影像数据动态投影。