地图开发中WebGL着色器32位浮点数精度损失问题

以下内容转载自木的树的文章《WebGL着色器32位浮点数精度损失问题》

作者：木的树

链接：https://www.cnblogs.com/dojo-lzz/p/11250327.html

来源：博客园

著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

前言

Javascript API GL是基于WebGL技术打造的3D版地图API，3D化的视野更为自由，交互更加流畅。
提供丰富的功能接口，包括点、线、面绘制，自定义图层、个性化样式及绘图、测距工具等，使开发者更加容易的实现产品构思。
充分发挥GPU的并行计算能力，同时结合WebWorker多线程技术，大幅度提升了大数据量的渲染性能。最高支持百万级点、线、面绘制，同时可以保持高帧率运行。

同步推出基于Javascript API GL的 位置数据可视化API库，欢迎体验。

问题

WebGL浮点数精度最大的问题是就是因为js是64位精度的，js往着色器里面穿的时候只能是32位浮点数，有效数是8位，精度丢失比较严重。

分析

在基础底图中，所有的要素拿到的都是瓦片里面的相对坐标，坐标范围在0-256之间。在每次渲染时都会重新实时计算瓦片相对中心点的一个偏移来计算瓦片自己的矩阵，这种情况下精度损失比较小，而且每个zoom级别都会加载新的瓦片，不会出现精度损失过大问题。

但是对于一些覆盖物，比如marker、polyline、label使用的都是经纬度，经纬度小数点后位数比较多，从js的数字传入到gl中使用的gl.FLOAT是32位浮点数，小数点只能保证到后4位或者5位。在18级会出现严重的抖动问题。

文章中提到了几种解决方案，像mapbox使用的是第二种方案，将覆盖物比如marker、polyline、polygon都按照瓦片切分，经纬都转换成瓦片网格里面的0-256数字。这种方法每次zoom变换都要按照新的网格来重新切分。尤其到了18级往后，比如室内图22级，网格非常小，导致切分时间特别长。
继续尝试发现mapbox中也有类似问题：https://github.com/mapbox/mapbox-gl-js/issues/7268

mapbox这里也是使用了转换到视空间。但这种方式并不适合我们。

继续思考，实际这个问题原因是32位浮点数有效位不够，我们要找一个相对坐标为基准，其他的覆盖物坐标都是以这个点为基准，这个相对原点的坐标保留大部分数字，剩下的相对坐标数字尽量小，这样有效位尽量留给更多的小数位。然后把这个相对坐标分为两部分Math.fround（lat），lat - Math.fround(lat)；然后两部分分别在着色器重进行计算结果在相加。

6.17号第一次按照这个逻辑执行了，搞到凌晨四点多，发现并不能解决浮点数精度问题。18号跟安哥讨论了下，首先这个高位和低位不能直接在着色器里相加后进行计算。尽管设置了highp类型的float还是不行，这里面可能是因为后面有做了一些大数的乘法计算导致精度被消磨掉了。而后有做了高位的低位分别计算最后在相加，结果也不行，猜测是因为里面做了瓦片坐标转换，有一部分256 x 2^n这种计算，导致精度损失。也有可能是在某些机型上即使设置了highp实际使用的浮点数也是32位的，按照这篇文章说法https://blog.csdn.net/abcdu1/article/details/75095781来看，下面这个确实是得到32位浮点数https://developer.mozilla.org/en-US/docs/Web/API/WebGL_API/WebGL_best_practices

map.renderEngin.gl.getShaderPrecisionFormat( map.renderEngin.gl.VERTEX_SHADER, map.renderEngin.gl.HIGH_FLOAT )

解决

最终从deck.gl中找到了一种解决方案，也是将传入的数据拆分成一个高位和低位。

project_uCoordinateOrigin使用的是地图中心点的经纬度坐标

其中着色器中的一部分关键是project_uCommonUnitsPerWorldUnit和project_uCommonUnitsPerWorldUnit2这两个uniform量。跟踪代码后发现在这里有计算：

