如何解决 WebGL 绘制地理信息的精度损失问题

引言：热力图为什么在抖动呢？

Deck.GL 是 Uber 开源的地理数据渲染框架，在使用 Deck.GL 绘制热力图的时候，发现不断放大地图时，地图层明显地抖动，且热力的聚合结果也有问题。下面的 demo 展示了这个现象，黄色图层是热力图图层，黑点代表原始数据，显然不断放大地图时，热力图的点并没有和原始数据点对应，且在不断抖动。

代码地址

但官方的 demo 却并没有出现这个现象，那么问题出在了哪里呢？

如果非要说两个 demo 之间有何不同，就是数据不一样。测试 demo 的地图数据是室内地图级别的，而官方 demo 数据是在城市级别的。室内地图数据到小数点后五六位才出现不同，城市级别的数据在小数点后一两位，那么很有可能是数据精度损失导致的。

基于这个猜想，上网一搜确实有不少文章介绍由于 WebGL 的精度损失带来的问题，看来是个共性的问题，那么是如何解决的呢？我们先从数据开始分析，了解不同地理渲染框架是如何解决这个问题的。

前置背景

Web 墨卡托投影

因为地球是圆的，要将地图展示在平面上，需要通过一定的投影变换绘制到平面上。墨卡托投影又称“等角正轴圆柱投影”，其等角的特性可以保证对象形状不变性，也可以保证方向和相互位置的正确性。具体的原理不在此赘述，有兴趣的可以自行了解。Web 墨卡托投影则将地球的椭圆球体简化为原型球体，坐标转换的公式如下：

从公式可以看出，纬度坐标到 Y 轴坐标的转换是非线性的，计算不仅依赖于三角函数和对数运算，且必须在每一帧的渲染中对每个坐标都进行计算，显然会带来大量的计算开销。

精度损失

地理渲染框架中都需要将元素的经纬度坐标转换为屏幕的像素坐标，随着地图不断放大，经纬度坐标转换到像素坐标的变换矩阵的位移值越来越大，即像素坐标值越来越大。因为 JS 数据是 64 位双精度浮点数，而着色器程序 GLSL 的数据只能是 32 位双精度浮点数，因此从 JS 往着色器内部传数据时，必然会出现精度丢失现象。如果不对精度丢失问题做处理，那么当放大到一定程度以后，移动地图时就会发现这些元素出现抖动现象。

Math.fround 方法可以将数据从 64 位转换到 32 位：

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

假设坐标是 [-122.4000588, 37.7900699] ，将其转换为 32 位浮点数是 [-122.40006256103516, 37.790069580078125]，这两个点之间的距离在真实世界的差值有 0.3325 米。

那么不同的地理渲染框架是如何处理大量计算和精度损失的问题呢？

Mapbox 的做法

Mapbox 采用了瓦片坐标系。地理信息的展示要素通常是静止的，因此可以预先将地图分成若干瓦片，每个瓦片包含了实际的地理信息要素。这样每次相机发生变化时，只需要以视口内的瓦片为单位渲染数据就行了。瓦片坐标系也很好理解，下图是缩放等级 z 为 2 下的瓦片坐标系：

Mapbox 所有的过程都是在平面坐标系下进行，因此首先通过墨卡托投影将要素的经纬度坐标投影至平面坐标，在每次渲染过程中都重新实时计算瓦片相对中心点的偏移矩阵，将数据变换到瓦片坐标系中：

function pixelsToTileUnits(tile: {tileID: OverscaledTileID, tileSize: number},
 pixelValue: number, z: number): number {
    return pixelValue * (EXTENT / (tile.tileSize * Math.pow(2,
 z - tile.tileID.overscaledZ)));
}

const translation = [
  inViewportPixelUnitsUnits ? translate[0] : pixelsToTileUnits(tile, translate[0], this.transform.zoom),
  inViewportPixelUnitsUnits ? translate[1] : pixelsToTileUnits(tile, translate[1], this.transform.zoom),
  0
];

