lisp xy轴不等比缩放_解决高缩放等级下的抖动问题

最近看到 deck.gl 在 medium 上的一篇文章，介绍了解决高缩放等级下「抖动」问题的思路。

How (sometimes) assuming the Earth is “flat” helps speed up rendering in deck.glmedium.com

我写了一个简单的 DEMO 演示这种现象：

高缩放等级下的“抖动”现象

GLSL 浮点数精度

出现抖动的原因是什么呢？

deck.gl 不同于 Mapbox，墨卡托坐标转换是在 vertex shader 中完成的。在 shader 中进行墨卡托投影，我们在本专栏上一篇「绘制地图」中已经介绍过了：

Wiki - Web Mercator projection

// deck.gl/shaders/project.glslvec2 project_mercator_(vec2 lnglat) {return vec2(radians(lnglat.x) + PI,PI - log(tan(PI * 0.25 + radians(lnglat.y) * 0.5)));
}

因此我们需要将经纬度坐标传入 shader 中，然而 GLSL 中使用的是 32-bit float，不同于 JS 默认的 64-bit 双精度浮点数。在转换过程中就会丢失精度，为此 JS 提供了 Math.fround 方法从 64-bit 转换到 32-bit ：

Math.fround(-122.4000588);
// -122.40006256103516

但是在转换过程中就会出现明显的精度丢失，在高缩放等级下就会造成肉眼可见的明显偏移，即「抖动」现象。并且随着缩放等级提升，误差将越来越大。下图来自 deck.gl medium 原文。

不同缩放等级在不同纬度下的 Y 轴像素误差

为了解决这个问题，deck.gl 早期在 GLSL 中使用两个 float 拼接成 64-bit 解决，但是对于 shader 编译和解析性能都有影响，需要在 JS 和 GLSL 中频繁 encode/decode，同时也会增加 CPU 向 GPU 传递的数据带宽。

而 v6.2 版本后 deck.gl使用了一种以屏幕中心作为动态坐标原点的「Offset Coords」方案，仅仅使用 32-bit 就解决了这个问题。

偏移坐标系

很自然想到，相近的两个坐标之差正好可以将高位抹去，只需要使用 32-bit 来存储差值，精度就完全足够了。

因此我们需要选取一个固定点，用来与其他顶点取差值。那么如何选择这个固定点呢？使用视口中心点是个很好的选择，而计算固定点和每个顶点坐标之差应该在 shader 中完成。那么能不能在 CPU 中算好差值传入呢？GPU 对于这种简单运算效率较高，而如果在 CPU 中，每一帧视口中心点的经纬度坐标都（可能）在改变，需要将原始数据都计算一遍，显然不行：

// deck.gl/shaders/project.glslvec4 project_position(vec4 position) {// 处理经纬度 offsetif (project_uCoordinateSystem == COORDINATE_SYSTEM_LNGLAT_AUTO_OFFSET) {// 与视口中心点的偏移，在经纬度坐标下因此可以保留低位，32-bit float 足够用float X = position.x - project_coordinate_origin.x;float Y = position.y - project_coordinate_origin.y;// 后续着重介绍return project_offset_(vec4(X, Y, position.z, position.w));} else {// 省略正常处理经纬度 -> 世界坐标return vec4(project_mercator(position.xy) * WORLD_SCALE * u_project_scale,project_scale(position.z),position.w);}
}

因此 deck.gl 依据当前缩放等级，在正常和 Offset 两种坐标系之间切换，如果缩放等级大于 12，则需要计算出视口中心在经纬度坐标系下的坐标：

const LNGLAT_AUTO_OFFSET_ZOOM_THRESHOLD = 12;
if (coordinateZoom < LNGLAT_AUTO_OFFSET_ZOOM_THRESHOLD) {} else {// 使用 Offset 坐标，传入经纬度坐标系下的视口中心点const lng = Math.fround(viewport.longitude);const lat = Math.fround(viewport.latitude);shaderCoordinateOrigin = [lng, lat];
}

回到 vertex shader 中，最终像素空间的计算结果可以拆解成两部分：世界坐标系偏移进行正常矩阵运算以及视口中心点的最终投影结果。

// 处理 offset 和正常经纬度到世界坐标系转换
vec4 project_pos = project_position(vec4(a_pos, 0.0, 1.0));
gl_Position = u_mvp_matrix * project_pos + u_viewport_center_projection;

视口中心点最终的投影结果可以在 CPU 中每一个渲染帧计算，例如在 mapbox custom layer 中可以使用相机参数创建视图和投影矩阵。值得一提的是 uber 也开源了 viewport-mercator-project ，能够很方便地完成视图和投影矩阵的计算：

