1分鐘教你完美解決地圖開發中WebGL着色器32位浮點數精度損失問題

如下內容轉載自木的樹的文章《WebGL着色器32位浮點數精度損失問題》

做者：木的樹

連接：https://www.cnblogs.com/dojo-...

來源：博客園

著做權歸做者全部。商業轉載請聯繫做者得到受權，非商業轉載請註明出處。

前言

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位的，按照這篇文章說法http://www.javashuo.com/article/p-pjfosbmo-hm.html來看，下面這個確實是獲得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));

    }

最終效果：