const translatedMatrix = new Float32Array(16);
mat4.translate(translatedMatrix, matrix, translation);

在得到平面坐标后，为了减少数据量，Mapbox 对某些要素进行简化，并根据瓦片信息剪裁要素，再获取当前视口内包含的瓦片，最后以瓦片为单位，渲染其包含的要素。具体的过程可以参考 Mapbox 的文章。

在渲染的这一步，每个瓦片的要素传入着色器中的坐标不是经纬度，也不是墨卡托的绝对坐标，而是相对于当前瓦片的坐标：

// 墨卡托坐标 -> 相对瓦片坐标 [0, 8192]
function transformPoint(x, y, extent, z2, tx, ty) {
  return [
    Math.round(extent * (x * z2 - tx)),
    Math.round(extent * (y * z2 - ty))];
}

由于使用的是相对瓦片坐标，GLSL 的 32 位精度足够用，因此精度问题也就不存在了。但是相应的，每次相机投影矩阵发生变化时，每个瓦片的投影矩阵也都需要重新计算。

一般缩放等级到 18 级以后，比如室内地图一般在 22 级，网格非常小，会导致切分时间非常长。考虑到用户体验，面对一组待渲染瓦片，Mapbox 按照距离屏幕中心的距离进行了排序，优先渲染中心瓦片：

// 屏幕中点坐标
const centerCoord = MercatorCoordinate.fromLngLat(this.center);
const centerPoint = new Point(numTiles * centerCoord.x - 0.5, numTiles * centerCoord.y - 0.5);

// 覆盖瓦片数组按屏幕中心距离排序
tiles.sort((a, b) => centerPoint.dist(a.canonical) - centerPoint.dist(b.canonical));

概括而言，Mapbox 将墨卡托投影计算放在 CPU 中，传入着色器程序中的坐标是相对瓦片坐标，避免了 GLSL 的精度损失，并且通过要素简化、分片剪裁等操作大幅减少数据，有效控制了在 CPU 中的运算量。

Deck.GL 的做法

Deck.GL 本身的定位是处理百万级别频繁变化的数据点，在 CPU 上进行墨卡托投影会严重影响性能。Deck.GL 将经纬度坐标直接传递给 GPU，在顶点着色器中进行转换。这样，必然会带来 JS 往 GLSL 中传数据时的精度损失问题。这些误差可能在地图范围大时还无法感知，在高缩放等级下就会造成肉眼可见的偏移，即“抖动”现象，并且随着缩放等级提升，误差将越来越大。

Deck.GL 曾测试过不同缩放等级在不同纬度下 Y 轴像素误差：

数据拆分为高位和低位

为了解决这个问题，Deck.GL v3 版本中，引入了一种在 GLSL 中模拟 64 位双精度浮点数的方法，将数据拆分为高位和低位，每个数字的高位和低位都将在 GPU 中计算：

• highPart = Math.fround(x)
• lowPart = x - highPart

然后通过使用 32 位浮点数的级联运算来模拟 64 位的浮点运算。显然代价是巨大的 GPU 消耗。例如一个 64 位除法运算需要映射到 11 个 32 位运算，64 位的矩阵运算 (mat4 to vec4) 需要 1952 个 32 位运算。

使用这种方案的确实解决了精度损失引起的抖动问题，但模拟 64 位矩阵运算严重影响了着色器编译和解析的性能，同时也会增加 CPU 向 GPU 传递的数据带宽。一些性能低的显卡驱动程序无法兼容，就算兼容可能也要几秒钟的时间来编译它，导致显示卡顿。

偏移坐标系

既想要保留高精度，又想避免过高的计算性能，在 v6.2 版本以后，Deck.GL 使用了一种以屏幕中心作为动态坐标原点的 "Offset Coords" 方案，解决了这个问题。