// 从 mapbox 上下文中获取相机参数
const center = this.map.getCenter();
const currentZoomLevel = this.map.getZoom();
const pitch = this.map.getPitch();
const center = this.map.getCenter();// https://github.com/uber-common/viewport-mercator-project
const viewport = new WebMercatorViewport({width: gl.drawingBufferWidth / 2,height: gl.drawingBufferHeight / 2,longitude: center.lng,latitude: center.lat,zoom: currentZoomLevel,pitch,bearing,
});// 当前缩放等级下墨卡托坐标 -> 世界坐标，在 JS 中使用 64-bit 浮点数计算
const positionPixels = viewport.projectFlat([ Math.fround(center.lng), Math.fround(center.lat) ],Math.pow(2, currentZoomLevel)
);// 视图投影后的最终结果
const projectionCenter = vec4.transformMat4([],[positionPixels[0], positionPixels[1], 0.0, 1.0],viewProjectionMatrix
);// 将视口中心点的最终投影结果传入 uniform
drawParams['u_viewport_center_projection'] = projectionCenter;

现在只剩下了最后一部分，我们已经得到了偏移坐标系下精确的经纬度差，如何使用这个差值来估计世界坐标系下的偏移量呢？

泰勒级数展开

为了使用墨卡托坐标系下的差值估计世界坐标系下的差值，可以使用泰勒级数展开。关于泰勒级数，可以阅读「怎样更好地理解并记忆泰勒展开式？」

有了

处的值，结合导数就可以对

在任意点进行估计。在偏移坐标系场景下，

就是动态的屏幕中心点的坐标，每个顶点与中心点的差就是

，而

就是偏移坐标系到世界坐标系的转换公式。

由于 X 轴方向是线性的，对于 x 就可以应用线性估计：

而 Y 轴方向的变化是非线性的，可以使用 2 级泰勒展开：

对于经纬度坐标下的偏移转换到世界坐标，其中：

在 GLSL 中使用向量运算可以快速实现上述转换公式：其中 u_pixels_per_degree 对应 [K11, K21] ，而 u_pixels_per_degree2 对应 [0, K22]：

// offset:[delta lng, delta lat]
vec4 project_offset(vec4 offset) {float dy = offset.y;dy = clamp(dy, -1., 1.);vec3 pixels_per_unit = u_pixels_per_degree + u_pixels_per_degree2 * dy;return vec4(offset.xyz * pixels_per_unit, offset.w);
}

对于u_pixels_per_degree 和 u_pixels_per_degree2 的计算过程如下。注释十分详细的介绍了求导过程：

export function getDistanceScales({latitude, longitude, zoom, scale, highPrecision = false}) {// Calculate scale from zoom if not providedscale = scale !== undefined ? scale : zoomToScale(zoom);const result = {};const worldSize = TILE_SIZE * scale;const latCosine = Math.cos(latitude * DEGREES_TO_RADIANS);/*** Number of pixels occupied by one degree longitude around current lat/lon:pixelsPerDegreeX = d(lngLatToWorld([lng, lat])[0])/d(lng)= scale * TILE_SIZE * DEGREES_TO_RADIANS / (2 * PI)pixelsPerDegreeY = d(lngLatToWorld([lng, lat])[1])/d(lat)= -scale * TILE_SIZE * DEGREES_TO_RADIANS / cos(lat * DEGREES_TO_RADIANS)  / (2 * PI)*/const pixelsPerDegreeX = worldSize / 360;const pixelsPerDegreeY = pixelsPerDegreeX / latCosine;/*** LngLat: longitude -> east and latitude -> north (bottom left)* UTM meter offset: x -> east and y -> north (bottom left)* World space: x -> east and y -> south (top left)** Y needs to be flipped when converting delta degree/meter to delta pixels*/result.pixelsPerDegree = [pixelsPerDegreeX, -pixelsPerDegreeY, altPixelsPerMeter];/*** Taylor series 2nd order for 1/latCosinef'(a) * (x - a)= d(1/cos(lat * DEGREES_TO_RADIANS))/d(lat) * dLat= DEGREES_TO_RADIANS * tan(lat * DEGREES_TO_RADIANS) / cos(lat * DEGREES_TO_RADIANS) * dLat*/if (highPrecision) {const latCosine2 = DEGREES_TO_RADIANS * Math.tan(latitude * DEGREES_TO_RADIANS) / latCosine;const pixelsPerDegreeY2 = pixelsPerDegreeX * latCosine2 / 2;result.pixelsPerDegree2 = [0, -pixelsPerDegreeY2, altPixelsPerDegree2];}// Main results, used for converting meters to latlng deltas and scaling offsetsreturn result;
}

至此对于偏移坐标系下的处理就完成了，高缩放等级下也不会抖动了，效果如下DEMO：

高缩放等级下无抖动

总结

Mapbox 是在 CPU 中进行墨卡托投影变换，因此并不会出现这个问题，可以查看官方的 DEMO。个人认为 Mapbox 对于数据进行了分片裁剪，有效控制了在 CPU 中进行的运算量，因此可以这么做，而 deck.gl 每次都需要渲染全量数据，所以将需要一次 log 和 tan 的投影运算放到 GPU 中就很有必要了。

另外关于 deck.gl，除了官方文档，知乎上 @hijiangtao 的系列文章 hijiangtao：从零开始学习时空数据可视化（序）也很值得关注。

本文使用到 DEMO 地址：

xiaoiver/custom-mapbox-layergithub.com

参考资料

「How (sometimes) assuming the Earth is “flat” helps speed up rendering in deck.gl」

「RFC: Improved 32-bit LNGLAT projection mode」

「Introduction to Taylor's theorem for multivariable functions」