getDistanceScales() {
        // {latitude, longitude, zoom, scale, highPrecision = false}

        let center = this.center;
        let latitude = center.lat;
        let longitude = center.lng;
        let scale = this.zoomScale(this.zoom);
        let highPrecision = true;
        // Calculate scale from zoom if not provided
        scale = scale !== undefined ? scale : this.zoomToScale(zoom);

        // assert(Number.isFinite(latitude) && Number.isFinite(longitude) && Number.isFinite(scale));
      
        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;
      
        /**
         * Number of pixels occupied by one meter around current lat/lon:
         */
        const altPixelsPerMeter = worldSize / EARTH_CIRCUMFERENCE / 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.pixelsPerMeter = [altPixelsPerMeter, altPixelsPerMeter, altPixelsPerMeter];
        result.metersPerPixel = [1 / altPixelsPerMeter, 1 / altPixelsPerMeter, 1 / altPixelsPerMeter];
      
        result.pixelsPerDegree = [pixelsPerDegreeX, pixelsPerDegreeY, altPixelsPerMeter];
        result.degreesPerPixel = [1 / pixelsPerDegreeX, 1 / pixelsPerDegreeY, 1 / altPixelsPerMeter];
      
        /**
         * Taylor series 2nd order for 1/latCosine
           f'(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;
          const altPixelsPerDegree2 = worldSize / EARTH_CIRCUMFERENCE * latCosine2;
          const altPixelsPerMeter2 = altPixelsPerDegree2 / pixelsPerDegreeY * altPixelsPerMeter;
      
          result.pixelsPerDegree2 = [0, pixelsPerDegreeY2, altPixelsPerDegree2];
          result.pixelsPerMeter2 = [altPixelsPerMeter2, 0, altPixelsPerMeter2];
        }
      
        // Main results, used for converting meters to latlng deltas and scaling offsets
        return result;
    }

对于project_uCommonUnitsPerWorldUnit来说就是计算在精度和纬度上，一度代表的瓦片像素数目。对于project_uCommonUnitsPerWorldUnit2来说这里面用了一个泰勒级数的二阶展开（咨询了下管戈，泰勒级数展开项越多代表模拟值误差越小，这里用到了第二级）主要是在着色器中在project_uCommonUnitsPerWorldUnit + project_uCommonUnitsPerWorldUnit2 * dy这里做精度补偿

这里也有一些疑点，这里数字也不小，有效位的保留也不多，难道是uniform这种能够保留的有效位多一些？（也可能是转化成了瓦片像素坐标不需要那么高的精度吧。只需要整数的瓦片位，个人猜测可能不对）

gl.uniform3f(this.project_uCommonUnitsPerWorldUnit,distanceScles.pixelsPerDegree[0],distanceScles.pixelsPerDegree[1],distanceScles.pixelsPerDegree[2]);

整体来说使用这种方案解决精度损失引起的抖动问题，为后续的点、线、面、seiya都做了精度基础。

vec2 project_offset(vec2 offset) {
      float dy = offset.y;
      // if (project_uCoordinateSystem == COORDINATE_SYSTEM_LNGLAT_AUTO_OFFSET) {
        dy = clamp(dy, -1., 1.);
      // }
      vec3 commonUnitsPerWorldUnit = project_uCommonUnitsPerWorldUnit + project_uCommonUnitsPerWorldUnit2 * dy;
      // return vec4(offset.xyz * commonUnitsPerWorldUnit, offset.w);
      return vec2(offset.xy * commonUnitsPerWorldUnit.xy);
    }

    // 返回在v3 api中的3d坐标系下的坐标, 采用高精度模式
    vec2 project_view_local_position3(vec2 latlngHigh, vec2 latlngLow) {
      vec2 centerCoordHigh = project_position(center.xy + center.zw, zoom);

      // Subtract high part of 64 bit value. Convert remainder to float32, preserving precision.
      float X = latlngHigh.x - center.x;
      float Y = latlngHigh.y - center.y;

      return project_offset(vec2(X + latlngLow.x, Y + latlngLow.y));

    }

最终效果：