偏移坐标系的基本想法是，相近的两个坐标之差正好可以将高位抹去，只需要使用 32 位来存储差值，精度就完全足够了。因此我们需要选取一个固定点，用来计算差值。固定点选择视口中心点，计算偏移部分的过程应该在着色器中完成。因为每一帧的视口中心点的经纬度坐标都可能在改变，如果在 CPU 中每帧都重复进行偏移部分的矩阵运算，性能无法支撑。

下面的代码是根据当前坐标系，选择是进行处理正常经纬度还是处理经纬度的差值。

// deck.gl/shaders/project.glsl
vec4 project_position(vec4 position) {
  // 处理经纬度 offset
  if (project_uCoordinateSystem == COORDINATE_SYSTEM_LNGLAT_AUTO_OFFSET) {
    // 与视口中心点的偏移，在经纬度坐标下保留低位，32 位足够用
    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];
}

因此在顶点着色器中，最终像素空间坐标的计算结果可以拆分为两部分：世界坐标系偏移部分的矩阵运算和视口中心的投影结果：

// 处理 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 中的每一渲染帧中计算，偏移部分的矩阵运算则在着色器中完成，因此每一帧的计算只需要更新少量的 uniform，几乎可以在 CPU 或 GPU 上以零成本完成。

测试结果如下图，新的混合坐标系（黄色）具有与 64 位模式（红色）相当的精度，即使只使用了 32 位，而原先 32 位模式（蓝色）在相同缩放级别下出现了抖动。

计算差值 —— 泰勒展开

上述计算过程中，还有一个细节值得注意。如何根据世界坐标系下（经纬度）的差值，来估计墨卡托投影坐标系下的差值呢？在偏移坐标系场景下， 是动态的屏幕中心坐标，其他点与中心点的差值是 ， 函数是世界坐标系到偏移坐标系的转换函数。泰勒展开做的就是根据 处的值，结合 函数的导数，可以对 函数任意点的值进行估计。泰勒展开的级数越多，代表模拟值的误差越小。

Web 墨卡托投影公式：

由于 X 轴方向是线性的，可以使用线性估计：

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

在 GLSL 中使用向量运算可以快速实现上述转换公式：其中 u_pixels_per_degree 对应 ，而 u_pixels_per_degree2 对应 ， 的值通过 导数得到。

// 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);
}

总结

那么回到本文最开始的问题，去看热力图的源码就会发现原因。Deck.GL 的热力图模块在传递坐标时没有转换，传入着色器中的坐标精度出现了损失。通过上一节可以知道，Deck.GL 在不同版本中对精度损失问题有不同的策略，所以可能是在策略迁移过程中没有测试覆盖，导致了热力图模块仍存在问题。所以我主动提了 issue，并很快得到了解决。

使用新版本的热力图之后，问题得到了解决：

代码地址

总结一下，WebGL 渲染高缩放等级下地理信息的抖动问题，有以下几个解决方法：

• 使用相对于瓦片的坐标系，可以有效解决精度问题。但是当缩放程度越来越大时，瓦片分割的时间越来越长，而且如果要解决非瓦片数据的精度问题，还需要将其换算到相应的瓦片中。

• 将数据拆分为高位和低位，将一个 Float64Array 拆分为两个 Float32Array，虽然可以解决精度问题，但是代价是显著增加了内存开销和 GPU 计算。

• 偏移坐标系，相近的两个坐标之差正好可以将高位抹去，只需要使用 32 位来存储差值，精度就完全足够了。

除了数据抖动之外，由于 WebGL 精度损失带来的现象还有 z-fighting，Z 缓冲区精度丢失问题，本文就不再赘述，有兴趣可以网上搜索相关资料了解。

参考资料

• How (sometimes) assuming the Earth is "flat" helps speed up rendering in deck.gl
• Rendering big geodata on the fly with GeoJSON-VT