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

ES2049
English

引言:热力图为什么在抖动呢?

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

代码地址

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

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

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

前置背景

Web 墨卡托投影

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

image.png

image.png

从公式可以看出,纬度坐标到 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 米。

img

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

Mapbox 的做法

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

img

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 轴像素误差:

img

数据拆分为高位和低位

为了解决这个问题,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 位模式(蓝色)在相同缩放级别下出现了抖动。

img

计算差值 —— 泰勒展开

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

Web 墨卡托投影公式:

image.png

image.png

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

image.png

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

image.png

在 GLSL 中使用向量运算可以快速实现上述转换公式:其中 u_pixels_per_degree 对应 [K11, K21],而 u_pixels_per_degree2 对应 [0, K22],K11K21K22 的值通过 f 导数得到。

// 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,并很快得到了解决。

img

img

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

代码地址

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

  • 使用相对于瓦片的坐标系,可以有效解决精度问题。但是当缩放程度越来越大时,瓦片分割的时间越来越长,而且如果要解决非瓦片数据的精度问题,还需要将其换算到相应的瓦片中。
  • 将数据拆分为高位和低位,将一个 Float64Array 拆分为两个 Float32Array,虽然可以解决精度问题,但是代价是显著增加了内存开销和 GPU 计算。
  • 偏移坐标系,相近的两个坐标之差正好可以将高位抹去,只需要使用 32 位来存储差值,精度就完全足够了。

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

参考资料

作者:ES2049 | timeless
文章可随意转载,但请保留此原文链接。
非常欢迎有激情的你加入 ES2049 Studio,简历请发送至 caijun.hcj@alibaba-inc.com
阅读 532

ES2049 Studio
阿里巴巴 - CRO 技术部 - 体验技术
3.2k 声望
2.8k 粉丝
0 条评论
你知道吗?

3.2k 声望
2.8k 粉丝
宣传栏